Editing Richardson extrapolation (section)

==Example pseudocode for Richardson extrapolation==

The following pseudocode in MATLAB style demonstrates Richardson extrapolation to help solve the ODE <math>y'(t) = -y^2</math>, <math>y(0) = 1</math> with the [[Trapezoidal method]]. In this example we halve the step size <math>h</math> each iteration and so in the discussion above we'd have that <math>t = 2</math>. 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 <math>t</math> to the second power and to take powers of <math>4 = 2^2 = t^2</math> in the pseudocode. We want to find the value of <math>y(5)</math>, which has the exact solution of <math>\frac{1}{5 + 1} = \frac{1}{6} = 0.1666...</math> since the exact solution of the ODE is <math>y(t) = \frac{1}{1 + t}</math>. This pseudocode assumes that a function called <code>Trapezoidal(f, tStart, tEnd, h, y0)</code> exists which attempts to compute <code>y(tEnd)</code> by performing the trapezoidal method on the function <code>f</code>, with starting point <code>y0</code> and <code>tStart</code> and step size <code>h</code>.

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.

<syntaxhighlight lang="matlab">
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

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

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.
% The first row of this (lower triangular) matrix has now been filled.
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
    % Halve the previous value of h since this is the start of a new row.
    h = h/2

    % Starting filling row i+1 from the left by calling
    % the Trapezoidal function with this new smaller step size
    A(i + 1, 1) = Trapezoidal(f, tStart, tEnd, h, y0)

    % Go across this current (i+1)-th row until the diagonal is reached
    for j = 1 : i
        % To compute A(i + 1, j + 1), which is the next Richardson extrapolate, 
        % use the most recently computed value (i.e. A(i + 1, j))
        % and the value from the row above it (i.e. A(i, j)).

        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
        % Display the result of the Richardson extrapolation
        print("y = ", A(i + 1, i + 1))
        haveWeFoundSolution = true
        % Done, so leave the loop
        break
    end
end

% If we were not able to find a solution to within the desired tolerance
if (not haveWeFoundSolution)
    print("Warning: Not able to find solution to within the desired tolerance of ", tolerance);
    print("The last computed extrapolate was ", A(maxRows, maxRows))
end
</syntaxhighlight>