Course Webpage

Contents

Gauss-Seidel method

Consider the system $A\mathbf{x}=\mathbf{b}$ with

$$ A=\left(\begin{array}{cccc} 20 & 1 & 0 & 1\\ 1 & 20 & 3 & 1\\ 0 & 3 & 20 & 1\\ 1 & 0 & 1 & 20 \end{array}\right)\quad\mathbf{b}=\left(\begin{array}{c} 10\\ 7\\ 4\\ 6 \end{array}\right) $$

Solve it using the Gauss-Seidel method.

Use as stopping condition $\left\Vert \mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\right\Vert _{\infty}\lt 0.01$

Let us use python

In [1]:
import numpy as np
np.set_printoptions(precision=3)

We write the data that are the coefficient matrix $A$ and $\mathbf{b}$

In [2]:
A = np.array([[20., 1, 0, 1],[1, 20, 3, 1],[0, 3, 20, 1],[1, 0, 1, 20]])
b = np.array([10., 7, 4, 6])

We initialize the two vectors where we will store the last two approximate solutions

In [3]:
x0 = np.zeros(4)
x1 = np.zeros(4)

If we have a system of $4$ equations and $4$ unknowns

$$ \begin{array}{cll} a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}+a_{14}x_{4} & = & b_{1}\\ a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}+a_{24}x_{4} & = & b_{2}\\ a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}+a_{34}x_{4} & = & b_{3}\\ a_{41}x_{1}+a_{42}x_{2}+a_{43}x_{3}+a_{44}x_{4} & = & b_{4} \end{array} $$

and we solve the first equation for the first unknown, the second equation for the second unknown and so on

$$ \begin{array}{cll} x_{1}& = &\dfrac{b_{1}-a_{12}x_{2}-a_{13}x_{3}-a_{14}x_{4}}{a_{11}} \\ x_{2} & = & \dfrac{b_{2}-a_{21}x_{1}-a_{23}x_{3}-a_{24}x_{4}}{a_{22}}\\ x_{3} & = & \dfrac{b_{3}-a_{31}x_{1}-a_{32}x_{2}-a_{34}x_{4}}{a_{33}}\\ x_{4} & = & \dfrac{b_{4}-a_{41}x_{1}-a_{42}x_{2}-a_{43}x_{3}}{a_{44}} \end{array} $$

And we get the values of the first iteration from the initial values, but as we get new results, we use them

$$ \begin{array}{cll} x_{1}^{\color{red}{(1)}} & = & \dfrac{b_{1}-a_{12}x_{2}^{(0)}-a_{13}x_{3}^{(0)}-a_{14}x_{4}^{(0)}}{a_{11}} \\ x_{2}^{\color{red}{(1)}} & = & \dfrac{b_{2}-a_{21}x_{1}^{\color{red}{(1)}}-a_{23}x_{3}^{(0)}-a_{24}x_{4}^{(0)}}{a_{22}}\\ x_{3}^{\color{red}{(1)}} & = & \dfrac{b_{3}-a_{31}x_{1}^{\color{red}{(1)}}-a_{32}x_{2}^{\color{red}{(1)}}-a_{34}x_{4}^{(0)}}{a_{33}}\\ x_{4}^{\color{red}{(1)}}& = & \dfrac{b_{4}-a_{41}x_{1}^{\color{red}{(1)}}-a_{42}x_{2}^{\color{red}{(1)}}-a_{43}x_{3}^{\color{red}{(1)}}}{a_{44}} \end{array} $$

Algorithm

  • Given an initial guess $\mathbf{x}^{(0)}$
  • For $k=1,2,\ldots,$MaxIter
    • For $i=1,2,\ldots,n$, compute $$x_{i}^{\color{red}{(k)}}=\dfrac{1}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{\color{red}{(k)}}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}\right)$$
    • If the stopping condition is met, take $\mathbf{x}^{(k)}$ as approximation of the solution.

Iteration 1

From $\mathbf{x^{(0)}}$ we compute $\mathbf{x^{(1)}}$, but to update the values of an unknown as soon as we have obtained it, we can use a single vector during the iteration (only x0). However, we need to save the previous iteration to check the stopping criterion.

In [4]:
k = 1
xant = np.copy(x0)
x0[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x0[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x0[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x0[3] = (b[3] - A[3,0] * x0[0] - A[3,1] * x0[1] - A[3,2] * x0[2]) / A[3,3]

We print the new iteration, the difference vector and the norm of said difference vector, which will be used to control the error

In [5]:
print ('x0                = ', x0)
print ('x0-xant           = ', x0-xant)
print ('max(abs(x0-xant)) = ', max(abs(x0-xant)))

x0                =  [0.5   0.325 0.151 0.267]
x0-xant           =  [0.5   0.325 0.151 0.267]
max(abs(x0-xant)) =  0.5

As the norm is greater than $0.01$, we go on with the iterations

Iteration 2

From $\mathbf{x^{(1)}}$ we compute $\mathbf{x^{(2)}}$

In [6]:
k = 2
xant = np.copy(x0)
x0[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x0[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x0[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x0[3] = (b[3] - A[3,0] * x0[0] - A[3,1] * x0[1] - A[3,2] * x0[2]) / A[3,3]

We print the new iteration, the difference vector and the norm of said difference vector, which will be used to control the error

In [7]:
print ('x0                = ', x0)
print ('x0-xant           = ', x0-xant)
print ('max(abs(x0-xant)) = ', max(abs(x0-xant)))

x0                =  [0.47  0.29  0.143 0.269]
x0-xant           =  [-0.03  -0.035 -0.008  0.002]
max(abs(x0-xant)) =  0.03457828125000001

As the norm is greater than $0.01$, we go on with the iterations

Iteration 3

From $\mathbf{x^{(2)}}$ we compute $\mathbf{x^{(3)}}$

In [8]:
k = 3
xant = np.copy(x0)
x0[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x0[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x0[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x0[3] = (b[3] - A[3,0] * x0[0] - A[3,1] * x0[1] - A[3,2] * x0[2]) / A[3,3]

We print the results of this iteration

In [9]:
print ('x0                = ', x0)
print ('x0-xant           = ', x0-xant)
print ('max(abs(x0-xant)) = ', max(abs(x0-xant)))

x0                =  [0.472 0.291 0.143 0.269]
x0-xant           =  [ 1.634e-03  1.052e-03 -2.522e-04 -6.911e-05]
max(abs(x0-xant)) =  0.0016343965429687635

Now that the norm of the difference vector is less than $0.01$ and we stop. Our approximate solution is

In [10]:
print ('x0 = ', x0)

x0 =  [0.472 0.291 0.143 0.269]

For comparison sake, the exact solution is

In [11]:
x = np.linalg.solve(A,b);
print(x)

[0.472 0.292 0.143 0.269]

And the error

In [12]:
print ('error = ', max(abs(x0-x)))

error =  5.146953503498697e-05