The Jacobi method

7.1. The Jacobi method#

The Jacobi method, named after German mathematician Carl Jacobi, computes each component of the new approximation using only values from the previous iteration. Since all updates use values from iteration \(k\), each component \(x_i^{(k+1)}\) can be computed independently.

Given a linear system of the form

\[ A \mathbf{x} = \mathbf{b}, \]

then the \(i\)th equation in the system is

\[ a_{i1} x_1 + a_{i2} x_2 + \ldots + a_{ii} x_i + \ldots + b_{in} x_n = b_i. \]

Rearraning to make \(x_i\) the subject gives

\[ x_i = \frac{1}{a_{ii}} (b_i - a_{i1}x_1 + \ldots + a_{i,i-1}x_{i-1} + a_{i,i+1}x_{i+1} + \ldots + a_{in}x_n). \]

Let the value on the left-hand side be the new approximation and the values on the right-hand side by the values from the previous iteration then we can write

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

where \(x_i^{(k+1)}\) and \(x_i^{(k)}\) denote the new and current values of \(x_i\).

Definition 7.1 (The Jacobi method)

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

(7.2)#\[\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}\]

where \(a_{ii} \ne 0\).

Note

The Jacobi method will fail if \(a_{ii} = 0\). If this is the case then row swaps can be performed to ensure the pivot elements are non-zero. The convergence criterion of iterative methods is explored later.

7.1.1. The residual#

The Jacobi method is applied by iterating equation (7.2) until the approximation \(\mathbf{x}^{(k)}\) is sufficiently accurate. Since the exact solution \(\mathbf{x}\) is unknown, we cannot compute the error directly and instead require a quantity that can be computed during the iteration.

The residual of the \(k\)-th iteration is defined as

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

The residual measures how closely the current approximation satisfies the original linear system \(A \mathbf{x} = \mathbf{b}\). If \(\mathbf{x}^{(k)}\) is the exact solution, then

\[ A \mathbf{x}^{(k)} = \mathbf{b}, \]

and therefore \(\mathbf{r}^{(k)} = \mathbf{0}\). To understand the relationship between the residual and the error let

\[ \mathbf{e}^{(k)} = \mathbf{x} - \mathbf{x}^{(k)}, \]

be the error at the \(k\)-th iteration. Since

\[ \mathbf{x} = \mathbf{x}^{(k)} + \mathbf{e}^{(k)}, \]

substituting into the linear system \(A \mathbf{x} = \mathbf{b}\) gives

\[ A ( \mathbf{x}^{(k)} + \mathbf{e}^{(k)} ) = \mathbf{b}. \]

Rearranging gives

\[ A \mathbf{e}^{(k)} = \mathbf{b} - A \mathbf{x}^{(k)}, \]

hence

\[ \mathbf{r}^{(k)} = A \mathbf{e}^{(k)}. \]

As the approximations converge to the exact solution, the error \(\mathbf{e}^{(k)}\) tends to zero, and consequently the residual \(\mathbf{r}^{(k)}\) also tends to zero. For this reason the residual is used to monitor convergence and determine when the iterations should cease.

Definition 7.2 (Residual)

The residual corresponding to the approximation \(\mathbf{x}^{(k)}\) is

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

A common convergence criterion is

\[ \| \mathbf{r}^{(k)} \|_\infty < \text{tol}, \]

where

\[ \| \mathbf{r}^{(k)} \|_\infty = \max_i |r_i^{(k)}|,\]

is the infinity norm.

The smaller the value of \(\text{tol}\) the closer \(\mathbf{x}^{(k)}\) is to the exact solutionm 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

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 \(\mathbf{x} = \mathbf{0}\). Computing the first iteration

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

Compute the residual

\[\begin{split} \begin{align*} \mathbf{r}^{(1)} = \mathbf{b} - A \mathbf{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 \(\| \mathbf{r}^{(1)} \|_\infty = 6.0 > 10^{-4}\) we continue iterating. Computing the second iteration

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

Compute the residual

\[\begin{split} \begin{align*} \mathbf{r}^{(2)} = \mathbf{b} - A \mathbf{x}^{(2)} = \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 \(\| \mathbf{r}^{(2)} \|_\infty = 5.0 > 10^{-4}\) we continue iterating.

The Jacobi method was iterated until \(\|\mathbf{r}\|_\infty < 10^{-4}\) and a selection of the iteration values are given in the table below.

\(k\)	\(x_{1}^{(k)}\)	\(x_{2}^{(k)}\)	\(x_{3}^{(k)}\)	\(\\| \mathbf{r}^{(k)} \\|_\infty\)
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.2. Code#

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

Python

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}")

MATLAB

% 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

The Jacobi method

Contents

7.1. The Jacobi method#

7.1.1. The residual#

7.1.2. Code#