Numerical Methods: Relaxation Techniques

The Big Picture: Iterative Averaging

Poisson's equation, $\nabla^2 V = -\rho/\varepsilon_0$, has a beautifully simple meaning: at every point, the value of $V$ differs from the average of its surroundings by exactly the amount of source sitting there. No source means "equal to the average"; positive source means "a bit above the average"; negative source means "a bit below." That single sentence is the entire physical content of the equation — and it is also, almost word-for-word, an algorithm.

Imagine a heated metal plate whose edges are held at fixed temperatures and whose interior has a few heat sources sprinkled around. To find the steady-state temperature distribution, you could solve a huge sparse linear system. Or you could do something far simpler: walk around the plate and adjust each interior point so it matches the average of its neighbours, plus a little kick from any local heat source. Repeat. The temperature field will relax into the steady state — hence the name relaxation methods.

The core idea in one line. A relaxation method takes the discretized Poisson equation, $V_{i,j} = \tfrac{1}{4}\big( V_{i-1,j} + V_{i+1,j} + V_{i,j-1} + V_{i,j+1} + h^2 f_{i,j} \big)$, and applies it as an assignment rule, sweep after sweep, until the values stop changing. The PDE's defining property — every cell averages its neighbours — becomes the literal source code.

In this section we will walk through the three classical relaxation methods in order of cleverness:

1. Jacobi — sweep the grid using only the old neighbour values. Maximum parallelism, slowest convergence.
2. Gauss-Seidel — sweep the grid using the freshest values available. Roughly twice as fast, harder to parallelize.
3. Successive Over-Relaxation (SOR) — do Gauss-Seidel, then deliberately overshoot by a factor $\omega \in (1, 2)$. With a well-chosen $\omega$, an order of magnitude faster than either of the previous two.

Then we will sketch how multigrid attacks the same problem from a completely different angle, achieving the asymptotically optimal $O(N^2)$ cost in 2D — meaning constant work per unknown, regardless of grid size.

From PDE to Grid: The 5-Point Stencil

Before we can iterate, we need a finite version of $\nabla^2 V$. Replace the continuous square $[0, L]^2$ with a uniform grid of $N \times N$ points and spacing $h = L/(N-1)$. Write $V_{i,j} = V(i h, j h)$. Now Taylor-expand:

Add them. Odd-order terms cancel:

Solving for $V_{xx}$ gives the famous centred second difference:

Do the same on the $y$-axis and add. The discrete Laplacian is:

This is the 5-point stencil. Set it equal to $-f_{i,j}$ (where $f = \rho/\varepsilon_0$), and solve for $V_{i,j}$:

This is the master formula every relaxation method below applies, over and over. The only thing the three methods disagree about is which values of the neighbours to use when computing the average.

Picture the stencil. The interior cell $(i,j)$ reaches out to its four nearest neighbours — up, down, left, right — and to a single source term. No diagonals. The pattern looks like a plus-sign:

           V_{i, j+1}
                |
V_{i-1, j} --- V_{i, j} --- V_{i+1, j}     +   (h² / 4) · f_{i,j}
                |
           V_{i, j-1}

Jacobi: The Cleanest Idea

Carl Gustav Jacob Jacobi's relaxation idea (1845) is the most literal possible reading of the discretized equation. Make two copies of the grid. In one sweep, read from the old copy, write to the new copy:

Then swap the copies and repeat. The superscript $(k)$ means "value at the start of iteration $k$." The defining property: every cell on iteration $k+1$ uses only iteration-$k$ values. The order in which we update cells does not matter, so we can do them all simultaneously — perfect for vectors and GPUs.

Why this works (intuition). Picture a hot point in the middle of a cold plate. After the first Jacobi sweep, only the four cells immediately touching the hot point feel anything; everyone else is still cold. After the second sweep, the next ring out feels it. Heat — or potential, or whatever $V$ represents — diffuses outward by exactly one cell per sweep. To reach a boundary that is $N$ cells away takes at least $N$ sweeps. That is the geometric reason Jacobi is slow.

