Module for Jacobi and Gauss-Seidel Iteration

Numerical Methods using Mathematica  
Complementary software supplements for   
Numerical Methods textbooks  
John H. Mathews and Kurtis Fink, 2004

Background

	Iterative schemes require time to achieve sufficient accuracy and are reserved for large systems of equations where there are a majority of zero elements in the matrix. Often times the algorithms are taylor-made to take advantage of the special structure such as band matrices.  Practical uses include applications in circuit analysis, boundary value problems and partial differential equations. 
	
	Iteration is a popular technique finding roots of equations.  Generalization of fixed point iteration can be applied to systems of linear equations to produce accurate results.  The method Jacobi iteration is attributed to Carl Jacobi (1804-1851) and Gauss-Seidel iteration is attributed to  Johann Carl Friedrich Gauss (1777-1855) and Philipp Ludwig von Seidel (1821-1896).

Consider that the n×n square matrix A is split into three parts, the main diagonal D, below diagonal L and above diagonal U.  We have  A = D - L - U. 

(({{(Subscript[a, 1, 1]), (Subscript[a, 1, 2]), (Subscript[a, 1, 3]), ⋯, (Subscript[a, 1, n - 2]), (Subscript[a, 1, 
                    n - 1]), (Subscript[a, 1, n])}, {(Subscript[a, 2, 1]), (Subscript[a, 2, 2]), (Subscript[a, 2, 3]), ⋯, (Subscript[a, 2, n - 2]), (Subscript[a, 2, 
                    n - 1]), (Subscript[a, 2, n])}, {(Subscript[a, 3, 1]), (Subscript[a, 3, 2]), (Subscript[a, 3, 3]), ⋯, (Subscript[a, 3, n - 2]), (Subscript[a, 3, 
                    n - 1]), (Subscript[a, 3, n])}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮, ⋮}, {(Subscript[a, n - 2, 1]), (Subscript[a, n - 2, 2]), (Subscript[a, n - 2, 3]), ⋯, (Subscript[a, n - 2, n - 2]), (Subscript[a, n - 2, 
                    n - 1]), (Subscript[a, n - 2, n])}, {(Subscript[a, n - 1, 1]), (Subscript[a, n - 1, 2]), (Subscript[a, n - 1, 3]), ⋯, (Subscript[a, n - 1, n - 2]), (Subscript[a, n - 1, 
                    n - 1]), (Subscript[a, n - 1, n])}, {(Subscript[a, n, 1]), (Subscript[a, n, 2]), (Subscript[a, n, 3]), ⋯, (Subscript[a, n, n - 2]), (Subscript[a, n, 
                    n - 1]), (Subscript[a, n, n])}})) = 

(({{(Subscript[a, 1, 1]), 0, 0, ⋯, 0, 0, 0}, {0, (Subscript[a, 2, 2]), 0, ⋯, 0, 0, 0}, {0, 0, (Subscript[a, 3, 3]), ⋯, 0, 0, 0}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮, ⋮}, {0, 0, 0, ⋯, (Subscript[a, n - 2, n - 2]), 0, 0}, {0, 0, 0, ⋯, 0, (Subscript[a, n - 1, n - 1]), 0}, {0, 0, 0, ⋯, 0, 0, (Subscript[a, n, n])}}))-(({{0, 0, 0, ⋯, 0, 0, 0}, {-(Subscript[a, 2, 1]), 0, 0, ⋯, 0, 0, 0}, {-(Subscript[a, 3, 1]), -(Subscript[a, 3, 2]), 0, ⋯, 0, 0, 0}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮, ⋮}, {-(Subscript[a, n - 2, 1]), -(Subscript[a, n - 2, 2]), \(-a\_\(n - 2, 3\)\), ⋯, 0, 0, 0}, {-(Subscript[a, n - 1, 1]), -(Subscript[a, n - 1, 2]), \(-a\_\(n - 1, 3\)\), ⋯, -(Subscript[a, n - 1, n - 2]), 0, 0}, {-(Subscript[a, n, 1]), -(Subscript[a, n, 2]), -(Subscript[a, n, 3]), ⋯, -(Subscript[a, n, n - 2]), \(-a\_\(n, n - 1\)\), 0}}))-(({{0, -(Subscript[a, 1, 2]), -(Subscript[a, 1, 3]), ⋯, -(Subscript[a, 1, n - 2]), \(-a\_\(1, n - 1\)\), -(Subscript[a, 1, n])}, {0, 0, -(Subscript[a, 2, 3]), ⋯, -(Subscript[a, 2, n - 2]), \(-a\_\(2, n - 1\)\), -(Subscript[a, 2, n])}, {0, 0, 0, ⋯, -(Subscript[a, 3, n - 2]), \(-a\_\(3, n - 1\)\), -(Subscript[a, 3, n])}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮, ⋮}, {0, 0, 0, ⋯, 0, -(Subscript[a, n - 2, n - 1]), -(Subscript[a, n - 2, n])}, {0, 0, 0, ⋯, 0, 0, -(Subscript[a, n - 1, n])}, {0, 0, 0, ⋯, 0, 0, 0}}))

