#
# WORKSHEET # 25
#
# ROMBERG INTEGRATION
#
#
--------------------------------------------------------------------------------
# In this worksheet we will explore the idea of Romberg integration. We start with the
# trapezoidal rule and apply a sequence of Richardson extrapolations.
--------------------------------------------------------------------------------
> 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 /
--------------------------------------------------------------------------------
# Let T[n] denote Trap(n). Then we know that the error term in T[n] has the following
# form, where the constants a[k] are independent of the step size h and the constant in the
# big O term is also independent of h.
--------------------------------------------------------------------------------
> 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 /
--------------------------------------------------------------------------------
# Richardson extrapolation cancels the first term in the error by taking a suitable linear
# combination of T[n] and T[2*n]. We will apply a sequence of such extrapolations to the
# example Int(x/(x^2+1/10),x=0..1). Specifically we will consider the ttrapezoidal
# approximations T[2^j] for j from 1 to 5. Thus we start with step size 1/2 and end with
# step size 1/32.
--------------------------------------------------------------------------------
> a:=0;b:=1;
a := 0
b := 1
--------------------------------------------------------------------------------
> J:=f->int(f(x),x=a..b);
b
/
|
J := f -> | f(x) dx
|
/
a
--------------------------------------------------------------------------------
> f:=x->x/(x^2+1/10);
x
f := x -> ---------
2
x + 1/10
--------------------------------------------------------------------------------
> J(f);
1/2 ln(11)
--------------------------------------------------------------------------------
> evalf(");
1.198947637
--------------------------------------------------------------------------------
> for j from 1 to 5 do R[j,1]:=evalf(Trap(2^j)) od;j:='j':
R[1, 1] := .9415584416
R[2, 1] := 1.138413473
R[3, 1] := 1.184736526
R[4, 1] := 1.195437378
R[5, 1] := 1.198072507
--------------------------------------------------------------------------------
# These are the trapezoidal approximations. In each case the error behaves like a constant
# times the square of the step size, i.e. a[2]*(1/2^j)^2. In the next calculation we do
# repeated Richardson extrapolation. For R[j,k] we expect the error to behave like a
# constant times the (2*k)^{th} power of the step size, i.e. const*(1/2^j)^(2*k).
--------------------------------------------------------------------------------
> for k from 2 to 5 do \
for j from k to 5 do \
R[j,k]:=(4^(k-1)*R[j,k-1]-R[j-1,k-1])/(4^(k-1)-1) od od;j:='j':k:='k':
> with(linalg):
Warning: new definition for norm
Warning: new definition for trace
--------------------------------------------------------------------------------
> for k from 2 to 5 do for j from 1 to k-1 do R[j,k]:=` ` od od; j:='j':k:='k':
--------------------------------------------------------------------------------
# The command R[j,k]=` ` uses the backquote key and means that the R[j.k] entries for
# j=1..k-1 are empty. For example the matrices below are upper triangular.
--------------------------------------------------------------------------------
> R[1,2];
--------------------------------------------------------------------------------
> R[1,1];
.9415584416
--------------------------------------------------------------------------------
> R[2,1];
1.138413473
--------------------------------------------------------------------------------
> R[2,2];
1.204031817
--------------------------------------------------------------------------------
> matrix([seq([seq(R[j,k],j=2..5)],k=1..5)]);
[ 1.138413473 1.184736526 1.195437378 1.198072507 ]
[ ]
[ 1.204031817 1.200177544 1.199004329 1.198950883 ]
[ ]
[ 1.199920592 1.198926115 1.198947320 ]
[ ]
[ 1.198910329 1.198947656 ]
[ ]
[ 1.198947802 ]
--------------------------------------------------------------------------------
# The entries in this matrix are produced by repeated extrapolation. The first row
# corresponds to k=1 the second to k=2, and so on. The first column corresponds to j=2, the
# second to j=3, etc. To produce a particular entry in the second row we look at the entry
# immediately above it, multiply it by 4, subtract the entry to the left, and then divide by 3.
# The procedure for the third row is analogous except that we now multiply by 16 and
# divide by 15.
--------------------------------------------------------------------------------
# Richardson extrapolation also allows us to estimate the error made in approximating J by
# R[j,k]. The approximate error is E[j,k]=(R[j,k]-R[j-1,k])/(4^k-1). In fact this is just the
# first term in the actual error.
--------------------------------------------------------------------------------
> for k from 1 to 4 do for j from k+1 to 5 do
> Er[j,k]:=(R[j,k]-R[j-1,k])/(4^k-1) od od;j:='j':k:='k':
--------------------------------------------------------------------------------
> for k from 2 to 4 do for j from 1 to k do\
Er[j,k]:=` ` od od;j:='j':k:='k':
> matrix([seq([seq(Er[j,k],j=2..5)],k=1..4)]);
[ .06561834379 .01544101767 .003566950666 .0008783763332 ]
[ ]
[ -5 ]
[ -.0002569515333 -.00007821433334 -.3563066667*10 ]
[ ]
[ -6 ]
[ -.00001578534920 .3365873015*10 ]
[ ]
[ -6 ]
[ .1463803921*10 ]
--------------------------------------------------------------------------------
# For comparison here are the true errors.
--------------------------------------------------------------------------------
> for k from 1 to 5 do for j from k to 5 do Ertrue[j,k]:=evalf(J(f))-R[j,k] od od;j:='j':k:='k':\
for k from 2 to 5 do for j from 1 to k-1 do\
Ertrue[j,k]:= `` od od; j:='j':k:='k':
--------------------------------------------------------------------------------
> matrix([seq([seq(Ertrue[j,k],j=2..5)],k=1..5)]);
[ .060534164 .014211111 .003510259 .000875130 ]
[ ]
[ -5 ]
[ -.005084180 -.001229907 -.000056692 -.3246*10 ]
[ ]
[ -6 ]
[ -.000972955 .000021522 .317*10 ]
[ ]
[ -7 ]
[ .000037308 -.19*10 ]
[ ]
[ -6 ]
[ -.165*10 ]
--------------------------------------------------------------------------------
