Engineering at Alberta Courses » Gauss-Seidel Method

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.

MATLAB file for the Gauss-Seidal Method Download

Lecture Video

The following video covers the Gauss-Seidel Method

Open Educational Resources

Iterative Methods: Gauss-Seidel Method

The Method

Example

Programming The Method

Lecture Video

Leave a Reply Cancel reply