7.2. Gauss-Seidel method

Fig. 7.2 Carl Friedrich Gauss (1840 - 1887)

Fig. 7.3 Philipp Ludwig von Seidel (1821 - 1896)

In the previous section on the Jacobi method we saw that the solution to the linear system of equations \(A \vec{x} = \vec{x}\) can be calculated by estimating the solution \(\vec{x}\) and calculate an improved estimation \(\vec{x}^{(k+1)}\) using values from the current estimation \(\vec{x}^{(k)}\). We continue to calculate the improved estimates until the largest absolute value in the {prf:ref}residual` is less than some convergence tolerance. When calculating the values in \(\vec{x}^{(k+1)}\) we only use values from the previous estimate \(\vec{x}^{k}\)

\[\begin{split} \begin{align*} x_i^{(k+1)} &= \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^n a_{ij} x_j^{(k)} \right), \\ &= \frac{1}{a_{ii}} \left( b_i - \underbrace{a_{i1} x_1^{(k)} - \cdots - a_{i,j-1} x_{j-1}^{(k)}}_{\mathsf{already\,calculated}} - \underbrace{a_{i,j+1} x_{j+1}^{(k)} - \cdots - a_{in} x_n^{(k)}}_{\mathsf{yet\, to\, be\, calculated}} \right) \end{align*} \end{split}\]

so when calculating \(x_i\), we have already calculated the next values of \(x_1, x_2, \ldots, x_{i-1}\). We can improve the speed of which the Jacobi method converges by using these new values to calculate \(x_i^{(k+1)}\). This gives us the Gauss-Seidel method named after German mathematicians Carl Gauss and Philipp von Seidel.

Definition 7.3 (The Gauss-Seidel method)

The Gauss-Seidel method for solving a system of linear equations of the form \(A \vec{x} = \vec{b}\) is

(7.6)\[ x_i^{(k+1)} = \frac{1}{a_{ii} }\left(b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} -\sum_{j=i+1}^n a_{ij} x_j^{(k)} \right), \qquad i=1, \ldots , n \]

7.2.1. The iteration matrix for the Gauss-Seidel method

The iteration matrix for the Gauss-Seidel method can be found by rearranging \((L+D+U)\vec{x}=\vec{b}\)

\[\begin{split} \begin{align*} (L+D+U)\vec{x}&=\vec{b}\\ (L+D)\vec{x}&=\vec{b}-U\vec{x}\\ \vec{x}&=-(L+D)^{-1} U\vec{x}+(L+D)^{-1} \vec{x}. \end{align*} \end{split}\]

So the matrix equation for the Gauss-Seidel method is

\[ \begin{align*} \vec{x}^{(k+1)} =-(L+D)^{-1} U\vec{x}^{(k)} +(L+D)^{-1} \vec{x}^{(k+1)}, \end{align*} \]

and the iteration matrix is

(7.7)\[ T_{GS} =-(L+D)^{-1} U. \]

Example 7.2

Calculate the solution to the system of linear equations from Example 7.1 (shown below) using the Gauss-Seidel method with an accuracy tolerance of \(tol = 10^{-4}\)

\[\begin{split} \begin{align*} 4x_1 +3x_2 &=-2,\\ 3x_1 +4x_2 -x_3 &=-8,\\ -x_2 +4x_3 &=14. \end{align*} \end{split}\]

Solution (click to show)

The Gauss-Seidel method for this system is

\[\begin{split} \begin{align*} x_{1}^{(k+1)} &= \frac{1}{4} \left( -2 - 3 x_{2}^{(k)} \right), \\ x_{2}^{(k+1)} &= \frac{1}{4} \left( -8 - 3 x_{1}^{(k+1)} + x_{3}^{(k)} \right), \\ x_{3}^{(k+1)} &= \frac{1}{4} \left( 14 + x_{2}^{(k+1)} \right). \end{align*} \end{split}\]

Using starting values of \(\vec{x} = \vec{0}\). Calculating the first iteration

\[\begin{split} \begin{align*} x_{1}^{(1)} &= \frac{1}{4} \left( -2 - 3 (-1.625) \right) = -0.5, \\ x_{2}^{(1)} &= \frac{1}{4} \left( -8 - 3 (-0.5) + -3.09375 \right) = -1.625, \\ x_{3}^{(1)} &= \frac{1}{4} \left( 14 + 1.625 \right) = 3.09375. \end{align*} \end{split}\]

Calculate the residual

\[\begin{split} \begin{align*} \vec{r}^{(1)} = \vec{b} - A \vec{x}^{(1)} = \begin{pmatrix} -2 \\ -8 \\ 14 \end{pmatrix} - \begin{pmatrix} 4 & 3 & 0 \\ 3 & 4 & -1 \\ 0 & -1 & 4 \end{pmatrix} \begin{pmatrix} -0.5 \\ -1.625 \\ 3.09375 \end{pmatrix} = \begin{pmatrix} 4.875 \\ 3.09375 \\ 0.0 \end{pmatrix}. \end{align*} \end{split}\]

Since \(\max(| \vec{r}^{(1)} |) = 4.875 > 10^{-4}\) we continue iterating. Calculating the second iteration

\[\begin{split} \begin{align*} x_{1}^{(2)} &= \frac{1}{4} \left( -2 - 3 (-1.765625) \right) = 0.71875, \\ x_{2}^{(2)} &= \frac{1}{4} \left( -8 - 3 (0.71875) + -3.058594 \right) = -1.765625, \\ x_{3}^{(2)} &= \frac{1}{4} \left( 14 + 1.765625 \right) = 3.058594. \end{align*} \end{split}\]

Calculate the residual

\[\begin{split} \begin{align*} \vec{r}^{(2)} = \vec{b} - A \vec{x}^{(1)} = \begin{pmatrix} -2 \\ -8 \\ 14 \end{pmatrix} - \begin{pmatrix} 4 & 3 & 0 \\ 3 & 4 & -1 \\ 0 & -1 & 4 \end{pmatrix} \begin{pmatrix} 0.71875 \\ -1.765625 \\ 3.058594 \end{pmatrix} = \begin{pmatrix} 0.421875 \\ -0.035156 \\ 0 \end{pmatrix}. \end{align*} \end{split}\]

Since \(\max(\| \vec{r}^{(2)} \|) = 0.421875 > 10^{-4}\) we continue iterating. The Gauss-Seidel method was iterated until \(\|\vec{r}\| < 10^{-4}\) and a selection of the iteration values are given in the table below. Note that the Gauss-Seidel method took 20 iterations to achieve convergence to \(tol=10^{-4}\) whereas the Jacobi method took 49 iterations to achieve the same accuracy.

\(k\)	\(x_{1}\)	\(x_{2}\)	\(x_{3}\)	\(\max( \| \vec{r} \|)\)
0	0.000000	0.000000	0.000000	14.000000
1	-0.500000	-1.625000	3.093750	4.88e+00
2	0.718750	-1.765625	3.058594	4.22e-01
3	0.824219	-1.853516	3.036621	2.64e-01
4	0.890137	-1.908447	3.022888	1.65e-01
5	0.931335	-1.942780	3.014305	1.03e-01
\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)
19	0.999905	-1.999921	3.000020	1.43e-04
20	0.999940	-1.999950	3.000012	8.93e-05

7.2.2. Code

The code below defines the function gauss_seidel() which solves a linear system of equations of the for \(A \vec{x} = \vec{b}\) using the Gauss-Seidel method.

Python

def gauss_seidel(A, b, tol=1e-6):
    n = len(b)
    x = np.zeros(n)
    maxiter = 100
    for k in range(maxiter):
        for i in range(n):
            sum_ = 0
            for j in range(n):
                if i != j:
                    sum_ += A[i,j] * x[j]
        
            x[i] = (b[i] - sum_) / A[i,i]
            
        r = b - np.dot(A, x)

        if max(abs(r)) < tol:
            break
    
    return x

MATLAB

function x = gauss_seidel(A, b, tol)

n = length(b);
x = zeros(n, 1);
maxiter = 100;
for k = 1 : 100
    for i = 1 : n
        sum_ = 0;
        for j = 1 : n
            if j ~= i
                sum_ = sum_ + A(i,j) * x(j);
            end
        end
        x(i) = (b(i) - sum_) / A(i,i);
    end
    r = b - A * x;
    if max(abs(r)) < tol
        break
    end
end

end