Formally, Jacobi's iteration matrix has spectral radius $\rho_J = \cos(\pi h)$ for the discretized Dirichlet Poisson problem on a unit square. For small $h$ that is $1 - \tfrac{1}{2}\pi^2 h^2 + \mathcal{O}(h^4)$, which means each iteration cuts the error by a factor that is very close to 1. Halving $h$ roughly quadruples the number of iterations needed. Jacobi is $O(N^2)$ iterations on an $N \times N$ grid.

Gauss-Seidel: Use the Fresh News

Imagine the same sweep, but as you walk from cell to cell, you immediately write the new value back. By the time you reach cell $(i, j)$, the neighbours $(i-1, j)$ and $(i, j-1)$ (assuming row-major sweep order) have already been updated this iteration. So when you compute the new $V_{i,j}$, you use:

Mixed superscripts — two updated, two not yet — because cells are visited in lexicographic order. This is the Gauss-Seidel method (Gauss 1823 in letters, Seidel 1874 in print).

Intuition. Why should using fresher data be faster? Picture the hot-point example again. In Jacobi, a single update at the hot spot only sends information to its immediate neighbours by the end of the sweep. In Gauss-Seidel, once the hot spot updates, every cell scanned after it in the same sweep already sees the change. Information propagates further per iteration. In fact, the spectral radius drops to $\rho_{GS} = \cos^2(\pi h) = \rho_J^2$, so Gauss-Seidel takes roughly half the iterations of Jacobi to reach the same accuracy.

The downside: the update order matters. Cell $(i, j)$ depends on the freshly-updated cell $(i-1, j)$, which we cannot compute in parallel with $(i, j)$. On a GPU this is a disaster — we would lose the very parallelism that made Jacobi attractive. The fix is a red-black ordering (chessboard colouring), which we will use in the PyTorch implementation below: update all red cells simultaneously, then all black cells, alternating sweeps. With this trick Gauss-Seidel and SOR vectorize cleanly.

SOR: Overshoot on Purpose

Here is a strange-sounding observation: when Gauss-Seidel updates a cell, the new value almost always lies between the old value and the eventual converged value. The iteration is, on average, under-correcting. What if we deliberately jumped a bit past the Gauss-Seidel update — would we get there faster? That gamble is Successive Over-Relaxation, due to David M. Young Jr. in his 1950 thesis.

Let $V^{\text{GS}}_{i,j}$ denote the value that Gauss-Seidel would assign. SOR replaces it with:

The factor $\omega$ is the relaxation parameter. Three regimes:

The miracle is that for the discretized Poisson problem there is an exact formula for the optimal value:

For $N = 33$ this is $\omega \approx 1.826$; for $N = 257$ it is $\omega \approx 1.976$. As the grid gets finer the optimal omega creeps toward 2 (but never reaches it).

The payoff. With optimal omega, the spectral radius becomes $\rho_{SOR} = \omega_{\text{opt}} - 1$, which behaves like $1 - 2\pi h$ for small $h$. Compare to Jacobi's $1 - \tfrac{1}{2}\pi^2 h^2$. Halving $h$ only doubles the iteration count for SOR, while Jacobi's quadruples. That is the difference between $O(N)$ and $O(N^2)$ iterations — a giant speedup for the cost of two extra multiplies per cell.

Why SOR Wins: Spectral Radius Intuition

Each relaxation method can be written as a linear iteration $V^{(k+1)} = G\, V^{(k)} + c$ for some iteration matrix $G$. The error $e^{(k)} = V^{(k)} - V^{*}$ evolves as $e^{(k+1)} = G\, e^{(k)}$, so after $k$ iterations the error is reduced by a factor of approximately $\rho(G)^k$, where $\rho(G)$ is the spectral radius (largest eigenvalue magnitude) of $G$.

For the discretized Poisson problem on an $N \times N$ unit square the three numbers are known in closed form:

For an example, set $N = 100$:

• Jacobi: ≈ 1400 iterations.
• Gauss-Seidel: ≈ 700 iterations.
• SOR(optimal): ≈ 70 iterations.

SOR is a full 20× faster than Jacobi at $N = 100$, and the gap grows linearly with $N$. The interactive plot below sweeps $\omega$ on a 32×32 grid and shows the characteristic basin — under-relaxation slows GS down, optimal omega crashes the residual, and over-shooting past 2 makes the iteration explode:

