Jacobi and gauss-seidel

Numerical methods to solve system of linear equations

1. ITERATIVE METHODS FOR SOLVING LINEAR EQUATIONS


There are other methods that can be used to solve a set of linear equations that are based
on iteration. In these cases, an initial estimate of the parameters is estimated and then the
equations are solved, yielding an updated version of the parameters. These new values
are then inserted back into the equations and the process continues until the desired
solution is reached. The two iterative methods discussed here are the Jacobi method and
the Gauss-Seidel method.

Jacobi Iteration Method

Given a set of linear equations,

nnn2211
2nn2211
1nn2211
sxdxdxd
sxbxbxb
sxaxaxa
=+++
=+++
=+++
L
M
L
L


the problem is one of solving for x1, x2, …, xn. The right hand side of these equations, si,
represents the solution. We begin by rearranging these equations i n the form of solving
for the unknown parameters one equation at a time. Thus,

() () ()
() () ()
() () ()0
1n
n
1n0
2
n
20
1
n
1
n
n
n
0
n
2
n0
3
2
30
1
2
1
12
2
2
0
n
1
n0
3
1
30
2
1
2
1
1
1
x
d
d
x
d
d
x
d
d
a
s
x
x
b
b
x
b
b
x
b
b
b
s
x
x
a
a
x
a
a
x
a
a
a
s
x
-
-
----=
----=
----=
L
M
L
L


The superscript (0) indicates the initial estimate of the parameters. For the first pass,
these parameters are given the value zero. The equations are then solved which results in
an updated value of the parameters. These current estimates are then inserted back into
the equations and a newer set of parameters is arrived at by solving these equations. The
process continues until the solution converges.

As an example, take the following linear equations:

6x12x4x
12x4x9x2
18xx3x7
321
321
321
=+-
=+-
=++

2. Iterative Methods for Solving Linear Equations Page 2

Rearrange these equations

213213
312312
321321
x333.0x083.0500.0xx
12
4
x
12
1
12
6
x
x444.0x222.0333.1xx
9
4
x
9
2
9
12
x
x143.0x429.0571.2xx
7
1
x
7
3
7
18
x
+-=Þ+-=
++-=Þ++-=
--=Þ--=


Use as the initial estimates: x 1
(0)
= x 2
(0)
= x 3
(0)
= 0. Insert these estimates into these
equations yielding new estimates of the parameters.

