Matrix-based formula

Then A can be decomposed into a diagonal component D, a lower triangular part L and an upper triangular part U:

$${\displaystyle A=D+L+U\qquad {\text{where}}\qquad D={\begin{bmatrix}a_{11}&0&\cdots &0\\0&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &a_{nn}\end{bmatrix}}{\text{ and }}L+U={\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\a_{21}&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &0\end{bmatrix}}.}$$

The solution is then obtained iteratively via

${\displaystyle \mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -(L+U)\mathbf {x} ^{(k)}).}$

Element-based formula

The element-based formula for each row ${\displaystyle i}$ is thus:

$${\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j\neq i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.}$$

The computation of ${\displaystyle x_{i}^{(k+1)}}$ requires each element in ${\displaystyle \mathbf {x} ^{(k)}}$ except itself. Unlike the Gauss–Seidel method, we can't overwrite ${\displaystyle x_{i}^{(k)}}$ with ${\displaystyle x_{i}^{(k+1)}}$, as that value will be needed by the rest of the computation. The minimum amount of storage is two vectors of size n.

Example 1

A linear system of the form ${\displaystyle Ax=b}$ with initial estimate ${\displaystyle x^{(0)}}$ is given by

${\displaystyle A={\begin{bmatrix}2&1\\5&7\\\end{bmatrix}},\ b={\begin{bmatrix}11\\13\\\end{bmatrix}}\quad {\text{and}}\quad x^{(0)}={\begin{bmatrix}1\\1\\\end{bmatrix}}.}$

We use the equation ${\displaystyle x^{(k+1)}=D^{-1}(b-(L+U)x^{(k)})}$, described above, to estimate ${\displaystyle x}$. First, we rewrite the equation in a more convenient form ${\displaystyle D^{-1}(b-(L+U)x^{(k)})=Tx^{(k)}+C}$, where ${\displaystyle T=-D^{-1}(L+U)}$ and ${\displaystyle C=D^{-1}b}$. From the known values

$${\displaystyle D^{-1}={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}},\ L={\begin{bmatrix}0&0\\5&0\\\end{bmatrix}}\quad {\text{and}}\quad U={\begin{bmatrix}0&1\\0&0\\\end{bmatrix}}.}$$

we determine ${\displaystyle T=-D^{-1}(L+U)}$ as

$${\displaystyle T={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}\left\{{\begin{bmatrix}0&0\\-5&0\\\end{bmatrix}}+{\begin{bmatrix}0&-1\\0&0\\\end{bmatrix}}\right\}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}.}$$

Further, ${\displaystyle C}$ is found as

$${\displaystyle C={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}.}$$

With ${\displaystyle T}$ and ${\displaystyle C}$ calculated, we estimate ${\displaystyle x}$ as ${\displaystyle x^{(1)}=Tx^{(0)}+C}$:

$${\displaystyle x^{(1)}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}{\begin{bmatrix}1\\1\\\end{bmatrix}}+{\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}={\begin{bmatrix}5.0\\8/7\\\end{bmatrix}}\approx {\begin{bmatrix}5\\1.143\\\end{bmatrix}}.}$$

The next iteration yields

$${\displaystyle x^{(2)}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}{\begin{bmatrix}5.0\\8/7\\\end{bmatrix}}+{\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}={\begin{bmatrix}69/14\\-12/7\\\end{bmatrix}}\approx {\begin{bmatrix}4.929\\-1.714\\\end{bmatrix}}.}$$

This process is repeated until convergence (i.e., until ${\displaystyle \|Ax^{(n)}-b\|}$ is small). The solution after 25 iterations is

${\displaystyle x={\begin{bmatrix}7.111\\-3.222\end{bmatrix}}.}$

Example 2

Suppose we are given the following linear system:

${\displaystyle {\begin{aligned}10x_{1}-x_{2}+2x_{3}&=6,\\-x_{1}+11x_{2}-x_{3}+3x_{4}&=25,\\2x_{1}-x_{2}+10x_{3}-x_{4}&=-11,\\3x_{2}-x_{3}+8x_{4}&=15.\end{aligned}}}$

If we choose (0, 0, 0, 0) as the initial approximation, then the first approximate solution is given by

$${\displaystyle {\begin{aligned}x_{1}&=(6+0-(2*0))/10=0.6,\\x_{2}&=(25+0+0-(3*0))/11=25/11=2.2727,\\x_{3}&=(-11-(2*0)+0+0)/10=-1.1,\\x_{4}&=(15-(3*0)+0)/8=1.875.\end{aligned}}}$$

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.

${\displaystyle x_{1}}$	${\displaystyle x_{2}}$	${\displaystyle x_{3}}$	${\displaystyle x_{4}}$
0.6	2.27272	-1.1	1.875
1.04727	1.7159	-0.80522	0.88522
0.93263	2.05330	-1.0493	1.13088
1.01519	1.95369	-0.9681	0.97384
0.98899	2.0114	-1.0102	1.02135

The exact solution of the system is (1, 2, −1, 1).

Python example

import numpy as np

ITERATION_LIMIT = 1000

# initialize the matrix
A = np.array([[10., -1., 2., 0.],
              [-1., 11., -1., 3.],
              [2., -1., 10., -1.],
              [0.0, 3., -1., 8.]])
# initialize the RHS vector
b = np.array([6., 25., -11., 15.])

# prints the system
print("System:")
for i in range(A.shape[0]):
    row = [f"{A[i, j]}*x{j + 1}" for j in range(A.shape[1])]
    print(f'{" + ".join(row)} = {b[i]}')
print()

x = np.zeros_like(b)
for it_count in range(ITERATION_LIMIT):
    if it_count != 0:
        print(f"Iteration {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)

Convergence in the symmetric positive definite case

In case that the system matrix ${\displaystyle A}$ is of symmetric positive-definite type one can show convergence.

Let ${\displaystyle C=C_{\omega }=I-\omega D^{-1}A}$ be the iteration matrix. Then, convergence is guaranteed for

${\displaystyle \rho (C_{\omega })<1\quad \Longleftrightarrow \quad 0<\omega <{\frac {2}{\lambda _{\text{max}}(D^{-1}A)}}\,,}$

where ${\displaystyle \lambda _{\text{max}}}$ is the maximal eigenvalue.

The spectral radius can be minimized for a particular choice of ${\displaystyle \omega =\omega _{\text{opt}}}$ as follows

$${\displaystyle \min _{\omega }\rho (C_{\omega })=\rho (C_{\omega _{\text{opt}}})=1-{\frac {2}{\kappa (D^{-1}A)+1}}\quad {\text{for}}\quad \omega _{\text{opt}}:={\frac {2}{\lambda _{\text{min}}(D^{-1}A)+\lambda _{\text{max}}(D^{-1}A)}}\,,}$$

where ${\displaystyle \kappa }$ is the matrix condition number.