Fourier intuition for what relaxation does. Any error on the grid can be decomposed into spatial frequencies (modes). Jacobi and GS damp high-frequency modes (the wiggly errors) very quickly, but low-frequency modes (the long-wavelength ones) decay agonisingly slowly — they need information to propagate across the whole grid. SOR's overshoot is engineered to cancel the slowest mode at every step, dramatically accelerating the low-frequency damping. Multigrid (see below) attacks the same problem in a different way: solve the low-frequency part on a coarser grid where it is no longer low-frequency.

Interactive: One-Step Relaxation Animator

Click Step ▶ repeatedly and watch a 7-cell 1D grid relax. The single red bar is a source at $i = 3$; the grey bars at the ends are the Dirichlet boundaries $V_0 = V_6 = 0$. The dashed green outlines show the exact converged solution. Switch methods between Jacobi, Gauss-Seidel, and SOR; with SOR enabled, slide $\omega$ and see how aggressively the cells overshoot.

Things to notice as you play:

1. Jacobi. Information spreads one cell per step. After one step, only the centre cell is non-zero. After two, the immediate neighbours light up. After three, the next ring out. Each cell takes many sweeps to reach its true value.
2. Gauss-Seidel. After a single sweep the cells right-of-centre already feel the source (because we sweep left-to-right and the centre updates first). After two or three sweeps the field is already close to the truth.
3. SOR. At $\omega = 1.5$ you will see cells overshoot the truth, then settle. At $\omega \approx 1.7$ the convergence is near-optimal. Push $\omega$ above 1.9 and the iteration may oscillate or refuse to settle.

Worked Example: 1D Poisson by Hand

Let's do the same problem on paper to feel the mechanics. Consider the 1D Poisson equation $-V''(x) = f(x)$ on $[0, 5]$ with $V(0) = V(5) = 0$, a single source $f_3 = 8$ at $x = 3$ and $h = 1$. The discrete equation is:

The unknowns are $V_1, V_2, V_3, V_4$. The exact converged solution comes from a 4×4 linear system — derived inside the collapsible below — and equals $(V_1, V_2, V_3, V_4) = (3.2,\, 6.4,\, 9.6,\, 4.8)$. Below we run Jacobi, Gauss-Seidel, and SOR for a few sweeps against this target and watch how quickly each one closes the gap.

V_1 = ½ (V_0 + V_2 + h² f_1) = ½ (0 + V_2 + 0)      = ½ V_2
V_2 = ½ (V_1 + V_3 + 0)                              = ½ (V_1 + V_3)
V_3 = ½ (V_2 + V_4 + h² f_3) = ½ (V_2 + V_4 + 8)
V_4 = ½ (V_3 + V_5 + 0)      = ½ (V_3 + 0)           = ½ V_3

V_2 = ½ (V_2/2 + V_3)   →   2 V_2 = V_2/2 + V_3   →   V_3 = (3/2) V_2

(3/2) V_2 = ½ (V_2 + (3/2) V_2 · ½ + 8)
(3/2) V_2 = ½ (V_2 + (3/4) V_2 + 8)
(3/2) V_2 = (V_2 + (3/4) V_2 + 8) / 2
3 V_2     =  V_2 + (3/4) V_2 + 8
3 V_2 - V_2 - (3/4) V_2 = 8
(5/4) V_2 = 8
V_2 = 32 / 5 = 6.4

V_1 = V_2 / 2     = 3.2
V_2               = 6.4
V_3 = (3/2) V_2   = 9.6
V_4 = V_3 / 2     = 4.8

Iter 1:
  V_1^GS = ½(0 + 0 + 0) = 0     →  V_1 = -0.5·0 + 1.5·0    = 0.000
  V_2^GS = ½(0 + 0 + 0) = 0     →  V_2 = -0.5·0 + 1.5·0    = 0.000
  V_3^GS = ½(0 + 0 + 8) = 4     →  V_3 = -0.5·0 + 1.5·4    = 6.000
  V_4^GS = ½(6 + 0 + 0) = 3     →  V_4 = -0.5·0 + 1.5·3    = 4.500

