= Multigrid method
{wiki}
A <multigrid method> solves a <linear system> by combining local relaxation with <coarse-grid correction> on a hierarchy of <meshes>. The key is to divide the work by scale: a suitable relaxation removes rapidly oscillating <approximation error>, while <coarse-grid correction> addresses slowly varying <approximation error> that relaxation barely changes. A broad error hump occupies many fine-grid points but fewer coarse-grid points; it becomes easier to resolve by relaxation after enough coarsening. This complementary action is why <multigrid methods> are useful for large <finite difference method> and <finite element method> systems.
For a concrete model, take the one-dimensional <Poisson equation> with homogeneous <Dirichlet boundary conditions>:
$$
-u''(x)=f(x),\qquad 0<x<1,\qquad u(0)=u(1)=0.
$$
Use $n$ interior points $x_j=jh$, $h=1/(n+1)$, and the centered <finite difference> system
$$
(A_hu_h)_j=\frac{2u_{h,j}-u_{h,j-1}-u_{h,j+1}}{h^2}=b_{h,j},\qquad b_{h,j}=f(x_j),\qquad u_{h,0}=u_{h,n+1}=0.
$$
Thus $A_h=h^{-2}\operatorname{tridiag}(-1,2,-1)$ is a <symmetric matrix> and a <positive-definite matrix>. The minus sign in $-u''$ fixes this positive convention; changing both $A_h$ and $b_h$ to their negatives leaves the relaxation update unchanged. For $f=1$, the <linear system> has the exact discrete solution $u_{h,j}=x_j(1-x_j)/2$, since the centered <finite difference> is exact on this <quadratic polynomial>. In general, the exact discrete solution need not equal the sampled continuous solution: solving the <linear system> removes algebraic <approximation error>, not the <global discretization error>.
Let $u_h^*$ solve $A_hu_h^*=b_h$, let $v_h$ be the current approximation, and define the signed <approximation error> $e_h=u_h^*-v_h$. The <linear-system residual> is $r_h=b_h-A_hv_h$, so $A_he_h=r_h$. In particular, a <linear-system residual> is available without knowing $u_h^*$; it supplies the right-hand side for a correction equation. With $D=2h^{-2}I$, a <weighted Jacobi method> sweep is
$$
v_h\leftarrow v_h+\omega D^{-1}r_h,\qquad e_h\leftarrow S_\omega e_h,\qquad S_\omega=I-\omega D^{-1}A_h.
$$
The <sine> <eigenvectors> shared by $A_h$ and the <iteration matrix> $S_\omega$ are $w_j^{(k)}=\sin(k\pi j/(n+1))$, $1\leq k\leq n$. Their <eigenvalues> are
$$
\lambda_k(A_h)=4h^{-2}\sin^2\frac{k\pi}{2(n+1)},\qquad \mu_k(S_\omega)=1-2\omega\sin^2\frac{k\pi}{2(n+1)}.
$$
For $n=63$, one <weighted Jacobi method> sweep with $\omega=2/3$ multiplies the lowest <Fourier frequency> mode $k=1$ by about $0.999197$, but the highest <Fourier frequency> mode $k=63$ by about $-0.332530$. Four sweeps retain about $99.68\%$ of the first amplitude and $1.22\%$ of the second. The sign reversal of the second amplitude is harmless when its magnitude shrinks. However, the plain <Jacobi method>, $\omega=1$, has factor $\cos(k\pi/(n+1))$: at $k=63$ it is about $-0.998795$. A nearly alternating mode near the grid's Nyquist limit therefore persists almost as strongly as the smoothest mode. Plain <Jacobi method> convergence on a fixed <Dirichlet boundary condition> grid does not imply effective smoothing.<Introduction to Multigrid Methods by Long Chen>{ref} The <high-frequency smoothing factor of weighted Jacobi> is uniformly at most $1/3$ for $\omega=2/3$ and $k\geq(n+1)/2$; the <exact finite-grid high-frequency Jacobi smoothing optimum> distinguishes this convenient uniform choice from the finite-grid optimum.<CS267: Supplemental Notes for Lecture 17 — How Multigrid Damps High Frequencies by James Demmel>{ref}
\Image[/multigrid-coarse-grid-smoothing.png]
{title=Analytic Poisson error modes: weighted Jacobi damps short wavelengths while plain Jacobi struggles near the Nyquist limit}
One two-grid cycle applies a few <weighted Jacobi method> sweeps, computes the <linear-system residual>, performs <coarse-grid correction>, and applies a few more sweeps. The first relaxation makes the remaining <approximation error> easier to represent on a coarse <mesh>; the last relaxation damps oscillatory components left after the transferred correction. In a <geometric multigrid V-cycle>, replace the expensive exact coarse solve by one recursive cycle on the coarse correction equation, starting its correction from zero. Repeat until a tiny coarsest <linear system> can be solved directly, then return through the levels by adding each interpolated correction and post-smoothing. The descent and return give the cycle its V shape.<A Multigrid Tutorial (lecture slides) by William L. Briggs>[slides 59–63]{ref}
For this one-dimensional model, fixed sweep counts and bounded local <finite difference> and transfer work give a <geometric multigrid V-cycle> cost proportional to $n+n/2+n/4+\cdots<2n$, hence $O(n)$ in <Big O notation>. A fixed error reduction also needs a contraction factor bounded below one independently of $h$; per-cycle work alone does not establish that. Starting from arbitrary <approximation error> and demanding a relative reduction $\varepsilon$ generally adds a factor proportional to $\log(1/\varepsilon)$ in the cycle count. <Full multigrid> uses coarse solutions as initial guesses on successively finer <meshes> to target the <global discretization error> efficiently. The <two-grid Poisson factor with weighted Jacobi> illustrates favorable convergence for this model. Other <partial differential equations>, strong <anisotropy>, or poorly chosen coarse spaces and relaxation can defeat these choices; a linear-work guarantee for arbitrary <linear systems> does not follow.
A historical example of the multilevel viewpoint is Achi Brandt's April 1977 paper, which analyzes adaptive <meshes> and reports numerical experiments on linear and nonlinear <boundary value problems>. This is a documented development of the approach, not a claim that the paper was its first invention.<Multi-Level Adaptive Solutions to Boundary-Value Problems by Achi Brandt>[pp. 333–390; abstract]{ref} The figure above is generated from the displayed <eigenvalues>, with lines joining discrete nodes; it contains no experimental measurements.
Back to article page
