Richardson Extrapolation

1 Richardson ExtrapolationThere are many approximation procedures in which one first picks a step sizehandthen generates an approximationA(h) to some desired quantityA. Often the order of theerror generated by the procedure is known. In other wordsA=A(h) +Khk+K hk+1+K hk+2+ withkbeing some known constant andK, K , K , being some other (usually unknown)constants. For example,Amight be the valuey(tf) at some final timetffor the solution toan initial value problemy =f(t, y), y(t0) =y0. ThenA(h) might be the approximationtoy(tf) produced by Euler s method with step sizeh. In this casek= 1. If the improvedEuler s method is usedk= 2. If Runge-Kutta is usedk= notationO(hk+1) is conventionally used to stand for a sum of terms of orderhk+1and higher . So the above equation may be writtenA=A(h) +Khk+O(hk+1)(1)If we were to drop the, hopefully tiny, termO(hk+1) from this equation, we would have onelinear equation,A=A(h) +Khk, in the two unknownsA, K.

2 But this is really a differentequation for each different value ofh. We can get a second such equation just by using adifferent step size. Then the two equations may be solved, yielding approximate values ofAandK. This approximate value ofAconstitutes a new improved approximation,B(h), forthe exactA. We do this now. Taking 2ktimesA=A(h/2) +K(h/2)k+O(hk+1)(2)(note that, in equations (1) and (2), the symbol O(hk+1) is used to stand for twodifferentsums of terms of orderhk+1and higher) and subtracting equation (1) gives(2k 1)A= 2kA(h/2) A(h) +O(hk+1)A=2kA(h/2) A(h)2k 1+O(hk+1)c Joel Feldman. 2000. All rights if we defineB(h) =2kA(h/2) A(h)2k 1(3)thenA=B(h) +O(hk+1)(4)and we have generated an approximation whose error is of orderk+1, one better thanA(h) widely used numerical integration algorithm, called Romberg integration, applies thisformula repeatedly to the trapezoidal (1) = wherey(t) obeysy(0) = 1, y = 1 t+ (h) =approximate value fory(1) given by improved Euler with step (h) =2kA(h/2) A(h)2k 1withk= (h) % #B(h) %#.

3 065 .62 .0083 .04 160 The % column gives the percentage error and the # column gives the number of evalu-ations off(t, y) , by subtracting equation (2) from equation (1), we can =A(h) A(h/2) +Khk(1 12k)+O(hk+1)K=A(h/2) A(h)hk(1 12k)+O(h)Once we knowKwe can estimate the error inA(h/2) byE(h/2) =A A(h/2)=K(h/2)k+O(hk+1)=A(h/2) A(h)2k 1+O(hk+1)If this error is unacceptably large, we can useE(h) =Khkc Joel Feldman. 2000. All rights determine a step sizehthat will give an acceptable error. This is the basis for a numberof algorithms that incorporate automatic step size thatA(h/2) A(h)2k 1=B(h) A(h/2). One cannot get a still better guess forAby combiningB(h) andE(h/2). that we wished to use improved Euler to find a numerical approximationtoA=y(1), whereyis the solution to the initial value problemy =y 2t y(0) = 3 Suppose further that we are aiming for an error of 10 6.

4 If we run improved Euler withstep size (5 steps) and again with step size (10 steps) we get the approximate valuesA( ) = andA( ) = Since improved Euler hask= 2, Theseapproximate values obeyA=A( ) +K( )2+ higher order = +K( )2+ higher orderA=A( ) +K( )2+ higher order = +K( )2+ higher orderSubtracting0 = +K( )2 K( )2+ higher order + thatK error for step sizehisKh2+O(h3), so to achieve an error of 10 6we needKh2+O(h3) = 10 6 10 6 h 10 = we run improved Euler with step size1617we get the approximate valueA(1617) = In this illustrative, and purely artifical, example, we can solve the intial valueproblem exactly. The solution isy(t) = 2+2t+et, so that the exact value ofy(1) = ,to eight decimal places, and the error inA(162) is Joel Feldman. 2000. All rights