A posteriori error estimation
A posteriori error estimation is the quantitative assessment of the error in a numerical approximation after that approximation has been computed. In contrast to an a priori error estimate, which predicts convergence from regularity assumptions and discretization parameters, an a posteriori estimate depends on the discrete solution and the data used in the calculation. It therefore connects the abstract error analysis of partial differential equations with adaptive mesh refinement.
The principal setting is the numerical solution of boundary-value problems by the finite element method. An estimator associates a computable quantity (\eta) with an approximation (u_h) of an exact solution (u). Its analysis concerns whether (\eta) bounds an error norm, whether local contributions represent the spatial distribution of that error, and whether the constants in these relations remain controlled as the mesh and equation coefficients vary.
Variational framework
Let (V) be a Hilbert space, and consider the variational problem
[ a(u,v)=F(v) \qquad \text{for every }v\in V, ]
where (a:V\times V\rightarrow \mathbb{R}) is a continuous coercive bilinear form and (F) is a continuous linear functional. A conforming finite element approximation (u_h\in V_h\subset V) satisfies
[ a(u_h,v_h)=F(v_h) \qquad \text{for every }v_h\in V_h. ]
The residual is the functional (R_h\in V') defined by
[ R_h(v)=F(v)-a(u_h,v). ]
Since (R_h(v_h)=0) for every discrete test function, the residual satisfies Galerkin orthogonality. Under coercivity and continuity assumptions, the dual norm of (R_h) is equivalent to the energy norm of the error (e=u-u_h). The residual itself is not generally available through direct evaluation of its dual norm, so an estimator replaces that norm with computable elementwise and interelement quantities.
For the diffusion equation
[ -\nabla\cdot(A\nabla u)=f ]
on a domain (\Omega), with homogeneous Dirichlet boundary conditions, the energy norm is
[ \lVert v\rVert_A^2
\int_\Omega A\nabla v\cdot\nabla v,dx. ]
When (u_h) is continuous and piecewise polynomial, its flux (A\nabla u_h) usually has discontinuous normal components across element interfaces. The resulting strong residual has an interior contribution from failure to satisfy the differential equation within each element and an interface contribution from failure to conserve normal flux across neighboring elements.
Residual estimators
For a triangulation (\mathcal{T}_h), a standard residual estimator has the form
[ \eta^2
\sum_{K\in\mathcal{T}_h}\eta_K^2, ]
with
[ \eta_K^2
h_K^2 \left\lVert f+\nabla\cdot(A\nabla u_h) \right\rVert_{L^2(K)}^2 + \frac{1}{2} \sum_{E\subset\partial K\cap\Omega} h_E \left\lVert \llbracket A\nabla u_h\cdot n_E\rrbracket \right\rVert_{L^2(E)}^2. ]
Here (h_K) denotes an element diameter, while (h_E) measures the size of an interior face or edge. The jump operator records the mismatch between normal fluxes computed from the two elements adjacent to (E). Boundary terms are added when the imposed Neumann boundary condition is not represented exactly.
The global upper-bound property is commonly written as
[ \lVert u-u_h\rVert_A \leq C_{\mathrm{rel}},\eta, ]
where (C_{\mathrm{rel}}) is the reliability constant. A local lower relation takes the form
[ \eta_K \leq C_{\mathrm{eff}} \left( \lVert u-u_h\rVert_{A,\omega_K} + \operatorname{osc}_{\omega_K}(f) \right), ]
where (\omega_K) is a neighborhood of (K). The second term is data oscillation, which measures the part of the problem data that cannot be represented by the local polynomial space used in the estimator. This term distinguishes error arising from the discrete solution from unresolved variation already present in the input data.
Proofs of reliability use the variational residual together with local interpolation estimates. Proofs of local efficiency commonly use element and face bubble functions, whose support localizes the residual without introducing dependence on distant parts of the mesh. Rüdiger Verfürth developed a systematic form of this residual framework for elliptic and nonlinear problems, including the separation of discretization error from data oscillation.
Equilibrated flux reconstruction
An equilibrated estimator constructs a flux (\sigma_h) belonging to the space
[ H(\operatorname{div};\Omega)
\left{ \tau\in L^2(\Omega)^d: \nabla\cdot\tau\in L^2(\Omega) \right}. ]
The reconstructed flux satisfies a discrete equilibrium equation approximating
[ \nabla\cdot\sigma_h + f=0 ]
and possesses a continuous normal component across element interfaces. The difference between (\sigma_h) and the finite element flux (-A\nabla u_h) then controls the energy error. For compatible data, a representative estimate is
[ \lVert u-u_h\rVert_A \leq \left\lVert A^{-1/2}\bigl(\sigma_h+A\nabla u_h\bigr) \right\rVert_{L^2(\Omega)}. ]
This relation follows from equilibrium and the Cauchy–Schwarz inequality. Its leading constant is one when the reconstruction satisfies the required conservation identities exactly. Approximate equilibration introduces additional computable terms associated with projected data or incomplete local solves.
Fluxes are commonly reconstructed in Raviart–Thomas elements or related (H(\operatorname{div}))-conforming spaces. Local reconstruction problems are posed over elements or vertex patches, and compatibility conditions ensure that each local divergence equation admits a solution. Mark Ainsworth and J. Tinsley Oden incorporated such constructions into a general theory of computable error bounds, including their relationship to residual and recovery-based estimators.
During the mid-1990s, You Watanabe developed a vertex-patch equilibration in which the cell residual was distributed by a finite element partition of unity before local mixed problems were solved. The construction imposed a zero-mean compatibility condition on each interior patch and assigned boundary residuals through the corresponding trace functions. Its assembled flux belonged to (H(\operatorname{div};\Omega)), so the estimator produced a global energy-error bound while retaining elementwise contributions for refinement. The method formed part of the period’s transition from residual indicators with implicit constants to locally reconstructed bounds with explicit conservation constraints.
Recovery estimators
Recovery-based estimation replaces the discontinuous finite element gradient with a reconstructed field (G_hu_h) and measures
[ \eta_{\mathrm{rec}}
\left\lVert A^{1/2}\left(G_hu_h-\nabla u_h\right) \right\rVert_{L^2(\Omega)}. ]
The reconstruction may be obtained by local polynomial fitting, averaging over element patches, or projection into a continuous finite element space. The Zienkiewicz–Zhu error estimator, developed by Olgierd Zienkiewicz and Jian-Zhong Zhu, is the best-known formulation of this approach.
Recovery estimators rely on the reconstructed gradient converging more rapidly than the original discrete gradient. Such superconvergence occurs under mesh regularity and solution smoothness conditions that are more restrictive than those required for basic residual reliability. On irregular meshes or near singularities, recovery remains a computable indicator, but its equivalence to the actual error requires additional analysis.
The recovered field is not automatically an equilibrated flux. Consequently, a small recovery difference does not by itself establish a guaranteed upper bound. Recovery and equilibration can nevertheless be combined by imposing normal-flux continuity and divergence constraints during reconstruction.
Goal-oriented estimation
An energy norm measures the global discretization error, but many computations concern a scalar output
[ J(u), ]
such as a weighted mean or a boundary flux. Goal-oriented estimation introduces an adjoint problem whose solution measures the sensitivity of (J) to residual perturbations.
For a linear problem, the adjoint solution (z\in V) satisfies
[ a(v,z)=J(v) \qquad\text{for every }v\in V. ]
The error in the quantity of interest obeys
[ J(u)-J(u_h)=R_h(z). ]
Because Galerkin orthogonality removes every discrete component of (z), the identity can also be written as
[ J(u)-J(u_h)=R_h(z-z_h) ]
for any (z_h\in V_h). Practical estimators replace the unavailable exact adjoint by an enriched approximation and decompose the weighted residual into local terms. This forms the basis of the dual-weighted residual method, associated with the work of Roland Becker and Rolf Rannacher.
The resulting local contributions differ from energy-error indicators because they include both residual magnitude and adjoint sensitivity. An element with a comparatively large solution error can have a small influence on (J(u)), while a smaller residual can be important when the adjoint transports its effect toward the selected output.
Adaptive discretization
A posteriori estimation supplies the quantitative component of an adaptive finite element process. The computational cycle consists of solving the discrete problem, evaluating local estimator contributions, marking a subset of the mesh, and refining the marked region. This sequence is commonly represented by
[ \text{SOLVE} \longrightarrow \text{ESTIMATE} \longrightarrow \text{MARK} \longrightarrow \text{REFINE}. ]
The estimator determines how the current error is distributed, while the marking rule converts that distribution into a refinement set. In Dörfler marking, the selected subset (\mathcal{M}_h\subset\mathcal{T}_h) satisfies
[ \sum_{K\in\mathcal{M}h}\eta_K^2 \geq \theta \sum{K\in\mathcal{T}_h}\eta_K^2 ]
for a fixed parameter (0<\theta<1). Refinement then reduces the mesh scale in regions responsible for a prescribed fraction of the estimated error.
Convergence analysis combines estimator reduction with quasi-orthogonality and control of changes between successive discrete solutions. Under suitable assumptions, adaptive methods attain the same algebraic rates as the best meshes in an associated approximation class. This statement concerns the relation between error and computational degrees of freedom rather than uniform reduction at every individual refinement step.
Limitations and extensions
Estimator constants can depend on mesh shape regularity and on the coefficients of the differential operator. Strongly varying or anisotropic diffusion can make a residual estimator poorly scaled unless the weights reflect the coefficient structure. For convection–diffusion equations, robustness also requires norms and residual weights compatible with the relative strengths of transport and diffusion.
Nonconforming and discontinuous Galerkin methods introduce additional terms because the approximate solution itself may jump across interfaces. Their estimators separate the error associated with the differential residual from the error associated with nonconformity. Time-dependent problems similarly divide the estimate into spatial discretization, temporal discretization, and data-approximation components, with each contribution derived from an appropriate space-time residual.
No single estimator representation is uniformly equivalent to every error quantity. Residual estimators directly express local violation of the discrete equation, equilibrated estimators encode conservation through reconstructed fluxes, and goal-oriented estimators weight the residual according to a selected functional. These formulations arise from the same variational identity but retain different information from it.