Richardson extrapolation

In numerical analysis, Richardson extrapolation is a method used to estimate some value A ∗ = lim h → 0 A ( h ) {\displaystyle A^{\ast }=\lim _{h\to 0}A(h)} if the truncation error is known to have a (possibly 1-sided) power-series expansion at h = 0 {\displaystyle h=0}. The method is most often applied as a sequence acceleration method to improve the rate of convergence of iterative methods. It is named after Lewis Fry Richardson, who introduced the technique in the early 20th century, though a form of the idea was already known to Christiaan Huygens in the 17th century and used his calculation of π {\displaystyle \pi }. Practical applications of Richardson extrapolation include Romberg integration, which applies Richardson extrapolation to the trapezoid rule, and the Bulirsch–Stoer algorithm for solving ordinary differential equations. In the words of Birkhoff and Rota, "its usefulness for practical computations can hardly be overestimated."

General formula

Notation

Let A 0 ( h ) {\displaystyle A_{0}(h)} be an approximation of A ∗ {\displaystyle A^{*}}(exact value) that depends on a step size h (where 0 < h < 1 {\textstyle 0<h<1}) with an error formula of the form A ∗ − A 0 ( h ) = a 0 h k 0 + a 1 h k 1 + a 2 h k 2 + ⋯ {\displaystyle A^{*}-A_{0}(h)=a_{0}h^{k_{0}}+a_{1}h^{k_{1}}+a_{2}h^{k_{2}}+\cdots } where the a i {\displaystyle a_{i}} are unknown constants and the k i {\displaystyle k_{i}} are known constants such that h k i > h k i + 1 {\displaystyle h^{k_{i}}>h^{k_{i+1}}}. Furthermore, O ( h k i ) {\displaystyle O(h^{k_{i}})} represents the truncation error of the A i ( h ) {\displaystyle A_{i}(h)} approximation such that A ∗ − A i ( h ) = O ( h k i ) . {\displaystyle A^{*}-A_{i}(h)=O(h^{k_{i}}).} Note that by simplifying with Big O notation, the following asymptotic formulae are implied: A ∗ − A 0 ( h ) = a 0 h k 0 + a 1 h k 1 + a 2 h k 2 + ⋯ = a 0 h k 0 + O ( h k 1 ) = O ( h k 0 ) {\displaystyle {\begin{aligned}A^{*}-A_{0}(h)&=a_{0}h^{k_{0}}+a_{1}h^{k_{1}}+a_{2}h^{k_{2}}+\cdots \\&=a_{0}h^{k_{0}}+O(h^{k_{1}})\\&=O(h^{k_{0}})\end{aligned}}}

Purpose

Richardson extrapolation is a process that finds a better approximation of A ∗ {\displaystyle A^{*}} by changing the error formula from A ∗ = A 0 ( h ) + O ( h k 0 ) {\displaystyle A^{*}=A_{0}(h)+O(h^{k_{0}})} to A ∗ = A 1 ( h ) + O ( h k 1 ) . {\displaystyle A^{*}=A_{1}(h)+O(h^{k_{1}}).} Therefore, by replacing A 0 ( h ) {\displaystyle A_{0}(h)} with A 1 ( h ) {\displaystyle A_{1}(h)} the truncation error has reduced from O ( h k 0 ) {\displaystyle O(h^{k_{0}})} to O ( h k 1 ) {\displaystyle O(h^{k_{1}})} for the same step size h {\displaystyle h}. The general pattern occurs in which A i ( h ) {\displaystyle A_{i}(h)} is a more accurate estimate than A j ( h ) {\displaystyle A_{j}(h)} when i > j {\displaystyle i>j}. By this process, we have achieved a better approximation of A ∗ {\displaystyle A^{*}} by subtracting the largest term in the error which was O ( h k 0 ) {\displaystyle O(h^{k_{0}})}. This process can be repeated to remove more error terms to get even better approximations.

Process

Using the step sizes h {\displaystyle h} and h / t {\displaystyle h/t} for some constant t {\displaystyle t}, the two formulas for A ∗ {\displaystyle A^{*}} are:

A ∗ = A 0 ( h ) + a 0 h k 0 + a 1 h k 1 + a 2 h k 2 + O ( h k 3 ) {\displaystyle A^{*}=A_{0}(h)+a_{0}h^{k_{0}}+a_{1}h^{k_{1}}+a_{2}h^{k_{2}}+O(h^{k_{3}})}

A ∗ = A 0 ( h t ) + a 0 ( h t ) k 0 + a 1 ( h t ) k 1 + a 2 ( h t ) k 2 + O ( h k 3 ) {\displaystyle A^{*}=A_{0}\!\left({\frac {h}{t}}\right)+a_{0}\left({\frac {h}{t}}\right)^{k_{0}}+a_{1}\left({\frac {h}{t}}\right)^{k_{1}}+a_{2}\left({\frac {h}{t}}\right)^{k_{2}}+O(h^{k_{3}})}

