Richardson Extrapolation in Numerical Analysis

• Richardson Extrapolation is a classical method that accelerates the convergence of numerical approximations by cancelling leading error terms using multiple discretization levels.
• It forms linear combinations of approximations based on their asymptotic error expansion, transforming an O(h^p) error into a higher-order O(h^q) error.
• Applied in finite difference and finite element methods, it enhances precision in complex domains, often achieving significantly more correct digits in eigenvalue problems.

Richardson extrapolation is a classical numerical technique for accelerating the convergence of a sequence of approximations to a limiting value, most commonly the result of a discretization method in the small-mesh or small-parameter asymptotic regime. Its outstanding feature is that, given an error expansion with a dominant but unknown leading coefficient, one can combine several approximations at different parameter values to cancel the principal error term and thus systematically increase the convergence rate, often with minimal intrusive modification to existing computational codes.

1. Asymptotic Error Expansion and Principle

The foundation of Richardson extrapolation is the existence of an asymptotic expansion for the numerical method’s error as the discretization parameter $h$ tends to zero, typically of the form

$$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$

where $A^*$ is the exact value, and $p$ is the leading order of the method. In this classical setting, Richardson's device is to form a linear combination of $A(h)$ and $A(\gamma h)$ for some scaling $\gamma$ (often $1/2$ or $b$) that precisely eliminates the $O(h^p)$ error and yields a new approximation whose error is $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$0. For instance, in the basic two-point form: $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$1 Extensions to more complex error structures (e.g., noninteger exponents, reentrant corner singularities, multiple parameters) are handled by solving appropriate linear systems for the extrapolation weights (Amore et al., 2015).

2. Classical and Generalized Formulas

The operation of Richardson extrapolation is algebraically determined by the asymptotic orders present in the error expansion. With two or more grid spacings $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$2 and respective approximate values $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$3, the extrapolated value is chosen such that the combination cancels the leading error terms: $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$4 where each $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$5 is an exponent in the series. A full $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$6-point Richardson extrapolation (with $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$7 grid spacings) can in principle cancel the first $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$8 error terms, and the resulting weights are determined by linear algebra (Vandermonde or generalized systems for rational singular exponents) (Amore et al., 2015). In problems with asymptotic expansions containing noninteger or problem-specific exponents (e.g., domains with reentrant corners), the method remains valid upon suitable generalization of the basis functions in the extrapolation system.

3. Implementation on Complicated Domains

Richardson extrapolation is especially useful for computational spectral problems on non-tensor and irregular domains, where high-order finite difference (FD) or finite element (FE) methods are challenged by the complexity of boundaries or the presence of singularities. The standard procedure involves:

• Generating a sequence of grids (meshes) $$A(h) = A^* + C_1 h^p + C_2 h^{q} + O(h^{r}), \qquad p < q < r,$$9 that exactly sample or resolve the geometric boundary or features of the domain.
• Computing the numerical approximation (e.g., lowest eigenvalue) on each grid, ensuring numerical precision to expose the asymptotic regime.
• Assembling and solving the extrapolation system, either directly (via matrix solve) or recursively via Richardson tables.
• Optionally incorporating Padé-Richardson generalizations, where the asymptotic series is rationally approximated, to further enhance convergence—particularly valuable in problems with limited convergence radius or asymptotic series (Amore et al., 2015).

4. Precision Assessment and Empirical Asymptotics

Validation of Richardson-extrapolated values is performed by comparison against high-precision reference computations, such as those obtained by the Method of Particular Solutions (MPS). A characteristic behavior is observed: error measures, such as the difference between extrapolated and reference values,

$A^*$0

decay rapidly with increasing number of extrapolation points, attain a minimum (optimal truncation), and then increase—a U-shaped curve symptomatic of an asymptotic series. The optimal number of grids or expansion terms is determined at this minimum, yielding the maximal attainable precision for the method. In practical computations, Richardson extrapolation can achieve correct digits far surpassing those available from the raw discretization, often matching or exceeding the precision of more sophisticated boundary-adapted spectral methods (Amore et al., 2015).

5. Numerical Results and Comparative Gains

The effectiveness of Richardson extrapolation in FD eigenvalue problems for planar domains is borne out by extensive numerical studies:

• For the L-shaped domain, two-point Richardson increased the number of correct digits from $A^*$1–$A^*$2 (plain FD) to $A^*$3, with higher-order and fractional-exponent schemes reaching beyond $A^*$4, and Padé enhancements achieving over $A^*$5 digits—on par with 70-digit MPS references.
• For the fundamental mode of the H-shaped domain, classical algorithms achieved 7–11 digits; Richardson extrapolation yielded 19 digits.
• For isospectral drums (Gordon–Webb–Wolpert), the Richardson approach extended the precision from 12 to 17 digits.
• On singular domains (e.g., cracked squares), two-point Richardson and Padé generalizations achieved 22–24 digits, dramatically exceeding the baseline FD accuracy (Amore et al., 2015).

A concise table (subset) illustrates the progression:

Domain	Standard FD Digits	Two-Point Richardson	Multi-Point Richardson	Padé-Richardson	Reference (MPS/FE)
L-shape fundamental	$A^*$6	$A^*$7	$A^*$8	$A^*$9	$p$0
H-shape fundamental	$p$1–$p$2	$p$3	—	—	$p$4
Isospectral drum	$p$5	$p$6	—	—	$p$7
Cracked square	$p$8–$p$9	$A(h)$0	—	$A(h)$1	$A(h)$2

In all cases, the Richardson/Padé-Richardson extrapolation applied to a small hierarchy of meshes delivers rapid and robust error suppression.

6. Empirical Observations and Limitations

Empirical findings indicate that the FD error series on domains with geometric singularities is asymptotic rather than convergent; thus, error minimization and the choice of the number of extrapolation points must be guided by the observed "U-shape" of error decay versus extrapolation level. Excessive truncation leads to divergence in the asymptotic regime, while insufficient points forgoes attainable precision (Amore et al., 2015). The method is fundamentally non-intrusive and black-box, requiring only standard FD computations on meshes sampling the domain—but nothing is gained if the leading error structure is not sufficiently regular (i.e., for nonsmooth or nonconforming grids).

7. Synthesis and Best Practices

The successful application of Richardson extrapolation in high-precision spectral calculations for planar domains proceeds as follows:

1. Analyze the boundary regularity to determine the dominant exponents in the error expansion.
2. Calculate FD eigenvalues on a suitable mesh sequence.
3. Construct and solve the Richardson (or Padé–Richardson) system for the extrapolated value.
4. Identify the error minimum to select the optimal number of points.
5. Validate the results against high-precision or reference solutions.

This approach systematically extends the utility of standard discretization schemes, doubling or more the effective number of significant digits accessible in challenging geometries, and provides an algorithmically simple but powerful tool for precision eigenvalue computation in complicated planar domains (Amore et al., 2015).