Iter 2:
  V_1^GS = ½(0 + 0 + 0) = 0     →  V_1 = -0.5·0 + 1.5·0    = 0.000
  V_2^GS = ½(0 + 6 + 0) = 3     →  V_2 = -0.5·0 + 1.5·3    = 4.500
  V_3^GS = ½(4.5 + 4.5 + 8) = 8.5 → V_3 = -0.5·6 + 1.5·8.5 = 9.750
  V_4^GS = ½(9.75 + 0 + 0) = 4.88 → V_4 = -0.5·4.5 + 1.5·4.88 = 5.06

Interactive: Jacobi vs Gauss-Seidel vs SOR

Now back to the 2D Poisson problem on an $N \times N$ grid. All three methods start from $V = 0$ and try to drive the residual $\|\nabla^2 V + f\|$ to zero. The curves below show that residual on a log scale per iteration. The thumbnails on the bottom show the final field each method produces.

Slide the grid resolution and watch what happens to the gap between the three curves. At $N = 16$ the gap is already noticeable; by $N = 64$ Jacobi has barely made progress in the same 400 iterations while SOR has bottomed out near machine precision. Try setting $\omega = 1.0$ — the SOR curve collapses onto Gauss-Seidel (as expected). Push $\omega$ above 1.95 and the SOR curve may bend back upward as you approach divergence.

Python: All Three Solvers in One File

Here is a complete, runnable Python implementation of Jacobi, Gauss-Seidel, and SOR. The test uses the manufactured-solution method: choose $V_{\text{exact}}(x, y) = \sin(\pi x) \sin(\pi y)$, compute $-\nabla^2 V_{\text{exact}} = 2\pi^2 \sin(\pi x) \sin(\pi y)$, feed that as the source, and check the solver recovers $\sin(\pi x) \sin(\pi y)$ back. With this gold-standard test you cannot fool yourself: the iteration count gap and the final-error gap separately tell you about solver speed and about discretization quality.

np.zeros((32, 32))  →  shape (32, 32) of float64 zeros

np.linspace(0,1,5) = array([0.0, 0.25, 0.5, 0.75, 1.0])

x=[0,1,2]
y=[10,20]
X = [[0,0],[1,1],[2,2]]
Y = [[10,20],[10,20],[10,20]]

f[16, 16] = 2π² · sin(π·0.5)² ≈ 2π² ≈ 19.74

f.shape = (33, 33)  →  N = 33

1.0 / (33 - 1) = 0.03125

V = np.zeros((33, 33))

id(V_new) != id(V)  →  True (separate memory)

V[1:-1, 1:-1] shape (31, 31) — the interior of a 33×33 grid

After (i=10, j=10) is updated, the very next call at (i=11, j=10) reads the FRESH V[10, 10].

If V_old = 5.0 and GS would give 6.0, then SOR(ω=1.5) gives:
  V_new = (1 - 1.5)·5 + 1.5·6 = -2.5 + 9 = 6.5
→ we move 50% further than GS would.

N=33 → ω_opt = 2 / (1 + sin(π/33)) = 2 / (1 + 0.0951) ≈ 1.8262

f[16, 16] = 2π²·sin(π·0.5)·sin(π·0.5) ≈ 19.74

len(hj) ≈ 1800

len(hg) ≈ 900

len(hs) ≈ 120

max|V_jacobi - V_exact| ≈ 2.5e-3   (this is the discretization error, not solver error)

