f16-0

Sample 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).Chapter 16. Integration of Ordinary Differential Equations 16.0 Introduction Problems involving ordinary differential equations (ODEs) can always be reduced to the study of sets of first-order differential equations. For example the second-order equation d2y dx2+q(x)dy dx=r(x)( 16.0.1 ) can be rewritten as two first-order equations dy dx=z(x) dz dx=r(x)−q(x)z(x)(16.0.2 ) where zisanewvariable. ThisexemplifiestheprocedureforanarbitraryODE.The usualchoiceforthenewvariablesistoletthembejustderivativesofeachother(andoftheoriginalvariable). Occasionally,itis usefultoincorporateintotheirdefinition some other factors in the equation, or some powers of the independent variable, for the purpose of mitigating singular behavior that could result in overflows orincreased roundoff error. Let common sense be your guide: If you find that the original variables are smooth in a solution, while your auxiliary variables are doing crazy things, then figure out why and choose different auxiliary variables. The generic problem in ordinary differential equations is thus reduced to the study of a set of Ncoupledfirst-order differential equations for the functions y i,i=1 ,2,...,N, having the general form dy i(x) dx=fi(x, y 1,...,y N),i =1 ,...,N (16.0.3 ) where the functions fion the right-hand side are known. A problem involving ODEs is not completely specified by its equations. Even more crucial in determining how to attack the problem numerically is the nature ofthe problem’s boundary conditions. Boundary conditions are algebraic conditions on the values of the functions y iin (16.0.3). In general they can be satisfied at 701 702 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).discretespecifiedpoints,butdonotholdbetweenthosepoints,i.e.,arenotpreserved automaticallybythedifferentialequations. Boundaryconditionscanbeassimpleasrequiring that certain variables have certain numerical values, or as complicated as a set of nonlinear algebraic equations among the variables. Usually, it is the nature of the boundary conditions that determines which numerical methods will be feasible. Boundary conditions divide into two broad categories. •Ininitialvalueproblems allthe y iaregivenatsomestartingvalue xs,and it is desired to find the yi’s at some final point xf, or at some discrete list of points (for example, at tabulated intervals). •Intwo-point boundary value problems , on the other hand, boundary conditions are specified at more than one x. Typically, some of the conditions will be specified at xsand the remainder at xf. This chapterwill consider exclusivelythe initial value problem,deferringtwo- pointboundaryvalueproblems,which aregenerallymoredifficult, to Chapter17. Theunderlyingideaofanyroutineforsolvingtheinitialvalueproblemisalways this: Rewritethe dy’sand dx’sin(16.0.3)asfinitesteps ∆yand∆x,andmultiplythe equationsby ∆x. Thisgivesalgebraicformulasforthechangeinthefunctionswhen theindependentvariable xis“stepped”byone“stepsize” ∆x. Inthelimitofmaking thestepsizeverysmall,agoodapproximationtotheunderlyingdifferentialequation is achieved. Literal implementation of this procedure results in Euler’s method (16.1.1,below), which is, however, notrecommendedfor any practical use. Euler’s methodisconceptuallyimportant,however;onewayoranother,practicalmethodsall comedowntothissameidea: Addsmallincrementstoyourfunctionscorrespondingto derivatives (right-hand sides of the equations) multiplied by stepsizes. In this chapter we consider three major types of practical numerical methods for solving initial value problems for ODEs: •Runge-Kutta methods •RichardsonextrapolationanditsparticularimplementationastheBulirsch- Stoer method •predictor-corrector methods. A brief description of each of these types follows.1.Runge-Kutta methods propagate a solution over an interval by combining the informationfrom several Euler-style steps (each involvingone evaluationof the right-hand f’s), and then using the information obtained to match a Taylor series expansion up to some higher order. 2.Richardsonextrapolation usesthepowerfulideaofextrapolatingacomputed result to the value that wouldhave been obtained if the stepsize had been very much smaller than it actually was. In particular, extrapolation to zero stepsize is the desired goal. The first practical ODE integrator that implemented this idea wasdeveloped by Bulirsch and Stoer, and so extrapolation methods are often called Bulirsch-Stoer methods. 3.Predictor-corrector methods store the solution along the way, and use those results to extrapolate the solution one step advanced; they then correct the extrapolation using derivative information at the new point. These are best forvery smooth functions. Runge-Kutta is what you use when (i) you don’t know any better, or (ii) you haveanintransigentproblemwhereBulirsch-Stoerisfailing,or(iii)youhaveatrivial 16.0Introduction 703Sample 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).problem where computational efficiency is of no concern. Runge-Kutta succeeds virtuallyalways; but it is notusuallyfastest, exceptwhenevaluating fiis cheapand moderate accuracy ( <∼10−5) is required. Predictor-corrector methods, since they use past information, are somewhat more difficult to start up, but, for many smooth problems,theyarecomputationallymoreefficientthanRunge-Kutta. InrecentyearsBulirsch-Stoer has been replacing predictor-corrector in many applications, but it is too soon to say that predictor-corrector is dominated in all cases. However, it appears that only rather sophisticated predictor-corrector routines are competitive. Accordingly, we have chosen notto give an implementation of predictor-corrector in this book. We discuss predictor-corrector further in §16.7, so that you can use a canned routine should you encounter a suitable problem. In our experience, the relatively simple Runge-Kutta and Bulirsch-Stoer routines we give are adequate for most problems. Each of the three types of methods can be organized to monitor internal consistency. This allows numerical errors which are inevitably introduced into the solution to be controlled by automatic, ( adaptive) changing of the fundamental stepsize. We always recommend that adaptive stepsize control be implemented, and we will do so below. In general, all three types of methods can be applied to any initial value problem. Eachcomes with its own set ofdebits and credits that must be understood before it is used. We have organized the routines in this chapter into three nested levels. The lowest or “nitty-gritty” level is the piece we call the algorithm routine. This implementsthebasicformulasofthemethod,startswithdependentvariables y iatx, andreturnsnewvaluesofthedependentvariablesat thevalue x+h. Thealgorithm routine also yields up some information about the quality of the solution after the step. The routine is dumb, however, and it is unable to make any adaptive decision about whether the solution is of acceptable quality or not. That quality-control decision we encode in a stepperroutine. The stepper routinecallsthealgorithmroutine. Itmayrejecttheresult,setasmallerstepsize,and call the algorithm routine again, until compatibility with a predetermined accuracycriterion has been achieved. The stepper’s fundamental task is to take the largest stepsizeconsistentwithspecifiedperformance. Onlywhenthisisaccomplisheddoes the true power of an algorithm come to light. Above the stepper is the driverroutine, which starts and stops the integration, stores intermediateresults, and generallyacts as an interfacewith the user. Thereisnothing at all canonical about our driver routines. You should consider them to be examples, and you can customize them for your particular application. Oftheroutinesthatfollow, rk4,rkck,mmid,stoerm,and simprarealgorithm routines; rkqs,bsstep,stiff, and stifbsare steppers; rkdumbandodeint are drivers. Section16.6ofthis chaptertreatsthesubjectof stiff equations ,relevantbothto ordinarydifferentialequationsandalsotopartialdifferentialequations(Chapter19). 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,yn)( 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,yn) 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