To improve our approximation from O ( h k 0 ) {\displaystyle O(h^{k_{0}})} to O ( h k 1 ) {\displaystyle O(h^{k_{1}})} by removing the first error term, we multiply equation 2 by t k 0 {\displaystyle t^{k_{0}}} and subtract equation 1 to give us ( t k 0 − 1 ) A ∗ = [ t k 0 A 0 ( h t ) − A 0 ( h ) ] + ( t k 0 a 1 ( h t ) k 1 − a 1 h k 1 ) + ( t k 0 a 2 ( h t ) k 2 − a 2 h k 2 ) + O ( h k 3 ) . {\displaystyle (t^{k_{0}}-1)A^{*}={\bigg [}t^{k_{0}}A_{0}\left({\frac {h}{t}}\right)-A_{0}(h){\bigg ]}+{\bigg (}t^{k_{0}}a_{1}{\bigg (}{\frac {h}{t}}{\bigg )}^{k_{1}}-a_{1}h^{k_{1}}{\bigg )}+{\bigg (}t^{k_{0}}a_{2}{\bigg (}{\frac {h}{t}}{\bigg )}^{k_{2}}-a_{2}h^{k_{2}}{\bigg )}+O(h^{k_{3}}).} This multiplication and subtraction was performed because [ t k 0 A 0 ( h t ) − A 0 ( h ) ] {\textstyle {\big [}t^{k_{0}}A_{0}\left({\frac {h}{t}}\right)-A_{0}(h){\big ]}} is an O ( h k 1 ) {\displaystyle O(h^{k_{1}})} approximation of ( t k 0 − 1 ) A ∗ {\displaystyle (t^{k_{0}}-1)A^{*}}. We can solve our current formula for A ∗ {\displaystyle A^{*}} to give A ∗ = [ t k 0 A 0 ( h t ) − A 0 ( h ) ] t k 0 − 1 + ( t k 0 a 1 ( h t ) k 1 − a 1 h k 1 ) t k 0 − 1 + ( t k 0 a 2 ( h t ) k 2 − a 2 h k 2 ) t k 0 − 1 + O ( h k 3 ) {\displaystyle A^{*}={\frac {{\bigg [}t^{k_{0}}A_{0}\left({\frac {h}{t}}\right)-A_{0}(h){\bigg ]}}{t^{k_{0}}-1}}+{\frac {{\bigg (}t^{k_{0}}a_{1}{\bigg (}{\frac {h}{t}}{\bigg )}^{k_{1}}-a_{1}h^{k_{1}}{\bigg )}}{t^{k_{0}}-1}}+{\frac {{\bigg (}t^{k_{0}}a_{2}{\bigg (}{\frac {h}{t}}{\bigg )}^{k_{2}}-a_{2}h^{k_{2}}{\bigg )}}{t^{k_{0}}-1}}+O(h^{k_{3}})} which can be written as A ∗ = A 1 ( h ) + O ( h k 1 ) {\displaystyle A^{*}=A_{1}(h)+O(h^{k_{1}})} by setting A 1 ( h ) = t k 0 A 0 ( h t ) − A 0 ( h ) t k 0 − 1 . {\displaystyle A_{1}(h)={\frac {t^{k_{0}}A_{0}\left({\frac {h}{t}}\right)-A_{0}(h)}{t^{k_{0}}-1}}.}

Recurrence relation

A general recurrence relation can be defined for the approximations by A i + 1 ( h ) = t k i A i ( h t ) − A i ( h ) t k i − 1 {\displaystyle A_{i+1}(h)={\frac {t^{k_{i}}A_{i}\left({\frac {h}{t}}\right)-A_{i}(h)}{t^{k_{i}}-1}}} where k i + 1 {\displaystyle k_{i+1}} satisfies A ∗ = A i + 1 ( h ) + O ( h k i + 1 ) . {\displaystyle A^{*}=A_{i+1}(h)+O(h^{k_{i+1}}).}

Properties

The Richardson extrapolation can be considered as a linear sequence transformation.

