10.4 Elliptic problems: Jacobi and Gauss–Seidel relaxation

Hyperbolic and parabolic PDEs (10.3) have a natural time-marching structure: discretise time, take one step at a time, and the algorithm is essentially a long sequence of small local updates. Elliptic PDEs have no time. There is nothing to march; the equation says the field’s value at every interior point is the average of its neighbours, and the boundary fixes everything. Discretising an elliptic problem produces a giant linear system, and solving it requires a fundamentally different algorithm.

This lesson covers the simplest such algorithm: Jacobi relaxation. The idea is brutally simple — keep replacing every interior cell by the average of its four neighbours — and the convergence is correspondingly slow. We then mention faster relaxation methods (Gauss–Seidel, SOR) and gesture at multigrid, which is the kind of algorithm that closes the practical performance gap with the time-marching schemes.

The discrete Laplacian

The second-order centred-difference approximation to $\partial_x^2$ is

\partial_x^2 u(x_j) \;\approx\; \frac{u_{j+1} - 2 u_j + u_{j-1}}{\Delta x^2}.

In 2-D, the same idea applied to both $x$ and $y$ gives the 5-point stencil for the Laplacian:

\nabla^2 u(x_i, y_j) \;\approx\; \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4 u_{i,j}}{\Delta x^2}.

(Assuming $\Delta x = \Delta y$ .) Setting this equal to zero for Laplace’s equation $\nabla^2 u = 0$ gives the discrete equation at every interior point:

\boxed{\;u_{i,j} \;=\; \tfrac{1}{4}\, \bigl(u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1}\bigr).\;}

The discrete value at any interior cell is the average of its four neighbours — the discrete mean-value property of harmonic functions, made explicit. On an $N \times N$ grid with $\sim N^2$ interior cells, this is $\sim N^2$ linear equations in $\sim N^2$ unknowns. A sparse linear system: each row has only 5 non-zero entries.

Jacobi iteration

The simplest way to solve this system: take the discrete-Laplace equation literally as an iteration rule. Start with any guess $u^{(0)}$ (boundary values pinned, interior arbitrary). At each step, replace every interior cell by the average of its current neighbours:

u_{i, j}^{(k+1)} \;=\; \tfrac{1}{4}\, \bigl(u_{i+1, j}^{(k)} + u_{i-1, j}^{(k)} + u_{i, j+1}^{(k)} + u_{i, j-1}^{(k)}\bigr).

Notice the superscript pattern: the new value at $(i, j)$ depends on the old values at the four neighbours. This is Jacobi iteration. It is trivially parallelisable: every interior cell can be updated independently using only the old grid.

preset:

Jacobi iteration replaces each interior cell with the average of its four neighbours, the discrete form of the mean-value property of harmonic functions. The interior starts at zero; with each pass, information about the boundary values diffuses one cell inward. After enough iterations, the field reaches the true Laplace solution. The residual decays exponentially — Jacobi is a *linear* iterative method with convergence rate ρ_J = cos(π / N) ≈ 1 − (π/N)²/2, so the residual takes roughly N²/π² iterations to fall by a factor of e. Faster relaxation methods (Gauss–Seidel, SOR, multigrid) accelerate this; Jacobi is the slow but simple workhorse.

The interactive shows the heatmap of $u$ on the left and the residual (the change per step) on the right, in log scale. Pick a boundary condition set and step the iteration. The residual decays exponentially — every step is multiplying it by a factor strictly less than 1. The interior temperature distribution gradually fills in from the boundary toward the centre.

Convergence rate

How fast does Jacobi converge? Slowly. The standard analysis: write $u^{(k)} = u^* + e^{(k)}$ where $u^*$ is the true solution. The error $e^{(k)}$ has the same Jacobi update rule (boundary values cancel), and decomposing it in the eigenmodes of the discrete Laplacian on an $N \times N$ grid produces:

e^{(k)} \;=\; \sum_{m, n} c_{mn}\, \rho_{mn}^k\, \phi_{mn},

where $\phi_{mn}$ is the $(m, n)$ -th eigenmode and $\rho_{mn}$ is the Jacobi spectral radius for that mode:

\rho_{m n} \;=\; \tfrac{1}{2}\, \bigl(\cos(m \pi / N) + \cos(n \pi / N)\bigr).

The slowest-decaying mode is $(m, n) = (1, 1)$ :