1import numpy as np
2
3def discretize_poisson(N, source_fn):
4    """Build the right-hand side f_{i,j} = rho_{i,j}/eps0 for a unit square."""
5    x = np.linspace(0, 1, N)
6    y = np.linspace(0, 1, N)
7    X, Y = np.meshgrid(x, y, indexing="ij")
8    return source_fn(X, Y)
9
10def jacobi(f, max_iter=4000, tol=1e-8):
11    """Jacobi relaxation: every new value uses OLD neighbours."""
12    N = f.shape[0]
13    h = 1.0 / (N - 1)
14    V = np.zeros((N, N))
15    history = []
16    for k in range(max_iter):
17        V_new = V.copy()
18        V_new[1:-1, 1:-1] = 0.25 * (
19            V[:-2, 1:-1] + V[2:, 1:-1]
20          + V[1:-1, :-2] + V[1:-1, 2:]
21          + h * h * f[1:-1, 1:-1]
22        )
23        diff = np.max(np.abs(V_new - V))
24        V = V_new
25        history.append(diff)
26        if diff < tol:
27            break
28    return V, history
29
30def gauss_seidel(f, max_iter=4000, tol=1e-8):
31    """Gauss-Seidel: update IN PLACE so the right side already sees fresh values."""
32    N = f.shape[0]
33    h = 1.0 / (N - 1)
34    V = np.zeros((N, N))
35    history = []
36    for k in range(max_iter):
37        diff = 0.0
38        for j in range(1, N - 1):
39            for i in range(1, N - 1):
40                new = 0.25 * (V[i-1, j] + V[i+1, j]
41                            + V[i, j-1] + V[i, j+1]
42                            + h * h * f[i, j])
43                diff = max(diff, abs(new - V[i, j]))
44                V[i, j] = new
45        history.append(diff)
46        if diff < tol:
47            break
48    return V, history
49
50def sor(f, omega=1.85, max_iter=4000, tol=1e-8):
51    """Successive Over-Relaxation: GS + a deliberate overshoot factor omega."""
52    N = f.shape[0]
53    h = 1.0 / (N - 1)
54    V = np.zeros((N, N))
55    history = []
56    for k in range(max_iter):
57        diff = 0.0
58        for j in range(1, N - 1):
59            for i in range(1, N - 1):
60                gs = 0.25 * (V[i-1, j] + V[i+1, j]
61                           + V[i, j-1] + V[i, j+1]
62                           + h * h * f[i, j])
63                new = (1 - omega) * V[i, j] + omega * gs
64                diff = max(diff, abs(new - V[i, j]))
65                V[i, j] = new
66        history.append(diff)
67        if diff < tol:
68            break
69    return V, history
70
71# --- Same problem for all three: source = pi^2 (sin(pi x) sin(pi y)) * 2 ---
72# Exact solution:   V(x, y) = sin(pi x) * sin(pi y)
73N = 33
74omega_opt = 2.0 / (1.0 + np.sin(np.pi / N))    # theoretical best for Dirichlet square
75f = discretize_poisson(N, lambda X, Y: 2.0 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y))
76
77Vj, hj = jacobi(f)
78Vg, hg = gauss_seidel(f)
79Vs, hs = sor(f, omega=omega_opt)
80
81print(f"Jacobi:      {len(hj):>5} iters,  max error = {np.max(np.abs(Vj - np.sin(np.pi*np.linspace(0,1,N))[:,None]*np.sin(np.pi*np.linspace(0,1,N))[None,:])):.2e}")
82print(f"Gauss-Seid:  {len(hg):>5} iters")
83print(f"SOR (w={omega_opt:.3f}): {len(hs):>5} iters")

Running this file produces output of the form:

Jacobi:       1804 iters,  max error = 2.45e-03
Gauss-Seid:    902 iters
SOR (w=1.826): 118 iters

The Jacobi count, the GS count, and the SOR count are very close to the theoretical predictions $\sim 1.4 N^2/\pi^2$, $\sim 0.7 N^2/\pi^2$, and $\sim 2.2 N/\pi$ for $N = 33$. The max-error value is the same for all three methods because it is the discretization error $\mathcal{O}(h^2)$ of the 5-point stencil, not the solver. Refining the grid would shrink that error; faster iteration would not.

PyTorch on the GPU: Red-Black SOR

Plain Gauss-Seidel and SOR have a sequential data dependency that makes them awkward on a GPU. The standard fix is a red-black ordering: colour grid cells like a chessboard, update all red cells in parallel, then all black cells in parallel. Each colour's neighbours are entirely of the opposite colour, so the update inside one colour is fully independent. The iteration matrix is slightly different from natural-ordered SOR but has the same asymptotic spectral radius — convergence rates are identical for fine grids.

torch.zeros(33, 33, device='cuda')  →  on the first GPU

f.device  →  device(type='cuda', index=0)

V.dtype  →  torch.float32