Additionally, the general formula can be used to estimate k 0 {\displaystyle k_{0}} (leading order step size behavior of Truncation error) when neither its value nor A ∗ {\displaystyle A^{*}} is known a priori. Such a technique can be useful for quantifying an unknown rate of convergence. Given approximations of A ∗ {\displaystyle A^{*}} from three distinct step sizes h {\displaystyle h}, h / t {\displaystyle h/t}, and h / s {\displaystyle h/s}, the exact relationship A ∗ = t k 0 A i ( h t ) − A i ( h ) t k 0 − 1 + O ( h k 1 ) = s k 0 A i ( h s ) − A i ( h ) s k 0 − 1 + O ( h k 1 ) {\displaystyle A^{*}={\frac {t^{k_{0}}A_{i}\left({\frac {h}{t}}\right)-A_{i}(h)}{t^{k_{0}}-1}}+O(h^{k_{1}})={\frac {s^{k_{0}}A_{i}\left({\frac {h}{s}}\right)-A_{i}(h)}{s^{k_{0}}-1}}+O(h^{k_{1}})} yields an approximate relationship (please note that the notation here may cause a bit of confusion, the two O appearing in the equation above only indicates the leading order step size behavior but their explicit forms are different and hence cancelling out of the two O terms is only approximately valid) A i ( h t ) + A i ( h t ) − A i ( h ) t k 0 − 1 ≈ A i ( h s ) + A i ( h s ) − A i ( h ) s k 0 − 1 {\displaystyle A_{i}\left({\frac {h}{t}}\right)+{\frac {A_{i}\left({\frac {h}{t}}\right)-A_{i}(h)}{t^{k_{0}}-1}}\approx A_{i}\left({\frac {h}{s}}\right)+{\frac {A_{i}\left({\frac {h}{s}}\right)-A_{i}(h)}{s^{k_{0}}-1}}} which can be solved numerically to estimate k 0 {\displaystyle k_{0}} for some arbitrary valid choices of h {\displaystyle h}, s {\displaystyle s}, and t {\displaystyle t}.

As t ≠ 1 {\displaystyle t\neq 1}, if t > 0 {\displaystyle t>0} and s {\displaystyle s} is chosen so that s = t 2 {\displaystyle s=t^{2}}, this approximate relation reduces to a quadratic equation in t k 0 {\displaystyle t^{k_{0}}}, which is readily solved for k 0 {\displaystyle k_{0}} in terms of h {\displaystyle h} and t {\displaystyle t}.

Example of Richardson extrapolation

Suppose that we wish to approximate A ∗ {\displaystyle A^{*}}, and we have a method A ( h ) {\displaystyle A(h)} that depends on a small parameter h {\displaystyle h} in such a way that A ( h ) = A ∗ + C h n + O ( h n + 1 ) . {\displaystyle A(h)=A^{\ast }+Ch^{n}+O(h^{n+1}).}

Let us define a new functionR ( h , t ) := t n A ( h / t ) − A ( h ) t n − 1 {\displaystyle R(h,t):={\frac {t^{n}A(h/t)-A(h)}{t^{n}-1}}}where h {\displaystyle h} and h t {\displaystyle {\frac {h}{t}}} are two distinct step sizes.

Then R ( h , t ) = t n ( A ∗ + C ( h t ) n + O ( h n + 1 ) ) − ( A ∗ + C h n + O ( h n + 1 ) ) t n − 1 = A ∗ + O ( h n + 1 ) . {\displaystyle R(h,t)={\frac {t^{n}(A^{*}+C\left({\frac {h}{t}}\right)^{n}+O(h^{n+1}))-(A^{*}+Ch^{n}+O(h^{n+1}))}{t^{n}-1}}=A^{*}+O(h^{n+1}).} R ( h , t ) {\displaystyle R(h,t)} is called the Richardson extrapolation of A(h), and has a higher-order error estimate O ( h n + 1 ) {\displaystyle O(h^{n+1})} compared to A ( h ) {\displaystyle A(h)}.

Very often, it is much easier to obtain a given precision by using R(h) rather than A(h′) with a much smaller h′. Where A(h′) can cause problems due to limited precision (rounding errors) and/or due to the increasing number of calculations needed (see examples below).

Example pseudocode for Richardson extrapolation

The following pseudocode in MATLAB style demonstrates Richardson extrapolation to help solve the ODE y ′ ( t ) = − y 2 {\displaystyle y'(t)=-y^{2}}, y ( 0 ) = 1 {\displaystyle y(0)=1} with the Trapezoidal method. In this example we halve the step size h {\displaystyle h} each iteration and so in the discussion above we'd have that t = 2 {\displaystyle t=2}. The error of the Trapezoidal method can be expressed in terms of odd powers so that the error over multiple steps can be expressed in even powers; this leads us to raise t {\displaystyle t} to the second power and to take powers of 4 = 2 2 = t 2 {\displaystyle 4=2^{2}=t^{2}} in the pseudocode. We want to find the value of y ( 5 ) {\displaystyle y(5)}, which has the exact solution of 1 5 + 1 = 1 6 = 0.1666... {\displaystyle {\frac {1}{5+1}}={\frac {1}{6}}=0.1666...} since the exact solution of the ODE is y ( t ) = 1 1 + t {\displaystyle y(t)={\frac {1}{1+t}}}. This pseudocode assumes that a function called Trapezoidal(f, tStart, tEnd, h, y0) exists which attempts to compute y(tEnd) by performing the trapezoidal method on the function f, with starting point y0 and tStart and step size h.

Note that starting with too small an initial step size can potentially introduce error into the final solution. Although there are methods designed to help pick the best initial step size, one option is to start with a large step size and then to allow the Richardson extrapolation to reduce the step size each iteration until the error reaches the desired tolerance.

External links

Feldman, Joel (2000). (PDF).
Israel, Robert. .
for generic Richardson extrapolation.
for generic Richardson extrapolation.