Open Educational Resources

Iterative Methods: Gauss-Seidel Method

The Method

The Gauss-Seidel method offers a slight modification to the Jacobi method which can cause it to converge faster. In the Gauss-Seidel method, the system is solved using forward substitution so that each component uses the most recent value obtained for the previous component. This is different from the Jacobi method where all the components in an iteration are calculated based on the previous iteration. First notice that a linear system Ax=b of size n can be written as:

    \[\left(\begin{matrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\vdots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\\vdots\\x_n\end{array}\right)=\left(\begin{array}{c}b_1\\b_2\\\vdots\\b_n\end{array}\right)\]


The left hand side can be decomposed as follows:

    \[\left(\begin{matrix}0&a_{12}&\cdots&a_{1n}\\0&0&\cdots&a_{2n}\\\vdots&\vdots&\vdots&\vdots\\0&0&\cdots&0\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\\vdots\\x_n\end{array}\right)+\left(\begin{matrix}a_{11}&0&\cdots&0\\a_{21}&a_{22}&\cdots&0\\\vdots&\vdots&\vdots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\\vdots\\x_n\end{array}\right)=\left(\begin{array}{c}b_1\\b_2\\\vdots\\b_n\end{array}\right)\]

Effectively, we have separated A into two additive matrices: A=U+L where U is an upper triangular matrix with zero diagonal components, while L is a lower triangular matrix with its diagonal components equal to those of the diagonal components of A. If we start with nonzero diagonal components for A, then L can be used to solve the system using forward substitution:

    \[Ax=b\Rightarrow (U+L)x=b\Rightarrow Lx=b-Ux\]

the same stopping criterion as in the Jacobi method can be used for the Gauss-Seidel method.

Example

Consider the system Ax=b defined as:

    \[ \left(\begin{matrix}3&-0.1&-0.2\\0.1&7&-0.3\\0.3&-0.2&10\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\end{array}\right)=\left(\begin{array}{c}7.85\\-19.3\\71.4\end{array}\right) \]

The Gauss-Seidel method with a stopping criterion of \varepsilon_s=0.0001 will be used. First the system is rearranged to the form:

    \[\left(\begin{matrix}3&0&0\\0.1&7&0\\0.3&-0.2&10\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\end{array}\right)=\left(\begin{array}{c}7.85\\-19.3\\71.4\end{array}\right)-\left(\begin{matrix}0&-0.1&-0.2\\0&0&-0.3\\0&0&0\end{matrix}\right)\left(\begin{array}{c}x_1\\x_2\\x_3\end{array}\right)\]

Then, the initial guesses for the components x_1^{(0)}=1, x_2^{(0)}=1, x_3^{(0)}=1 are used to calculate the new estimates using forward substitution. Notice that in forward substitution, the value of x_2^{(1)} uses the value for x_1^{(1)} and that x_3^{(1)} uses the values of x_1^{(1)} and x_2^{(1)}:

    \[\begin{split}3x_1^{(1)}=7.85+0.1(1)+0.2(1)\Rightarrow x_1^{(1)}=2.716667\\0.1x_1^{(1)}+7x_2^{(1)}=-19.3+0.3(1)\Rightarrow x_2^{(1)}=-2.7531\\0.3x_1^{(1)}-0.2x_2^{(1)}+10x_3^{(1)}=71.4\Rightarrow x_3^{(1)}=7.00344\end{split}\]


The relative approximate error in this case is

    \[\varepsilon_r=\frac{\sqrt{(2.716667-1)^2+(-2.7531-1)^2+(7.00344-1)^2}}{\sqrt{2.716667^2+(-2.7531)^2+7.00344^2}}=0.91>\varepsilon_s\]

For iteration 2:

    \[\begin{split}3x_1^{(2)}=7.85+0.1(-2.7531)+0.2(7.00344)\Rightarrow x_1^{(2)}=2.99179\\0.1x_1^{(2)}+7x_2^{(2)}=-19.3+0.3(7.00344)\Rightarrow x_2^{(2)}=-2.49974\\0.3x_1^{(2)}-0.2x_2^{(2)}+10x_3^{(2)}=71.4\Rightarrow x_3^{(2)}=7.00025\end{split}\]


The relative approximate error in this case is

    \[\varepsilon_r=\frac{\sqrt{(2.99179-2.716667)^2+(-2.49974+2.7531)^2+(7.00025-7.00344)^2}}{\sqrt{2.99179^2+(-2.49974)^2+7.00025^2}}=0.047>\varepsilon_s\]

For iteration 3:

    \[\begin{split}3x_1^{(3)}=7.85+0.1(-2.49974)+0.2(7.00025)\Rightarrow x_1^{(3)}=3.00003\\0.1x_1^{(3)}+7x_2^{(3)}=-19.3+0.3(7.00025)\Rightarrow x_2^{(3)}=-2.49999\\0.3x_1^{(3)}-0.2x_2^{(3)}+10x_3^{(3)}=71.4\Rightarrow x_3^{(3)}=7.\end{split}\]


The relative approximate error in this case is

    \[\varepsilon_r=\frac{\sqrt{(3.00003-2.99179)^2+(-2.49999+2.49974)^2+(7.-7.00025)^2}}{\sqrt{3.00003^2+(-2.49999)^2+7.^2}}=0.00103>\varepsilon_s\]


And for iteration 4:

    \[\begin{split}3x_1^{(4)}=7.85+0.1(-2.49999)+0.2(7)\Rightarrow x_1^{(4)}=3.\\0.1x_1^{(4)}+7x_2^{(4)}=-19.3+0.3(7.)\Rightarrow x_2^{(4)}=-2.5\\0.3x_1^{(4)}-0.2x_2^{(4)}+10x_3^{(4)}=71.4\Rightarrow x_3^{(4)}=7.\end{split}\]


The relative approximate error in this case is

    \[\varepsilon_r=\frac{\sqrt{(3.-3.00003)^2+(-2.5+2.49999)^2+(7.-7.)^2}}{\sqrt{3^2+(-2.5)^2+7.^2}}=3.4\times 10^{-6}<\varepsilon_s\]

Therefore, the system converges after iteration 4 similar to the Jacobi method. It is not obvious in this example, but the Gauss-Seidel method can converge faster than the Jacobi method.

Programming The Method

Using forward substitution in the equation Lx^{(k+1)}=b-Ux^{(k)}, then, the i^{th} component x_i^{(k+1)} in iteration k+1 can be computed as follows:

(1)   \begin{equation*}x_i^{(k+1)}=\frac{b_i-\sum_{j=i+1}^nA_{ij}x_j^{(k)}-\sum_{j=1}^{i-1}A_{ij}x_j^{(k+1)}}{A_{ii}}\end{equation*}

The following code utilizes the procedure “GS” that takes a matrix A, a vector b, and a vector x of some guesses to return the new guess for the components of the vector x.

View Mathematica Code
GS[A_, b_, x_] :=
 (
  n = Length[A];
  xnew = Table[0, {i, 1, n}];
  Do[xnew[[i]] = (b[[i]] - Sum[A[[i, j]] x[[j]], {j, i + 1, n}] - Sum[A[[i, j]] xnew[[j]], {j, 1, i - 1}])/A[[i, i]], {i, 1, n}];
  xnew
  )
A = {{3, -0.1, -0.2}, {0.1, 7, -0.3}, {0.3, -0.2, 10}}
b = {7.85, -19.3, 71.4}
x = {{1, 1, 1.}};
MaxIter = 100;
ErrorTable = {1};
eps = 0.001;
i = 2;

While[And[i <= MaxIter, Abs[ErrorTable[[i - 1]]] > eps], 
 xi = GS[A, b, x[[i - 1]]]; x = Append[x, xi]; 
 ei = Norm[x[[i]] - x[[i - 1]]]/Norm[x[[i]]]; 
 ErrorTable = Append[ErrorTable, ei]; i++]
x // MatrixForm
ErrorTable // MatrixForm
View Python Code
import numpy as np
def GS(A, b, x):
  A = np.array(A)
  n = len(A)
  xnew = np.zeros(n)
  for i in range(n):
    xnew[i] = (b[i] - sum([A[i, j]*x[j] for j in range(i+1, n)]) - sum([A[i, j]*xnew[j] for j in range(i)]))/A[i, i]
  return xnew

A = [[3, -0.1, -0.2], [0.1, 7, -0.3], [0.3, -0.2, 10]]
b = [7.85, -19.3, 71.4]
x = [[1, 1, 1.]]
MaxIter = 100
ErrorTable = [1]
eps = 0.001
i = 1

while i <= MaxIter and abs(ErrorTable[i - 1]) > eps:
  xi = GS(A, b, x[i - 1])
  x.append(xi) 
  ei = np.linalg.norm(x[i] - x[i - 1])/np.linalg.norm(x[i]) 
  ErrorTable.append(ei)
  i+=1
print("x:",np.array(x))
print("ErrorTable:",np.vstack(ErrorTable))

The example MATLAB code below is written to emphasize the approach and explicitly works with L and U. However, the algorithm can be programmed more compactly, for example using the forward substitution form presented in Equation 1.

Lecture Video

The following video covers the Gauss-Seidel Method

Leave a Reply

Your email address will not be published.