## Jacobi iteration

## Jacobi iteration

(OP)

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

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

## RE: Jacobi iteration

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.

## RE: Jacobi iteration