f16-7

740 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).16.7 Multistep, Multivalue, and Predictor-Corrector Methods Thetermsmultistepandmultivaluedescribetwodifferentwaysofimplementing essentially the same integrationtechniqueforODEs. Predictor-correctoris a partic- ular subcategrory of these methods — in fact, the most widely used. Accordingly,the name predictor-correctoris often loosely used to denote all these methods. We suspectthatpredictor-correctorintegratorshavehadtheirday,andthatthey are no longerthe methodof choice for most problems in ODEs. For high-precision applications,orapplicationswhereevaluationsoftheright-handsidesareexpensive, Bulirsch-Stoer dominates. For convenience, or for low precision, adaptive-stepsizeRunge-Kuttadominates. Predictor-correctormethodshavebeen,wethink,squeezed out in the middle. There is possibly only one exceptional case: high-precision solution of very smooth equations with very complicated right-hand sides, as wewill describe later. Nevertheless, these methods have had a long historical run. Textbooks are full of information on them, and there are a lot of standard ODE programs around that are based on predictor-corrector methods. Many capable researchers have a lot of experience with predictor-corrector routines, and they see no reason to makea precipitous change of habit. It is not a bad idea for you to be familiar with the principlesinvolved,andevenwith the sorts of bookkeepingdetails that are the bane ofthesemethods. Otherwisetherewillbeabigsurpriseinstorewhenyoufirst haveto fix a problem in a predictor-corrector routine. Let us first consider the multistep approach. Think about how integrating an ODEisdifferentfromfindingtheintegralofafunction: Forafunction,theintegrand has a known dependence on the independent variable x, and can be evaluated at will. For an ODE, the “integrand” is the right-hand side, which depends both onxand on the dependent variables y. Thus to advance the solution of y /prime=f(x, y ) from xntox,w eh a v e y(x)=yn+/integraldisplayx xnf(x/prime,y)dx/prime(16.7.1 ) In a single-step method like Runge-Kuttaor Bulirsch-Stoer,the value yn+1atxn+1 dependsonlyon yn. Inamultistepmethod,weapproximate f(x, y )byapolynomial passing through severalprevious points xn,x n−1,...and possibly also through xn+1. Theresult ofevaluatingthe integral(16.7.1)at x=xn+1is thenofthe form yn+1=yn+h(β0y/prime n+1+β1y/prime n+β2y/prime n−1+β3y/prime n−2+···)( 16.7.2 ) where y/prime ndenotes f(xn,yn), andso on. If β0=0, the methodis explicit; otherwise it is implicit. The order of the method depends on how many previous steps we use to get each new value of y. Considerhowwemightsolveanimplicitformulaoftheform(16.7.2)for yn+1. Two methods suggest themselves: functional iteration andNewton’s method .I n functionaliteration,we takesome initialguess for yn+1,insert it intothe right-hand side of (16.7.2)to get an updated value of yn+1, insert this updated value back into theright-handside,andcontinueiterating. Buthowarewetogetaninitialguessfor 16.7Multistep,Multivalue,andPredictor-CorrectorMethods 741Sample 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).yn+1? Easy! Just use some explicitformula of the same form as (16.7.2). This is called the predictor step . In the predictor step we are essentially extrapolating the polynomial fit to the derivative from the previous points to the new point xn+1and then doing the integral (16.7.1) in a Simpson-like manner from xntoxn+1. The subsequent Simpson-like integration, using the prediction step’s value of yn+1to interpolate the derivative, is called the corrector step . The difference between the predictedandcorrectedfunctionvalues supplies informationon the local truncation error that can be used to control accuracy and to adjust stepsize. If one corrector step is good, aren’t many better? Why not use each corrector as an improvedpredictorand iterate to convergenceon each step? Answer: Even ifyou had a perfectpredictor, the step would still be accurate only to the finite order of the corrector. This incurable error term is on the same order as that which your iterationissupposedtocure,soyouareatbestchangingonlythecoefficientinfrontoftheerrortermbyafractionalamount. Sodubiousanimprovementiscertainlynot worththeeffort. Yourextraeffortwouldbebetterspentintakingasmallerstepsize. As describedso far,youmightthinkit desirableornecessaryto predictseveral intervals ahead at each step, then to use all these intervals, with various weights, in a Simpson-like corrector step. That is not a good idea. Extrapolation is the least stable partof the procedure,andit is desirableto minimizeits effect. Therefore,the integrationstepsofapredictor-correctormethodareoverlapping,eachoneinvolving several stepsize intervals h, but extending just one such interval farther than the previousones. Onlythatoneextendedintervalisextrapolatedbyeachpredictorstep. The most popular predictor-corrector methods are probably the Adams- Bashforth-Moulton schemes, which have good stability properties. The Adams-Bashforth part is the predictor. For example, the third-order case is predictor: y n+1=yn+h 12(23y/prime n−16y/prime n−1+5y/prime n−2)+O(h4)(16.7.3 ) Hereinformationatthecurrentpoint xn,togetherwiththetwopreviouspoints xn−1 andxn−2(assumed equally spaced), is used to predict the value yn+1at the next point, xn+1. The Adams-Moultonpart is the corrector. The third-ordercase is corrector: yn+1=yn+h 12(5y/prime n+1+8y/prime n−y/prime n−1)+O(h4)(16.7.4 ) Without the trial value of yn+1from the predictor step to insert on the right-hand side, the corrector would be a nasty implicit equation for yn+1. There are actually three separate processes occurring in a predictor-corrector method: the predictor step, which we call P, the evaluation of the derivative y/prime n+1 from the latest value of y, which we call E, and the corrector step, which we call C. In this notation, iterating mtimes with the corrector (a practice we inveighed against earlier)would be written P(EC)m. One also has the choiceof finishing with a C or an E step. The lore is that a final E is superior, so the strategy usually recommended is PECE. Notice that a PC method with a fixed number of iterations (say, one) is an explicit method! When we fix the number of iterations in advance, then the final valueof yn+1canbewrittenassomecomplicatedfunctionofknownquantities. Thus fixed iteration PC methods lose the strong stability properties of implicit methods andshould only be used for nonstiff problems . 742 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).For stiff problems we mustuse an implicit method if we want to avoid having tiny stepsizes. (Not all implicit methods are goodfor stiff problems,but fortunatelysomegoodonessuchas the Gearformulasare known.) We then appearto havetwo choices for solving the implicit equations: functional iteration to convergence, or Newtoniteration. However,itturnsoutthatforstiffproblemsfunctionaliterationwillnot even converge unless we use tiny stepsizes, no matter how close our prediction is! Thus Newton iteration is usually an essential part of a multistep stiff solver. For convergence,Newton’smethoddoesn’tparticularlycarewhatthestepsizeis,aslong as the prediction is accurate enough. Multistep methods, as we have described them so far, suffer from two serious difficulties when one tries to implement them: •Since the formulas require results from equally spaced steps, adjusting the stepsize is difficult. •Starting and stopping present problems. For starting, we need the initial values plus several previous steps to prime the pump. Stopping is a problem because equal steps are unlikely to land directly on the desiredtermination point. Older implementations of PC methods have various cumbersome ways of dealing with these problems. For example, they might use Runge-Kutta to start and stop. Changing the stepsize requires considerable bookkeeping to do some kind of interpolation procedure. Fortunately both these drawbacks disappear withthe multivalue approach. For multivalue methods the basic data available to the integrator are the first fewterms ofthe Taylorseries expansionofthesolutionat thecurrentpoint x n. The aimis toadvancethesolutionandobtaintheexpansioncoefficientsat thenextpointx n+1. This is in contrast to multistep methods, where the data are the values of the solution at xn,x n−1,.... We’ll illustrate the idea by considering a four-value method, for which the basic data are yn≡ yn hy/prime n (h2/2)y/prime/prime n (h3/6)y/prime/prime/prime n  (16.7.5 ) It is also conventionalto scale the derivativeswith the powersof h=xn+1−xnas shown. Note that here we use the vector notation yto denote the solution and its first few derivativesat a point,notthe fact that we aresolvinga system ofequationswith many components y. In terms of the data in (16.7.5), we can approximate the value of the solution yat some point x: y(x)=y n+(x−xn)y/prime n+(x−xn)2 2y/prime/prime n+(x−xn)3 6y/prime/prime/prime n (16.7.6 ) Setx=xn+1in equation (16.7.6) to get an approximation to yn+1. Differentiate equation(16.7.6)andset x=xn+1togetanapproximationto y/prime n+1,andsimilarlyfor y/prime/prime n+1andy/prime/prime/prime n+1. Calltheresultingapproximation /tildewideyn+1,wherethetildeisareminder 16.7Multistep,Multivalue,andPredictor-CorrectorMethods 743Sample 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).that all we have done so far is a polynomial extrapolation of the solution and its derivatives;wehavenotyetusedthedifferentialequation. Youcaneasilyverifythat /tildewideyn+1=B·yn (16.7.7 ) where the matrix Bis B= 1111 01230013 0001  (16.7.8 ) We now write the actual approximation to y n+1that we will use by adding a correction to /tildewideyn+1: yn+1=/tildewideyn+1+αr (16.7.9 ) Hererwill be a fixed vector of numbers, in the same way that Bis a fixed matrix. We fix αby requiring that the differential equation y/prime n+1=f(xn+1,yn+1)( 16.7.10 ) be satisfied. The second of the equations in (16.7.9) is hy/prime n+1=h/tildewidey/prime n+1+αr 2 (16.7.11 ) and this will be consistent with (16.7.10) provided r2=1,α =hf(xn+1,yn+1)−h/tildewidey/prime n+1 (16.7.12 ) Thevaluesof r1,r3, and r4arefreefortheinventorofa givenfour-valuemethodto choose. Different choices give different orders of method (i.e., through what order inhthe final expression 16.7.9 actually approximates the solution), and different stability properties. An interesting result, not obviousfrom our presentation, is that multivalue and multistep methods are entirely equivalent. In other words, the value yn+1given by a multivalue method with given Bandris exactly the same value given by some multistep method with given β’s in equation (16.7.2). For example, it turns out thattheAdams-Bashforthformula(16.7.3)correspondstoa four-valuemethodwithr 1=0,r3=3/4, and r4=1/6. The method is explicit because r1=0. The Adams-Moultonmethod(16.7.4)correspondstotheimplicitfour-valuemethodwith r1=5/12,r3=3/4, and r4=1/6. Implicit multivalue methods are solved the same way as implicit multistep methods: either by a predictor-corrector approach usingan explicitmethodforthepredictor,orbyNewtoniterationforstiff systems. Why go to all the trouble of introducing a whole new method that turns out to be equivalent to a method you already knew? The reason is that multivalue methodsallowaneasysolutiontothetwodifficultieswementionedaboveinactuallyimplementing multistep methods. Consider first the question of stepsize adjustment. To change stepsize from h toh /primeat some point xn, simply multiply the components of ynin (16.7.5) by the appropriate powers of h/prime/h, and you are ready to continue to xn+h/prime. 744 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).Multivalue methods also allow a relatively easy change in the orderof the method: Simply change r. The usual strategy for this is first to determine the new stepsize with the current order from the error estimate. Then check what stepsize would be predicted using an order one greater and one smaller than the current order. Choose the order that allows you to take the biggest next step. Being able tochangeorderalso allows an easy solutionto the starting problem: Simplystart with afirst-ordermethodandlettheorderautomaticallyincreasetotheappropriatelevel. For low accuracy requirements, a Runge-Kutta routine like rkqsis almost always the most efficient choice. For high accuracy, bsstepis both robust and efficient. For very smooth functions, a variable-order PC method can invoke veryhigh orders. If the right-hand side of the equation is relatively complicated, so that the expense of evaluating it outweighs the bookkeeping expense, then the best PC packages can outperform Bulirsch-Stoer on such problems. As you can imagine,however,suchavariable-stepsize,variable-ordermethodisnottrivialtoprogram. If yoususpectthat yourproblemis suitableforthis treatment,we recommenduse of a cannedPCpackage. ForfurtherdetailsconsultGear [1]orShampineandGordon [2]. Our prediction, nevertheless, is that, as extrapolation methods like Bulirsch- Stoer continue to gain sophistication, they will eventually beat out PC methods in all applications. We are willing, however, to be corrected. CITED REFERENCES AND FURTHER READING: Gear,C.W.1971, NumericalInitialValueProblemsinOrdinaryDifferentialEquations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 9. [1] Shampine, L.F., and Gordon, M.K. 1975, Computer Solution of Ordinary Differential Equations. The Initial Value Problem. (San Francisco: W.H Freeman). [2] Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 5. Kahaner,D.,Moler,C.,andNash,S.1989, NumericalMethods andSoftware (EnglewoodCliffs, NJ: Prentice Hall), Chapter 8. Hamming, R.W. 1962, Numerical Methods for Engineers and Scientists ; reprinted 1986 (New York: Dover), Chapters 14–15. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 7.