# 2.14 ROMBERG INTEGRATION

```#
# 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   ]
--------------------------------------------------------------------------------

```