Richardson extrapolation

Suppose that we wish to approximate ${\displaystyle A^{*}}$, and we have a method ${\displaystyle A(h)}$ that depends on a small parameter ${\displaystyle h}$ in such a way that

${\displaystyle A(h)=A^{\ast }+Ch^{n}+O(h^{n+1}).}$

Let us define a new function

${\displaystyle R(h,t):={\frac {t^{n}A(h/t)-A(h)}{t^{n}-1}}}$ where ${\displaystyle h}$ and ${\displaystyle {\frac {h}{t}}}$ are two distinct step sizes.

Then

${\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}).}$

${\displaystyle R(h,t)}$ is called the Richardson extrapolation of A(h), and has a higher-order error estimate ${\displaystyle O(h^{n+1})}$ compared to ${\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).

Let ${\displaystyle A(h)}$ be an approximation of ${\displaystyle A^{*}}$(exact value) that depends on a positive step size h with an error formula of the form

${\displaystyle A^{*}=A(h)+a_{0}h^{k_{0}}+a_{1}h^{k_{1}}+a_{2}h^{k_{2}}+\cdots }$

where the a_i are unknown constants and the k_i are known constants such that h^k_i > h^k_i+1.

k₀ is the leading order step size behavior of Truncation error as ${\displaystyle A^{*}=A(h)+O(h^{k_{0}})}$

The exact value sought can be given by

${\displaystyle A^{*}=A(h)+a_{0}h^{k_{0}}+a_{1}h^{k_{1}}+a_{2}h^{k_{2}}+\cdots }$

which can be simplified with Big O notation to be

${\displaystyle A^{*}=A(h)+a_{0}h^{k_{0}}+O(h^{k_{1}}).\,\!}$

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

${\displaystyle A^{*}=A(h)+a_{0}h^{k_{0}}+O(h^{k_{1}})\,\!}$

${\displaystyle A^{*}=A\!\left({\frac {h}{t}}\right)+a_{0}\left({\frac {h}{t}}\right)^{k_{0}}+O(h^{k_{1}}).}$

Multiplying the second equation by t^k₀ and subtracting the first equation gives

${\displaystyle (t^{k_{0}}-1)A^{*}=t^{k_{0}}A\left({\frac {h}{t}}\right)-A(h)+(t^{k_{0}}-1)O(h^{k_{1}})}$

which can be solved for ${\displaystyle A^{*}}$ to give

${\displaystyle A^{*}={\frac {t^{k_{0}}A\left({\frac {h}{t}}\right)-A(h)}{t^{k_{0}}-1}}+O(h^{k_{1}}).}$

Therefore, using ${\displaystyle R(h,t)={\frac {t^{k_{0}}A\left({\frac {h}{t}}\right)-A(h)}{t^{k_{0}}-1}}}$ the truncation error has been reduced to ${\displaystyle O(h^{k_{1}})}$. This is in contrast to ${\displaystyle A(h)}$ where the truncation error is ${\displaystyle O(h^{k_{0}})}$ for the same step size ${\displaystyle h}$

By this process, we have achieved a better approximation of A by subtracting the largest term in the error which was O(h^k₀). This process can be repeated to remove more error terms to get even better approximations.

A general recurrence relation beginning with ${\displaystyle A_{0}=A(h)}$ can be defined for the approximations by

${\displaystyle A_{i+1}(h)={\frac {t^{k_{i}}A_{i}\left({\frac {h}{t}}\right)-A_{i}(h)}{t^{k_{i}}-1}}}$

where ${\displaystyle k_{i+1}}$ satisfies

${\displaystyle A^{*}=A_{i+1}(h)+O(h^{k_{i+1}})}$.

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

Additionally, the general formula can be used to estimate k₀ (leading order step size behavior of Truncation error) when neither its value nor A^* (exact value) is known a priori. Such a technique can be useful for quantifying an unknown rate of convergence. Given approximations of A from three distinct step sizes h, h / t, and h / s, the exact relationship

