f16-1

704 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). Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 5. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 7. Lambert, J. 1973, Computational Methods inOrdinaryDifferential Equations (NewYork: Wiley). Lapidus, L., and Seinfeld, J. 1971, Numerical Solution of Ordinary Differential Equations (New York: Academic Press). 16.1 Runge-Kutta Method The formula for the Euler method is yn+1=yn+hf(xn,y n)( 16.1.1 ) whichadvancesasolutionfrom xntoxn+1≡xn+h. Theformulaisunsymmetrical: It advances the solution through an interval h, but uses derivative information only atthebeginningofthatinterval(seeFigure16.1.1). Thatmeans(andyoucanverifyby expansion in power series) that the step’s error is only one power of hsmaller than the correction, i.e O(h 2)added to (16.1.1). ThereareseveralreasonsthatEuler’smethodis notrecommendedforpractical use, among them, (i) the method is not very accurate when compared to other, fancier, methods run at the equivalent stepsize, and (ii) neither is it very stable(see§16.6 below). Consider, however, the use of a step like (16.1.1) to take a “trial” step to the midpoint of the interval. Then use the value of both xand yat that midpoint to compute the “real” step across the whole interval. Figure 16.1.2 illustrates the idea. In equations, k 1=hf(xn,y n) k2=hf/parenleftbig xn+1 2h, y n+1 2k1/parenrightbig yn+1=yn+k2+O(h3)(16.1.2 ) As indicated in the error term, this symmetrization cancels out the first-order error term, making the method second order . [A method is conventionally called nth order if its error term is O(hn+1).] In fact, (16.1.2) is called the second-order Runge-Kutta ormidpoint method. We needn’t stop there. There are many ways to evaluate the right-hand side f(x, y )thatallagreetofirstorder,butthathavedifferentcoefficientsofhigher-order error terms. Adding up the right combination of these, we can eliminate the error termsorderbyorder. ThatisthebasicideaoftheRunge-Kuttamethod. Abramowitz andStegun [1],andGear [2],givevariousspecificformulasthatderivefromthisbasic 16.1Runge-KuttaMethod 705Sample 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).y(x) 12 x1 x2 x3 x Figure 16.1.1. Euler ’s method. In this simplest (and least accurate) method for integrating an ODE, the derivative at the starting point of each interval is extrapolated to find the next function value. The method has first-order accuracy. y(x) 12 x1 x2 x3 x34 5 Figure 16.1.2. Midpoint method. Second-order accuracy is obtained by using the initial derivative at each step to find a point halfway across the interval, then using the midpoint derivative across the full width of the interval. In the figure,filled dots represent final function values, while open dots represent function values that are discarded once their derivatives have been calculated and used. idea. By far the most often used is the classical fourth-order Runge-Kuttaformula , which has a certain sleekness of organization about it: k1=hf(xn,y n) k2=hf(xn+h 2,y n+k1 2) k3=hf(xn+h 2,y n+k2 2) k4=hf(xn+h, y n+k3) yn+1=yn+k1 6+k2 3+k3 3+k4 6+O(h5)( 16.1.3 ) The fourth-order Runge-Kutta method requires four evaluations of the right- hand side per step h(see Figure 16.1.3). This will be superior to the midpoint method(16.1.2) ifatleast twiceas largeastep ispossiblewith(16.1.3)forthesame accuracy. Is that so? The answer is: often, perhaps even usually, but surely not always! Thistakesusbacktoacentraltheme,namelythat highorder doesnotalways meanhighaccuracy . Thestatement “fourth-orderRunge-Kuttaisgenerallysuperior to second-order ”is a true one, but you should recognize it as a statement about the 706 Chapter16. Integrationof OrdinaryDifferentialEquationsSample 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).1 2 3 4yn + 1yn Figure 16.1.3. Fourth-order Runge-Kutta method. In each step the derivative is evaluated four times: once at the initial point, twice at trial midpoints, and once at a trial endpoint. From these derivatives thefinal function value (shown as a filled dot) is calculated. (See text for details.) contemporarypracticeofscienceratherthanasastatementaboutstrictmathematics. Thatis,itre flectsthenatureoftheproblemsthatcontemporaryscientistsliketosolve. Formanyscienti ficusers,fourth-orderRunge-Kuttaisnotjustthe firstwordon ODE integrators,but thelast wordas well. In fact,youcan getprettyfar onthis old workhorse, especially if you combine it with an adaptive stepsize algorithm. Keep in mind, however, that the old workhorse ’s last trip may well be to take you to the poorhouse: Bulirsch-Stoer or predictor-corrector methods can be very much more efficient for problems where very high accuracy is a requirement. Those methods are the high-strung racehorses. Runge-Kutta is for ploughing the fields. However, even the old workhorse is more nimble with new horseshoes. In §16.2 we will give amodernimplementationofaRunge-Kuttamethodthatis quitecompetitiveas longas very high accuracy is not required. An excellent discussion of the pitfalls in constructing a good Runge-Kutta code is given in [3]. Here is the routine for carrying out one classical Runge-Kutta step on a set ofndifferential equations. You input the values of the independent variables, and you get out new values which are stepped by a stepsize h(which can be positive or negative). You will notice that the routine requires you to supply not only function derivsfor calculating the right-hand side, but also values of the derivatives at the starting point. Why not let the routine call derivsfor thisfirst value? The answer will become clear only in the next section, but in brief is this: This call may not be your only one with these starting conditions. You may have taken a previous step with too large a stepsize, and this is your replacement. In that case, you do not want to call derivsunnecessarily at the start. Note that the routine that follows has, therefore, only three calls to derivs. SUBROUTINE rk4(y,dydx,n,x,h,yout,derivs) INTEGER n,NMAXREAL h,x,dydx(n),y(n),yout(n) EXTERNAL derivs PARAMETER (NMAX=50) Set to the maximum number of functions. Given values for the variables y(1:n) and their derivatives dydx(1:n) known at x,u s e the fourth-order Runge-Kutta method to advance the solution over an interval hand return the incremented variables as yout(1:n) , which need not be a distinct array from y.T h e user supplies the subroutine derivs(x,y,dydx) , which returns derivatives dydx atx. INTEGER i REAL h6,hh,xh,dym(NMAX),dyt(NMAX),yt(NMAX) hh=h*0.5h6=h/6. xh=x+hh 16.1Runge-KuttaMethod 707Sample 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).do11i=1,n First step. yt(i)=y(i)+hh*dydx(i) enddo 11 call derivs(xh,yt,dyt) Second step. do12i=1,n yt(i)=y(i)+hh*dyt(i) enddo 12 call derivs(xh,yt,dym) Third step. do13i=1,n yt(i)=y(i)+h*dym(i) dym(i)=dyt(i)+dym(i) enddo 13 call derivs(x+h,yt,dyt) Fourth step. do14i=1,n Accumulate increments with proper weights. yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i)) enddo 14 returnEND The Runge-Kutta method treats every step in a sequence of steps in identical manner. Prior behavior of a solution is not used in its propagation. This is mathematicallyproper,sinceanypointalongthetrajectoryofanordinarydifferential equationcanserveasaninitialpoint. Thefactthatallstepsaretreatedidenticallyalso makes it easy to incorporateRunge-Kuttaintorelativelysimple “driver”schemes. Weconsideradaptivestepsizecontrol,discussedinthenextsection,anessential forseriouscomputing. Occasionally,however,youjustwanttotabulateafunctionat equallyspacedintervals,andwithoutparticularlyhighaccuracy. Inthemostcommoncase, you want to produce a graph of the function. Then all you need may be a simpledriverprogramthatgoesfromaninitial x stoafinal xfinaspecifiednumber of steps. To check accuracy,double the numberof steps, repeat the integration,and compareresults. This approachsurely does not minimize computertime, and it can failforproblemswhosenature requiresavariablestepsize,butitmaywellminimize user effort. On small problems, this may be the paramount consideration. Here is sucha driver,self-explanatory,which tabulates theintegratedfunctions in a common block path. SUBROUTINE rkdumb(vstart,nvar,x1,x2,nstep,derivs) INTEGER nstep,nvar,NMAX,NSTPMXPARAMETER (NMAX=50,NSTPMX=200) Maximum number of functions and maximum number of values to be stored.REAL x1,x2,vstart(nvar),xx(NSTPMX),y(NMAX,NSTPMX) EXTERNAL derivs COMMON /path/ xx,y Storage of results. C USES rk4 Starting from initial values vstart(1:nvar) known at x1use fourth-order Runge-Kutta to advance nstep equal increments to x2. The user-supplied subroutine derivs(x,v,dvdx) evaluates derivatives. Results are stored in the common block path .B e s u r e t o d i m e n s i o n the common block appropriately. INTEGER i,k REAL h,x,dv(NMAX),v(NMAX)do 11i=1,nvar Load starting values. v(i)=vstart(i) y(i,1)=v(i) enddo 11 xx(1)=x1 x=x1 h=(x2-x1)/nstepdo 13k=1,nstep Take nstep steps. call derivs(x,v,dv) 708 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).call rk4(v,dv,nvar,x,h,v,derivs) if(x+h.eq.x)pause ’stepsize not significant in rkdumb’ x=x+h xx(k+1)=x Store intermediate steps. do12i=1,nvar y(i,k+1)=v(i) enddo 12 enddo 13 return END CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.5. [1] Gear,C.W.1971, NumericalInitialValueProblemsinOrdinaryDifferentialEquations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 2. [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. [3] Rice, J.R. 1983, Numerical Methods, Software, andAnalysis (New York: McGraw-Hill), §9.2. 16.2 AdaptiveStepsizeControlforRunge-Kutta AgoodODEintegratorshouldexertsomeadaptivecontroloveritsownprogress, makingfrequentchangesinitsstepsize. Usuallythepurposeofthisadaptivestepsize control is to achieve some predetermined accuracy in the solution with minimumcomputational effort. Many small steps should tiptoe through treacherous terrain, while a few great strides should speed through smooth uninteresting countryside. The resulting gains in ef ficiency are not mere tens of percents or factors of two; they can sometimes be factors of ten, a hundred, or more. Sometimes accuracy may be demanded not directly in the solution itself, but in some related conservedquantity that can be monitored. Implementationofadaptivestepsizecontrolrequiresthatthesteppingalgorithm returninformationaboutitsperformance,mostimportant,anestimateofitstruncationerror. Inthissectionwewilllearnhowsuchinformationcanbeobtained. Obviously, the calculation of this information will add to the computational overhead, but the investment will generally be repaid handsomely. With fourth-order Runge-Kutta, the most straightforward technique by far is step doubling (see, e.g., [1]). We take each step twice, once as a full step, then, independently, as two half steps (see Figure 16.2.1). How much overhead is this, say in terms of the numberof evaluationsof the right-handsides? Each of the three separate Runge-Kutta steps in the procedure requires 4 evaluations, but the singleanddoublesequencesshareastartingpoint,sothetotalis11. Thisistobecompared not to 4, but to 8 (the two half-steps), since —stepsize control aside —we are achievingthe accuracyof the smaller (half)stepsize. The overheadcost is therefore a factor 1.375. What does it buy us?