Carl Jacobi
The simplest iterative method that was first applied for solving systems of linear equations is known as the Jacobi’s method,
Carl Gustav Jacob Jacobi (1804--1851) was a German mathematician who made fundamental contributions to elliptic functions, dynamics, differential equations, determinants, and number theory.
Although much of his
work was in applied mathematics,
Jacobi believed in the importance of
doing mathematics for its own sake.
A fine teacher, he held positions
at the Universities of Berlin and
Königsberg (now it is city of Kaliningrad) and was one of the most
famous mathematicians in Europe.
In numerical linear algebra, the Jacobi method (a.k.a. the
Jacobi iteration method) is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The method was discovered by Carl Gustav Jacob Jacobi in 1845.
The idea of this iteration method is to extract the main diagonal of the matrix of the given vector equation
\[
{\bf A}\,{\bf x} = {\bf b} ,
\]
and consider the rest of matrix A as perturbation:
\[
{\bf A} = \Lambda + {\bf M} .
\]
Here A is the given square n × n matrix, Λ is the diagonal matrix built from A (it is assumed that the diagonal entries of A do not vanish), and M is the matrix obtained from A by annihilating its diagonal entries. So
Since the new equation x = d + B x is of the form x = f(x), it is suitable for iteration procedure: x(k+1) = d + B x(k). Jacobi’s method uses diagonal part of matrix A for converting the given system of equation into fixed point form.
Here is the algorithm for Jacobi iteration method:
Input: initial guess x(0) to the solution, (diagonal dominant) matrix A, right-hand side
vector b, convergence criterion
Output: solution when convergence is reached
Comments: pseudocode based on the element-based formula above
k = 0
while convergence not reached do
for i := 1 step until n do
σ = 0
for j := 1 step until n do
if j ≠ i then
σ = σ + aij xj(k)
end
end
xi(k+1) = (bi − σ) / aii
end
increment k
end
We now need an initial approximation to the solution. It turns out that it does not
matter what this initial approximation is, so we might as well take x₁ = 0, x₂ = 0. We
use these values in Equations (1.1) to get new values of x₁ and x₂:
We repeat this process (using the old values of x₁
and x₂ to get the new values of x₁ and x₂), producing the sequence of approximations
given in the following Table
The successive vectors \( \displaystyle \quad \begin{bmatrix} x_1^{(k)} \\ x_2^{(k)} \end{bmatrix} \quad \) are called iterates. For instance, when k = 5, the fifth iterate is
\( \displaystyle \quad \begin{bmatrix} 1.8110 \\ 0.740849 \end{bmatrix} \quad \)
Here is Mathematica code to produce the result you see above, first using a module, which is used in other examples also.
In code above, parameter n_ is for the number of iterations you want - the display is always the last 5 iterations; fp_ is for floating point, the default being "True" - if you input "False" as the last argument, you will get an explicit solution, which produces fractions when input coefficients are from ℚ.
Clear[A, M, B, b, diag, d, \[CapitalLambda], x, l, u, x1, x2, x3];
Example 1.
Although exact calculations involve rational numbers, their usage in actual computations is entirely impractical because the number of bits required
to maintain exactness increases exponentially with the number of iterations. Therefore, next we demonstrate Jacobi's iterative method with floating point arithmetic.
Now we repeat calculations by making 50 iterative steps.
It took less than a second for Mathematica to execute the code with 50 steps. We do not present all output because it is expressed with huge numbers not suitable for human's eyes. For instance,
Expressing the rational numbers above as floating point reals, we get the same result, which is more suitable for human eyes.
So, we repeat calculations using floating point representation.
Here are the component parts of the code above, producing the same output in steps.
A = ({
{6, -5},
{4, -7}
});
d = DiagonalMatrix[Diagonal[A]];
l = LowerTriangularize[A, -1];
u = UpperTriangularize[A, 1];
b = ({
{7},
{1}
});
x[1] = ({
{0},
{0}Sdemo3
});
Grid[{{"A", "d", "l", "u", "b", "x0"},
MatrixForm[#] & /@ {A, d, l, u, b, x[1]}}, Frame -> All]
Starting point does not matter for final destination of Jacobi's ietrating scheme. Let's test that with calculations using different starting points, random floating point numbers and integers.
Below, a function is provided that works with a list of starting points.
The eigenvalues of matrix −Λ−1M are λ₁ ≅ 1.44914 and λ₂ ≅ −1.44914. Since the magnitudes of these eigenvalues exceed 1, we conclude that the Jacobi's iteration scheme diveges for this matrix.
Unlike direct methods that converge to the theoretical solution in a finite time, iterative methods provide approximations since they converge,
under some conditions, to the exact solution of the system of linear equations in infinite steps.
The Jacobi method to be applied, matrix A should
not display zeros on its main diagonal.
It is also convenient to organize the unknown values x1, x2, … , xn into a column vector, denoted by x.
To start an iteration method, you need an initial guess vector
Since convergence of Jacobi's scheme does not depend on your initial choice x(0), it could be any vector, for instance, zero.
At the (k + 1)st iteration x(k+1) = d + Bx(k), the vector x(k+1) is calculated by
The vector x(k) is used in calculating x(k+1). Consequently, these two vectors must be
stored separately.
If the diagonal elements of matrix A are much larger than the off-diagonal elements, the
entries of B should all be small and the Jacobi iteration should converge. We introduce a very important definition that reflects diagonally dominant property of a matrix.
A square matrix is said to be diagonally dominant if, for every row of the matrix, the magnitude of the diagonal entry in a row is larger than or equal to the sum of the magnitudes of all the other (non-diagonal) entries in that row. More precisely, the n×n matrix A = [𝑎ij] is diagonally dominant if
\[
\left\vert a_{ii} \right\vert \ge \sum_{j\ne i} \left\vert a_{ij} \right\vert \qquad \forall i.
\]
If A is diagonally dominant, the matrix B = −Λ−1M of the Jacobi iteration will have the property
It follows, then, that if A is diagonally dominant, the Jacobi iteration will converge to
the solution of A x = b.
Theorem 2:
If a system of n linear equations in n variables has a strictly diagonally dominant coefficient matrix, then it has a unique solution and both the Jacobi and the
Gauss--Seidel method converge to it.
Although we don't need to prove this theorem because it is a special (finite-dimensional) case of more general fixed point theorem, we present its proof..
For j = 1, 2, … , n, let
\[
\alpha_j = \sum_{i=1}^{j-1} \left\vert a_{j,i} \right\vert , \qquad \beta_j = \sum_{i=j+1}^{n} \left\vert a_{j,i} \right\vert \qquad\mbox{and} \qquad M_j = \frac{\beta_j}{| a_{j,i} | - \alpha_j} .
\]
Since A is diagonally dominant, it follows that
\[
\left\vert a_{j,j} \right\vert > \alpha_j + \beta_j , \qquad j=1,2,\ldots , n.
\]
Consequently, Mj < 1 for j = 1, 2, … , n. Thus,
\[
M = \max_{1 \le j \le n} M_j < 1.
\]
We will show that
\[
\| {\bf B} \|_{\infty} = \max_{{\bf x} \ne 0} \frac{\| {\bf B}\,{\bf x} \|_{\infty}}{\| {\bf x} \|_{\infty}} \le M < 1 .
\]
Let x be a nonzero vector in ℝn and let y =
B x. Choose k so that
\[
\| {\bf y} \|_{\infty} = \max_{1 \le i \le n} \, | y_i | = | y_k | .
\]
It follows from the definition of B that
\[
{\bf y} = {\bf B} \,{\bf x} = {\bf S}^{-1} {\bf U}\, {\bf x} ,
\]
and hence
\[
{\bf y} = {\bf B}\,{\bf x} = \Lambda^{-1} \left( {\bf L}\,{\bf y} + {\bf U}\, {\bf x} \right) ,
\]
where Λ is the diagonal matrix extracted from A and
L = S − Λ is strictly lower triangular matrix. Comparing the kth coordinates of each side, we see that
\[
y_k = -\frac{1}{a_{k,k}} \left( \sum_{i=1}^{k-1} a_{k,i} y_i + \sum_{i=k+1}^n a_{k,i} x_i \right) ,
\]
and hence
\[
\| {\bf y} \|_{\infty} = \left\vert y_k \right\vert \le \frac{1}{| a_{k,k} |} \left( \alpha_k \| {\bf y} \|_{\infty} + \beta_k \| {\bf x} \|_{\infty} \right) .
\]
From this inequality, it follows
\[
\frac{\| {\bf B}\,{\bf x} \|_{\infty}}{\| {\bf x} \|_{\infty}} =
\frac{\| {\bf y} \|_{\infty}}{\| {\bf x} \|_{\infty}} \le M_k \le M .
\]
So
\[
\| {\bf B} \|_{\infty} = \max_{{\bf x} \ne 0} \frac{\| {\bf B}\,{\bf x} \|_{\infty}}{\| {\bf x} \|_{\infty}} \le M < 1 ,
\]
and hence the iteration will converge to the solution of A x = b.
Example 3:
We consider system of linear equations
\[
\begin{split}
5\,x_1 - 4\, x_2 + x_3 &= 5, \\
5\,x_1 + 7\, x_2 -4\, x_3 &= 21, \\
3\, x_1 - 7\, x_2 + 8\, x_3 &= -9 ,
\end{split}
\tag{3.1}
\]
with matrix not being diagonally dominant:
\[
{\bf A} = \begin{bmatrix} 5 & -4 & \phantom{-}1 \\
5 & \phantom{-}7 & -4 \\ 3 & -7 & \phantom{-}8 \end{bmatrix} .
\]
Indeed,
\[
5 = |-4| +1, \qquad 7 < 5 + |-4| = 9, \qquad 8 < 3 + |-7| = 10 .
\]
We can check that the given system has the unique solution x₁ = 2, x₂ = 1, x₃ = −1.
A = {{5, -4, 1}, {5, 7, -4}, {3, -7, 8}};
b = {5, 21, -9};
LinearSolve[A, b]
{2, 1, -1}
However, the code above will not work if you define vector b as column vector. To overcome this, use the following modification.
b = ({
{5},
{21},
{-9}
});
LinearSolve[A, b][[All,1]]
Upon extracting the diagonal matrix from A, we split the given system (3.1) into the following form:
\[
\left( \Lambda + {\bf M} \right) {\bf x} = {\bf b} \qquad \iff \qquad {\bf x} =
\Lambda^{-1} {\bf b} - \Lambda^{-1} {\bf M} \,{\bf x} ,
\]
where
\[
{\bf d} = \Lambda^{-1} {\bf b} = \begin{bmatrix} \frac{1}{5} & 0 & 0 \\ 0 & \frac{1}{7} & 0 \\ 0 & 0 & \frac{1}{8} \end{bmatrix} \begin{pmatrix} \phantom{-}5 \\ 21 \\ -9 \end{pmatrix} = \begin{bmatrix} 1 \\ 3 \\ - \frac{9}{8} \end{bmatrix} , \qquad {\bf B} = - \Lambda^{-1} {\bf M} = \begin{bmatrix} 0 & \frac{4}{5} & - \frac{1}{5} \\ - \frac{5}{7} & 0 & \frac{4}{7} \\ - \frac{3}{8} & \frac{7}{8} & 0 \end{bmatrix} .
\]
Matrix B has one real eigenvalue λ₁ ≈ -0.362726,
and two complex conjugate eigenvalues λ₂ ≈ 0.181363 +
0.308393j, λ₃ ≈ 0.181363 - 0.308393j, where j is the imaginary unit of the complex plane ℂ, so j² = −1.
Since all eigenvalues have absolute value less than 1, we expect Jacobi's method to be convergent,
Sqrt[0.181363^2 + 0.308393^2]
0.357769
The code above shows that the absolute value of complex number λ = x + jy, with x = 0.181363 and y = 0.308393, is
\( \displaystyle |\lambda | = | x + {\bf j}\,y | = \sqrt{x^2 + y^2} \approx 0.357769 < 1 . \)
Although we know exact solution, x₁ = 2, x₂ = 1,
x₃ = −1, of the given system of linear equations (3.1), we know that convergence of Jacobi's iteration scheme does not depend on the initial guess. So we start with zero:
\[
{\bf x}^{(0)} = \left[ 0, 0, 0 \right] .
\]
d = {1,3,-9/8};
B= {{0, 4/5, -1/5},{-5/7,0,4/7},{-3/8,7/8,0}};
x3 = d + B.x2
{117/56, 59/56, -(67/64)}
N[%]
{2.08929, 1.05357, -1.04688}
\[
{\bf x}^{(3)} = {\bf d} + {\bf B}\, {\bf x}^{(2)} = \begin{bmatrix} \frac{117}{56} \\ \frac{59}{56} \\ -\frac{67}{64} \end{bmatrix} =
\begin{bmatrix} 2.08929 \\ 1.05357 \\ -1.04688 \end{bmatrix} .
\]
Now it is right time to organize our calculations in a loop.
Theorem 3:
If the Jacobi or the Gauss-Seidel method converges for a system of n linear
equations in n unknowns, then it must converge to the solution of the system.
The idea of the proof is quite simple: to apply an iterative method, we first transferred the linear algebra problem A x = b to a fixed point problem x = f(x), and only then applied the iterative method to the latter. To prove Theorem 3, we go backward: from the fixed point equation we go to the linear algebra one.
Convergence means that "as iterations increase, the values of the iterates get
closer and closer to a limiting value:' This means that x(k) converge to r, As usual, we denote by x = c the true solution of the system A x = b.
We must show that r = c. Recall that the Jacobi's (GS is similar) iterates satisfy the equation
\[
\Lambda\,{\bf x}^{(k+1)} = {\bf b} + \left( \Lambda - {\bf A} \right) {\bf x}^{(k)} , \qquad k=0,,2,\ldots .
\tag{P3.1}
\]
The Gauss-Seidel method leads to a similar equation where A is further decomposed into upper and lower triangular parts. Application of limit to both sides of Eq.(P3.1), we obtain
\[
\Lambda\,{\bf r} = {\bf b} + \left( \Lambda - {\bf A} \right) {\bf r} \qquad \iff \qquad {\bf 0} = {\bf b} - {\bf A}\,{\bf r} .
\]
Therefore, vector \( \displaystyle {\bf r} = \lim_{k\to \infty} {\bf x}^{(k)} \) satisfies the equation A x = b, so it must be equal to c.
Example 4:
Let us reconsider the system four equations with four unknowns:
Although we know one of solutions, x₁ = 2, x₂ = 1, x₃ = 3, x₄ = 4, we demonstrate the iterative method---Jacobi's one. Since the nullspace of A is spanned on the vector [12, 31, 124, 120], the general solution of Eq.(4.1) depends on real parameter t ∈ ℝ:
\[
{\bf x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} 2 \\ 1 \\ 3 \\ 4 \end{bmatrix} + t \begin{bmatrix} 12 \\ 31 \\ 124 \\ 120 \end{bmatrix} \in \mathbb{R}^{4\times 1} .
\]
The Jacobi scheme generates approximants according to the formula
\[
{\bf x}^{(k+1)} = {\bf d} + {\bf B}\,{\bf x}^{(k)} , \qquad k=0, 1, 2, \ldots .
\tag{4.3}
\]
Eigenvalues of matrix B can be determined with Mathematica
N[Eigenvalues[B]]
{1., -0.576902 + 0.664232 I, -0.576902 - 0.664232 I, 0.153804}
Since one of the eigenvalues λ₁ = 1, we do not expect convergence of Jacobi's iteration procedure. Other eigenvalues λ₂ = 0.153804, λ₃ = -0.576902 + 0.664232j, λ₄ = -0.576902 − 0.664232j, have magnitudes less than 1.
Sqrt[0.5769020858621225^2 + 0.6642324328522395^2]
0.879784
Starting with zero vector x(0) = 0, we obtain the first approximant to be x(1) = d. The second and third approximations become
\[
{\bf x}^{(2)} = {\bf d} + {\bf B}\,{\bf x}^{(1)} = {\bf d} + {\bf B}\,{\bf d} = \begin{pmatrix} \frac{209}{63} \\ -\frac{89}{112} \\ \frac{1}{2} \\ \frac{37}{84} \end{pmatrix}
\]
and
\[
{\bf x}^{(3)} = {\bf d} + {\bf B}\,{\bf x}^{(2)} = {\bf d} + {\bf B} \left( {\bf d} + {\bf B}\,{\bf b} \right) = {\bf b} + {\bf B}\,{\bf b} + {\bf B}^2 {\bf b} = \begin{pmatrix} \frac{40}{189} \\ \frac{175}{289} \\ -\frac{1235}{1764} \\ \frac{425}{1512} \end{pmatrix} .
\]
diag = {{9, 0, 0, 0}, {0, 8, 0, 0}, {0, 0, 7, 0}, {0, 0, 0, 6}};
M = A - diag ;
d = {1, 7/8, -(11/7), 3/2};
B = -Inverse[diag] . M ;
x1 = d + B . d
{209/63, -(89/112), 1/2, 37/84}
x2 = d + B . d + B . B . d
{40/189, 175/288, -(1235/1764), 425/1512}
Now it is right time to organize our calculations in a loop.
The numerical values in the table above confirm that Jacobi's iteration diverges for this system of equations.
We are going to make a numerical experiment and check whether Jacobi's iteration depends on the initial guess, which we take closer to the true values:
x(0) = [1.5, 0.5, 2.5, 3.5]. Using iteration (4.3), we calculate Jacobi's approximants.
This table of iterative values shows that the first coordinate plays around true value x₁ = 2, but we cannot claim convergence of other coordinates, except probably x₂ = 1.
Input: square matrix A and vector bn is the size of A and vector b
whle precision condition do
for i=1 : n do
s=0
for j=1 : n do
if j ≠ i then
s = s + 𝑎i,jxj
endif
endfor
yi = (bi−s)/𝑎i,i
endfor
x=y
endwhile
Anton, Howard (2005), Elementary
Linear Algebra (Applications Version) (9th ed.), Wiley International
