Coding: Richardson Extrapolation in Mathematica

When we have a complicated algebraic procedure like Richardson extrapolation, it’s tempting to todo the algebra by computer (instead of by hand) in the hopes of avoiding calculation errors. This notebook is an exercise in computer algebra which shows you how to work out Richardson extrapolation formulae in general using Mathematica. 

Credits: This started as a homework project by 2010 student Nona Dowling, was continued for the Fall 2013 class, and was rewritten in Fall 2017.

Theory.

The basic idea of Richardson extrapolation is that we start with an estimator for a quantity of interest whose power series contains only even powers of h. (Note: This is not unusual! Any estimator which is symmetric in h and -h has this property.) We’ll call this

ϕ(h) = L + (∑_(k=1)^∞)(Subscript[a, 2k])(h^(2k)).

We want to improve this estimate by taking linear combinations of ϕ(h), ϕ((h/2)), ϕ((h/4)) and so forth. 

Definition. We define D(n,0) := ϕ((h/(2^n))) = L + (∑_(k=1)^∞)A(k,0) ((((h/(2^n))))^(2k)), where the A(k,0) = (Subscript[a, 2k]).

We’ll call these the “level 0” estimators. Of course, each of these still has error O((h^2)). 

We will now recursively define higher level estimators D(n,m), and prove that they have better error.

Definition. We define D(n,m) recursively by the relation

	D(n,m) = ((4^m)/((4^m)-1 ))D(n,m-1) - (1/((4^m)-1))D(n-1,m-1).
	
Proposition. D(n,m) = L + (∑_(k=m+1)^∞)A(k,m) ((((h/(2^n))))^(2k)).

Proof. By induction on m. The base case (m=0) follows immediately from the definition. Suppose it holds for m-1. Then we should have

D(n,m) = ((4^m)/((4^m)-1 ))(L+(∑_(k=m)^∞)A(k,m-1) ((((h/(2^n))))^(2k))) - (1/((4^m)-1))(L+ (∑_(k=m)^∞)A(k,m-1) ((((h/(2^(n-1)))))^(2k)))

On the right hand side, we can split things into two parts:

(multiples of L): ((4^m)/((4^m)-1 ))L - (1/((4^m)-1))L = L.

(the (h^(2k)) term): A(k,m) ((((h/(2^n))))^(2k))= ((4^m)/((4^m)-1 ))A(k,m-1) ((((h/(2^n))))^(2k)) - (1/((4^m)-1))A(k,m-1) ((((h/(2^(n-1)))))^(2k)) 
				= A(k,m-1) ((((h/(2^n))))^(2k))(((4^m)/((4^m)-1 ))-(1/((4^m)-1))× (1/(2^(-2k)))) .

Simplifying the right hand side further, we see that 

A(k,m) = ((((4^m)-(4^k))/((4^m)-1)))A(k,m-1).

Notice that A(k,k) = 0 (the numerator is zero), so we might as well start the sum at k=m+1. ▮

Implementation.

The best way to define something recursively in Mathematica is to use transformation rules. 
However, you have to be careful: the recursive formula above is not supposed to apply unless m > 0. We can make sure Mathematica knows this using a conditional pattern with the /; (such that) operator.

We can define a “D”[n,m] (the quotes are to prevent Mathematica from confusing this with the built-in D function, which differentiates) by recursion as follows:

RecursionRule := "D"[n_,m_ /; m > 0] :> (((4^m)/((4^m)-1)))"D"[n,m-1] - ((1/((4^m)-1))) "D"[n-1,m-1];
BaseRule := "D"[n_,0] :> ϕ[(h/(2^n))];

We can now apply these rules to define the higher order Richardson extrapolations. Notice that we’re using the “Row” construct to display the expressions in a nice format.

Richardson[n_] := "D"[n,n] //. {RecursionRule,BaseRule};

Manipulate[
Row[{"D"[n,n], "=", Richardson[n]}," "],{n,1,5,1}]

These are certainly ugly, but they look a lot nicer if we simplify them:

Manipulate[Row[{"D"[n,n], "=", Simplify[Richardson[n]]}," "],{n,1,5,1}]

Richardson estimators for the first derivative.

In the estimates above, we haven’t actually used any particular definition for ϕ[h]. But for the first derivative, we know that 

ϕ(h) = ((f(x+h) - f(x-h))/(2h)) = f’(x) + (∑_(k=1)^∞)(Subscript[a, 2k])(h^(2k)),

as we need to get started. Plugging this in, we see:

Manipulate[
Row[{"D"[n,n], "=", Simplify[(Richardson[n] //. ϕ[h_]:> ((f[x+h]-f[x-h])/(2h)))]}," "],{n,1,5,1}]

We can see that our estimators really are just weighted sums of function values!

We can double check our work by expanding each of the functions as a Taylor series and simplifying. Recall that the Series command is very happy to work with an arbitrary function.

Series[f[x+h],{h,0,4}]

If we now substitute in the Taylor series for each of the f terms, we can see that the Richardson estimators really do reduce to f’[x] plus a term of high order in h.

Manipulate[Row[{"D"[n,n], "=",(Richardson[n] //. {ϕ[h_]:> ((f[x+h]-f[x-h])/(2h)),f[a_] :> Series[f[a],{h,0,2n+3}]})," "}],{n,1,5,1}]

Richardson Extrapolation for the Second Derivative.

To continue this process for the second derivative, we need to start with an estimator for the second derivative which has all even terms in the power series (in particular, it just needs to be symmetric in h and -h). It’s pretty clear at this point that the plan is to write out the Taylor series for f(x+h) and f(x-h) and combine them to eliminate f’(x), then solve for f’’(x).

Mathematica makes this easy!

TaylorSeriesEquations = 
{f[x+h] == Normal[Series[f[x+h],{h,0,4}]],
f[x-h] == Normal[Series[f[x-h],{h,0,4}]]}

Eliminate[TaylorSeriesEquations,{f'[x]}]

Solve[Eliminate[TaylorSeriesEquations,{f'[x]}],f''[x]]

We see that for f’’[x], 

ϕ[h] = ((f(x-h) - 2 f(x) + f(x+h))/(h^2)).

is exactly what we need to get Richardson extrapolation started. The resulting estimators look sort of mysterious:

Manipulate[
Row[{"D"[n,n], "=", Simplify[(Richardson[n] //. ϕ[h_]:> ((f[x-h]-2f[x] + f[x+h])/(h^2)))]}," "],{n,0,5,1}]

But again, you can see that if we substitute in the Taylor expansions for the function and let Mathematica simplify the expression, we do indeed reduce to f’’[x] + some high order term in h.

Manipulate[Row[{"D"[n,n], "=",(Richardson[n] //. {ϕ[h_]:> ((f[x+h]-2f[x]+f[x-h])/(h^2)),f[a_] :> Series[f[a],{h,0,2n+5}]})," "}],{n,1,5,1}]

Richardson Extrapolation for the third derivative?

Challenge: Can you follow along with the ideas above to find a starting point for Richardson extrapolation for the third derivative? You’ll have to evaluate the function at more points, and eliminate both f’[x] and f’’[x] from the resulting system of equations.