What Jacobi Cannot Kill in 10,000 Sweeps, a V-Cycle Kills in 10 — The Real Secret of Multigrid

Why 10,000 Jacobi sweeps refuse to finish#

Engineer B was solving a pressure Poisson on a 1024×1024 grid with point Jacobi. Two hundred seconds in, the residual was still flickering at the same order of magnitude. The CFL was tiny. The grid was orthogonal. The code was correct. Some errors Jacobi never kills.

The same error changes frequency when the grid changes. The V-cycle is that sentence turned into an algorithm.

Move Pre/Post between 1 and 5 and, on the same 64-cell grid with a zero right-hand side, only the V-cycle residual steepens, then the extra smoothing stops paying.

The orange trace is Jacobi. The wiggles die and the broad hump stays. The cyan log-residual falls almost as a straight line, while orange flattens on the low-frequency shelf.

Error frequencies — who dies fast and who lives forever#

Suppose the linear system is $A u = f$. The error in the current guess $u^{(k)}$ is

one vector. The residual is $r^{(k)} = f - A u^{(k)} = A e^{(k)}$.

A spectral look at weighted Jacobi says a Fourier mode $\sin(j \pi x)$ in the error decays per sweep as

where $\omega$ is the relaxation weight, $j$ the mode number, and $h$ the spacing. The absolute damping factors:

Short waves shrink to a third each sweep. Long waves barely move. What Jacobi cannot kill is the smooth hill spread over the whole grid.

On a coarser grid that low frequency looks high#

This is where multigrid splits from a plain smoother. The longest wave on $N$ cells is $\lambda \sim N h$. Move it to $N/2$ cells and the same wave occupies half as many cells. It looks shorter. More precisely, the lowest fine mode is a middle mode on the coarse grid, and Jacobi there kills middle modes.

Coarsen again and the wave shortens again. After $\log_2 N$ levels every mode is high-frequency somewhere. A few Jacobi sweeps on every level, and in theory every error dies. Leaving the low modes to a Krylov comparison is the point of convergence from Jacobi through BiCGSTAB.

V-cycle — the path that goes down and comes back up#

Turn the split into steps and the shape is a V.

1. Pre-smoothing — $\nu_1$ Jacobi sweeps on the fine grid.
2. Residual — $r_h = f_h - A_h u_h$.
3. Restriction — $r_{2h} = R\, r_h$. Carry the residual down.
4. Recurse — solve $A_{2h} e_{2h} = r_{2h}$ with another V-cycle, or directly if the grid is small enough.
5. Prolongation — $e_h = P\, e_{2h}$. Lift the correction.
6. Correct — $u_h \leftarrow u_h + e_h$.
7. Post-smoothing — $\nu_2$ Jacobi sweeps on the fine grid.

Downward, only the residual moves. Upward, only the correction. One cycle costs $O(N)$: one fine sweep plus a geometric tail. One cycle drops about one digit of residual.

Restriction and prolongation — the exchange rate between two grids#

The two transfer operators have to be a pair. On a 1D vertex-centered grid the plain pair is this.

Full-weighting restriction: coarse point $i$ is a weighted average of fine points $2i-1, 2i, 2i+1$.

Linear prolongation: copy onto coarse-aligned nodes, average the two neighbors in between.

They are transposes up to a grid-ratio constant ($P = c R^T$). That pair is part of the convergence argument. An asymmetric pair can diverge.

A two-grid V-cycle in NumPy#

Solve $-u''(x) = f(x)$ on $[0,1]$ with $u(0)=u(1)=0$. Start at 64 cells.

import numpy as np
 
class GridLevel:
    """State at one grid level: solution, RHS, spacing."""
    def __init__(self, n):
        self.n = n              # number of cells
        self.h = 1.0 / n
        self.u = np.zeros(n + 1)
        self.f = np.zeros(n + 1)
 
def jacobi_relax(u, f, h, iters):
    """Weighted Jacobi for 1D Poisson. Dirichlet zero at both ends."""
    n = len(u) - 1
    nxt = np.zeros_like(u)
    for _ in range(iters):
        nxt[1:n] = 0.5 * (u[0:n-1] + u[2:n+1] + h * h * f[1:n])
        u[1:n] = nxt[1:n]
    return u
 
def residual_norm(u, f, h):
    n = len(u) - 1
    r = np.zeros_like(u)
    inv = 1.0 / (h * h)
    r[1:n] = f[1:n] - (-u[0:n-1] + 2 * u[1:n] - u[2:n+1]) * inv
    return float(np.sqrt(np.mean(r[1:n] ** 2)))
 
def restrict_full_weighting(rh):
    n = len(rh) - 1
    n2 = n // 2
    r2 = np.zeros(n2 + 1)
    r2[1:n2] = 0.25 * rh[1:n-1:2] + 0.5 * rh[2:n:2] + 0.25 * rh[3:n+1:2]
    return r2
 
def prolong_linear(e2h):
    n2 = len(e2h) - 1
    n = 2 * n2
    eh = np.zeros(n + 1)
    eh[2:n:2] = e2h[1:n2]
    eh[1:n+1:2] = 0.5 * (e2h[0:n2] + e2h[1:n2+1])
    return eh
 
