f16-3

716 Chapter16. IntegrationofOrdinaryDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).CITED REFERENCES AND FURTHER READING: Gear,C.W.1971, NumericalInitialValueProblemsinOrdinaryDifferentialEquations (Englewood Cliffs, NJ: Prentice-Hall). [1] Cash,J.R.,andKarp,A.H.1990, ACMTransactionsonMathematicalSoftware ,vol.16,pp.201– 222. [2] Shampine,L.F.,andWatts,H.A.1977,in MathematicalSoftwareIII ,J.R.Rice,ed.(NewYork:Aca- demic Press), pp. 257–275; 1979, Applied Mathematics and Computation , vol. 5, pp. 93– 121. Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall). 16.3 Modified Midpoint Method This section discusses the modified midpoint method , which advances a vector of dependent variables y(x)from a point xto a point x+Hby a sequence of n substeps each of size h, h=H/n (16.3.1 ) Inprinciple,onecouldusethemodifiedmidpointmethodinitsownrightasanODE integrator. In practice, the method finds its most important application as a part of the more powerful Bulirsch-Stoer technique, treated in §16.4. You can therefore consider this section as a preamble to §16.4. The number of right-hand side evaluations required by the modified midpoint method is n+1. The formulas for the method are z0≡y(x) z1=z0+hf(x, z 0) zm+1=zm−1+2hf(x+mh, z m)for m=1,2,...,n −1 y(x+H)≈yn≡1 2[zn+zn−1+hf(x+H, z n)] (16.3.2 ) Herethe z’sareintermediateapproximationswhichmarchalonginstepsof h,while ynis the final approximation to y(x+H). The method is basically a “centered difference”or “midpoint”method(compareequation16.1.2),exceptat the first and last points. Those give the qualifier “modified.” Themodifiedmidpointmethodisasecond-ordermethod,like(16.1.2),butwith theadvantageofrequiring(asymptoticallyforlarge n)onlyonederivativeevaluation per step hinstead of the two required by second-orderRunge-Kutta. Perhaps there are applications where the simplicity of (16.3.2), easily coded in-line in some other program,recommendsit. Ingeneral,however,useofthe modifiedmidpointmethodby itself will be dominated by the embedded Runge-Kutta method with adaptive stepsize control, as implemented in the preceding section. TheusefulnessofthemodifiedmidpointmethodtotheBulirsch-Stoertechnique (§16.4)derivesfroma “deep”result aboutequations(16.3.2),dueto Gragg. Itturns 16.3ModifiedMidpointMethod 717Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).out thatthe errorof (16.3.2),expressedas a powerseries in h, the stepsize, contains onlyevenpowers of h, yn−y(x+H)=∞/summationdisplay i=1αih2i(16.3.3 ) where His held constant, but hchanges by varying nin (16.3.1). The importance of this even power series is that, if we play our usual tricks of combining steps to knock out higher-order error terms, we can gain twoorders at a time! For example, suppose nis even, and let yn/2denote the result of applying (16.3.1)and (16.3.2)with half as many steps, n→n/2. Then the estimate y(x+H)≈4yn−yn/2 3(16.3.4 ) isfourth-order accurate, the same as fourth-order Runge-Kutta, but requires only about 1.5 derivative evaluations per step hinstead of Runge-Kutta’s 4 evaluations. Don’tbe too anxious to implement(16.3.4),since we will soon do even better. Now would be a good time to look back at the routine qsimpin§4.2, and especially to compare equation (4.2.4) with equation (16.3.4) above. You will see that thetransitionin Chapter4to the ideaof Richardsonextrapolation,as embodied in Romberg integrationof §4.3, is exactly analogousto the transition in going from this section to the next one. Here is the routine that implements the modified midpoint method, which will be used below. SUBROUTINE mmid(y,dydx,nvar,xs,htot,nstep,yout,derivs) INTEGER nstep,nvar,NMAX REAL htot,xs,dydx(nvar),y(nvar),yout(nvar) EXTERNAL derivsPARAMETER (NMAX=50) Modified midpoint step. Dependent variable vector y(1:nvar) and its derivative vector dydx(1:nvar) are input at xs. Also input is htot , the total step to be made, and nstep , the number of substeps to be used. The output is returned as yout(1:nvar) , which need not be a distinct array from y; if it is distinct, however, then yanddydx are returned undamaged. INTEGER i,n REAL h,h2,swap,x,ym(NMAX),yn(NMAX)h=htot/nstep Stepsize this trip. do 11i=1,nvar ym(i)=y(i)yn(i)=y(i)+h*dydx(i) First step. enddo 11 x=xs+hcall derivs(x,yn,yout) Will use yout for temporary storage of derivatives. h2=2.*h do 13n=2,nstep General step. do12i=1,nvar swap=ym(i)+h2*yout(i)ym(i)=yn(i) yn(i)=swap enddo 12 x=x+h call derivs(x,yn,yout) 718 Chapter16. IntegrationofOrdinaryDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).enddo 13 do14i=1,nvar Last step. yout(i)=0.5*(ym(i)+yn(i)+h*yout(i)) enddo 14 return END CITED REFERENCES AND FURTHER READING: Gear,C.W.1971, NumericalInitialValueProblemsinOrdinaryDifferentialEquations (Englewood Cliffs, NJ: Prentice-Hall), §6.1.4. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §7.2.12. 16.4 Richardson Extrapolation and the Bulirsch-Stoer Method The techniques described in this section are not for differential equations containing nonsmooth functions. For example, you might have a differential equationwhoseright-handsideinvolvesafunctionthatisevaluatedbytablelook-up and interpolation. If so, go back to Runge-Kutta with adaptive stepsize choice: Thatmethoddoesanexcellentjoboffeelingitswaythroughrockyordiscontinuous terrain. It is also an excellent choice for quick-and-dirty, low-accuracy solutionof a set of equations. A second warning is that the techniques in this section are not particularly good for differential equations that have singular points insidethe interval of integration. A regular solution must tiptoe very carefully across suchpoints. Runge-Kuttawithadaptivestepsizecansometimeseffectthis;moregenerally, there are special techniquesavailable forsuch problems,beyondour scopehere. Apart from those two caveats, we believe that the Bulirsch-Stoer method, discussed in this section, is the best known way to obtain high-accuracy solutions to ordinary differential equations with minimal computational effort. (A possibleexception, infrequently encountered in practice, is discussed in §16.7.) Three key ideas are involved. The first is Richardson’s deferred approach to the limit , which we already met in §4.3 on Romberg integration. The idea is to consider the final answer of a numerical calculation as itself being an analytic function (if a complicated one) of an adjustable parameter like the stepsize h. That analytic function can be probed by performing the calculation with various values ofh,noneof them being necessarily small enough to yield the accuracy that we desire. When we know enough about the function, we fitit to some analytic form, and then evaluateit at that mythical and golden point h=0(see Figure 16.4.1). Richardson extrapolation is a method for turning straw into gold! (Lead into gold for alchemist readers.) Thesecondideahastodowithwhatkindoffittingfunctionisused. Bulirschand Stoer first recognizedthe strength of rational functionextrapolation in Richardson- type applications. That strength is to break the shackles of the power series and its limitedradiusofconvergence,outonlytothedistanceofthefirstpoleinthecomplex