# 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). #