def v_cycle(u, f, h, n_pre=2, n_post=2):
    n = len(u) - 1
    if n <= 4:
        return jacobi_relax(u, f, h, 50)
    jacobi_relax(u, f, h, n_pre)
    inv = 1.0 / (h * h)
    r = np.zeros_like(u)
    r[1:n] = f[1:n] - (-u[0:n-1] + 2 * u[1:n] - u[2:n+1]) * inv
    r2 = restrict_full_weighting(r)
    e2 = np.zeros_like(r2)
    v_cycle(e2, r2, 2 * h, n_pre, n_post)
    e = prolong_linear(e2)
    u += e
    jacobi_relax(u, f, h, n_post)
    return u
 
# --- run ---
n = 64
h = 1.0 / n
x = np.linspace(0, 1, n + 1)
f = (np.pi ** 2) * np.sin(np.pi * x)   # exact u = sin(pi x)
 
u_jac = np.zeros(n + 1)
u_mg = np.zeros(n + 1)
 
for k in range(20):
    jacobi_relax(u_jac, f, h, 4)
    v_cycle(u_mg, f, h, n_pre=2, n_post=2)
    print(f"step {k:3d} | jacobi r={residual_norm(u_jac, f, h):.3e} "
          f"v-cycle r={residual_norm(u_mg, f, h):.3e}")

One printed step spends 4 Jacobi sweeps, and the V-cycle spends 4 fine sweeps (2 pre + 2 post). Similar work, different residuals. After 20 steps the V-cycle residual is usually below $10^{-10}$ and Jacobi is stuck near $10^{-2}$.

The table closes at $\omega = 2/3$. The fastest mode, $jh = 1$, has $\cos(\pi) = -1$, so the factor is $1 - 2\omega = -1/3$. The absolute value matches the table entry $1/3$. Ten sweeps on that mode alone leave $(1/3)^{10} \approx 1.7 \times 10^{-5}$. The middle entry is $\approx 0.67$, and $0.67^{20} \approx 3.3 \times 10^{-4}$. The slowest factor is $1 - O(h^2)$. At $N = 64$, $h = 1/64$, so $h^2 = 1/4096 \approx 2.4 \times 10^{-4}$. That is the size of one chop on the long wave.

The printout lines up with the table. Four Jacobi sweeps per step, twenty steps, eighty sweeps, and the residual is still near $10^{-2}$. The fast mode is already dead at a factor $1/3$. The leftover $10^{-2}$ is the low mode in the table. The same twenty V-cycles, 2 pre and 2 post on the fine grid, usually sit under $10^{-10}$. Eight orders separate $10^{-2}$ and $10^{-10}$.

The coarsest level is 50 Jacobi sweeps once $n \le 4$. Sixty-four cells restrict as 64, 32, 16, 8, 4 — four transfers. Full weighting is $1/4$, $1/2$, $1/4$. The coarse work sums as $N/2 + N/4 + \cdots$, which does not exceed $N$. Ten V-cycles repeat that sum ten times. Ten thousand Jacobi sweeps are a different count before the constant even matters.

Put the same factor on the slowest mode. At $N = 64$, $jh = 1/N$ and $\omega = 2/3$, the factor is about $0.9992$. The listing's eighty sweeps multiply to about $0.94$. One digit costs about $2900$ sweeps. At $N = 1024$ the factor is about $0.999997$, and one digit costs about $7.3 \times 10^{5}$ sweeps. Ten thousand sweeps still leave that mode near $0.97$. The 1024×1024 run that sat in one order of magnitude for 200 seconds is the same stall. The formula is one-dimensional weighted Jacobi. The engineer's grid was two-dimensional. A V-cycle that drops one digit per cycle does in about ten cycles what those sweeps do not do to the slow mode.

The factor $0.9992$ is not a measured residual. It is the table formula on the slowest mode. Eighty sweeps move it by about six percent. The printed Jacobi residual, near $10^{-2}$ after twenty steps, is what a factor that close to 1 looks like once the fast mode — damped by $1/3$ — is already gone.

Why the V-cycle is not always faster — two traps#

A V-cycle is not magic. Two places break it in practice.

Rough coefficients. On anisotropic or jumping-coefficient Poisson, Jacobi is a poor smoother. In $-\nabla \cdot (k(x) \nabla u) = f$, if $k$ jumps by orders of magnitude the standard V-cycle parks the residual at some digit. The fix is line-Jacobi, a line-aligned grid, or algebraic multigrid, which builds the hierarchy from the matrix.

Non-nested grids. On unstructured meshes or AMR the fine-to-coarse map is not an index halving. Agglomeration multigrid bundles fine cells into coarse cells, restricts by volume weights, and prolongs with a cell-centered interpolant. OpenFOAM's GAMG is that pattern.

The rate is mesh-independent only when the smoother, the operator, and the coarsening ratio all match. Miss one and the iteration count creeps up with $N$. That is when a W-cycle or an F-cycle is the next try.

Where to put the standard V-cycle down#

If $k$ jumps by orders of magnitude, plain Jacobi is not a smoother and the residual sticks at some digit. Switch to line-Jacobi or AMG. The stretched-grid version of that switch is a linelet preconditioner.

Unstructured meshes and AMR do not build a coarse level by halving an index. OpenFOAM GAMG agglomerates cells. If the iteration count starts tracking $N$, try a W-cycle or an F-cycle. The 64-cell run, 50 Jacobi sweeps at $n \le 4$, and the 20-step gap of $10^{-10}$ versus $10^{-2}$ were read on aligned 1D Poisson only.