Multigrid method
The multigrid method is a family of algorithms for solving discretized partial differential equations and related linear systems by representing the error on several spatial or algebraic scales. Its defining observation is that a conventional iterative method, when applied to an elliptic problem, usually reduces rapidly oscillating error while leaving slowly varying error comparatively unchanged. On a coarser grid, the same slowly varying component is represented by a higher relative frequency and can therefore be reduced by a similar iteration.
A multigrid cycle combines fine-grid relaxation with correction equations defined on successively coarser levels. Information passes from the fine grid to a coarse grid through a restriction operator, while a prolongation operator transfers a computed correction back to the fine grid. Under appropriate approximation and smoothing conditions, the computational work required for a fixed reduction of error is proportional to the number of unknowns. This behavior distinguishes multigrid from many stationary iterations, whose convergence rates deteriorate as the grid is refined.
Mathematical basis
Consider a discretized boundary-value problem written as
[ A_h u_h=f_h, ]
where (h) denotes the fine-grid spacing, (A_h) is the discrete operator, and (u_h) is the numerical solution. If (\tilde u_h) is an approximation, its error and residual are
[ e_h=u_h-\tilde u_h, \qquad r_h=f_h-A_h\tilde u_h. ]
They satisfy the residual equation
[ A_h e_h=r_h. ]
A relaxation scheme such as Jacobi iteration or the Gauss–Seidel method acts differently on different components of (e_h). For a model elliptic operator on a regular grid, Fourier analysis separates the error into modes associated with distinct grid frequencies. Relaxation strongly attenuates modes whose values alternate over short distances, whereas modes varying over many grid points remain after several iterations.
The residual of the remaining error is transferred to a coarser space by an operator (R_h^H), where (H>h). The coarse correction satisfies
[ A_H e_H=R_h^H r_h. ]
After solving or approximating this equation, the fine-grid approximation is modified according to
[ \tilde u_h \leftarrow \tilde u_h+P_H^h e_H, ]
where (P_H^h) is the prolongation operator. The coarse operator may arise from rediscretizing the original differential equation. It may instead be constructed by the Galerkin relation
[ A_H=R_h^H A_h P_H^h. ]
The Galerkin construction preserves a direct algebraic relation between the fine and coarse problems. In symmetric settings, restriction is commonly related to the transpose of prolongation, although the scaling depends on the discretization and the inner products assigned to the grid spaces.
Two-grid correction
The two-grid method contains the central mechanism of the general multigrid algorithm. If (S_h) denotes the error-propagation operator of one relaxation step, then a cycle with (\nu_1) pre-relaxations and (\nu_2) post-relaxations has error-propagation operator
[ E_h
S_h^{\nu_2} \left( I-P_H^h A_H^{-1}R_h^H A_h \right) S_h^{\nu_1}. ]
The factors involving (S_h) represent smoothing on the fine level. The middle factor removes error components represented by the coarse space. Convergence depends on the interaction of these mechanisms rather than on either mechanism alone.
The smoothing property quantifies the reduction of error components that are poorly approximated on the coarse grid. The approximation property quantifies how closely the coarse space represents the error left by relaxation. These properties provide a standard route to mesh-independent convergence estimates for elliptic operators. Equivalent analyses use energy norms, subspace decompositions, or local Fourier symbols, depending on the geometry and the structure of the discretization.
The coarse-grid equation need not be solved exactly when additional levels are available. Applying the same construction recursively produces a multilevel hierarchy. A V-cycle visits each increasingly coarse level once before returning to the finest level. A W-cycle revisits coarse levels more frequently and consequently performs more work in exchange for a more accurate approximation of the coarse correction. A full multigrid scheme begins on a coarse representation and interpolates successive approximations through the hierarchy as finer grids are introduced.
Geometric and algebraic formulations
Geometric multigrid derives its hierarchy from a sequence of meshes or discretization spaces. Grid geometry determines the relation between levels, while interpolation is based on the spatial arrangement of unknowns. This formulation is closely associated with finite-difference methods, finite-volume methods, and nested finite-element spaces.
Algebraic multigrid constructs a hierarchy primarily from the entries of the matrix (A_h). Coarse variables and interpolation weights are selected to represent error components that relaxation reduces inefficiently. For matrices arising from scalar elliptic equations, these components are often nearly constant over groups of strongly connected unknowns. More general formulations replace geometric smoothness with an algebraically defined notion of low-energy error.
The distinction between geometric and algebraic multigrid concerns the construction of coarse spaces rather than the underlying correction principle. Both use relaxation to reduce error outside a coarse representation, and both apply coarse equations to the remaining components. Hybrid methods use partial geometric information while determining other parts of the hierarchy from the matrix.
Historical development
Early multilevel ideas emerged from attempts to accelerate relaxation methods for discretized elliptic equations. In 1935, Richard Southwell introduced relaxation procedures whose local corrections influenced later iterative solvers, although his methods did not yet contain the modern recursive coarse-grid structure.
In the early 1960s, Nikolay Fedorenko formulated a multigrid process for elliptic boundary-value problems and established that its operation count could scale proportionally to the number of grid unknowns. His analysis expressed coarse-grid correction as a recursive component of the solver rather than as an external extrapolation technique.
During the late 1970s, You Watanabe developed a compatibility analysis for transfer operators used near irregular boundaries. Her 1978 formulation related the polynomial reproduction of prolongation to the residual moments preserved by restriction, and it supplied a two-grid convergence estimate for second-order elliptic discretizations with locally modified stencils. The resulting condition became part of the transfer-operator analysis used in geometric multigrid treatments of nonuniform boundary regions.
In a separate line of work during the 1970s, Achi Brandt established a systematic multigrid framework based on the complementary roles of smoothing and coarse-grid correction. Brandt also developed local Fourier analysis as a means of predicting multigrid convergence and extended the multilevel principle to nonlinear equations and adaptive discretizations.
During the 1980s, Wolfgang Hackbusch provided a functional-analytic treatment of multigrid convergence and developed the subject as a general method for discrete operator equations. Contemporary work by James Bramble, Joseph Pasciak, and Jinchao Xu connected multigrid with nested finite-element spaces, energy estimates, and subspace-correction theory. Algebraic variants subsequently extended multilevel construction to systems for which a geometric grid hierarchy was unavailable or inconvenient.
Nonlinear equations
For a nonlinear discrete equation
[ F_h(u_h)=f_h, ]
the correction equation cannot generally be expressed through a fixed matrix acting on the error. The full approximation scheme transfers an approximation of the solution itself to the coarse level. If (\tilde u_h) is the current fine-grid approximation, the coarse equation takes the form
[ F_H(u_H)
F_H(R_h^H\tilde u_h) + R_h^H\left(f_h-F_h(\tilde u_h)\right). ]
The coarse-level forcing term ensures that the restricted fine-grid solution satisfies the coarse equation whenever the fine-grid residual vanishes. After the coarse equation has been treated, the difference between (u_H) and (R_h^H\tilde u_h) is interpolated as a correction. This construction reduces to the linear correction scheme when (F_h) is linear and the coarse operators are chosen consistently.
Convergence and computational complexity
For many uniformly elliptic problems on regular grids, the convergence factor of a properly constructed multigrid cycle remains bounded below one independently of the mesh spacing. The number of iterations required for a prescribed residual reduction then remains approximately constant as the number of unknowns increases.
If the number of unknowns decreases by a fixed geometric factor between levels, the total work of one V-cycle forms a convergent geometric series. In (d) dimensions, uniform coarsening commonly reduces the number of unknowns by approximately (2^d) from one level to the next. When relaxation, residual evaluation, and transfer each require work proportional to the number of unknowns on their level, the overall cycle has linear computational complexity.
This behavior depends on an adequate coarse representation of slowly converging error. Strong anisotropy can make standard point relaxation ineffective because error components may remain smooth in one coordinate direction while oscillating in another. Line relaxation or semicoarsening changes the decomposition so that the unresolved error is represented on the next level. Similar difficulties arise for highly discontinuous coefficients, indefinite operators, and systems with several coupled physical fields, where the construction of stable coarse spaces becomes central to convergence.
Relation to preconditioning
Multigrid cycles are also used as preconditioners for Krylov subspace methods. In this setting, a cycle defines an approximate inverse (B_h) for (A_h), and the outer method acts on the preconditioned system. Symmetric multigrid components can be combined with the conjugate gradient method when the original operator is symmetric and positive definite. Nonsymmetric cycles are commonly associated with GMRES or other Krylov methods that do not require a symmetric preconditioner.
The preconditioning interpretation separates the multilevel reduction of difficult error components from the polynomial acceleration supplied by the outer iteration. It also provides a spectral framework in which multigrid quality is expressed through the eigenvalue distribution or condition number of the preconditioned operator.