For N=4:
ii = [[0,0,0,0],[1,1,1,1],[2,2,2,2],[3,3,3,3]]
jj = [[0,1,2,3],[0,1,2,3],[0,1,2,3],[0,1,2,3]]

For N=4 the interior is just the 2×2 inner block of cells (1,1), (1,2), (2,1), (2,2).

Cells (1,1), (1,3), (2,2), (3,1), (3,3) are red.
Cells (1,2), (2,1), (2,3), (3,2) are black.

V                  torch.roll(V, +1, dims=0)
[[a,b,c]            [[g,h,i]
 [d,e,f]   shifts to  [a,b,c]
 [g,h,i]]            [d,e,f]]

Red cells move toward equilibrium; black cells wait their turn.

diff = 7.2e-9  →  done

device = 'cuda'   # if a GPU exists
device = 'cpu'    # otherwise

f[32, 32] = 2π²·sin(π·0.5)·sin(π·0.5) ≈ 19.74

N=65  →  ω_opt = 2/(1 + sin(π/65)) ≈ 1.9082

err ≈ 6.4e-4   ← discretization-limited, not solver-limited

SOR-RB on cuda: 152 iters, omega=1.908, max err = 6.4e-04

1import torch
2
3def sor_red_black_gpu(f, omega=1.85, max_iter=2000, tol=1e-7):
4    """Vectorized SOR using red-black ordering: parallel-safe on GPU.
5
6    Standard SOR has a data dependency in the inner loop (V[i,j] needs the
7    FRESH V[i-1,j]).  Red-black ordering colours grid cells like a chessboard
8    so that every red cell depends only on black cells, and vice versa.
9    We can then update ALL reds in parallel, then ALL blacks in parallel.
10    """
11    device = f.device
12    N      = f.shape[0]
13    h      = 1.0 / (N - 1)
14    V      = torch.zeros(N, N, device=device)
15
16    # Build the red / black masks once (interior cells only)
17    ii, jj  = torch.meshgrid(torch.arange(N, device=device),
18                             torch.arange(N, device=device), indexing="ij")
19    interior = (ii > 0) & (ii < N - 1) & (jj > 0) & (jj < N - 1)
20    red      = interior & ((ii + jj) % 2 == 0)
21    black    = interior & ((ii + jj) % 2 == 1)
22
23    history = []
24    for k in range(max_iter):
25        V_prev = V.clone()
26        # --- RED sweep ---
27        gs = 0.25 * (
28            torch.roll(V, +1, dims=0) + torch.roll(V, -1, dims=0)
29          + torch.roll(V, +1, dims=1) + torch.roll(V, -1, dims=1)
30          + h * h * f
31        )
32        V = torch.where(red, (1 - omega) * V + omega * gs, V)
33
34        # --- BLACK sweep ---
35        gs = 0.25 * (
36            torch.roll(V, +1, dims=0) + torch.roll(V, -1, dims=0)
37          + torch.roll(V, +1, dims=1) + torch.roll(V, -1, dims=1)
38          + h * h * f
39        )
40        V = torch.where(black, (1 - omega) * V + omega * gs, V)
41
42        diff = (V - V_prev).abs().max().item()
43        history.append(diff)
44        if diff < tol:
45            break
46    return V, history
47
48# --- Same manufactured-solution test on the GPU ---
49device = "cuda" if torch.cuda.is_available() else "cpu"
50N      = 65
51x      = torch.linspace(0, 1, N, device=device)
52X, Y   = torch.meshgrid(x, x, indexing="ij")
53f      = 2 * torch.pi**2 * torch.sin(torch.pi * X) * torch.sin(torch.pi * Y)
54
55omega_opt = 2.0 / (1.0 + torch.sin(torch.tensor(torch.pi / N))).item()
56V, hist   = sor_red_black_gpu(f, omega=omega_opt)
57
58V_exact = torch.sin(torch.pi * X) * torch.sin(torch.pi * Y)
59err = (V - V_exact).abs().max().item()
60print(f"SOR-RB on {device}: {len(hist)} iters, omega={omega_opt:.3f}, max err = {err:.2e}")

