Jacobi & Gauss-Seidel iteration

These are the oldest iterative solvers, and the cleanest way to meet the idea. Take Ax = b and rewrite the i-th equation to isolate the i-th unknown: x_i equals (b_i minus the other terms) divided by A_ii. That formula begs to be iterated — plug in your current guess for the other unknowns to get a new guess for x_i, and repeat. Jacobi uses only old values within a sweep; Gauss-Seidel reuses each new value the moment it is computed.

Both fit a single template called matrix splitting. Write A = M - N where M is easy to invert. The iteration is M x_{k+1} = N x_k + b, equivalently x_{k+1} = M^-1 (N x_k + b). For Jacobi, M is the diagonal of A; for Gauss-Seidel, M is the lower-triangular part (diagonal included). The map's behavior is governed entirely by the iteration matrix G = M^-1 N.

Convergence has a clean criterion: the iteration converges for every starting guess if and only if the spectral radius rho(G) — the largest absolute eigenvalue of G — is strictly less than 1. The smaller rho(G), the faster: the error shrinks by roughly a factor of rho(G) per sweep. Gauss-Seidel typically beats Jacobi because using fresh values usually gives a smaller spectral radius, and unlike Jacobi it needs no second copy of the vector.

In modern practice these are rarely used as standalone solvers — their convergence is too slow on hard problems. But they live on as smoothers inside multigrid methods and as building blocks for preconditioners, where a few cheap sweeps damp the high-frequency error components beautifully.