Convergence of Gauss-Seidel Method

Jamie Trahan, Autar Kaw, Kevin Martin 

University of South Florida

United States of America

kaw@eng.usf.edu

Introduction

Gauss-Seidel method is an advantageous approach to solving a system of simultaneous linear equations because it allows the user to control round-off error that is inherent in elimination methods such as Guassian Elimination. However, this method is not without its pitfalls. Gauss-Seidel method is an iterative technique whose solution may or may not converge. Convergence is only ensured is the coefficient matrix, (Subscript[[A], nxn]),is diagonally dominant, otherwise the method may or may not converge.

A diagonally dominant square matrix [A] is defined by the following:

\(\(|\)\(a\_\(i, 
                i\)\)\(|\)\(\(≥\)\(\ \)\(\((∑\_\(j = 1\)\%n\(\(|\)\(a\_\(i, j\)\)\(|\)\))\)\_\(i ≠ j\)\)\)\)      for all i, and

|(Subscript[a, i, 
                i])|>(Subscript[(((∑_(j = 1)^n)|(Subscript[a, i, 
                            i])|)), i ≠ j])       for at least one i.

Fortunately, many physical systems that result in simultaneous linear equations have diagonally dominant coefficient matrices, or with the exchange of a few equations, the coefficient matrix can become diagonally dominant. To learn more about diagonally dominant matrices as well as how to perform Gauss-Seidel method, click here.

The following simulation illustrates the convergence of Gauss-Seidel method.

Section 1: Input

The following are the input parameters for the simulation. The user may change values only in this section. Once entered, Mathematica will produce plots that demonstrate the convergence of each solution (Subscript[X, i]) as a function of the iteration number.

Number of equations, n

n = 3

nxn coefficient matrix, [A]

A = 
        Table[{{12, 7, 3, 1}, {1, 5, 1, 2}, {2, 7, -11, 1}, {9, 2, 1, 
              13}}];
A // MatrixForm

nx1 right hand side array, [RHS]

RHS = Table[{22, 7, -2, 3}];
RHS // MatrixForm

nx1 initial guess of the solution vector

X = Table[{1, 2, 1, 3.001}];
X // MatrixForm

maximum number of iterations

maxit = 8

Section 2: Gauss-Seidel Procedure

The following procedure will use Gauss-Seidel method to calculate the value of the solution for the above system of equations using maxit iterations. It will then store each approximate solution, (Subscript[X, i]), from each iteration in a matrix with maxit columns. Thereafter, Mathematica will plot the solutions as a function of the iteration number.

Variable Parameters:
A = nxn coefficient matrix
RHS = nx1 right hand side array
n = number of equations
Xinitial = nx1 initial guess solution vector
maxit = maximum number of iterations

gausssiedel[A_, RHS_, n_, Xinitial_, maxit_] := 
        Module[{}, 
epsa = 
            Array[x, n]; 
Xnew = 
            Array[x, n]; 
epsmax = 
            Array[x, maxit]; 
Xstore = 
            Array[s, {n, maxit}]; 
Xprev = 
            Xinitial; 
For[k = 1, 
            k ≤ 
              maxit, k++, 
epsmax[[k]] = 
              0; 
For[i = 1, 
              i ≤ n, i++, 
summ = 
                0; 
For[j = 1, 
                j ≤ n, j++, 
If[
                  i ≠ j, 
                  summ = summ + 
                      A[[i, j]]*
                        Xprev[[j]]]]; 
Xnew[[i]] = 
                N[(RHS[[i]] - summ)/
                    A[[i, i]]]; 
epsa[[i]] = 
                Abs[(Xnew[[i]] - Xprev[[i]])/
                      Xnew[[i]]]*100; 
If[
                epsmax[[k]] ≤ epsa[[i]], 
                epsmax[[k]] = 
                  epsa[[i]]]; 
Xprev[[i]] = 
                Xnew[[i]]; 
Xstore[[i, k]] = 
                Xnew[[i]];]];];

Section 3: Results

gausssiedel[A, RHS, n, X, maxit];

The following table displays the value of the solution for (Subscript[X, i]) in the (i^th) row after each given iteration.

Xvar = Table[(Subscript["X", J]), {J, 1, n, 1}];
Iter = 
        Table["iteration"[i], {i, 1, maxit, 
            1}];
TableForm[Xstore, TableHeadings -> {Xvar, Iter}]

The following matrix stores the maximum absolute relative approximate error after each iteration.

Iter = 
        Table["iteration"[i], {i, 1, maxit, 
            1}];
Print["Maximum absolute relative approximate error for each \
iteration:"]
TableForm[epsmax, TableDirections -> {Row}, 
      TableHeadings -> {Iter}]

Section 4: Convergence Graphs

The graphs below plot the value of (Subscript[X, i]) as a function of the iteration number in order to display the convergence of each solution.

\(For[i = 1, 
      i ≤ 
        n, \(i++\), 
Print["  "]; 
ttl = 
        Print["Value of", " ", "X"\_i, " ", "as a function of the iteration number"]; 
data = 
        Table[{k, Xstore[\([i, k]\)]}, {k, 1, 
            maxit}]; 
ListPlot[
        data, 
PlotJoined -> 
          True, 
PlotRange -> 
          All, 
Frame -> 
          True, 
FrameLabel -> "X"\_i, 
PlotStyle -> {Thickness[0.012], 
            RGBColor[0, 0, 1]}]; 
Print[
        data]; 
Print["  "]; 
      Print["____________________________________"]; 
      Print["  "]]\)

The graph below plots the value of the maximum absolute relative approximate error as a function of the iteration number.

data2 = 
      Table[{k, epsmax[[k]]}, {k, 1, maxit, 
          1}]
\(plot2 = 
      ListPlot[data2, 
PlotLabel -> "Maximum absolute relative approximate 
error as a function of iteration number", 
PlotJoined -> True, 
Frame -> 
          True, 
FrameLabel -> "epsmax", 
PlotStyle -> {Thickness[0.012], 
            RGBColor[1, 0, 1]}]\)

Conclusion

Mathematica helped us to study the convergence of Gauss-Seidel method.

Question: Solve a set of equations for which the coefficient matrix is not diagonally dominant. For example,

5x + 6y + 7z = 18
6x + 3y + 9z = 18
7x + 9y + 10z = 26

Choose an initial solution vector guess of [2, 5, 7]. Does the solution converge? Now choose [0.99, 0.995, 0.997] as the initial guess of the solution vector. Does the solution converge now?

References

[1] Autar Kaw, Holistic Numerical Methods Institute, http://numericalmethods.eng.usf.edu/mws, See
How does Gauss-Seidel method work?