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, cases, an initial estimate of the parameters parameters is estimated and and then the equations are solved, yielding an updated updated version version of the parameters. parameters. These new values values are then inserted back into the equations and the process continues until the desired solution is reached. reached. The two iterative methods discussed here are the Jacobi method and and the Gauss-Seidel method.
Jacobi Iteration Method Given a set of linear equations, a 1 x 1 + a 2 x 2 + L + a n x n = s1 b1 x1 + b 2 x 2 + L + bn x n = s 2 M
d1 x1 + d 2 x 2 +L + d n x n = s n 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 rearranging these equations in in the form of solving solving for the unknown parameters parameters one equation equation at a time. Thus, x1 = x2 =
s1
−
a2
a1 s2
x (20 ) −
a1
−
b 12
b1 b2
a3 a1
x 3(0 ) − L −
x 1(0 ) −
b3
x 1(0 ) −
d2
an a1
x n(0 )
x 3(0 ) − L −
bn
x 2(0 ) − L −
d n −1
b2
b2
x n(0 )
M
xn =
sn an
−
d1 dn
dn
dn
x n(0−)1
The superscript superscript (0) indicates indicates the initial initial estimate of of the parameters. parameters. For the first pass, these parameter parameterss are given the value value zero. The equations equations are then solved which results in an updated value value of the parameters. These current estimates estimates are then inserted back into the equations and a newer set of parameters parameters is arrived at by solving solving these equations. equations. The process continues until the solution converges. As an example, take the following linear equations: 7x 1 + 3x 2 + x 3 = 18 2 x 1 − 9x 2 + 4 x 3 = 12 x 1 − 4 x 2 + 12 x 3 = 6
Rearrange these equations x1 =
18 7
x2 = −
x3 =
3
1
7
7
− x2 − x3
12
2
+
9
6 12
−
9 1 12
x1 +
x1 +
4 9 4
12
⇒
x 1 = 2.571 − 0.429x 2 − 0.143x 3
x3
⇒
x 2 = −1.333 + 0.222x 1 + 0.444 x 3
x2
⇒
x 3 = 0.500 − 0.083x 1 + 0.333x 2
Use as the initial estimates: x 1(0) = x2 (0) = x3(0) = 0. equations yielding new estimates of the parameters.
Insert these estimates into these
x 1 = 2.571 − 0.429 (0) − 0.143 (0) = 2.571 (1)
x 2(1) = −1.333 + 0.222 (0 ) + 0.444 (0 ) = −1.333 x (31 ) = 0.500 − 0.083 (0 ) + 0.333 (0) = 0.500 Insert these updated estimates back into original equation again, yielding ( 2)
x1
= 2.571 − 0.429 (− 1.333) − 0.143 (0.500) = 3.071
x 2( 2 ) = −1.333 + 0.222 (2.571) + 0.444 (0.500) = −0.540 x (32 ) = 0.500 − 0.083 (2.571) + 0.333 (−1.333) = −0.159 This process is continued until the desired results are obtained. The following table shows the solutions 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. SOLVING LINEAR EQUATION USING THE JACOBI ITERATION METHOD The estimated results after each iteration are shown as: Iteration 1 2 3 4 5 6
x(1) 2.57143 3.07143 2.82540 2.87141 2.85811 2.85979
x(2) -1.33333 -.53968 -.72134 -.67695 -.68453 -.68261
x(3) .50000 -.15873 .06415 .02410 .03506 .03365
dx(1) 2.57143 .50000 -.24603 .04601 -.01330 .00168
dx(2) -1.33333 .79365 -.18166 .04439 -.00757 .00192
dx(3) .50000 -.65873 .22288 -.04005 .01096 -.00142
How good are these results? result s? Lets take our equations and put them them into into an an augmented augmented matrix and solve using Gauss-Jordan elimination. Iterative Methods for Solving Linear Equations
Page 2
6 R + 31R 7 3 1 18 R ↔ R 1 − 4 12 6 R − 2R 1 − 4 12 1 3 2 1 3 2 2 − 9 4 12 2 − 9 4 12 0 −1 − 20 0 → R→ → − −R2 7 R 3 1 1 − 4 12 6 7 3 1 18 0 31 − 83 − 24 6 R +4 R 1 0 92 6 R −92 R 1 0 0 2.859 1 − 4 12 1 2 1 3 20 0 0 1 20 0 0 1 0 − 0.683 → 0 1 R→ r2 − 20 R 3 3 0 0 − 703 − 24 703 0 0 1 0.034 0 0 1 0.034 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 c c c
Program Jacobi Program to solve a linear equation using the Jacobi Iteration method 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 c c c c c c c
The data are entered into the program using a data file called jacobi.dat. It has the following row values number of equations number of variables x(1) x(2) x(3) solution for the first equation x(1) x(2) x(3) solution for the second equation x(1) x(2) x(3) solution for the third equation OPEN (4, file = 'jacobi.dat') OPEN (6, file = 'results')
c c c
5
c c c c
no is the number of equations and nv is the number of variables read(4,*) no do 5 i=1,no xn(i) = 0.d0 continue read(4,*) nv write(6,901) The coefficients for the variables are read in the matrix x with the solution to the equations being the last column do 10 i=1,no read(4,*)(x(i,j),j=1,no+1)
c c c c c
d is the coefficient for the variable that is being solved for it forms the denominator to compute the real number for the remaining coefficients d = x(i,i) do 7 j=1,no+1
Iterative Methods for Solving Linear Equations
Page 3
7 c c c c c c
10
13 15 c c c
coef(i,j) = x(i,j)/d end do Because the Jacobi method solves for the respect to the current estimates of the coefficient for the variable is made to use in the loop to compute the adjusted
unknown variable with other variables, the be zero for subsequent estimates
coef(i,i) = 0.d0 write(6,900)(x(i,j),j=1,nv+1) end do write(6,902) do 13 i=1,no write(6,900)(coef(i,j),j=1,nv+1) end do write(6,903) iter = 0 iterate is just a counter to keep track of the number of iterations iterate = iterate+1
c c c
18 20 c c c c
Solve for the estimate of the unknown parameters 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) end do end do dx is a vector showing the change in the estimate of the variable with respect to the estimated value used in the previous iteration do 50 i=1,no dx(i) = xnp(i) - xn(i)
c c c c c c
Test to see if the change is greater than the threshold If it is, then the variable iter is made equal to 1 At the beginning of each loop, this value is made equal to 0 If iter is 1 then this means to iterate again if (dabs(dx(i)).gt.0.01d0) iter = 1
c c c 50
900 901
902 903
904
Update the estimated parameter value xn(i) = xnp(i) continue write(6,904)iterate,(xn(i),i=1,nv),(dx(i),i=1,nv) if (iter.gt.0) go to 15 FORMAT(5x,4(f10.4,5x)) FORMAT(15x,'SOLVING FORMAT(15x,'SOLVIN G LINEAR EQUATION USING USING THE JACOBI JACOBI ITERATION ITERATION MET 1HOD',//,'The coefficients to the equations with the solution at th 1e end are:',/) FORMAT(//,'Rearranging FORMAT(//,'Rearran ging the equations to solve for the unknown vari 1ables yields',/,5x,'the following coefficients: ',/) FORMAT(//,'The estimated results results after each iteration iteration are shown as 1:'//,'Iteration',2x,'x(1)',9x,'x(2)', 1:'//,'Iteration',2 x,'x(1)',9x,'x(2)',9x,'x(3)',7x,'dx(1) 9x,'x(3)',7x,'dx(1)',7x,'dx(2 ',7x,'dx(2 2)',7x,'dx(3)') FORMAT(i3,4x,6(f10.5,2x)) stop end
Iterative Methods for Solving Linear Equations
Page 4
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. method. When solving solving for x2, insert the just computed value for x 1. In other words, for each calculation, the most current estimate of the parameter parameter value is used. 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: x1(1) = 2.571 − 0.429 x (20) − 0.143 x 3(0 ) = 2.571 − 0.429 (0 ) − 0.143 (0 ) = 2.571 x 2(1) = −1.333 + 0.222 x (11 ) + 0.444 x 3(0 ) = −1.333 + 0.222 (2.571) + 0.444 (0 ) = −0.762 x (31) = 0.500 − 0.083 x1(1) + 0.333 x 2(1 ) = 0.500 − 0.083 (2.571) + 0.333 (− 0.762 ) = 0.033 The next iteration is performed in a similar fashion. It can be shown as: (2)
(1 )
(1 )
x 1 = 2.571 − 0.429x 2 − 0.143x 3 = 2.571 − 0.429 (− 0.762) − 0.143 (0.033) = 2.893 x (22 ) = −1.333 + 0.222 x 1(2 ) + 0.444 x 3(1) = −1.333 + 0.222 (2.893) + 0.444 (0.033) = −0.676 x (32 ) = 0.500 − 0.083x1(2 ) + 0.333x 2(2 ) = 0.500 − 0.083 (2.893) + 0.333 (− 0.676 ) = 0.034 A Fortran program was written to solve this problem. problem. The results are shown in in the next next table. SOLVING LINEAR EQUATION USING THE GAUSS-SEIDEL ITERATION METHOD The estimated results after each iteration are shown as: Iteration 1 2 3 4
x(1) 2.57143 2.89342 2.85646 2.85957
x(2) -.76190 -.67624 -.68369 -.68273
x(3) .03175 .03347 .03406 .03412
dx(1) -2.57143 -.32200 .03696 -.00311
dx(2) .76190 -.08566 .00745 -.00096
dx(3) -.03175 -.00172 -.00060 -.00006
As with the Jacobi method, the results from the Gauss-Seidel method are essentially correct. The Fortran Fortran program used to compute compute the Jacobi iteration method was modified to solve for the Gauss-Seidel iterative method. The program program is is shown as follows: follows: c c c c
Program Gaus_sdl.for Program to solve a linear equation using the Gauss-Seidel Iteration method IMPLICIT none REAL*8 coef(3,4), d, dx(3), x(3,4), xn(3), xnp(3)
Iterative Methods for Solving Linear Equations
Page 5
INTEGER i, iter, iterate,j, no, nv DATA iterate /0/ c c c c c c c c c
The data are entered into the program using a data file called gauss.dat. It has the following row values number of equations number of variables x(1) x(2) x(3) solution for the first equation x(1) x(2) x(3) solution for the second equation x(1) x(2) x(3) solution for the third equation OPEN (4, file = 'gauss.dat') OPEN (6, file = 'results')
c c c
5
c c c c
no is the number of equations and nv is the number of variables read(4,*) no do 5 i=1,no xn(i) = 0.d0 xnp(i) = 0.d0 continue read(4,*) nv write(6,901) The coefficients for the variables are read in the matrix x with the solution to the equations being the last column do 10 i=1,no read(4,*)(x(i,j),j=1,no+1)
c c c c c
7 c c c c c c
10
13 15 c c c
d is the coefficient for the variable that is being solved for it forms the denominator to compute the real number for the remaining coefficients d = x(i,i) do 7 j=1,no+1 coef(i,j) = x(i,j)/d end do Because the Jacobi method solves for the respect to the current estimates of the coefficient for the variable is made to use in the loop to compute the adjusted
unknown variable with other variables, the be zero for subsequent estimates
coef(i,i) = 0.d0 write(6,900)(x(i,j),j=1,nv+1) end do write(6,902) do 13 i=1,no write(6,900)(coef(i,j),j=1,nv+1) end do write(6,903) iter = 0 iterate is just a counter to keep track of the number of iterations iterate = iterate+1
c c c
Solve for the estimate of the unknown parameters do 20 i=1,no xn(i) = coef(i,nv+1) do 18 j=1,nv
Iterative Methods for Solving Linear Equations
Page 6
18 20 c c c c
xn(i) = xn(i) - coef(i,j)*xn(j) end do end do dx is a vector showing the change in the estimate of the variable with respect to the estimated value used in the previous iteration do 50 i=1,no dx(i) = xnp(i) - xn(i) xnp(i) = xn(i)
c c c c c c
Test to see if the change is greater than the threshold If it is, then the variable iter is made equal to 1 At the beginning of each loop, this value is made equal to 0 If iter is 1 then this means to iterate again if (dabs(dx(i)).gt.0.01d0) iter = 1
c c c 50
900 901
902 903
904
Update the estimated parameter value xn(i) = xnp(i) continue write(6,904)iterate,(xn(i),i=1,nv),(dx(i),i=1,nv) if (iter.gt.0) go to 15 FORMAT(5x,4(f10.4,5x)) FORMAT(15x,'SOLVING FORMAT(15x,'SOLVIN G LINEAR EQUATION USING USING THE GAUSS-SEIDEL GAUSS-SEIDEL ITERATI ITERATI 1ON METHOD',//,'The coefficients to the equations with the solution 1 at the end are:',/) FORMAT(//,'Rearranging FORMAT(//,'Rearran ging the equations to solve for the unknown vari 1ables yields',/,5x,'the following coefficients: ',/) FORMAT(//,'The estimated results results after each iteration iteration are shown as 1:'//,'Iteration',2x,'x(1)',9x,'x(2)', 1:'//,'Iteration',2 x,'x(1)',9x,'x(2)',9x,'x(3)',7x,'dx(1) 9x,'x(3)',7x,'dx(1)',7x,'dx(2 ',7x,'dx(2 2)',7x,'dx(3)') 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 Gaussc Seidel iteration. An initial 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 accomplished 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
Iterative Methods for Solving Linear Equations
Page 7
c c
a(n,n) b(n)
coefficient matrix with largest values on diagonal 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 (Note: 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 Number 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 to the calling routine: c ierr = 1 return end
Iterative Methods for Solving Linear Equations
Page 8
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 process of copying copying and pasting to add more lines lines to solve the equations. 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.
Iterative Methods for Solving Linear Equations
Page 9
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 that the solution from from either of these these methods methods could be in error. For example, lets lets look at the the following equations: equations: 12x1 + 3x 2 + 4x 3 = 48 6 x1 + 15x 2 − 4 x 3 = 54.6 9x 1 − 4 x 2 + 6x 3 = 22.3 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) 1 4.00000 2 1.85111 !
13 14
2.65235 2.61850
x(2) 3.64000 3.03111
x(3) 3.71667 .14333
dx(1) 4.00000 -2.14889
dx(2) 3.64000 -.60889
dx(3) 3.71667 -3.57333
3.05960 3.07234
1.84979 1.77788
.04225 -.03384
-.01281 .01274
.11114 -.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) 1 4.00000 2 3.79778 3 3.77474
x(2) 2.04000 1.87467 1.93538
x(3) -.92333 -.73022 -.65519
dx(1) -4.00000 .20222 .02304
dx(2) -2.04000 .16533 -.06071
dx(3) .92333 -.19311 -.07503
These are markedly different results.
Iterative Methods for Solving Linear Equations
Page 10
