How to stop iteration

[proteus01527]

Hello,
I have some problems with the code in which I would like for the following formula:

g(i,j)=0.25*(p(i,j+1)+p(i,j-1)+p(i+1,j)+p(i-1,j))

to keep repeating itself until there is some very small difference between previous and current result. Then the iteration would stop and print the results. I have set a small example 4 x 4 array in which the inner values are unknown and are computed as averages of surrounding values. The only known values are those at the borders:

1 1 1 1
4 ? ? 9
4 ? ? 9
3 3 3 3

Here's the code I wrote I tought that would work (in fortran77):

real p(4,4), g(4,4), eps
integer i, j
eps=0.0001

do i=2,3
do j=1,4
p(1,j)=1
p(i,4)=9
p(4,j)=3
p(i,1)=4
enddo
enddo

do while ((g(i,j)-p(i,j)).GT.eps)
do i=2,3
do j=2,3
g(i,j)=0.25*(p(i,j+1)+p(i,j-1)+p(i+1,j)+p(i-1,j))
enddo
enddo
do i=2,3
do j=2,3
p(i,j)=g(i,j)
enddo
enddo
if ((g(i,j)-p(i,j)).LT.eps) then
goto 1
else
continue
endif
enddo

1 do i=1,4
do j=1,4
print 2, p(i,j)
2 format (1X,'',F5.2)
enddo
enddo

end

However, the results I get are wrong. It looks like the it doesn't do enough iterations, even at such small eps value. Since I've just began to learn programming I guess I did something wrong in above code. These are the results from that code (those should be the unknown "?" inner values of that small example array above):

1.35
2.54
1.79
3.10

But they are wrong. It should be:

3.37
4.62
3.87
5.12

Thanks for any help with the above code

AndreaS

 

[xwb]

I'm really surprised you got any answers at all.

1) G is not initialized.
2) In the statement

        do while ((g(i,j)-p(i,j)).GT.eps)
        ...
        if ((g(i,j)-p(i,j)).LT.eps) then

i and j are undefined: they could be 3 ad 4 or 4 and 5, depending on the compiler
3) In the statements

        do i=2,3
        do j=2,3
        p(i,j)=g(i,j)
        enddo
        enddo

the values are clobbered so where you do the test after that, whichever index you are checking, the result will be zero so it will drop out of the loop straight away.

 

[ArkM]

And don't forget the wonderful abs function to compare deltas with eps...

 

[proteus01527]

ok, I somehow got it working now