Solving a Linear System by Jacobi Iteration

Solve Ax = b without factoring the matrix by repeatedly updating each unknown from the previous sweep's values.

The splitting

Split A = D + R, where D is the diagonal and R the rest. Rearranging Ax = b gives x = D^-1 (b - R x). Jacobi turns this into an iteration: x_new[i] = (b[i] - sum of A[i,j] x_old[j] for j != i) / A[i,i]. Every new component uses only old values.

import numpy as np
A=np.array([[4.0,1.0,0.0],[1.0,5.0,1.0],[0.0,1.0,3.0]])
b=np.array([5.0,7.0,4.0])
x=np.zeros(3); D=np.diag(A)
for it in range(100):
    xnew=(b-(A@x-D*x))/D
    if np.max(np.abs(xnew-x))<1e-10: break
    x=xnew
print(np.round(x,5),'in',it,'iters')

Convergence condition

Jacobi converges if A is strictly diagonally dominant - each diagonal entry exceeds the sum of the absolute values of the others in its row. More generally it converges when the spectral radius of the iteration matrix D^-1 R is below 1. The smaller that radius, the faster it converges.

Why iterate instead of factor

Direct methods like LU cost order n^3 and fill in sparse matrices. Jacobi uses only matrix-vector products, preserves sparsity, and parallelizes trivially because every component updates independently from the same old vector. For the huge sparse systems from discretized PDEs, that scalability wins.

Relation to Gauss-Seidel

Using each freshly computed component immediately, rather than waiting for the next sweep, gives Gauss-Seidel, which typically converges about twice as fast but is harder to parallelize. Adding over-relaxation on top yields SOR, faster still when its parameter is tuned.