Jacobi iteration

[Bunnycalculator]

Hello, i'm studying geology and i have problems with making simple program in fortran.

we have a system of n equations
ai1x1 + ai2x2 + · · · + aiixi + · · · + ainxn = bi
Solve this equation for xi.
xi = [bi − (ai1x1 + ai2x2 + · · · + ai,i−1xi−1 + ai,i+1xi+1 + · · · + ainxn)]/aii
Summation notation greatly simplifies this statement.
xi = ( bi - suma (i/=j)aij*xj)/aii

program jacobi
real A(50,50), B(10), xold, xnew,suma
integer i,j,n

open (5, file='input.txt')
open (6, file='output.txt')

read (5,*) n
read (5,*) ((A(i,j), j=1,n), i=1,n)
read (5,*) (B(i), i=1,n)

do 70 k=1,100
xold=0
do 10 i=1,n
suma=0
do 20 j=1,n
if (j.ne.i) go to 20
suma=suma+A(i,j)*xold(j)
20 continue
xnew(i)=(b(i)-suma)/(A(i,i))
10 continue
if (A(i,i).eq.0) then
write (6,100)
go to 40
else go to 10
if (abs(xnew-xold).LT.10**(-5)) then
go to 40
else xold(j) = xnew(i)
end if
write (6,200) i, (xnew(i), i=1,n)
70 continue
100 format ('nije moguce izracunati xnovi zbog dijeljenja s nulom')
200 format ('x',I1,'=',10F8.4)
40 close (5)
close (6)
end

 

[salgerman]

A few things

You have declared xnew and xold as scalars; yet, you are using them as if they were arrays.

I don't think "if (A(i,i).eq.0) then" is going to work well...the chances of a real number being exactly zero are very low...you'd better use some tolerance as in the other "if" (and 'abs' if it can be negative).

Regarding the other "if", I wouldn't bother evaluating a power inside a loop, let alone a known result: instead of evaluating 10**(-5), simply write 0.00001

Lastly, please stay away from unit numbers 5 and 6 to attach your own files; these are predefined for standard input and output and you may need for that. There are many other number to choose from, pick two.

 

[Bunnycalculator]

Thank you very much for help!