# Lesson 19. Richardson Extrapolation and Romberg Integration
# ------------------------------------------------------------
#
# Last time we defined the Newton-Cotes rules for numerical
# integration of int(f(x), x=a..b).
#
> restart:
> h:= n ->(b-a)/n:
> X:= (i,n) -> a + i*(b-a)/n:
> a:= 0: b:= 1:
> J:= f -> int(f(x), x=a..b):
> rule:= (n,k) ->
> h(n) * sum(sum(d[j]*f(X(k*i+j,n)), j=0..k), i=0..n/k-1):
> eqns:= k -> eval({ seq(subs(f=unapply(x^ii,x), rule(k,k)=J(f)), ii=0..k) }):
Warning, `ii` is implicitly declared local
> for kk from 2 to 12 do
> ncrule[kk]:= unapply(subs(solve(eqns(kk)), rule(n,kk)), n) od:
# We saw that the Newton-Cotes rule ncrule[k](n) gives the
# correct answer when f is a polynomial of degree q = k (if k is
# odd) or k+1 (if k is even). The error is a constant times
# (D@@(q+1))(f)(t)/n^(q+1).
#
# If we compare a Newton-Cotes rule with 1/n^p in the error to
# one with 1/n^q in the error, where q > p, the second one is
# more accurate when n is large. However, as we saw last time
# this is not necessarily true for any given n, especially when
# the higher order derivatives of f are large.
> f:= x -> 1/(x^2 + 1/20):
> map(k -> evalf(sqrt(20)*arctan(sqrt(20))-ncrule[k](36)), [2,3,4,6,9,12]);
-7 -6 -5 -5 -5 -6
[.62*10 , .223*10 , -.1006*10 , .3563*10 , .1121*10 , -.503*10 ]
# Here's an even more striking example. This time we take n=k.
# For each n we find what would appear to be the best possible
# approximation to the integral, given the values at the points
# x[0],...,x[n] - an approximation that is exact for polynomials
# of the greatest possible degree. But in this case it turns
# out that the errors actually get worse as n increases.
> f:= x -> 1/((8*x-4)^2+1):
> map(k -> evalf(J(f)-ncrule[k](k)), [2,4,6,8,10,12]);
[-.3548200938, .0467485336, -.0846453499, .0888176280, -.1179906340, .1646286874
]
# Richardson extrapolation is a very widely applicable idea,
# which can often be used to improve the results of a numerical
# method.
# Suppose you want to approximate a quantity J, and you have
# available approximations Q(n). Typically, you know that this
# approximation is of a certain order, say J = Q(n) + O(1/n^p).
# But often, more can be said:
> eq1:= J = Q(n) + A/n^p + O(1/n^(p+1)):\
#
# The idea of Richardson extrapolation is to take two different
# values of n, and eliminate the A term. Typically, we take n
# and n/2. Now
> eq2:= J = Q(n/2) + 2^p*A/n^p + O(1/n^(p+1)):
# We can take a combination of this equation and the last one
# and eliminate the A/n^p term which is the main contribution
# to the error.
> solve({eq1,eq2},{J,A});
p 1 1 p
Q(n) 2 - Q(1/2 n) - O(--------) + O(--------) 2
(p + 1) (p + 1)
n n
{J = -------------------------------------------------,
p
- 1 + 2
p
n (- Q(n) + Q(1/2 n))
A = - ----------------------}
p
- 1 + 2
# Thus (ignoring the higher-order terms)
> R:=subs(O(1/n^(p+1))=0, subs(",J));
p
Q(n) 2 - Q(1/2 n)
R := ------------------
p
- 1 + 2
# should be a better approximation to the true value. Probably
# more important is the fact that
> subs("", A/n^p);
- Q(n) + Q(1/2 n)
- -----------------
p
- 1 + 2
# is an approximation for the error in Q(n). The error in our
# improved approximation R is not known (though it should be
# O(1/n^(p+1))), but to be conservative we can take the value
# above as an indication of its error as well.
#
# Suppose we apply this idea to the Trapezoid Rule.
> traprule:= n ->\
h(n) * sum((f(X(i,n))+f(X(i+1,n)))/2, i=0..n-1):
# It turns out that the error in traprule(n) for smooth functions
# f is of the form
# sum(A[2*k]/n^(2*k), k=1..N) + O(1/n^(2N+1))
# The leading term in the error is A[2]/n^2. Richardson
# extrapolation removes this:
> T[1]:= n -> (4*traprule(n)-traprule(n/2))/3;
T[1] := n -> 4/3 traprule(n) - 1/3 traprule(1/2 n)
> f:= 'f': T[1](2);
1/6 f(0) + 2/3 f(1/2) + 1/6 f(1)
# Look familiar? This is the Simpson's Rule for n=2. In fact,
# T1 is just the same as Simpson's Rule. The error in T1 should
# be O(1/n^4), which is true for Simpson's Rule. What is new
# is the error estimate
> Er[1]:= n -> (traprule(n)-traprule(n/2))/3;
Er[1] := n -> 1/3 traprule(n) - 1/3 traprule(1/2 n)
# But there's no reason to stop here. We can apply Richardson
# Extrapolation to Simpson's Rule as well.
> T[2]:= n -> (16*T[1](n)-T[1](n/2))/15;
T[2] := n -> 16/15 T[1](n) - 1/15 T[1](1/2 n)
> T[2](4);
16 16
7/90 f(0) + ---- f(1/4) + 2/15 f(1/2) + ---- f(3/4) + 7/90 f(1)
45 45
# This turns out to be the same as the Newton-Cotes rule
# ncrule[4]. But after this the Richardson Extrapolation results
# are not Newton-Cotes rules.
> T[3]:= n -> (64*T[2](n)-T[2](n/2))/63;
T[3] := n -> 64/63 T[2](n) - 1/63 T[2](1/2 n)
> T[3](8);
31 512 176 512 218 512
--- f(0) + ---- f(1/8) + ---- f(1/4) + ---- f(3/8) + ---- f(1/2) + ---- f(5/8)
810 2835 2835 2835 2835 2835
176 512 31
+ ---- f(3/4) + ---- f(7/8) + --- f(1)
2835 2835 810
> ncrule[8](8);
989 2944 464 5248 454
----- f(0) + ----- f(1/8) - ----- f(1/4) + ----- f(3/8) - ---- f(1/2)
28350 14175 14175 14175 2835
5248 464 2944 989
+ ----- f(5/8) - ----- f(3/4) + ----- f(7/8) + ----- f(1)
14175 14175 14175 28350
# T[3] has error O(1/n^8) while ncrule[8] has error O(1/n^10), so
# ncrule[8] should be better for large n. Let's try the last
# example.
> f:= x -> 1/((8*x-4)^2+1):\
[seq(evalf(J(f)-T[3](8*jj)), jj= 1 .. 10)];
[.0085039025, -.0002856241, -.0002048834, -.0000297271, -.0000143341,
-6 -5 -6 -6 -6
.6357*10 , -.15203*10 , .4373*10 , -.2525*10 , .1042*10 ]
> [seq(evalf(J(f)-ncrule[8](8*jj)), jj= 1 .. 10)];
-5 -6
[.0888176280, -.0006237649, .0005724945, .93610*10 , .0000112229, .9989*10 ,
-6 -8 -7 -7
.5013*10 , -.15*10 , .452*10 , -.110*10 ]
# We can carry the Richardson extrapolation idea as far as we
# want. It's usually done in a table. Note that only at the
# first level (the trapezoid rule) do you actually evaluate f
# (and the values from traprule(n/2) can be reused in
# traprule(n)), so the following will only require the same
# number of function evaluations as traprule(32). I'm going to
# use a different example, that works out better than the last
# one.
> f:= x -> x/(x^2 + 1/10):
> for jj from 1 to 5 do R[jj,1]:= evalf(traprule(2^jj)) od:
> for kk from 2 to 5 do \
for jj from kk to 5 do \
R[jj,kk]:= (4^(kk-1)*R[jj,kk-1] - R[jj-1,kk-1])/(4^(kk-1)-1)\
od od;
> with(linalg, matrix):\
for kk from 2 to 5 do\
for jj from 1 to kk-1 do\
R[jj,kk]:= ` `\
od od;
> R[1,1];
.9415584416
> matrix([seq([seq(R[jj,kk], jj=2..5)], kk=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 true value, for comparison, is
> evalf(J(f));
1.198947637
# And here are the error estimates.
> for kk from 1 to 4 do for jj from kk+1 to 5 do \
Ee[jj,kk]:= (R[jj,kk]-R[jj-1,kk])/(4^kk-1)\
od od;
> for kk from 2 to 4 do\
for jj from 1 to kk do\
Ee[jj,kk]:= ` `\
od od;
> matrix([seq([seq(Ee[jj,kk], jj=2..5)], kk=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 kk from 1 to 5 do for jj from kk to 5 do \
Et[jj,kk]:= evalf(J(f))-R[jj,kk]\
od od;
> for kk from 2 to 5 do\
for jj from 1 to kk-1 do\
Et[jj,kk]:= ` `\
od od;
> matrix([seq([seq(Et[jj,kk], jj=2..5)], kk=1..5)]);
[ .060534164 .014211111 .003510259 .000875130 ]
[ ]
[ -5 ]
[ -.005084180 -.001229907 -.000056692 -.3246*10 ]
[ ]
[ -6 ]
[ -.000972955 .000021522 .317*10 ]
[ ]
[ -7 ]
[ .000037308 -.19*10 ]
[ ]
[ -6 ]
[ -.165*10 ]
# One indication of whether everything is going according to
# plan is whether the error estimates behave the right way as
# you go to across the row. In the k=1 row, where R[j,1] =
# traprule(2^j) which should have an error proportional to 1/n^2
# = 1/4^j, each error should be about 1/4 of the previous one.
> seq(Ee[jj,1]/Ee[jj+1,1], jj = 2 .. 4);
4.249612635, 4.328912597, 4.060845598
# In the k=2 row, the error should be proportional to 1/n^4, so
# each error should be about 1/16 of the previous one.
> seq(Ee[jj,2]/Ee[jj+1,2], jj = 3 .. 4);
3.285223083, 21.95140890
# The last one's not too bad, though the first one's way off.
# In the k=3 row, we should have about 1/64.
> Ee[4,3]/Ee[5,3];
-46.89823154
# Not only is it not 64, it's negative. It appears that the
# error estimates in the third row are not trustworthy (and
# indeed Ee[4,3] is quite bad).
#