() ()
() ()
()
() ()500.00333.00083.0500.0x
333.10444.00222.0333.1x
571.20143.00429.0571.2x
1
3
)1(
2
)1(
1
=+-=
-=++-=
=--=


Insert these updated estimates back into original equation again, yielding

( ) ( )
( ) ( )
()
( ) ( ) 159.0333.1333.0571.2083.0500.0x
540.0500.0444.0571.2222.0333.1x
071.3500.0143.0333.1429.0571.2x
2
3
)2(
2
)2(
1
-=-+-=
-=++-=
=---=


This process is continued until the desired results are obtained. The following table
shows the solutions arrived at after each iteration. These results are from the attached
Fortran program. The output shows x(i) that are the parameters x i. Also output is the
change in the parameters, dx(i), between each iteration.



How good are these results? Lets take our equations and put them into an augmented
matrix and solve using Gauss-Jordan elimination.
SOLVING LINEAR EQUATION USING THE JACOBI ITERATION METHOD

The estimated results after each iteration are shown as :

Iteration x(1) x(2) x(3) dx(1) dx(2) dx(3)
1 2.57143 -1.33333 .50000 2.57143 -1.33333 .50000
2 3.07143 -.53968 -.15873 .50000 .79365 -.65873
3 2.82540 -.72134 .06415 -.24603 -.18166 .22288
4 2.87141 -.67695 .02410 .04601 .04439 -.04005
5 2.85811 -.68453 .03506 -.01330 -.00757 .01096
6 2.85979 -.68261 .03365 .00168 .00192 -.00142

3. Iterative Methods for Solving Linear Equations Page 3


ú
ú
ú
û
ù
ê
ê
ê
ë
é
-
ú
ú
ú
û
ù
ê
ê
ê
ë
é
ú
ú
ú
û
ù
ê
ê
ê
ë
é
--
-
ú
ú
ú
û
ù
ê
ê
ê
ë
é
--
--
-
ú
ú
ú
û
ù
ê
ê
ê
ë
é
-
-
ú
ú
ú
û
ù
ê
ê
ê
ë
é
-
-
®®
®®®
-
-
+
+
-
-
-
«
034.0100
683.0010
859.2001
034.0100
02010
69201
2470300
02010
61241
2483310
02010
61241
18137
12492
61241
61241
12492
18137
31
32
21
3
23
2
12
13
31
R92R
R20r
R4R
703
R
R31R
R
R2R
R7R
RR


As one can see, the values using the Jacobi iterative method are very close. Following is
a Fortran program that can be used to use the Jacobi iteration to solve a set of equations.
The limitation now is that it is restricted to only a 3 x 3 matrix, due to formatting
procedures currently used in the program.

c Program Jacobi
c Program to solve a linear equation using the Jacobi It eration
c method
c
IMPLICIT none
REAL*8 coef(3,4), d, dx(3), x(3,4), xn(3), xnp(3)
INTEGER i, iter, iterate,j, no, nv
DATA iterate /0/
c
c The data are entered into the program using a data file called
c jacobi.dat. It has the f ollowing row values
c number of equations
c number of variables
c x(1) x(2) x(3) solution for the first equation
c x(1) x(2) x(3) solution for the second equation
c x(1) x(2) x(3) solution for the third equation
c
OPEN (4, file = 'jacobi.dat')
OPEN (6, file = 'results')
c
c no is the number of equations and nv is the number of variables
c
read(4,*) no
do 5 i=1,no
xn(i) = 0.d0
5 continue
read(4,*) nv
write(6,901)
c
c The coefficients for th e variables are read in the matrix x with
c the solution to the equations being the last column
c
do 10 i=1,no
read(4,*)(x(i,j),j=1,no+1)
c
c d is the coefficient for the variable that is being solved for
c it forms the denominator to compute th e real number for the
c remaining coefficients
c
d = x(i,i)
do 7 j=1,no+1

4. Iterative Methods for Solving Linear Equations Page 4

coef(i,j) = x(i,j)/d
7 end do
c
c Because the Jacobi method solves for the unknown variable with
c respect to the current estimates of the other variables, the
c coefficient for the variable is made to be zero for subsequent
c use in the loop to compute the adjusted estimates
c
coef(i,i) = 0.d0
write(6,900)(x(i,j),j=1,nv+1)
10 end do
write(6,902)
do 13 i=1,no
write(6,900)(coef(i,j),j=1,nv+1)
13 end do
write(6,903)
15 iter = 0
c
c iterate is just a counter to keep track of the number of iterations
c
iterate = iterate+1
c
c Solve for the estimate of the unknown parameters
c
do 20 i=1,no
xnp(i) = coef(i,nv+1)
do 18 j=1,nv
xnp(i) = xnp(i) - coef(i,j)*xn(j)
18 end do
20 end do
c
c dx is a vector showing the change in the estimate of the variable
c with respect to the estimated value used in the previous iteration
c
do 50 i=1,no
dx(i) = xnp(i) - xn(i)
c
c Test to see if the change is greater than the threshold
c If it is, then the variable iter is made equal to 1
c At the beginning of each loop, this value is made equal to 0
c If iter is 1 then this means to iterate again
c
if (dabs(dx(i)).gt.0.01d0) iter = 1
c
c Update the estimated parameter value
c
xn(i) = xnp(i)
50 continue
write(6,904)iterate,(xn(i),i=1,nv),(dx(i),i=1,nv)
if (iter.gt.0) go to 15
900 FORMAT(5x,4(f10.4,5x))
901 FORMAT(15x,'SOLVING LINEAR EQUATION USING THE JACOBI ITER ATION MET
1HOD',//,'The coefficients to the equations with the solution at th
1e end are:',/)
902 FORMAT(//,'Rearranging the equations to solve for the unknown vari
1ables yields',/,5x,'the following coefficients: ',/)
903 FORMAT(//,'The estimated results after each iteration are shown as
1:'//,'Iteration',2x,'x(1)',9x,'x(2)',9x,'x(3)',7x,'dx(1)',7x,'dx(2
2)',7x,'dx(3)')
904 FORMAT(i3,4x,6(f10.5,2x))
stop
end

5. Iterative Methods for Solving Linear Equations Page 5



Gauss-Seidel Iterative Method

The Gauss-Seidel iterative method of solving for a set of linear equations can be thought
of as just an extension of the Jacobi method. Start out using an initial value of zero for
each of the parameters. Then, solve for x1 as in the Jacobi method. When solving for x2,
insert the just computed value for x 1. In other words, for each calculation, the most
current estimate of the parameter value is used. To see how the Gauss-Seidel method
works, lets use the values in the last example. The initial estimates are set to zero. Then,
the results from the first iteration are shown as:

() ()
() ()
() ()
( ) ()
() () ()
( ) ( )033.0762.0333.0571.2083.0500.0x333.0x083.0500.0x
762.00444.0571.2222.0333.1x444.0x222.0333.1x
571.20143.00429.0571.2x143.0x429.0571.2x
1
2
1
1
1
3
0
3
1
1
)1(
2
0
3
0
2
)1(
1
=-+-=+-=
-=++-=++-=
=--=--=


The next iteration is performed in a similar fashion. It can be shown as:

() () ()
( ) ( )
() () ()
( ) ( )
() () ()
( ) ( )034.0676.0333.0893.2083.0500.0x333.0x083.0500.0x
676.0033.0444.0893.2222.0333.1x444.0x222.0333.1x
893.2033.0143.0762.0429.0571.2x143.0x429.0571.2x
2
2
2
1
2
3
1
3
2
1
2
2
1
3
1
2
2
1
=-+-=+-=
-=++-=++-=
=---=--=


A Fortran program was written to solve this problem. The results are shown in the next
table.



As with the Jacobi method, the results from the Gauss-Seidel method are essentially
correct. The Fortran program used to compute the Jacobi iteration method was modified
to solve for the Gauss-Seidel iterative method. The program is shown as follows:

c Program Gaus_sdl.for
c Program to solve a linear equation using the Gauss -Seidel
c Iteration method
c
IMPLICIT none
REAL*8 coef(3,4), d, dx(3), x(3,4), xn(3), xnp(3)
SOLVING LINEAR EQUATION USING THE GAUSS -SEIDEL ITERATION
METHOD

The estimated results after each iteration are shown as:

Iteration x(1) x(2) x(3) dx(1) dx(2) dx(3)
1 2.57143 -.76190 .03175 -2.57143 .76190 -.03175
2 2.89342 -.67624 .03347 -.32200 -.08566 -.00172
3 2.85646 -.68369 .03406 .03696 .00745 -.00060
4 2.85957 -.68273 .03412 -.00311 -.00096 -.00006

6. Iterative Methods for Solving Linear Equations Page 6

INTEGER i, iter, iterate,j, no, nv
DATA iterate /0/
c
c The data are entered into the program using a data file called
c gauss.dat. It has the following row values
c number of equations
c number of variables
c x(1) x(2) x(3) solution for the first equation
c x(1) x(2) x(3) solution for the second equation
c x(1) x(2) x(3) solution for the third equation
c
OPEN (4, file = 'gauss.dat')
OPEN (6, file = 'results')
c
c no is the number of equations and nv is the number of variables
c
read(4,*) no
do 5 i=1,no
xn(i) = 0.d0
xnp(i) = 0.d0
5 continue
read(4,*) nv
write(6,901)
c
c The coefficients for the variables are read in the matrix x with
c the solution to the equations being the last column
c
do 10 i=1,no
read(4,*)(x(i,j),j=1,no+1)
c
c d is the coefficient for the variable that is being solved for
c it forms the denominator to compute the real number for the
c remaining coefficients
c
d = x(i,i)
do 7 j=1,no+1
coef(i,j) = x(i,j)/d
7 end do
c
c Because the Jacobi method solves for the unknown variable with
c respect to the current estimates of the other variables, the
c coefficient for the variable is made to be zero for subsequent
c use in the loop to compute the adjusted estimates
c
coef(i,i) = 0.d0
write(6,900)(x(i,j),j=1,nv+1)
10 end do
write(6,902)
do 13 i=1,no
write(6,900)(coef(i,j),j=1,nv+1)
13 end do
write(6,903)
15 iter = 0
c
c iterate is just a counter to keep track of the number of iterations
c
iterate = iterate+1
c
c Solve for the estimate of the unknown parameters
c
do 20 i=1,no
xn(i) = coef(i,nv+1)
do 18 j=1,nv

7. Iterative Methods for Solving Linear Equations Page 7

xn(i) = xn(i) - coef(i,j)*xn(j)
18 end do
20 end do
c
c dx is a vector showing the change in the estimate of the variable
c with respect to the estimated value used in the previous iteration
c
do 50 i=1,no
dx(i) = xnp(i) - xn(i)
xnp(i) = xn(i)
c
c Test to see if the change is greater than the threshold
c If it is, then the variable iter is made equal to 1
c At the beginning of each loop, this value is made equal to 0
c If iter is 1 then this means to iterate again
c
if (dabs(dx(i)).gt.0.01d0) iter = 1
c
c Update the estimated parameter value
c
xn(i) = xnp(i)
50 continue
write(6,904)iterate,(xn(i),i=1, nv),(dx(i),i=1,nv)
if (iter.gt.0) go to 15
900 FORMAT(5x,4(f10.4,5x))
901 FORMAT(15x,'SOLVING LINEAR EQUATION USING THE GAUSS -SEIDEL ITERATI
1ON METHOD',//,'The coefficients to the equations with the solution
1 at the end are:',/)
902 FORMAT(//,'Rearranging the equations to solve for the unknown vari
1ables yields',/,5x,'the following coefficients: ',/)
903 FORMAT(//,'The estimated results after each iteration are shown as
1:'//,'Iteration',2x,'x(1)',9x,'x(2)',9x,'x(3)',7x, 'dx(1)',7x,'dx(2
2)',7x,'dx(3)')
904 FORMAT(i3,4x,6(f10.5,2x))
stop
end

A more sophisticated subroutine for solving the Gauss-Seidel is shown as [source
unknown]:

subroutine gsitrn(a,b,x,n,ndim,niter,tol,ierr)
c------------------ -----------------------------------------------------
c
c Gauss -Seidel Iterative Method
c *****************************
c
c This subroutine obtaines the solution to n linear equations by Gauss -
c Seidel iteration. An in itial approximation is sent to the subroutine
c in the vector x. The solution, as approximated by the subroutine is
c returned in x. The iterations are continued until the maximum change
c in any x component is less than tol. If this cannot be accomplis hed
c in niter iterations, a non -zero error flag is returned. The matrix
c is to be arranged so as to have the largest values on the diagonal.
c (from "Applied Numerical Analysis," C.F. Gerald, p 138)
c
c
c
c INPUT/OUTPUT VARIABLES:
c

8. Iterative Methods for Solving Linear Equations Page 8

c a(n,n) coefficient matrix with largest values on diagonal
c b(n) right hand side vector
c x(n) solution vector (initial guess)
c n Dimension of the system you're solving
c ndim Dimension of matrix a (No te: In the main program,
c matrix a may have been dimensioned larger than
c necessary, i.e. n, the size of the system you're
c decomposing, may be smaller than ndim.)
c niter Numbe r of iterations
c tol tolerance for solution
c ierr Error code: ierr=0 - no errors; ierr=1 - the
c solution was not obtained in maximum iterations
c
c------------------------------------------------------ -----------------
dimension a(ndim,ndim),b(ndim),x(ndim)
ierr = 0
c
c We can save some divisions by making all the diagonal
c elements equal to unity:
c
do 1 i=1,n
temp = 1.0 / a(i,i)
b(i) = b(i) * temp
do 2 j=1,n
a(i,j) = a(i,j) * temp
2 continue
1 continue
c
c Now we perform the iterations. Store max change in x values for
c testing against tol. Outer loop limits iterations to niter:
c
do 3 iter=1,niter
xmax = 0.0
do 4 i=1,n
temp = x(i)
x(i) = b(i)
do 5 j=1,n
if(i.ne.j) then
x(i) = x(i) - a(i,j)*x(j)
endif
5 continue
if(abs(x(i) -temp).gt.xmax) xmax = abs(x(i) -temp)
4 continue
if(xmax.le.tol) return
3 continue
c
c Normal exit from the loop means non -convergent solution. Flag the
c error and pass control back t o the calling routine:
c
ierr = 1
return
end

9. Iterative Methods for Solving Linear Equations Page 9

These iterative methods can also be very effectively programmed in a spreadsheet like
Excel. For example, the Jacobi method of solving for linear equations can be shown as:



The formulas for computing x1 and x2 within the spreadsheet are shown as:



It is a simple process of copying and pasting to add more lines to solve the equations. In
a similar fashion, the Gauss-Seidel method can also be programmed within Excel to
arrive at the same results, as shown in the following figure.

10. Iterative Methods for Solving Linear Equations Page 10



While the Gauss-Seidel method appears to be the best in this example, this is not always
the case. In fact, it is very possible that the solution from either of these methods could
be in error. For example, lets look at the following equations:

3.22x6x4x9
6.54x4x15x6
48x4x3x12
321
321
321
=+-
=-+
=++


Then, using a criteria of 0.1 (this is the significant figures for the input values) to stop the
iterations, the solution using the Jacobi method can be shown to be:

SOLVING LINEAR EQUATION USING THE JACOBI ITERATION METHOD

The estimated results after each iteration are shown as:

Iteration x(1) x(2) x(3) dx(1) dx(2) dx(3)
1 4.00000 3.64000 3.71667 4.00000 3.64000 3.71667
2 1.85111 3.03111 .14333 -2.14889 -.60889 -3.57333
!
13 2.65235 3.05960 1.84979 .04225 -.01281 .11114
14 2.61850 3.07234 1.77788 -.03384 .01274 -.07191


The same results using the Gauss-Seidel method, same criteria, are:


SOLVING LINEAR EQUATION USING THE GAUSS -SEIDEL ITERATION METHOD

The estimated results after each iteration are shown as:

Iteration x(1) x(2) x(3) dx(1) dx(2) dx(3)
1 4.00000 2.04000 -.92333 -4.00000 -2.04000 .92333
2 3.79778 1.87467 -.73022 .20222 .16533 -.19311
3 3.77474 1.93538 -.65519 .02304 -.06071 -.07503

These are markedly different results.