Definition (Diagonally Dominant).  The matrix  A  is strictly diagonally dominant if 

		|(Subscript[a, i, i])| > (∑_(j = 1)^(i - 1))|(Subscript[a, i, j])|  +  (∑_(j = i + 1)^n)|(Subscript[a, i, j])|   for   i = 1,2,...,n.

Theorem (Jacobi Iteration).  The solution to the linear system  AX = B  can be obtained starting with  (Subscript[P, 0]), and using iteration scheme 
	
		(Subscript[P, k + 1]) = (Subscript[M, J]) (Subscript[P, k]) + (Subscript[C, J])  
where  

		(Subscript[M, J])= (D^(-1))(L+U)  and  (Subscript[C, J])=(D^(-1))B.  

If  (Subscript[P, 0]) is carefully chosen a sequence {(Subscript[P, k])} is generated which converges to the solution  P,  i.e.  AP = B.  
A sufficient condition for the method to be applicable is that A is strictly diagonally dominant or diagonally dominant and irreducible.

In equation form the above equations are (Subscript[x, i]^((k)))= \(TraditionalForm\`\(b\_i - ∑\+\(j = 1, \ j ≠ i\)\%n\( a\_ij\) x\_j\%\((k - 1)\)\)\/a\_ii\).  For example, consider the following system of equations:

4  (Subscript[x, 1]) - (Subscript[x, 2]) + (Subscript[x, 3]) = 7
4  (Subscript[x, 1]) - 8  (Subscript[x, 2]) + (Subscript[x, 3]) = -21
-2 (Subscript[x, 1]) + (Subscript[x, 2]) + 5  (Subscript[x, 3]) = 15

Which has the following solution as determined by Mathematica

N[Solve[A . {x1, x2, x3} == {1, 2, 1}, {x1, x2, x3}]]

The matrices associated with this system are

A=(({{4, -1, 1}, {4, -8, 1}, {-2, 1, 5}})); X=(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}}));B=(({{7}, {-21}, {15}}));

The system of equations (Subscript[x, i]^((k)))= \(TraditionalForm\`\(b\_i - ∑\+\(j = 1, \ j ≠ i\)\%n\( a\_ij\) x\_j\%\((k - 1)\)\)\/a\_ii\) becomes

({{(Subscript[x, 1]^((k)))=((7  +  
                        (Subscript[x, 2]^((k - 1))) -   (Subscript[x, 3]^((k - 1))))/4).}, {(Subscript[x, 2]^((k)))=((-21  -  
                        4  (Subscript[x, 1]^((k - 1))) -   (Subscript[x, 3]^((k - 1))))/(-8)).}, {(Subscript[x, 3]^((k)))=((15  +  2 (Subscript[x, 1]^((k - 1))) -  
                        (Subscript[x, 2]^((k - 1))))/5).}})

The next 3 commands define functions that split a square matrix mat into a strictly upper triangular matrix  called U (0's on the diagonal and below the diagonal), a strictly lower triangular matrix called L (0's on and above the diagonal), and a diagonal matrix Diag (0's off the diagonal) so that mat=Diag-L-U.

U[mat_?MatrixQ] := 
      Table[If[i < j, -mat[[i, j]], 0], {i, Length[mat]}, {j, 
          Length[mat]}]

L[mat_?MatrixQ] := 
      Table[If[i > j, -mat[[i, j]], 0], {i, Length[mat]}, {j, 
          Length[mat]}]

MatrixForm[U[A]]

MatrixForm[L[A]]

Diag[mat_?MatrixQ] := mat + U[mat] + L[mat]

MatrixForm[Diag[A]]

The following computation shows that A = Diag[A]-L[A]-U[A]

MatrixForm[Diag[A] - L[A] - U[A]]

The result on Jacobi iteration shows that the (k+1)st approximation to the solution of AX=B can be found by (Subscript[P, k + 1]) = (Subscript[M, J]) (Subscript[P, k]) + (Subscript[C, J]) where  (Subscript[M, J])= (D^(-1))(L+U),   (Subscript[C, J])=(D^(-1))B , and  (Subscript[P, k]) is the kth approximation to the solution of AX=B. 
The next statement shows (Subscript[M, J]) for our example.

MatrixForm[mj = Inverse[Diag[A]] . (L[A] + U[A])]

and the next statement shows (Subscript[C, J]) for our example.

MatrixForm[cj = Inverse[Diag[A]] . B]

In terms of the variables, the new approximation would be computed as follows--replace (Subscript[x, 1]), (Subscript[x, 2]), and (Subscript[x, 3]) by some initial guesses and compute the new approximations to the solution as follows.  Iterate this as many times as you wish.

MatrixForm[new = mj . X + cj]

The next statement shows what happens if you replace X in the last statement by the current guess (old) for the solution and get a new approximation (new).  It starts with the guess {1.,1.,1.}.

old = {1. , 2. , 2. }

MatrixForm[new = mj . old + cj]

old = new

Now, repeat the sequence by clicking on the MatrixForm[new=mj.old+cj] followed by clicking on the old=new statement as many times as you wish.  This process is automated by the following Mathematica program.

jacobimat[amat_, bmat_, initialguess_, 
        niterations_] := 
Module[{umat, dmat, lmat, mj, cj, 
          old, new}, 
umat = U[amat]; lmat = L[amat]; 
        dmat = Diag[amat]; 
mj = 
          Inverse[dmat] . (lmat + umat); 
cj = 
          Inverse[dmat] . bmat; 
old = 
          initialguess; 
Print[(Subscript["P", 1]), " = ", 
          Flatten[old]]; 
For[ i = 1, 
          i ≤ niterations, i++, 
          new = mj . old + cj; 
old = new; 
          Print[(Subscript["P", i + 1]), " = ", 
            Flatten[new]];]; 
Flatten[
          new]
]

jacobimat[A, B, {1. , 2. , 2. }, 3]

Proof  Jacobi and Gauss-Seidel Iteration  Jacobi and Gauss-Seidel Iteration

Theorem (Gauss-Seidel Iteration).  The solution to the linear system  AX = B  can be obtained starting with  (Subscript[P, 0]), and using iteration scheme 
	
		(Subscript[P, k + 1]) = (Subscript[M, S]) (Subscript[P, k]) + (Subscript[C, S])  
where  

		(Subscript[M, S])=(((D - L))^(-1))U  and  (Subscript[C, S])=(((D - L))^(-1))B.  

If  (Subscript[P, 0]) is carefully chosen a sequence {(Subscript[P, k])} is generated which converges to the solution  P,  i.e.  AP = B.  
A sufficient condition for the method to be applicable is that A is strictly diagonally dominant or diagonally dominant and irreducible.

The Gauss-Seidel Iteration is slightly different than Jacobi iteration.  In Gauss-Seidel iteration the (k+1)st approximation to the solution of AX=B can be found by (Subscript[P, k + 1]) = (Subscript[M, s]) (Subscript[P, k]) + (Subscript[C, s]) where now (Subscript[M, s])=(((D - L))^(-1))U and (Subscript[C, s])=(((D - L))^(-1))B , and  (Subscript[P, k]) is the kth approximation to the solution of AX=B. 
For the example that we are studying, this gives:

MatrixForm[Inverse[Diag[A] - L[A]]]

MatrixForm[U[A]]

MatrixForm[ms = (Inverse[Diag[A] - L[A]] . U[A])]

MatrixForm[cs = Inverse[Diag[A] - L[A]] . B]

MatrixForm[new = ms . X + cs]

old = {1. , 2. , 2. }

MatrixForm[new = ms . old + cs]

old = new

gaussseidelmat[amat_, bmat_, initialguess_, 
        niterations_] := 
Module[{umat, dmat, lmat, ms, cs, 
          old, new}, 
umat = U[amat]; lmat = L[amat]; 
        dmat = Diag[amat]; 
ms = 
          Inverse[dmat - lmat] . umat; 
cs = 
          Inverse[dmat - lmat] . bmat; 
old = 
          initialguess; 
Print[(Subscript["P", 1]), " = ", 
          Flatten[old]]; 
For[ i = 1, 
          i ≤ niterations, i++, 
          new = ms . old + cs; 
old = new; 
          Print[(Subscript["P", i + 1]), " = ", 
            Flatten[new]];]; 
Flatten[
          new]
]

gaussseidelmat[A, B, {1. , 2. , 2. }, 3]

Proof  Jacobi and Gauss-Seidel Iteration  Jacobi and Gauss-Seidel Iteration

Computer Programs  Jacobi and Gauss-Seidel Iteration  Jacobi and Gauss-Seidel Iteration

Mathematica Subroutine (Jacobi Iteration).

\(\(\(Jacobi[A0_, B0_, P0_, 
          max_] := 
Module[{A = N[A0], B = N[B0], i, j, 
            k = 0, n = Length[P0], P = P0, 
            Pold = P0}, \ 
Print["P"\_0, " = ", 
            P]; \ 
While[\ 
            k < max, 
For[\ i = 1, 
              i ≤ 
                n, \(i++\), 
P\_\(\([[\)\(i\)\(]]\)\)\  = \ \(1\/A\_\(\([[\)\(i, i\)\(]]\)\)\) \((B\_\(\([[\)\(i\)\(]]\)\)\  + \ 
                      A\_\(\([[\)\(i, i\)\(]]\)\)\ Pold\_\(\([[\)\(i\)\(]]\)\)\  - \ ∑\+\(j = 1\)\%n A\_\(\([[\)\(i, j\)\(]]\)\)\ Pold\_\(\([[\)\(j\)\(]]\)\))\)\ ]; \ 
Print["P"\_\(k + 1\), " = ", 
              P]; \ 
Pold\  = \ 
              P; \ 
k\  = \ 
              k + 1;\ ]; \ 
Return[
            P];\ ];\)\(\ \ \)\)\)

Mathematica Subroutine (Gauss-Seidel Iteration).

\(\(\(GaussSeidel[A0_, B0_, P0_, 
          max_] := 
Module[{A = N[A0], B = N[B0], i, j, 
            k = 0, n = Length[P0], 
            P = P0}, \ 
Print["P"\_0, " = ", 
            P]; \ 
While[\ 
            k < max, 
For[\ i = 1, 
              i ≤ 
                n, \(i++\), 
P\_\(\([[\)\(i\)\(]]\)\)\  = \ \(1\/A\_\(\([[\)\(i, i\)\(]]\)\)\) \((B\_\(\([[\)\(i\)\(]]\)\)\  + \ 
                      A\_\(\([[\)\(i, i\)\(]]\)\)\ P\_\(\([[\)\(i\)\(]]\)\)\  - \ ∑\+\(j = 1\)\%n A\_\(\([[\)\(i, j\)\(]]\)\)\ P\_\(\([[\)\(j\)\(]]\)\))\)\ ]; \ 
Print["P"\_\(k + 1\), " = ", 
              P]; \ 
k\  = \ 
              k + 1;\ ]; \ 
Return[
            P];];\)\(\ \ \)\)\)

Example 1.  Use Jacobi iteration to solve the linear system  (({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{3}, {-2}, {5}, {4}})).  
Try 10, 20 and 30 iterations.

Solution 1.

Enter the matrix A, vector B and starting vector P.

A =(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}})); 
B  =  {3, -2, 5, 4};  
\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])"", \*""(Subscript[x, 4])""};\)
Print["Solve the system"];
Print[MatrixForm[A], MatrixForm[vars], " = ", 
        MatrixForm[B]];

Use 10 iterations.

P  = {1, 1, 1};  
X = Jacobi[A, B, P, 3];

Determine if the method has converged.

Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

Use 20 iterations.

P  = {0, -1, 1, 1};  
X = Jacobi[A, B, P, 20];

Determine if the method has converged.

Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

Use 30 iterations.

P  = {0, -1, 1, 1};  
X = Jacobi[A, B, P, 30];

Determine if the method has converged.

Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

Example 2.  Use Jacobi iteration to attempt solving the linear system  (({{2, 8, 3, 1}, {0, 2, -1, 4}, {7, -2, 1, 2}, {-1, 0, 5, 2}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{-2}, {4}, {3}, {5}})).
Try 10 iterations.  
Observe that something is not working.  In example 5 we will check to see if this matrix is diagonally dominant.

Solution 2.

Enter the matrix A, vector B and starting vector P.

A =(({{2, 8, 3, 1}, {0, 2, -1, 4}, {7, -2, 1, 2}, {-1, 0, 5, 2}}));  
B  = {-2, 4, 3, 5}; 
\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])"", \*""(Subscript[x, 4])""};\)
Print["Solve the system"];
Print[MatrixForm[A], MatrixForm[vars], " = ", 
        MatrixForm[B]];

Use 10 iterations.

P  = {-1, 1, 0, 1}; 
X = Jacobi[A, B, P, 10];

Determine if the method has converged.

Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

Was a solution found ?  Why ?

Print["A = ", MatrixForm[A]];

The matrix A is not diagonally dominant !

Example 3.  Use Gauss-Seidel iteration to solve the linear system  (({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{3}, {-2}, {5}, {4}})).  
Try 10, 20 iterations.  Compare the speed of convergence with Jacobi iteration.

Solution 3.

Enter the matrix A, vector B and starting vector P.

A =(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}})); 
B  =  {3, -2, 5, 4}; 
\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])"", \*""(Subscript[x, 4])""};\)
Print["Solve the system"];
Print[
        MatrixForm[A], MatrixForm[vars], " = ", 
        MatrixForm[B]];

Use 10 iterations.

P  = {0, -1, 1, 1}; 
X = Jacobi[A, B, P, 10]; 
Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

P  = {0, -1, 1, 1}; 
X = GaussSeidel[A, B, P, 10]; 
Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

Use 20 iterations.

P  = {0, -1, 1, 1}; 
X = Jacobi[A, B, P, 20]; 
Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

P  = {0, -1, 1, 1}; 
X = GaussSeidel[A, B, P, 20]; 
Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

Which method converges faster ?  Were all the iterations necessary ?  Why ?

Example 4  Use Gauss-Seidel iteration to attempt solving the linear system  (({{2, 8, 3, 1}, {0, 2, -1, 4}, {7, -2, 1, 2}, {-1, 0, 5, 2}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{-2}, {4}, {3}, {5}})).
Try 10 iterations.  
Observe that something is not working.  In example 5 we will check to see if this matrix is diagonally dominant.

Solution 4.

Enter the matrix A, vector B and starting vector P.

A =(({{2, 8, 3, 1}, {0, 2, -1, 4}, {7, -2, 1, 2}, {-1, 0, 5, 2}})); 
B  = {-2, 4, 3, 5}; 
\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])"", \*""(Subscript[x, 4])""};\)
Print["Solve the system"];
Print[
        MatrixForm[A], MatrixForm[vars], " = ", 
        MatrixForm[B]];

Use 10 iterations.

P  = {-1, 1, 0, 1};  
    X = GaussSeidel[A, B, P, 10]; 
Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

Was a solution found ?  Why ?

Print["A = ", MatrixForm[A]];

The matrix A is not diagonally dominant !

Warning.  

	Iteration does not always converge.  A sufficient condition for iteration to Jacobi iteration to converge is that A is strictly diagonally  dominant. The following subroutine will check to see if a matrix is strictly diagonally dominant.  It should be used before any call to Jacobi iteration or Gauss-Seidel iteration is made.  There exists a counter-example for which Jacobi iteration converges and Gauss-Seidel iteration does not converge.  The "true" sufficient condition for Jacobi iteration to converge is that the "spectral radius" of   M=(D^(-1))(A-D)  is less than 1, where  D is the diagonal of  A.  That is, the magnitude of the largest eigenvalue of M must be less than 1.  This condition seems harsh because numerical computation of eigenvalues is an advanced topic compared to solution of a linear system.

\(\(\(Dominant[
          A_] := 
Module[{cond = 1, 
            n = Length[A]}, 
For[\ i = 1, \ 
            i ≤ 
              n, \ \(i++\), 
\(If[\ ∑\+\(j = 1\)\%n 
                      Abs[A\_\(\([[\)\(i, j\)\(]]\)\)] > 
                  2\ Abs[A\_\(\([[\)\(i, i\)\(]]\)\)], \(cond = 0;\)\ ];\)\ ]; 
If[\ 
            cond == 
              1, 
Print["The matrix IS strictly diagonally dominant."]; 
Print["Jacobi iteration should converge."];\ , 
Print["The matrix is NOT strictly diagonally dominant."]; 
Print["Jacobi iteration will FAIL !"];\ \ ];\ \ ];\)\(\ \ \)\)\)

Example 5  Test the matrix A of Examples 2 and 4 to see if A is strictly diagonally dominant.

Solution 5.

A =(({{2, 8, 3, 1}, {0, 2, -1, 4}, {7, -2, 1, 2}, {-1, 0, 5, 2}})); 
Dominant[A];

We are done !

Aside. We can use Mathematica's Eigensystem subroutine and actually find the spectral radius of  M=(D^(-1))(A-D)  is less than 1.    This is just for fun.

d = 
        DiagonalMatrix[
          Table[(Subscript[A, [[i, i]]]), {i, Length[A]}]]; 
M = Abs[Inverse[d] . (A - d)]; 
Eigenvalues[N[M]]  
Abs[Eigenvalues[N[M]]]  
Max[Abs[Eigenvalues[N[M]]]]

Since the spectral radius is  5.37399 > 1,  iteration in Example 4 will not converge.

Aside.  Test the matrix of Example 1.

This is just for fun.

A =(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}})); 
Dominant[A];

Aside.  Test the matrix of Example 1to see if the spectral radius of  M=(D^(-1))(A-D)  is less than 1.

d = 
        DiagonalMatrix[
          Table[(Subscript[A, [[i, i]]]), {i, Length[A]}]]; 
M = Abs[Inverse[d] . (A - d)]; 
Eigenvalues[N[M]]  
Abs[Eigenvalues[N[M]]]  
Max[Abs[Eigenvalues[N[M]]]]

Since the spectral radius is  0.702014 < 1,  iteration in Example 1 will converge.

More efficient subroutines

	A tolerance can be supplied to either the Jacobi or Gauss-Seidel method which will permit it to exit the loop if convergence has been achieved.

Computer Programs  Jacobi and Gauss-Seidel Iteration  Jacobi and Gauss-Seidel Iteration

Mathematica Subroutine (Jacobi Iteration).

\(\(\(Jacobi[A0_, B0_, P0_, ϵ_, 
          max_] := 
Module[{A = N[A0], B = N[B0], e, i, j, 
            k = 0, n = Length[P0], P = P0, 
            Pold = P0}, \ 
Print["P"\_0, " = ", 
            P]; \ 
e = 1; \ 
While[\ 
            And[\ k < max, e > ϵ\ ], 
For[\ i = 1, 
              i ≤ 
                n, \(i++\), 
P\_\(\([[\)\(i\)\(]]\)\)\  = \ \(1\/A\_\(\([[\)\(i, i\)\(]]\)\)\) \((B\_\(\([[\)\(i\)\(]]\)\)\  + \ 
                      A\_\(\([[\)\(i, i\)\(]]\)\)\ Pold\_\(\([[\)\(i\)\(]]\)\)\  - \ ∑\+\(j = 1\)\%n A\_\(\([[\)\(i, j\)\(]]\)\)\ Pold\_\(\([[\)\(j\)\(]]\)\))\)\ ]; \ 
e = \@\(\((P - Pold)\) . \((P - Pold)\)\); \ 
Print["P"\_\(k + 1\), " = ", 
              P]; \ 
Pold = P; \ 
k = 
              k + 1]; \ 
Return[
            P];\ ];\)\(\ \ \)\)\)

Mathematica Subroutine (Gauss-Seidel Iteration).

\(\(\(GaussSeidel[A0_, B0_, P0_, ϵ_, 
          max_] := 
Module[{A = N[A0], B = N[B0], e, i, j, 
            k = 0, n = Length[P0], P = P0, 
            Pold = P0}, \ 
Print["P"\_0, " = ", 
            P]; \ 
e = 1; \ 
While[\ 
            And[\ k < max, e > ϵ\ ], 
Pold = 
              P; \ 
For[\ i = 1, 
              i ≤ 
                n, \(i++\), 
P\_\(\([[\)\(i\)\(]]\)\)\  = \ \(1\/A\_\(\([[\)\(i, i\)\(]]\)\)\) \((B\_\(\([[\)\(i\)\(]]\)\)\  + \ 
                      A\_\(\([[\)\(i, i\)\(]]\)\)\ P\_\(\([[\)\(i\)\(]]\)\)\  - \ ∑\+\(j = 1\)\%n A\_\(\([[\)\(i, j\)\(]]\)\)\ P\_\(\([[\)\(j\)\(]]\)\))\)\ ]; \ 
e = \@\(\((P - Pold)\) . \((P - Pold)\)\); \ 
Print["P"\_\(k + 1\), " = ", 
              P]; \ 
k = 
              k + 1]; \ 
Return[
            P];\ ];\)\(\ \ \)\)\)

Example 6.  Use Jacobi and Gauss-Seidel iteration to solve the linear system  (({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{3}, {-2}, {5}, {4}})).  
Use a tolerance of  (10^(-6))  and a maximum of  50  iterations.

Solution 6.

Enter the matrix A, vector B and starting vector P.

A =(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}})); 
B  =  {3, -2, 5, 4}; 
\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])"", \*""(Subscript[x, 4])""};\)
Print["Solve the system"];
Print[
        MatrixForm[A], MatrixForm[vars], " = ", 
        MatrixForm[B]];

P  = {0, -1, 1, 1}; 
X = Jacobi[A, B, P, 0.000001, 50]; 
Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

How many iterations did Jacobi iteration use ?

P  = {0, -1, 1, 1}; 
X = 
        GaussSeidel[A, B, P, 0.000001, 50]; 
Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

How many iterations did Gauss-Seidel iteration use ?

Subroutines using matrix commands

	In the Jacobi subroutine we can use fix point iteration as suggested by the theory.

Mathematica Subroutine (Jacobi Iteration).

\(\(\(JacobiFP[A0_, B0_, P0_, 
          max_] := 
Module[{A = N[A0], B = N[B0], i, k = 0, 
            n = Length[P0], 
            P = P0}, \ 
Print["P"\_0, " = ", 
            P]; \ 
d = 
            DiagonalMatrix[
              Table[A\_\(\([[\)\(i, i\)\(]]\)\), {i, Length[A]}]]; \ \ 
M = \(-\ 
              Inverse[d] . \((A - d)\)\); \ \ 
c = 
            Inverse[d] . \((B)\); \ \ \ 
For[k = 1, 
            k ≤ max, \(k++\), 
P = 
              c + M . P; \ \ 
Print["P"\_\(k + 1\), " = ", P];\ \ ]; \ 
Return[P];\ ];\)\(\ \ \)\)\)

Example 7.  Use the fixed point version of Jacobi iteration to solve the linear system  (({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}}))(({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}})) = (({{3}, {-2}, {5}, {4}})).

Solution 7.

A =(({{7, -2, 1, 2}, {2, 8, 3, 1}, {-1, 0, 5, 2}, {0, 2, -1, 4}})); 
B  =  {3, -2, 5, 4}; 
\(vars = {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])"", \*""(Subscript[x, 3])"", \*""(Subscript[x, 4])""};\)
Print["Solve the system"];
Print[
        MatrixForm[A], MatrixForm[vars], " = ", 
        MatrixForm[B]];

P  = {0, -1, 1, 1}; 
X = JacobiFP[A, B, P, 20]; 
Print["A X = ", MatrixForm[A . X], " ≈ ", 
        MatrixForm[B], " = B"];

Research Experience for Undergraduates

Jacobi and Gauss-Seidel Iteration  Jacobi and Gauss-Seidel Iteration  Internet hyperlinks to web sites and a bibliography of articles.