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. But this is really a differentequation for each different value ofh.

2 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) %#..065 .62 .0083.

3 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. 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.

4 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