= Richardson extrapolation
{c}
{title2=$Y_{\rm ext}=(2^pY_{h/2}-Y_h)/(2^p-1)$}
{wiki}
<Richardson extrapolation> improves an approximation by combining two resolutions so that their leading <truncation errors> cancel. The useful surprise is that the size of that leading error need not be known: its scaling under refinement is enough. This lets an existing <finite difference method> or <quadrature rule> produce a more accurate answer from two runs, without first deriving a new method from scratch.
Suppose $A(h)$ approximates a fixed <limit> $L$ as the positive resolution parameter $h$ tends to zero, with an <asymptotic expansion>
$$
A(h)=L+c h^p+O(h^{p+q}),\qquad p>0,\quad q>0.
$$
Here $p$ is known, $c$ is independent of $h$, and $O$ denotes <Big O notation>. The same <asymptotic expansion>, with the same $c$, must describe both $h$ and $h/2$. Nested grids make it natural to reuse samples, although nesting itself does not establish this error model. Define the <linear combination>
$$
R(h)=\frac{2^p A(h/2)-A(h)}{2^p-1}
=A(h/2)+\frac{A(h/2)-A(h)}{2^p-1}.
$$
The weights sum to one, preserving $L$, while $2^p c(h/2)^p-c h^p=0$. Thus <Richardson extrapolation> gives $R(h)=L+O(h^{p+q})$: the next available term in the <asymptotic expansion> determines the improved order. Justin Solomon's numerical analysis text develops this cancellation for <finite differences> and applies it to <numerical integration>.https://people.csail.mit.edu/jsolomon/share/book/numerical_book.pdf[Justin Solomon, Numerical Algorithms, §14.3.3 and Exercise 14.9]{ref}
Geometrically, regard the approximation as a <function> of $t=h^p$. Its first two terms are a straight line $L+ct$. The two computed points lie at positive $t$; extending the line through them to $t=0$ gives $R(h)$. This is why the operation is extrapolation: it estimates a value outside the sampled parameter interval. The remaining curvature corresponds to the uncancelled terms of the <asymptotic expansion>.
\b[Worked example: a centered derivative of sine]. For the <smooth function> $f(x)=\sin x$, the exact <derivative> at $x=0$ is $L=\cos0=1$. The centered <difference quotient>, a <central finite difference>, is
$$
A(h)=\frac{f(h)-f(-h)}{2h}=\frac{\sin h}{h}.
$$
The <Taylor series> of <sine> gives the signed <truncation error>
$$
A(h)-1=-\frac{h^2}{6}+\frac{h^4}{120}+O(h^6).
$$
Thus $p=2$, $q=2$ and $c=-1/6$. Halving $h$ divides the leading <truncation error> by four. Combining the two centered <difference quotients> with <Richardson extrapolation> gives
$$
R(h)=\frac{4A(h/2)-A(h)}{3},\qquad
R(h)-1=-\frac{h^4}{480}+O(h^6).
$$
Indeed, the $h^4$ coefficient becomes $(4/16-1)/(3\cdot120)=-1/480$; the $h^2$ coefficient vanishes. This is the even-power case of <Richardson extrapolation> discussed in Persson's authored lecture notes.https://persson.berkeley.edu/math128a/lecture_notes/math128a_lecture_notes_chapter_4.html[Per-Olof Persson, Math 128A, Chapter 4, §§4.1–4.2]{ref}
For $h=0.4$ and $h/2=0.2$, the actual computations are
$$
\begin{aligned}
A(0.4)&\approx0.973545855771626,\\
A(0.2)&\approx0.993346653975306,\\
R(0.4)&\approx0.999946920043199.
\end{aligned}
$$
The absolute error drops from about $2.6454\times10^{-2}$ for $A(0.4)$ to $5.3080\times10^{-5}$ for $R(0.4)$, an improvement by roughly a factor of 498. Refining again gives $R(0.2)\approx0.999996670632607$, whose absolute error is about $3.3294\times10^{-6}$. The refined <truncation errors> decrease by nearly four for $A$ and nearly sixteen for $R$, reflecting their second and fourth orders.
\Image[/richardson-derivative-error-orders.png]
{title=Centered sine derivative errors: second order before Richardson extrapolation and fourth order after cancellation}
{source=richardson-derivative-error-orders.py}
The <logarithmic scale> on each axis turns power laws into straight lines. The dashed reference curves are $h^2/6$ and $h^4/480$, the magnitudes of the leading <truncation errors>. The plotted errors are computed with 70-digit decimal arithmetic, rather than clipped to a display floor; this isolates the <asymptotic expansion> from <roundoff error>.
\b[Worked example: integrating an exponential]. Approximate the <definite integral> of the <exponential function> on $[0,1]$:
$$
L=\int_0^1 e^x\,dx=e-1\approx1.718281828459045.
$$
For <numerical integration>, the composite <trapezoidal rule> adds the areas of trapezoids between successive samples; this is the quadrature use of the endpoint average that also occurs in the ODE <trapezoidal rule>. On $n$ equal intervals, with $h=1/n$, define
$$
T_n=h\left[\frac{1+e}{2}+\sum_{j=1}^{n-1}e^{jh}\right].
$$
Summing the <geometric series> gives $T_n=(e-1)\frac h2\frac{e^h+1}{e^h-1}$. Expanding the <exponential function> near zero yields
$$
T_n-L=L\left(\frac{h^2}{12}-\frac{h^4}{720}+O(h^6)\right).
$$
The leading <truncation error> again has $p=2$, $q=2$. On the uniform nested grids with $n$ and $2n$ intervals, <Richardson extrapolation> therefore gives
$$
S_{2n}=\frac{4T_{2n}-T_n}{3},\qquad
S_{2n}-L=\frac{Lh^4}{2880}+O(h^6).
$$
This <linear combination> is exactly composite <Simpson's rule> on the finer grid, whose spacing is $h/2$. On each coarse panel, an old endpoint receives weight $h/6$ and the new <midpoint> weight $2h/3$; shared old endpoints receive weight $h/3$. These are the usual <Simpson's rule> weights $1,4,2,\ldots,4,1$ multiplied by $h/6$. The equivalence is explained directly in LeMesurier's university notes.https://lemesurierb.people.charleston.edu/elementary-numerical-analysis-python/notebooks/definite-integrals-3-simpson-richardson.html[Brenton LeMesurier, MATH 375, §19.1: The Basic Simpson's Rule by Richardson Extrapolation]{ref}
For two coarse panels and four fine panels,
$$
\begin{aligned}
T_2&\approx1.753931092464825,\\
T_4&\approx1.727221904557517,\\
S_4&\approx1.718318841921747.
\end{aligned}
$$
The absolute errors are about $0.0356493$, $0.00894008$ and $0.0000370135$, respectively. Reusing the fine-grid samples makes <Richardson extrapolation> especially attractive here: the resulting <Simpson's rule> needs no samples beyond those already used by $T_{2n}$.
\Image[/richardson-trapezoid-extrapolation.png]
{title=Exponential quadrature on nested grids: trapezoidal areas, reused sample nodes, and the fourth-order Simpson combination}
{source=richardson-trapezoid-extrapolation.py}
\b[When the improvement is justified]. <Richardson extrapolation> requires the stated <asymptotic expansion> and the correct positive $p$. If the actual leading exponent is $r\ne p$, its coefficient after combination is $c(2^{p-r}-1)/(2^p-1)$, so the leading <truncation error> generally remains. Merely knowing $A(h)-L=O(h^p)$ does not supply a common coefficient $c$ or a higher-order remainder. A <discontinuity of a function>, a singular endpoint, or unrelated changes between grids can invalidate the assumed <asymptotic expansion>; refinement alone does not guarantee improvement for every numerical problem.
<Roundoff error> and errors in the sampled values survive the cancellation of <truncation error>. If the coarse and fine values have absolute perturbations bounded by $\varepsilon$, the perturbation in $R(h)$ is bounded by $(2^p+1)\varepsilon/(2^p-1)$; for $p=2$ this is $5\varepsilon/3$. In <floating-point arithmetic>, a centered <difference quotient> can also amplify sample errors through division by $h$, so an indefinitely smaller $h$ need not be more accurate. Persson's notes describe this balance between <truncation error> and <roundoff error>.https://persson.berkeley.edu/math128a/lecture_notes/math128a_lecture_notes_chapter_4.html[Per-Olof Persson, Math 128A, Chapter 4, §4.1: Roundoff Error Instability]{ref}
The same cancellation idea appears in <step-doubling error estimation>. For an order-$p$ one-step method for an <ordinary differential equation>, two half steps each contribute <local truncation error> of order $h^{p+1}$, so their accumulated leading local error has ratio $2^{-p}$ to one full step. This explains the denominator $2^p-1$ there; the appropriate exponent follows from that local error model, rather than automatically from a generic use of <Richardson extrapolation>.
Back to article page