${\displaystyle A^{*}={\frac {t^{k_{0}}A\left({\frac {h}{t}}\right)-A(h)}{t^{k_{0}}-1}}+O(h^{k_{1}})={\frac {s^{k_{0}}A\left({\frac {h}{s}}\right)-A(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 approximately valid)

${\displaystyle A\left({\frac {h}{t}}\right)+{\frac {A\left({\frac {h}{t}}\right)-A(h)}{t^{k_{0}}-1}}\approx A\left({\frac {h}{s}}\right)+{\frac {A\left({\frac {h}{s}}\right)-A(h)}{s^{k_{0}}-1}}}$

which can be solved numerically to estimate k₀ for some arbitrary choices of h, s and t.

The following pseudocode in MATLAB style demonstrates Richardson extrapolation to help solve the ODE ${\displaystyle y'(t)=-y^{2}}$, ${\displaystyle y(0)=1}$ with the Trapezoidal method. In this example we halve the step size ${\displaystyle h}$ each iteration and so in the discussion above we'd have that ${\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 ${\displaystyle t}$ to the second power and to take powers of ${\displaystyle 4=2^{2}=t^{2}}$ in the pseudocode. We want to find the value of ${\displaystyle y(5)}$, which has the exact solution of ${\displaystyle {\frac {1}{5+1}}={\frac {1}{6}}=0.1666...}$ since the exact solution of the ODE is ${\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.

tStart = 0          % Starting time
tEnd = 5            % Ending time
f = -y^2            % The derivative of y, so y' = f(t, y(t)) = -y^2
                    % The solution to this ODE is y = 1/(1 + t)
y0 = 1              % The initial position (i.e. y0 = y(tStart) = y(0) = 1)
tolerance = 10^-11  % 10 digit accuracy is desired

maxRows = 20                % Don't allow the iteration to continue indefinitely
initialH = tStart - tEnd    % Pick an initial step size
haveWeFoundSolution = false % Were we able to find the solution to within the desired tolerance? not yet.

h = initialH

% Create a 2D matrix of size maxRows by maxRows to hold the Richardson extrapolates
% Note that this will be a lower triangular matrix and that at most two rows are actually
% needed at any time in the computation.
A = zeroMatrix(maxRows, maxRows)

%Compute the top left element of the matrix
A(1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

% Each row of the matrix requires one call to Trapezoidal
% This loops starts by filling the second row of the matrix, since the first row was computed above
for i = 1 : maxRows - 1 % Starting at i = 1, iterate at most maxRows - 1 times
    h = h/2             % Half the previous value of h since this is the start of a new row
    
    % Call the Trapezoidal function with this new smaller step size
    A(i + 1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)
    
    for j = 1 : i         % Go across the row until the diagonal is reached
        % Use the value just computed (i.e. A(i + 1, j)) and the element from the
        % row above it (i.e. A(i, j)) to compute the next Richardson extrapolate
     
        A(i + 1, j + 1) = ((4^j).*A(i + 1, j) - A(i, j))/(4^j - 1);
    end
    
    % After leaving the above inner loop, the diagonal element of row i + 1 has been computed
    % This diagonal element is the latest Richardson extrapolate to be computed.
    % The difference between this extrapolate and the last extrapolate of row i is a good
    % indication of the error.
    if (absoluteValue(A(i + 1, i + 1) - A(i, i)) < tolerance)  % If the result is within tolerance
        print("y(5) = ", A(i + 1, i + 1))                      % Display the result of the Richardson extrapolation
        haveWeFoundSolution = true
        break                                                  % Done, so leave the loop
    end
end

if (haveWeFoundSolution == false)   % If we weren't able to find a solution to within the desired tolerance
    print("Warning: Not able to find solution to within the desired tolerance of ", tolerance);
    print("The last computed extrapolate was ", A(maxRows, maxRows))
end

• Aitken's delta-squared process
• Takebe Kenko
• Richardson iteration

• Extrapolation Methods. Theory and Practice by C. Brezinski and M. Redivo Zaglia, North-Holland, 1991.
• Ivan Dimov, Zahari Zlatev, Istvan Farago, Agnes Havasi: Richardson Extrapolation: Practical Aspects and Applications, Walter de Gruyter GmbH & Co KG, ISBN 9783110533002 (2017).

• Fundamental Methods of Numerical Extrapolation With Applications, mit.edu
• Richardson-Extrapolation
• Richardson extrapolation on a website of Robert Israel (University of British Columbia)
• Matlab code for generic Richardson extrapolation.
• Julia code for generic Richardson extrapolation.