Example

Let us derive the centered three-point formula again, keeping explicitly terms up to order :

> restart;

> tayl:=convert(taylor(f(x),x=x0,9),polynom);

The error term is given by

> err:=(D@@9)(f)(xi(x))*(x-x0)^9/9!;

so that,

> f(x)=tayl+err;

We can then evaluate this expansion at the two neighbouring points:

> expr1:=f(xo+h)=subs(x=x0+h,tayl+err);

> expr2:=f(x0-h)=subs(x=x0-h,tayl+err);

We can eliminate the second derivative terms by subtracting both expressions :

> lhs(expr1)-lhs(expr2)=rhs(expr1)-rhs(expr2);

which gives the centered three-point difference formula. When is continuous over then the Intermediate Value Theorem says that

,

with in the interval .

Therefore, the three-point difference formula becomes:

> expr:=D(f)(x[0])=(f(x[0]+h)-f(x[0]-h))/(2*h)+C[1]*h^2+C[2]*h^4+c[3]*h^6 + O(h^8);

So we should be able to apply Richardson Extrapolation up to :

> restart:f:=x->ln(cos(x));

> df:=D(f):

> h:=0.8:df1[1]:=1/(2*h)*(f(0.4+h)-f(0.4-h)):acterr:=abs(df1[1]-df(0.4)): print(h,df1[1],acterr);

> for i from 2 to 8 do h:=0.8/(2^(i-1)): df1[i]:=1/(2*h)*(f(0.4+h)-f(0.4-h)): N[2,i]:=(4*df1[i]-df1[i-1])/3: newerr:=abs(N[2,i]-df(0.4)): print(h,df1[i],N[2,i],newerr);od:

The next sequence generated by Richardson extrapolation is

> for i from 3 to 8 do h:=0.8/(2^(i-1)):N[3,i]:=(16*N[2,i]-N[2,i-1])/15: newerr:=abs(N[3,i]-df(0.4)): print(h,N[3,i],newerr);od:

and finally:

> for i from 4 to 8 do h:=0.8/(2^(i-1)):N[4,i]:=(64*N[3,i]-N[3,i-1])/63: newerr:=abs(N[4,i]-df(0.4)): print(h,N[4,i],newerr);od:

Always take into account, as illustrated in the numbers above, that round-off errors will be present.