Jacobi and Gauss-Seidel

These classic stationary iterations split the matrix and sweep repeatedly; Gauss-Seidel reuses updated values and usually converges faster than Jacobi.

Splitting the matrix

Stationary iterative methods solve A x = b by splitting A into parts and iterating. Jacobi's method solves each equation for its diagonal unknown using the previous iterate's values for the rest, updating all components from old data. It is naturally parallel because every update is independent.

Gauss-Seidel makes one change: it uses the newest available values as it sweeps through the unknowns, so updates within a sweep see each other. This coupling usually roughly doubles the convergence rate compared with Jacobi, at the cost of the sweep being inherently sequential.

import numpy as np
def gauss_seidel(A, b, x, iters):
    n = len(b)
    for _ in range(iters):
        for i in range(n):
            s = A[i,:] @ x - A[i,i]*x[i]
            x[i] = (b[i] - s) / A[i,i]
    return x

When they converge

Both methods converge if A is strictly diagonally dominant or symmetric positive definite. Convergence is guaranteed when the spectral radius of the iteration matrix is below 1, and its size sets the rate. For many discretized PDEs the rate is close to 1, meaning slow convergence that worsens as the grid refines.

Role today

On their own, Jacobi and Gauss-Seidel are too slow for large modern problems. Their lasting value is as smoothers inside multigrid methods, where they efficiently damp high-frequency error, and as simple preconditioners for Krylov solvers. Understanding them is the gateway to those more powerful methods.

As smoothers within multigrid, these iterations help solve the large sparse elliptic systems that arise in breeder Hyperion field computations.