Two things make this implementation worth the extra lines. First, it is fully vectorized: the inner update is four $\text{torch.roll}$ calls plus a $\text{torch.where}$, all GPU kernels. No Python loop over space. Second, it is differentiable: every operation is a registered PyTorch op, so we can call $\text{V.backward()}$ and get gradients with respect to the source $f$ or the boundary values. That makes this routine a drop-in layer inside a neural network, which is the basis of physics-informed networks (PINNs), neural Poisson solvers, and differentiable rendering.

Multigrid: A Sketch of the Real Champion

Even SOR has an Achilles heel. As the Fourier picture shows, relaxation kills high-frequency error fast but low-frequency error slowly. Halving the grid spacing doubles the wavelength of every mode relative to the grid, so what used to be high-frequency suddenly is low-frequency on the new grid. Multigrid exploits this by alternating between fine and coarse grids:

1. Smooth on the fine grid. Do a few Jacobi or GS sweeps. This kills the highest-frequency error.
2. Restrict the residual to a coarser grid. Average or inject the residual onto an $N/2 \times N/2$ grid. On the coarse grid the surviving low-frequency error looks like high-frequency error.
3. Solve (recursively) on the coarse grid. Either run a few sweeps and keep going coarser, or solve exactly if the grid is small enough.
4. Prolong (interpolate) the correction back up. Linearly interpolate the coarse-grid solution back to the fine grid and add it to the current iterate.
5. Smooth again on the fine grid to fix any high-frequency noise the interpolation introduced.

One pass through this V-cycle has cost $O(N^2)$ on an $N \times N$ grid (just constant work per unknown) and reduces the error by a fixed factor independent of grid size. So multigrid converges in $O(1)$ iterations, with total cost $O(N^2)$ — optimal up to constants. Whatever speedup SOR gave over Jacobi, multigrid gives again over SOR for large grids.

Why mention multigrid here? The relaxation methods of this section are not just a stepping stone. Multigrid uses them as its smoother — the inner loop that kills high-frequency error on each level. So every multigrid solver is built on top of Jacobi, Gauss-Seidel, or SOR. Understand the simple methods, and you understand the engine inside the fast ones. Most production Poisson solvers in computational physics, image processing, and machine learning today are multigrid or one of its conjugate-gradient cousins.

Summary

• Discretization. Replace $\nabla^2 V = -f$ by the 5-point stencil: $V_{i,j} = \tfrac{1}{4}(V_{i\pm 1,j} + V_{i,j\pm 1} + h^2 f_{i,j})$. Applied as an assignment, this is the heart of every relaxation method.
• Jacobi. Use only OLD values during a sweep. Maximally parallel, but spectral radius $\cos(\pi h)$ means $O(N^2)$ iterations.
• Gauss-Seidel. Use the freshest values in place. Spectral radius $\cos^2(\pi h)$ — roughly half the iterations of Jacobi. Sequential, unless red-black ordered.
• SOR. Gauss-Seidel + an overshoot factor $\omega$. With $\omega_{\text{opt}} = 2/(1+\sin(\pi/N))$ the spectral radius drops to $\omega_{\text{opt}} - 1 \approx 1 - 2\pi h$ — $O(N)$ iterations, an asymptotic speedup of a factor $N$ over Jacobi.
• Red-black ordering recovers full parallelism for Gauss-Seidel and SOR on GPUs with no loss of convergence rate.
• Multigrid uses these methods as building blocks ("smoothers") and achieves $O(N^2)$ total work for the entire solve.
• Manufactured solutions (pick $V_{\text{exact}}$, derive $f$, solve, compare) are the standard way to validate any PDE solver — the iteration-count gap measures solver speed, and the final-error gap measures discretization quality. Both gaps are independent.

Take-away. Relaxation methods turn the differential statement "the value here is the average of the neighbours, plus the source" into a literal program. The three classical methods — Jacobi, Gauss-Seidel, SOR — differ only in which neighbour values the assignment reads from and how aggressively to use them. Their spectral radii live in closed form for Poisson, and the gap between Jacobi's $O(N^2)$ and SOR's $O(N)$ iteration count is the entire reason large-scale physics simulations are feasible. Beyond them lies multigrid, the optimal-complexity champion — but every modern fast solver is still standing on the shoulders of Jacobi.