2.13 RICHARDSON EXTRAPOLATION

```#
#
# WORKSHEET # 24
#
# RICHARDSON EXTRAPOLATION
#
#
--------------------------------------------------------------------------------
# In this worksheet we describe a quite elementary method, called Romberg integration, for
# improving on the accuracy of numerical methods of integration.   We apply it first to the
# trapezoidal rule.
--------------------------------------------------------------------------------
> restart;
--------------------------------------------------------------------------------
> h:=n->(b-a)/n;

b - a
h := n -> -----
n
--------------------------------------------------------------------------------
> X:=(i,n)->a+i*(b-a)/n;

i (b - a)
X := (i,n) -> a + ---------
n
--------------------------------------------------------------------------------
> Trap:=n->(h(n)/2)*sum(f(X(i,n))+f(X(i+1,n)),i=0..n-1);

/n - 1                              \
|-----                              |
| \                                 |
Trap := n -> 1/2 h(n) |  )   (f(X(i, n)) + f(X(i + 1, n)))|
| /                                 |
|-----                              |
\i = 0                              /
--------------------------------------------------------------------------------
# Writing T[n] for Trap[n] we can prove the following formula, under reasonable
# differentiability  assumptions on the function f(x).  The last term in the formula, that is
# O(h^(2*m+2)), indicates any function of h that is bounded by K*h^(2*m+2) as h tends to
# 0, where K is a constant.  The important point here is that the constant K and the
# coefficients a[k] depend only on f(x), a and b, and not on the step size h.  We refer to the
# trapezoidal rule as an order 2 method since the first term in the error is a[2]*h^2, i.e.
# degree 2 in the step size h.  Doubling n changes the step size from h to h/2 and
# consequently the first term in the error changes from a[2]*h^2 to (1/4)*a[2]*h^2, that is
# the leading term in the error has been quartered.
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=T[n]+Sum(a[k]*h^(2*k),k=1..m)+O(h^(2*m+2));

b                  /  m              \
/                  |-----            |
|                   | \          (2 k)|      (2 m + 2)
|  f(x) dx = T[n] + |  )   a[k] h     | + O(h         )
|                   | /               |
/                    |-----            |
a                    \k = 1            /
--------------------------------------------------------------------------------
# For m=2 this formula becomes
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=T[n]+a[2]*h^2+a[4]*h^4+O(h^6);

b
/
|                         2         4      6
|  f(x) dx = T[n] + a[2] h  + a[4] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
# Now repeat this for step size h/2, that is double the value of n to 2*n.
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=T[2*n]+a[2]*(h/2)^2+a[4]*(h/2)^4+O(h^6);

b
/
|                               2              4      6
|  f(x) dx = T[2 n] + 1/4 a[2] h  + 1/16 a[4] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
# By performing elementary arithmetic on these last 2 equations we can eliminate the h^2
# terms from the error.
--------------------------------------------------------------------------------
> (1/3)*(4*(")-"");

b
/
|                                   4      6
|  f(x) dx = 4/3 T[2 n] - 1/4 a[4] h  + O(h ) - 1/3 T[n]
|
/
a
--------------------------------------------------------------------------------
# That is (4*T[2n]-T[n])/3 is a "better approximation" since the first term in the error now
# has order 4.  This is a "new" rule of numerical integration.
--------------------------------------------------------------------------------
> T[1]:=n->(4*Trap(2*n)-Trap(n))/3;

T[1] := n -> 4/3 Trap(2 n) - 1/3 Trap(n)
--------------------------------------------------------------------------------
> a:=0;b:=1;

a := 0

b := 1
--------------------------------------------------------------------------------
>  f:=x->1/(x+1);\

1
f := x -> -----
x + 1
--------------------------------------------------------------------------------
> J:=int(f(x),x=a..b);\

J := ln(2)
--------------------------------------------------------------------------------
> evalf(ln(2));

.6931471806
--------------------------------------------------------------------------------
> seq([n,evalf(J-Trap(n)),evalf(J-T[1](n))],n=1..5);n:='n':

[1, -.0568528194, -.0012972638], [2, -.0151861527, -.0001067877],

-5
[3, -.0068528194, -.0000226126], [4, -.0038766289, -.73501*10  ],

-5
[5, -.0024877400, -.30501*10  ]
--------------------------------------------------------------------------------
# The error terms J-T[1](n) are a lot smaller than the error terms J-Trap(n).
--------------------------------------------------------------------------------
# Let's compare T[1] with simpson's rule.
--------------------------------------------------------------------------------
> Simpson:=n->(h(n)/3)*sum(f(X(2*i,n))+4*f(X(2*i+1,n))+f(X(2*i+2,n)),i=0..n/2-1);

Simpson :=

/1/2 n - 1                                                       \
|  -----                                                         |
|   \                                                            |
n -> 1/3 h(n) |    )     (f(X(2 i, n)) + 4 f(X(2 i + 1, n)) + f(X(2 i + 2, n)))|
|   /                                                            |
|  -----                                                         |
\  i = 0                                                         /
--------------------------------------------------------------------------------
> seq([evalf(J-T[1](n)),evalf(J-Simpson(2*n))],n=1..5);n:='n':

[-.0012972638, -.0012972638], [-.0001067877, -.0001067877],

-5            -5
[-.0000226126, -.0000226126], [-.73501*10  , -.73501*10  ],

-5            -5
[-.30501*10  , -.30501*10  ]
--------------------------------------------------------------------------------
# These numbers are the same, and in fact T[1](n)=Simpson(2*n).
--------------------------------------------------------------------------------
> f:='f':
--------------------------------------------------------------------------------
> T[1](1);

1/6 f(0) + 2/3 f(1/2) + 1/6 f(1)
--------------------------------------------------------------------------------
> Simpson(2);

1/6 f(0) + 2/3 f(1/2) + 1/6 f(1)
--------------------------------------------------------------------------------
> for k from 1 to 5 do print(T[1](k)-Simpson(2*k)) od;k:='k':

0

0

0

0

0
--------------------------------------------------------------------------------
> a:='a':b:='b':
--------------------------------------------------------------------------------
# We can use the trapezoidal rules T[n], T[2*n] to estimate the error made in
# approximating Int(f(x),x=a..b) by T[n] or T[2*n].
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=T[n]+a[2]*h^2+a[4]*h^4+O(h^6);

b
/
|                         2         4      6
|  f(x) dx = T[n] + a[2] h  + a[4] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=T[2*n]+a[2]*(h/2)^2+a[4]*(h/2)^4+O(h^6);

b
/
|                               2              4      6
|  f(x) dx = T[2 n] + 1/4 a[2] h  + 1/16 a[4] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
> ""-";

2    15        4
0 = T[n] + 3/4 a[2] h  + ---- a[4] h  - T[2 n]
16
--------------------------------------------------------------------------------
# There is a mistake here, due to the ambiguity in the big O notation  The big O terms are
# not the same and therefore should not be cancelled.  If we only consider the terms out to
# h^4 then the correct formula should read:
--------------------------------------------------------------------------------
> 0=T[n]-T[2*n]+(3/4)*a[2]*h^2+O(h^4);

2      4
0 = T[n] - T[2 n] + 3/4 a[2] h  + O(h )
--------------------------------------------------------------------------------
> (1/4)*a[2]*h^2=(T[2*n]-T[n])/3+O(h^4);

2                              4
1/4 a[2] h  = 1/3 T[2 n] - 1/3 T[n] + O(h )
--------------------------------------------------------------------------------
# Now recall that the lowest order term in the error for T[2*n] is (1/4)*a[2]*h^2.  In other
# words (T[2*n]-T[n])/3 gives a good indication of the error made in T[2*n].  In a similar
# fashion we see that 4*(T[2*n]-T[n])/3 is a good indication of the error made in T[n].
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)-T[2*n]=(T[2*n]-T[n])/3+O(h^4);

b
/
|                                                4
|  f(x) dx - T[2 n] = 1/3 T[2 n] - 1/3 T[n] + O(h )
|
/
a
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)-T[n]=4*(T[2*n]-T[n])/3+O(h^4);

b
/
|                                              4
|  f(x) dx - T[n] = 4/3 T[2 n] - 4/3 T[n] + O(h )
|
/
a
--------------------------------------------------------------------------------
# We can also "solve" for the coefficient a[2].
--------------------------------------------------------------------------------
> a[2]=4*(T[2*n]-T[n])/(3*h^2)+O(h^2);

T[2 n] - T[n]      2
a[2] = 4/3 ------------- + O(h )
2
h
--------------------------------------------------------------------------------
# If we do several calculations and see that the values of 4*(T[2*n]-T[n])/(3*h^2) are
# stabilizing to a certain value then we may reasonably say that a[2] is this value.  Once we
# have an idea of the value of a[2] we can control the first term in the error of the
# trapezoidal rule.
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=T[n]+a[2]*h^2+O(h^4);

b
/
|                         2      4
|  f(x) dx = T[n] + a[2] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
# For example, suppose we want to use the trapezoidal rule to estimate Int(f(x),x=a..b) to
# within 10^(-8).  We do several calculations of the trapezoidal rule first to get an idea of
# the value of a[2], and then we choose the step size h so that abs(a[2]*h^2)<10^(-8).
# Assuming that the higher order terms O(h^4) can be ignored this step size should be
# sufficient.
--------------------------------------------------------------------------------
# We can apply this idea of extrapolating rules to Simpson's rule.    Let S[n] denote
# Simpson(n).  Recall that n must be even for Simpson's rule.
--------------------------------------------------------------------------------
> Simpson:=n->(h(n)/3)*sum(f(X(2*i,n))+4*f(X(2*i+1,n))+f(X(2*i+2,n)),i=0..n/2-1);

Simpson :=

/1/2 n - 1                                                       \
|  -----                                                         |
|   \                                                            |
n -> 1/3 h(n) |    )     (f(X(2 i, n)) + 4 f(X(2 i + 1, n)) + f(X(2 i + 2, n)))|
|   /                                                            |
|  -----                                                         |
\  i = 0                                                         /
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=S[n]+b[4]*h^4+ O(h^6);

b
/
|                         4      6
|  f(x) dx = S[n] + b[4] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=S[2*n]+(1/16)*b[4]*h^4+ O(h^6);

b
/
|                                4      6
|  f(x) dx = S[2 n] + 1/16 b[4] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
# Again we do elementary arithmetic on these last 2 formulas.
--------------------------------------------------------------------------------
> (16*"-"")/15;

b
/
|             16              6
|  f(x) dx = ---- S[2 n] + O(h ) - 1/15 S[n]
|             15
/
a
--------------------------------------------------------------------------------
# Thus we have the "improved" approximation, improved in the sense that the first term in
# the error is now of order 6.
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=(16*S[2*n]-S[n])/15;

b
/
|             16
|  f(x) dx = ---- S[2 n] - 1/15 S[n]
|             15
/
a
--------------------------------------------------------------------------------
# Of course we can also solve for b[4].
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=S[n]+b[4]*h^4+ O(h^6);

b
/
|                         4      6
|  f(x) dx = S[n] + b[4] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
> Int(f(x),x=a..b)=S[2*n]+(1/16)*b[4]*h^4+ O(h^6);

b
/
|                                4      6
|  f(x) dx = S[2 n] + 1/16 b[4] h  + O(h )
|
/
a
--------------------------------------------------------------------------------
> "-"";

15        4
0 = S[2 n] - ---- b[4] h  - S[n]
16
--------------------------------------------------------------------------------
# We have to add the big O term in.
--------------------------------------------------------------------------------
> b[4]*h^4=(16/15)*(S[2*n]-S[n])+O(h^6);

4    16            16            6
b[4] h  = ---- S[2 n] - ---- S[n] + O(h )
15            15
--------------------------------------------------------------------------------
# This gives us the first order term in the error for S[n].
--------------------------------------------------------------------------------
> b[4]=(16/15)*(S[2*n]-S[n])/(h^4)+O(h^6);

16  S[2 n] - S[n]      6
b[4] = ---- ------------- + O(h )
15         4
h
--------------------------------------------------------------------------------
> Simpson[1]:=n->(16*Simpson(2*n)-Simpson(n))/15;

Simpson[1] := n -> 16/15 Simpson(2 n) - 1/15 Simpson(n)
--------------------------------------------------------------------------------
> a:=0:b:=1:n:='n':
--------------------------------------------------------------------------------

```