Jacobi method

In numerical linear algebra, the Jacobi 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 process is then iterated until it converges. This algorithm is a stripped-down version of the Jacobi transformation method of matrix diagonalization. The method is named after Carl Gustav Jacob Jacobi.
Description.
Let
be a square system of "n" linear equations, where:
formula_2
Then "A" can be decomposed into a diagonal component "D", a lower triangular part "L" and an upper triangular part "U":
The solution is then obtained iteratively via
where formula_5 is the "k"th approximation or iteration of formula_6 and formula_7 is the next or "k" + 1 iteration of formula_6. The element-based formula is thus:
The computation of formula_10 requires each element in x("k") except itself. Unlike the Gauss–Seidel method, we can't overwrite formula_11 with formula_10, as that value will be needed by the rest of the computation. The minimum amount of storage is two vectors of size "n".
Algorithm.
 Input: , (diagonal dominant) matrix formula_13, right-hand side vector formula_14, convergence criterion
 Output: 
 Comments: 
 
 while convergence not reached do
 for i := 1 step until n do
 for j := 1 step until n do
 if j ≠ i then
 end
 end
 end
 end
Convergence.
The standard convergence condition (for any iterative method) is when the spectral radius of the iteration matrix is less than 1:
A sufficient (but not necessary) condition for the method to converge is that the matrix "A" is strictly or irreducibly diagonally dominant. Strict row diagonal dominance means that for each row, the absolute value of the diagonal term is greater than the sum of absolute values of other terms:
The Jacobi method sometimes converges even if these conditions are not satisfied.
Note that the Jacobi method does not converge for every symmetric positive-definite matrix.
For example
Examples.
Example 1.
A linear system of the form formula_18 with initial estimate formula_19 is given by
We use the equation formula_21, described above, to estimate formula_22. First, we rewrite the equation in a more convenient form formula_23, where formula_24 and formula_25. From the known values
we determine formula_27 as
Further, formula_29 is found as
With formula_31 and formula_29 calculated, we estimate formula_22 as formula_34:
The next iteration yields
This process is repeated until convergence (i.e., until formula_37 is small). The solution after 25 iterations is
Example 2.
Suppose we are given the following linear system:
If we choose as the initial approximation, then the first approximate solution is given by
Using the approximations obtained, the iterative procedure is repeated until the desired accuracy has been reached. The following are the approximated solutions after five iterations.
The exact solution of the system is .
Python example.
import numpy as np
ITERATION_LIMIT = 1000
A = np.array(10., -1., 2., 0.],
 [-1., 11., -1., 3.],
 [2., -1., 10., -1.],
 [0.0, 3., -1., 8.)
b = np.array([6., 25., -11., 15.])
print("System:")
for i in range(A.shape[0]):
 row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
 print(" + ".join(row), "=", b[i])
print()
x = np.zeros_like(b)
for it_count in range(ITERATION_LIMIT):
 if it_count != 0:
 print("Iteration {0}: {1}".format(it_count, x))
 x_new = np.zeros_like(x)
 for i in range(A.shape[0]):
 s1 = np.dot(A[i, :i], x[:i])
 s2 = np.dot(A[i, i + 1:], x[i + 1:])
 x_new[i] = (b[i] - s1 - s2) / A[i, i]
 if x_new[i] == x_new[i-1]:
 break
 if np.allclose(x, x_new, atol=1e-10, rtol=0.):
 break
 x = x_new
print("Solution:")
print(x)
error = np.dot(A, x) - b
print("Error:")
print(error)
Weighted Jacobi method.
The weighted Jacobi iteration uses a parameter formula_41 to compute the iteration as
with formula_43 being the usual choice.
From the relation formula_44, this may also be expressed as
Convergence in the symmetric positive definite case.
In case that the system matrix formula_46 is of symmetric positive-definite type one can show convergence. 
Let formula_47 be the iteration matrix.
Then, convergence is guaranteed for
where formula_49 is the maximal eigenvalue.
The spectral radius can be minimized for a particular choice of formula_50 as follows
where formula_52 is the matrix condition number.