#
#
# 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':
--------------------------------------------------------------------------------