\rho_{1 1} \;=\; \cos(\pi / N) \;\approx\; 1 - \tfrac{\pi^2}{2 N^2}.

The number of iterations needed to reduce the residual by a factor $e$ is therefore $\sim 2 N^2 / \pi^2$ . The cost scales as $N^2$ iterations on an $N \times N$ grid. Combined with the $\mathcal{O}(N^2)$ cost per iteration, total cost is $\mathcal{O}(N^4)$ — terrible for fine grids.

Gauss–Seidel

A small modification gives a roughly $2\times$ speedup. Gauss–Seidel iteration updates each cell using the most recent available values for its neighbours, including those updated earlier in the same sweep:

u_{i, j}^{(k+1)} \;=\; \tfrac{1}{4}\, \bigl(u_{i-1, j}^{(k+1)} + u_{i+1, j}^{(k)} + u_{i, j-1}^{(k+1)} + u_{i, j+1}^{(k)}\bigr).

If you sweep left-to-right, top-to-bottom, the “above” and “left” neighbours have just been updated; the “below” and “right” neighbours are still from the previous iteration. The spectral radius is roughly $\rho_J^2$ instead of $\rho_J$ , so the iteration count drops by a factor of 2 — modest improvement but free.

Gauss–Seidel is not parallelisable in its natural form (it depends on the sweep order), unlike Jacobi. The standard workaround is red–black ordering: colour the grid like a checkerboard, update all red cells in parallel (each red cell depends only on black cells), then all black cells in parallel. Asymptotic convergence is the same as natural-order Gauss–Seidel.

Successive over-relaxation (SOR)

A more clever modification gives a much bigger speedup. SOR updates with an over-relaxation factor $\omega$ :

u_{i, j}^{(k+1)} \;=\; (1 - \omega)\, u_{i, j}^{(k)} \;+\; \omega \cdot \tfrac{1}{4}\, \bigl(\cdots\bigr),

where the parenthesised expression is the Gauss–Seidel update. Optimal $\omega \approx 2/(1 + \pi/N)$ approaches 2 as the grid grows; the spectral radius becomes $\omega - 1 \approx 1 - 2\pi/N$ , so the iteration count drops to $\mathcal{O}(N)$ instead of $\mathcal{O}(N^2)$ . Total cost becomes $\mathcal{O}(N^3)$ . Better than Jacobi by a factor of $N$ .

Multigrid

The asymptotically optimal method for elliptic problems is multigrid. The insight: Jacobi/Gauss–Seidel kill high-frequency error modes quickly (large $\rho$ close to zero) but low-frequency modes slowly ( $\rho$ close to 1). Multigrid alternates between fine and coarse grids — running a few smoothing iterations on the fine grid to kill high-frequency error, then projecting the residual onto a coarser grid (where the residual’s “low-frequency” content becomes “high-frequency” in the coarser representation) and smoothing there.

The result: each multigrid V-cycle reduces the residual by a constant factor (independent of $N$ ), and each cycle costs $\mathcal{O}(N^2)$ . Total cost is $\mathcal{O}(N^2)$ — optimal, because just reading the grid once costs $\mathcal{O}(N^2)$ .

We won’t develop multigrid here; the LaplaceSolver interactive uses plain Jacobi for simplicity. But multigrid is what you use when the problem is large.

Beyond Laplace

The same techniques extend to other linear elliptic equations:

Poisson’s equation $\nabla^2 u = \rho$ — same discretisation, the right-hand side becomes a $\rho$ contribution to each cell update.
Helmholtz equation $\nabla^2 u + k^2 u = 0$ — discrete version is non-symmetric in a particular sense and converges much more poorly. Iterative Krylov-subspace methods (GMRES, BiCGStab) are typically used for indefinite operators.
General sparse linear systems — most numerical libraries (PETSc, Trilinos) bundle a zoo of iterative solvers tuned for different elliptic operators.

What we use this for

LaplaceSolver (MATH-22) in Foundations 6.6 — Jacobi relaxation, exactly as described above.
Implicit time-stepping for the heat equation, when the FTCS stability constraint becomes painful, reduces to solving a sparse linear system at each step — usually by a Krylov method with multigrid preconditioning.

The next lesson, 10.5, takes us into a different kind of algorithm — the FFT — that turns out to be a workhorse for many numerical-PDE solvers too, via spectral methods that represent the field in a Fourier basis rather than on a grid.