7.1. Jacobi method#

Carl Jacobi

Fig. 7.1 Carl Gustav Jacob Jacobi (1804 - 1851)#

The Jacobi method named after German mathematician Carl Jacobi is the simplest indirect method. Given a system of linear equations of the form \(A \vec{x} = \vec{b}\), splitting the coefficient matrix \(A\) into the of elements from the lower triangular, diagonal and upper triangular parts of \(A\) to form matrices \(L\), \(D\) and \(U\) such that \(A = L + D + U\), e.g.,

\[\begin{split} \begin{align*} A \qquad \quad &= \qquad \quad L \qquad \quad + \qquad \quad D \qquad \qquad + \quad \qquad U \\ \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} &= \begin{pmatrix} 0 & 0 & 0 \\ a_{21} & 0 & 0 \\ a_{31} & a_{32} & 0 \end{pmatrix} + \begin{pmatrix} a_{11} & 0 & 0 \\ 0 & a_{22} & 0 \\ 0 & 0 & a_{33} \end{pmatrix} + \begin{pmatrix} 0 & a_{12} & a_{13} \\ 0 & 0 & a_{23} \\ 0 & 0 & 0 \end{pmatrix}. \end{align*} \end{split}\]

Rewriting the linear system \(A\vec{x}=\vec{b}\) using \(L\), \(D\) and \(U\) gives

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

Let the \(\vec{x}\) on the left-hand side be \(\vec{x}^{(k+1)}\) and the \(\vec{x}\) on the right-hand side be \(\vec{x}^{(k)}\) then

(7.2)#\[ \vec{x}^{(k+1)} = D^{-1} (\vec{b} - (L + U)\vec{x}^{(k)}), \]

and writing this out for each element gives the Jacobi method gives the following definition of the Jacobi method.

Definition 7.1 (The Jacobi method)

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

(7.3)#\[\begin{split} x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j = 1,\\j \neq i}^n a_{ij} x_j^{(k)} \right), \qquad i = 1, \ldots ,n. \end{split}\]

7.1.1. The iteration matrix for the Jacobi method#

The iteration matrix for the Jacobi method can be determined by rearranging equation (7.2)

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

and comparing to equation (7.1) we get the iteration matrix for the Jacobi method

(7.4)#\[ T_J = - D^{-1}(L + U). \]

7.1.2. The residual#

The Jacobi method is applied by iterating equation (7.3) until the estimate \(\vec{x}^{(k+1)}\) is accurate enough for our needs. Since we do not know what the exact solution is, we need a way to estimate the error in our approximations. Since \(\vec{x}^{(k)}\) is an approximation of the exact solution \(\vec{x}\) then if \(\vec{e}^{(k)}\) is the error of the \(k\)th iteration we have \(\vec{x} = \vec{x}^{(k)} + \vec{e}^{(k)}\). Substituting this into the linear system \(A\vec{x} = \vec{b}\) and rearranging gives

\[\begin{split} \begin{align*} A (\vec{x}^{(k)} +\vec{e}^{(k)}) &= \vec{b} \\ A\vec{e}^{(k)} &= \vec{b} - A\vec{x}^{(k)}. \end{align*} \end{split}\]

Let \(\vec{r}^{(k)} = A\vec{e}^{(k)}\) be an \(n \times 1\) column vector known as the residual.

Definition 7.2 (The residual)

The residual for the system of linear equations \(A \vec{x}^{(k)} = \vec{b}\) is

(7.5)#\[ \vec{r} = \vec{b} - A \vec{x}^{(k)}.\]

As \(\vec{x}^{(k)} \to \vec{x}\), \(\vec{r} \to 0\) so we can use the following convergence criteria

\[ \max(\|\vec{r}\|) < tol, \]

where \(tol\) is some small number. The smaller the value of \(tol\) the closer \(\vec{x}^{(k)}\) is to the exact solution but of course this will require more iterations. In practice a compromise is made between the accuracy required and the computational resources available. Typical values of \(tol\) are around \(10^{-4}\) or \(10^{-6}\).

Example 7.1

Calculate the solution to the following system of linear equations using the Jacobi 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 Jacobi 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)} + x_{3}^{(k)} \right), \\ x_{3}^{(k+1)} &= \frac{1}{4} \left( 14 + x_{2}^{(k)} \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 (-2.0) \right) = -0.5, \\ x_{2}^{(1)} &= \frac{1}{4} \left( -8 - 3 (-0.5) + -3.5 \right) = -2.0, \\ x_{3}^{(1)} &= \frac{1}{4} \left( 14 + 2.0 \right) = 3.5. \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 \\ -2.0 \\ 3.5 \end{pmatrix} = \begin{pmatrix} 6.0 \\ 5.0 \\ -2.0 \end{pmatrix}. \end{align*} \end{split}\]

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

\[\begin{split} \begin{align*} x_{1}^{(2)} &= \frac{1}{4} \left( -2 - 3 (-0.75) \right) = 1.0, \\ x_{2}^{(2)} &= \frac{1}{4} \left( -8 - 3 (1.0) + -3.0 \right) = -0.75, \\ x_{3}^{(2)} &= \frac{1}{4} \left( 14 + 0.75 \right) = 3.0. \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} 1.0 \\ -0.75 \\ 3.0 \end{pmatrix} = \begin{pmatrix} -3.75 \\ -5.0 \\ 1.25 \end{pmatrix}. \end{align*} \end{split}\]

Since \(\max(\| \vec{r}^{(2)} \|) = 5.0 > 10^{-4}\) we continue iterating. The Jacobi method was iterated until \(\max(\|\vec{r}\|) < 10^{-4}\) and a selection of the iteration values are given in the table below.

\(k\)

\(x_{1}\)

\(x_{2}\)

\(x_{3}\)

\(\max(\|\vec{r}\|)\)

0

0.000000

0.000000

0.000000

14.000000

1

-0.500000

-2.000000

3.500000

6.00e+00

2

1.000000

-0.750000

3.000000

5.00e+00

3

0.062500

-2.000000

3.312500

3.75e+00

4

1.000000

-1.218750

3.000000

3.12e+00

5

0.414062

-2.000000

3.195312

2.34e+00

\(\vdots\)

\(\vdots\)

\(\vdots\)

\(\vdots\)

\(\vdots\)

48

1.000000

-1.999975

3.000000

1.01e-04

49

0.999981

-2.000000

3.000006

7.57e-05

So the Jacobi method took 49 iterations to converge to the solution \(x_1 =1\), \(x_2 =-2\) and \(x_3 = 3\).

7.1.3. Code#

The code below defines a function called jacobi() which solves a linear system of equations of the form \(A \vec{x} = \vec{b}\) using the Jacobi method.

import numpy as np

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

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


 # Define system
A = np.array([[4, 3, 0], [3, 4, -1], [0, -1, 4]])
b = np.array([-2, -8, 14])

# Solve system
x = jacobi(A, b, tol=1e-4)
print(f"x = {x}")

% Define system
A = [4, 3, 0 ;
     3, 4, -1 ;
     0, -1, 4 ];
b = [-2 ; -8 ; 14];

% Solve system
tol = 1e-4;
x = jacobi(A, b, tol)

% --------------------------------------------------------------
function x = jacobi(A, b, tol)

n = length(b);
x = zeros(n, 1);
maxiter = 100;
for k = 1 : 100
    xold = x;
    for i = 1 : n
        sum_ = 0;
        for j = 1 : n
            if j ~= i
                sum_ = sum_ + A(i,j) * xold(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