Richardson extrapolation
Richardson extrapolation is a technique in numerical analysis used to enhance the accuracy of approximate solutions to mathematical problems by systematically eliminating leading-order error terms through the combination of results computed at multiple step sizes. Developed by British mathematician and physicist Lewis Fry Richardson in his 1911 paper on solving differential equations via finite differences, the method assumes that an approximation T(h) to a true value T(0) can be expressed as T(h) = T(0) + c_p h^p + c_{p+1} h^{p+1} + \cdots, where p is the order of the leading error term and c_p is a constant.[1][2] The core idea involves computing approximations at step sizes h and kh (often k=2), then forming a linear combination that cancels the h^[p](/page/P′′) term, resulting in an improved estimate of order p+1 or higher. For instance, the extrapolated value is given by T^*(h) = \frac{k^[p](/page/P′′) T(h) - T(kh)}{k^[p](/page/P′′) - 1}, which removes the dominant error assuming the expansion holds. This process can be iterated recursively to generate even higher-order approximations, forming a table of values where each entry refines the previous ones. The technique is particularly effective when the underlying approximation method has a known asymptotic error expansion, and it requires no additional assumptions about the problem beyond the existence of such an expansion.[3][2] Richardson extrapolation finds wide application in various numerical methods, including differentiation, where it upgrades central difference formulas from second-order to fourth-order accuracy; integration, serving as the foundation for Romberg integration, which achieves exponential convergence for smooth integrands; and the solution of ordinary and partial differential equations, improving finite difference schemes. Its efficiency stems from reusing prior computations in the recursive scheme, though it is sensitive to round-off errors in floating-point arithmetic, necessitating careful choice of step sizes to balance truncation and rounding effects. Despite these limitations, the method remains a cornerstone of computational mathematics due to its simplicity and power in accelerating convergence without modifying the base algorithm.[3][4]Fundamentals
Definition and Purpose
Richardson extrapolation is a sequence acceleration method that combines numerical approximations obtained at different step sizes to eliminate leading-order error terms and attain higher-order accuracy in estimating a limit value.[5] It is named after Lewis Fry Richardson, who introduced the technique in 1911 while developing methods for the approximate arithmetical solution of physical problems involving differential equations, with an application to the stresses in a masonry dam.[1] The primary purpose of Richardson extrapolation is to enhance the precision of numerical approximations in methods like finite differences, where errors decrease asymptotically with smaller step sizes, by effectively extrapolating toward the exact solution as the step size approaches zero.[6] This approach achieves improved accuracy without necessitating additional evaluations of the underlying function beyond the initial approximations at selected step sizes, making it computationally efficient for refining results in numerical analysis.[7] At its core, Richardson extrapolation relies on the assumption that the approximation error admits an asymptotic expansion in powers of the step size h. For instance, consider an approximation A(h) to a true value L, expressed as A(h) = L + c h^k + O(h^{k+1}), where c and k > 0 are constants, and higher-order terms are of order O(h^{k+1}). By computing approximations at step sizes h and h/2 (or multiples thereof) and linearly combining them, the leading h^k error term can be canceled, yielding an improved estimate with error dominated by the next higher-order term.[5] This process systematically boosts the order of convergence, provided the asymptotic regime is reached where higher-order terms are negligible.[6]Historical Background
Richardson extrapolation originated in the early 20th century through the work of British mathematician and physicist Lewis Fry Richardson, who developed the technique to improve the accuracy of finite difference approximations in solving partial differential equations. In 1911, Richardson introduced the method in a seminal paper addressing geophysical computations, particularly the numerical solution of Laplace's equation for determining stresses in a masonry dam using irregular grids.[1] This approach allowed him to refine coarse approximations by extrapolating toward the limit as the grid size approached zero, effectively eliminating leading-order truncation errors without requiring finer meshes.[1] During the 1920s, Richardson extended these ideas to practical applications in numerical weather prediction and structural stress analysis. In his 1922 book Weather Prediction by Numerical Process, he employed finite difference methods to compute atmospheric pressure tendencies, aiming to forecast weather patterns through direct numerical integration of hydrodynamic equations. This work represented an early attempt at systematic numerical simulation of complex geophysical phenomena, where such methods were applied on limited computational resources. Concurrently, Richardson applied similar refinement techniques to stress analysis in engineering structures, building on his 1911 foundations to handle irregular boundaries and non-uniform grids.[8] The method gained formal recognition in 1927 when Richardson, collaborating with J. Arthur Gaunt, published "The Deferred Approach to the Limit," which generalized the extrapolation process for sequences derived from lattice-based discretizations of differential equations.[9] This paper articulated the deferred correction strategy, emphasizing its utility in iteratively approaching exact solutions by postponing higher-order terms.[9] Following a period of limited adoption, Richardson extrapolation experienced a revival and formalization in the mid-20th century as numerical analysis matured alongside electronic computing. In their 1957 book Methods of Numerical Integration, Philip J. Davis and Philip Rabinowitz popularized the technique by integrating it into broader frameworks for sequence acceleration and error estimation in quadrature rules.[10] They linked it explicitly to extrapolation of asymptotic expansions, making it accessible for general numerical integration problems and highlighting its role in improving convergence rates.[10] A key milestone in the 1950s came with Werner Romberg's 1955 development of an iterative scheme for numerical integration based on repeated Richardson extrapolation applied to the trapezoidal rule, which achieved higher-order accuracy efficiently.[11] This connection, further disseminated in the 1960s through computational implementations and textbooks, embedded the method within standard numerical toolkits, influencing adaptive algorithms for solving integral equations.[12] Throughout its evolution, Richardson extrapolation has played a pivotal role in accelerating the convergence of iterative methods in numerical analysis, predating and inspiring modern adaptive techniques by providing a systematic way to combine multiple approximations for enhanced precision without proportional increases in computational cost.[13]Theoretical Framework
Notation
In Richardson extrapolation, the approximation obtained from a numerical method with step size h is denoted by A(h), which converges to the true value L as h approaches zero.[14] The sequence of step sizes is typically defined as h_m = h / b^m for nonnegative integers m \geq 0, where b > 1 is the extrapolation base, commonly taken as b = 2 to halve the step size iteratively.[14] The error in the approximation admits an asymptotic expansion of the form A(h) = L + \sum_{k=1}^\infty c_k h^{p k}, where p > 0 is the order of the underlying method and the coefficients c_k are constants independent of h. This form assumes the error terms appear at orders that are integer multiples of p, as is typical in certain numerical methods such as the trapezoidal rule for integration.[14] The following table summarizes the principal symbols employed:| Symbol | Description |
|---|---|
| h | Step size |
| T_{m,k} | Extrapolated approximation at level m and order k |
| b | Extrapolation base (typically b=2) |
General Formula
The general formula for Richardson extrapolation arises from the asymptotic error expansion of a numerical approximation A(h) to a limiting value L = \lim_{h \to 0} A(h), typically expressed as A(h) = L + c h^p + O(h^{2p}), where p > 0 is the known order of the leading error term and c is a constant.[3] To eliminate the h^p term, consider approximations at step sizes h and h/b (with b > 1), yielding A(h) = L + c h^p + O(h^{2p}) and A(h/b) = L + c (h/b)^p + O(h^{2p}).[15] Multiplying the second equation by b^p gives b^p A(h/b) = b^p L + c h^p + O(h^{2p}). Subtracting the first equation from this scaled version cancels the c h^p terms: b^p A(h/b) - A(h) = (b^p - 1) L + O(h^{2p}). Solving for L produces the extrapolated approximation T_0 = \frac{b^p A(h/b) - A(h)}{b^p - 1} = L + O(h^{2p}), which achieves an error of order $2p.[3] This formula, introduced by Lewis Fry Richardson in his 1911 paper on finite difference solutions to differential equations, assumes prior knowledge of the method's error order p and the step ratio b, often chosen as an integer like 2 for computational convenience.[1][15] For higher-order extrapolations, the process iterates on successively refined approximations. Define T_{m,k} as the k-th order extrapolation at level m. The recursive formula is T_{m,k} = \frac{b^{p k} T_{m,k-1} - T_{m-1,k-1}}{b^{p k} - 1}, with base cases T_{m,0} = A(h / b^m) for m = 0, 1, \dots.[3] This yields T_{m,m} = L + O(h^{p(m+1)}), progressively eliminating higher-order error terms through repeated application, provided the error expansion holds up to the desired order.[4] The assumption of known p and fixed b ensures the coefficients align correctly for cancellation, though variations exist for unknown or variable orders in advanced extensions.[15]Recurrence Relation
The recurrence relation for Richardson extrapolation provides a systematic way to construct higher-order approximations by iteratively combining prior estimates in a tableau. In standard notation, the entries T_{m,k} of the tableau satisfy the relation T_{m,k} = T_{m,k-1} + \frac{T_{m,k-1} - T_{m-1,k-1}}{b^{p k} - 1}, where m indexes the row (corresponding to step size level), k indexes the column (extrapolation order), b > 1 is the base refinement factor (often b = 2), and p is the order of the base approximation method.[16] This formula is derived from differencing consecutive levels to cancel the leading error term, assuming the error expansion A(h) = A + c_p h^p + c_{p+1} h^{p+1} + \cdots.[17] The tableau is constructed starting with the first column, where each entry is the base approximation at successively refined step sizes: T_{m,0} = A(h / b^m) for m = 0, 1, \dots, n-1. Subsequent columns are then filled recursively from left to right and top to bottom using the recurrence relation, with each T_{m,k} depending only on the immediately preceding entries in the same row and the diagonal above.[16] This process builds a lower triangular array, where the diagonal entries T_{m,m} provide the highest-order extrapolations at each level.[17] For n distinct step sizes, the recurrence enables computation of the full tableau in O(n^2) operations, yielding approximations up to order p + n - 1 along the anti-diagonal.[16] Regarding error propagation, the recurrence eliminates the h^p term at the first extrapolation level (k=1), the h^{p+1} term at the second level (k=2), and so on, with the k-th column entry T_{m,k} having error O(h^{p+k}) under suitable smoothness assumptions on the underlying function.[17] This stepwise cancellation leverages the known powers in the asymptotic error expansion to progressively refine accuracy without additional base evaluations beyond the initial column.[16]Properties
Richardson extrapolation converges to the true limit L under the assumption that the approximation possesses a complete asymptotic expansion in powers h^p, h^{2p}, h^{3p}, \dots of the step size h, as the step sizes are successively refined toward zero; the order of accuracy increases along the diagonals of the extrapolation tableau, achieving superlinear convergence rates under suitable conditions on the expansion coefficients.[18] The recurrence relation facilitates this diagonal order enhancement by systematically eliminating leading error terms.[18] The method exhibits bounded stability for the sequences in its tableau columns and diagonals when the underlying approximations satisfy the asymptotic assumptions, ensuring that perturbations do not amplify uncontrollably in the limit processes.[19] However, practical stability is compromised at higher extrapolation levels, where sensitivity to rounding errors intensifies due to the ill-conditioned nature of the computations, arising from subtractions of closely valued terms that magnify floating-point inaccuracies.[20] Richardson extrapolation is particularly effective and optimal when the leading error order p is known a priori, allowing precise cancellation of the dominant term; it underperforms or fails entirely if the error lacks the required asymptotic expansion, as occurs with non-smooth functions where higher-order terms vanish or do not conform to the power series form.[21] Key limitations stem from the necessity of step sizes forming a geometric sequence on the logarithmic scale—typically h_i = h_0 b^{-i} for integer b > 1—to align with the powers in the expansion; deviations disrupt the cancellation process.[18] Computationally, building the full tableau incurs a quadratic cost in the number of refinement levels, as each new level requires reevaluating all prior approximations.[21] Furthermore, the error magnification factor for perturbations at level k approximates b^{p k}, resulting in exponential growth of rounding error effects with deeper extrapolation.[20]Computational Aspects
Step-by-Step Process
The application of Richardson extrapolation follows a structured algorithm that begins with selecting an appropriate base numerical method and proceeds through iterative refinement to produce higher-order approximations. First, identify the base method, such as a finite difference scheme, and determine its error order p, which represents the leading term in the asymptotic error expansion; this order is typically known from the method's theoretical analysis.[22][23] Next, choose a refinement factor b, commonly 2 for simplicity and computational efficiency, and decide on the number of levels n based on desired accuracy and available resources; larger n allows for more extrapolation steps but increases computational cost, so it is selected by estimating convergence rates or monitoring error reduction in preliminary runs.[24][15] Compute the initial approximations by evaluating the base method at successively refined parameters: start with step size h to obtain A(h), then A(h/b), A(h/b^2), up to A(h/b^{n-1}); these values form the first column of the extrapolation tableau, a triangular array where rows correspond to the refinement levels.[22][23] To build the tableau, apply the recurrence relation row by row, using pairs of entries from the previous column to compute entries in the next column, thereby eliminating successive error terms and achieving higher-order accuracy in each subsequent column.[15][24] The process fills the tableau diagonally, with each new row incorporating a finer approximation and each new column providing an improved estimate. If the error order p is unknown, estimate it iteratively by comparing differences between consecutive approximations and adjusting the refinement weights accordingly, often starting with an assumed value and refining based on observed convergence behavior. Once the tableau is complete, select the improved estimate from the diagonal elements (for the highest order at the finest level) or the bottom-right corner, depending on the desired balance of order and refinement.[23][22] To ensure reliability, verify the underlying assumptions of the asymptotic expansion throughout the process, such as by plotting residuals between base and extrapolated approximations to confirm that errors diminish as expected with refinement; deviations may indicate invalid assumptions like non-dominant higher-order terms or instability. Stability checks, drawing from the method's properties, can be integrated by monitoring for oscillatory behavior in the tableau entries during construction. The overall workflow can be visualized as a flowchart: begin with input of the base method, order p, factor b, and level count n; branch to compute the sequence of approximations A(h/b^k) for k = 0 to n-1; proceed to initialize and populate the tableau row by row via recurrence; converge to output the selected estimate, followed by a verification loop for residual analysis and potential iteration adjustment.[22][23]Pseudocode Example
A pseudocode implementation of Richardson extrapolation typically constructs a tableau to iteratively refine a sequence of base approximations using the recurrence relation, assuming the approximations are ordered from coarsest to finest step size. The input consists of a list of base approximations A = [A_0, A_1, \dots, A_{n-1}], where A_i is computed with step size h / 2^i, and the known error order p (e.g., p=2 for quadratic convergence). The output is the highest-order extrapolated estimate T_{n-1, n-1}.[25] The algorithm initializes the first column of the tableau with the input approximations and then fills subsequent columns using the general recurrence. This process builds higher-order accurate estimates diagonally across the tableau.This pseudocode translates the step-by-step process into a programmatic form, suitable for languages like Python or MATLAB. The time and space complexity is O(n^2), as it involves n(n+1)/2 recurrence computations. For efficiency in array-based languages, vectorized implementations can compute entire columns simultaneously using broadcasting, reducing loop overhead.[26]function richardson_extrapolation(A_list, p): n = length(A_list) if n == 1: return A_list[0] // No extrapolation possible if p is unknown or invalid: error("Error order p must be specified") T = 2D array of size n x n // Initialize tableau for i = 0 to n-1: T[i][0] = A_list[i] // Base approximations in first column for k = 1 to n-1: for m = k to n-1: factor = 2^(p * k) T[m][k] = (factor * T[m][k-1] - T[m-1][k-1]) / (factor - 1) return T[n-1][n-1] // Highest-order estimatefunction richardson_extrapolation(A_list, p): n = length(A_list) if n == 1: return A_list[0] // No extrapolation possible if p is unknown or invalid: error("Error order p must be specified") T = 2D array of size n x n // Initialize tableau for i = 0 to n-1: T[i][0] = A_list[i] // Base approximations in first column for k = 1 to n-1: for m = k to n-1: factor = 2^(p * k) T[m][k] = (factor * T[m][k-1] - T[m-1][k-1]) / (factor - 1) return T[n-1][n-1] // Highest-order estimate
Applications and Examples
Simple Numerical Example
A simple numerical example of Richardson extrapolation involves approximating \sin(1) \approx 0.84147098 using the identity \sin(1) = \frac{\sin(1+h) + \sin(1-h)}{2 \cos h}, where \cos h is approximated by its Taylor series truncated after the h^2 term: \cos h \approx 1 - \frac{h^2}{2}. This truncation yields an approximation A(h) = \frac{\sin(1+h) + \sin(1-h)}{2 \left(1 - \frac{h^2}{2}\right)}, with leading error O(h^4) arising from the omitted h^4/24 term in the cosine expansion.[15] Computations are performed for h = 0.1, h = 0.05, and h = 0.025, assuming exact evaluation of the sine terms:- A(0.1) = 0.84147449
- A(0.05) = 0.84147120
- A(0.025) = 0.84147099
| j \setminus k | 0 | 1 | 2 |
|---|---|---|---|
| 0 | 0.84147449 | 0.84147098 | 0.84147098 |
| 1 | 0.84147120 | 0.84147098 | |
| 2 | 0.84147099 |
Use in Numerical Differentiation
Richardson extrapolation enhances the accuracy of numerical differentiation by combining finite difference approximations computed at successively halved step sizes to cancel lower-order error terms. The central difference formula provides a second-order approximation to the first derivative, D(h) = \frac{f(x + h) - f(x - h)}{2h} = f'(x) + \frac{h^2}{6} f'''(\xi) + O(h^4) for some \xi \in (x - h, x + h), where the leading error is O(h^2).[27] To apply Richardson extrapolation, compute D(h), D(h/2), and D(h/4), then form the first-level extrapolations as T^{(1)}(h) = \frac{4 D(h/2) - D(h)}{3}, \quad T^{(1)}(h/2) = \frac{4 D(h/4) - D(h/2)}{3}, which eliminate the O(h^2) term, yielding O(h^4) accuracy. A second-level extrapolation further improves the order: T^{(2)}(h) = \frac{16 T^{(1)}(h/2) - T^{(1)}(h)}{15}, removing the O(h^4) term for O(h^6) or higher accuracy. This process systematically increases the order of convergence without requiring higher-order finite difference stencils directly.[3] For illustration, consider approximating f'(1) where f(x) = \sin x, with exact value \cos 1 \approx 0.540302305868. Using initial step size h = 0.01, the approximations and extrapolations form the following tableau (values rounded to 7 decimal places for clarity; errors shown relative to the true value):| Level \ Step | h = 0.01 | h = 0.005 | h = 0.0025 |
|---|---|---|---|
| 0 (D(h)) | 0.5402933 (error: -9.0 \times 10^{-6}) | 0.5403001 (error: -2.2 \times 10^{-6}) | 0.5403017 (error: -5.6 \times 10^{-7}) |
| 1 (T^{(1)}) | 0.5403023 (error: $1.4 \times 10^{-7}) | 0.5403023 (error: $2.4 \times 10^{-8}) | |
| 2 (T^{(2)}) | 0.540302306 (error: < 10^{-9}) |