f16-2

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 efficiency 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? 16.2AdaptiveStepsizeControlforRunge-Kutta 709Sample 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).two small stepsbig step x Figure 16.2.1. Step-doubling as a means for adaptive stepsize control in fourth-order Runge-Kutta. Points where the derivative is evaluated are shown as filled circles. The open circle represents the same derivatives asthe filled circle immediately above it, sothe total numberofevaluations is11per twosteps. Comparing theaccuracy ofthebigstepwiththetwosmallstepsgives acriterion foradjusting thestepsizeon the next step, or for rejecting the current step as inaccurate. Letus denotetheexactsolutionforanadvancefrom xtox+2hbyy(x+2h) and the two approximate solutions by y1(one step 2h) and y2(2 steps each of size h). Since the basic method is fourth order, the true solution and the two numerical approximations are related by y(x+2h)=y1+( 2h)5φ+O(h6)+... y(x+2h)=y2+2 (h5)φ+O(h6)+...(16.2.1 ) where, to order h5, the value φremains constant over the step. [Taylor series expansion tells us the φis a number whose order of magnitude is y(5)(x)/5!.] The first expressionin (16.2.1)involves (2h)5since the stepsize is 2h, while the second expressioninvolves 2(h5)sincetheerroroneachstepis h5φ. Thedifferencebetween the two numerical estimates is a convenientindicator of truncation error ∆≡y2−y1 (16.2.2 ) It is this difference that we shall endeavor to keep to a desired degree of accuracy, neither too large nor too small. We do this by adjusting h. It might also occur to you that, ignoring terms of order h6and higher, we can solve the two equations in (16.2.1) to improve our numerical estimate of the true solution y(x+2h), namely, y(x+2h)=y2+∆ 15+O(h6)( 16.2.3 ) This estimate is accurate to fifth order , one order higher than the original Runge- Kuttasteps. However,wecan ’thaveourcakeandeat it: (16.2.3)maybe fifth-order accurate, but we have no way of monitoring itstruncation error. Higher order is not always higher accuracy! Use of (16.2.3) rarely does harm, but we have noway of directly knowing whether it is doing any good. Therefore we should use ∆as the error estimate and take as “gravy”any additional accuracy gain derived from (16.2.3). In the technical literature, use of a procedure like (16.2.3) is called“local extrapolation. ” An alternativestepsize adjustmentalgorithmis based on the embeddedRunge- Kutta formulas , originally invented by Fehlberg. An interesting fact about Runge- Kutta formulas is that for orders Mhigher than four, more than Mfunction 710 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).evaluations (though never more than M+2) are required. This accounts for the popularity of the classical fourth-order method: It seems to give the most bangfor the buck. However, Fehlberg discovered a fifth-order method with six function evaluations where another combination of the six functions gives a fourth-order method. The difference between the two estimates of y(x+h)can then be used as an estimate of the truncation error to adjust the stepsize. Since Fehlberg ’s original formula, several other embedded Runge-Kutta formulas have been found. Many practitioners were at one time wary of the robustness of Runge-Kutta- Fehlbergmethods. Thefeelingwas thatusingthesameevaluationpointstoadvance thefunctionandtoestimatetheerrorwasriskierthanstep-doubling,wheretheerrorestimate is based on independent function evaluations. However, experience has shownthatthisconcernisnotaprobleminpractice. Accordingly,embeddedRunge- Kutta formulas, which are roughly a factor of two more ef ficient, have superseded algorithms based on step-doubling. The general form of a fifth-order Runge-Kutta formula is k 1=hf(xn,y n) k2=hf(xn+a2h, y n+b21k1) ··· k6=hf(xn+a6h, y n+b61k1+···+b65k5) yn+1=yn+c1k1+c2k2+c3k3+c4k4+c5k5+c6k6+O(h6)(16.2.4 ) The embedded fourth-order formula is y∗ n+1=yn+c∗ 1k1+c∗ 2k2+c∗ 3k3+c∗ 4k4+c∗ 5k5+c∗ 6k6+O(h5)(16.2.5 ) and so the error estimate is ∆≡yn+1−y∗ n+1=6/summationdisplay i=1(ci−c∗ i)ki (16.2.6 ) Theparticularvaluesof the variousconstants thatwe favorare thosefoundbyCash and Karp [2], and given in the accompanying table. These give a more ef ficient methodthan Fehlberg ’soriginal values, with somewhat better error properties. Now that we know, at least approximately, what our error is, we need to consider how to keep it within desired bounds. What is the relation between ∆ andh? According to (16.2.4) –(16.2.5), ∆scales as h5. If we take a step h1 and produce an error ∆1, therefore, the step h0thatwould have given some other value ∆0is readily estimated as h0=h1/vextendsingle/vextendsingle/vextendsingle/vextendsingle∆0 ∆1/vextendsingle/vextendsingle/vextendsingle/vextendsingle0.2 (16.2.7 ) Henceforth we will let ∆0denote the desiredaccuracy. Then equation (16.2.7) is used in two ways: If ∆1is larger than ∆0in magnitude, the equation tells how much to decrease the stepsize when we retry the present (failed) step .I f ∆1is 16.2AdaptiveStepsizeControlforRunge-Kutta 711Sample 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).Cash-KarpParametersforEmbeddedRunga-KuttaMethod i ai bij ci c∗ i 137 3782825 27648 21 51 50 0 33 103 409 40250 62118575 48384 43 53 10−9 106 5125 59413525 55296 5 1 −11 545 2−70 2735 270277 14336 67 81631 55296175 512575 1382444275 110592253 4096512 17711 4 j=1 2 345 smaller than ∆0, on the other hand, then the equationtells how much we can safely increase the stepsize for the next step . Local extrapolation consists in accepting thefifth order value yn+1, even though the error estimate actually applies to the fourth order value y∗ n+1. Our notation hides the fact that ∆0is actually a vector of desired accuracies, one foreach equationin the set of ODEs. In general,ouraccuracyrequirementwill be that all equations are within their respective allowed errors. In other words, wewill rescale the stepsize accordingto the needs of the “worst-offender ”equation. Howis ∆ 0, thedesiredaccuracy,relatedtosomelooserprescriptionlike “geta solution good to one part in 106”? That can be a subtle question, and it depends on exactlywhat yourapplicationis! You maybe dealingwith a set of equationswhose dependent variables differ enormously in magnitude. In that case, you probablywanttousefractionalerrors, ∆ 0=/epsilon1y,where /epsilon1is thenumberlike 10−6orwhatever. On the other hand, you may have oscillatory functions that pass through zero but are bounded by some maximum values. In that case you probably want to set ∆0 equal to /epsilon1times those maximum values. A convenient way to fold these considerations into a generally useful stepper routine is this: One of the arguments of the routine will of course be the vector of dependent variables at the beginning of a proposed step. Call that y(1:n). Let us require the user to specify for each step another, corresponding, vector argumentyscal(1:n) , and also an overall tolerance level eps. Then the desired accuracy for the ith equation will be taken to be ∆ 0=eps×yscal(i) (16.2.8 ) If you desire constant fractional errors, plug yinto the yscalcalling slot (no need to copy the values into a different array). If you desire constant absolute errorsrelativetosomemaximumvalues,settheelementsof yscalequaltothosemaximum values. A useful “trick”for getting constant fractional errors except“very”near zero crossings is to set yscal(i) equal to |y(i)|+|h×dydx(i) |. (The routine odeint, below, does this.) Here is a more technical point. We have to consider one additional possibility foryscal. The error criteria mentioned thus far are “local,”in that they bound the errorofeachstep individually. Insomeapplicationsyoumaybeunusuallysensitive 712 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).about a“global”accumulation of errors, from beginning to end of the integration and in the worst possible case where the errors all are presumed to add with thesame sign. Then, the smaller the stepsize h, the smaller the value ∆ 0that you will need to impose. Why? Because there will be more steps between your starting and ending values of x. In such cases you will want to set yscalproportional to h, typically to something like ∆0=/epsilon1h×dydx(i) (16.2.9 ) Thisenforcesfractionalaccuracy /epsilon1notonthevaluesof ybut(muchmorestringently) ontheincrements to thosevaluesat eachstep. But nowlookbackat (16.2.7). If ∆0 has an implicit scaling with h, then the exponent 0.20is no longer correct: When thestepsizeisreducedfromatoo-largevalue,thenewpredictedvalue h1willfailto meet the desired accuracywhen yscalis also altered to this new h1value. Instead of0.20 = 1 /5,we must scale bythe exponent 0.25 = 1 /4for thingsto workout. The exponents 0.20and0.25are not really very different. This motivates us to adopt the following pragmatic approach, one that frees us from having to know in advance whether or not you, the user, plan to scale your yscal’s with stepsize. Wheneverwedecreaseastepsize,letususethelargervalueoftheexponent(whetherwe need it or not!), and whenever we increase a stepsize, let us use the smaller exponent. Furthermore, because our estimates of error are not exact, but only accuratetotheleadingorderin h,we areadvisedtoputinasafetyfactor Swhichis a few percent smaller than unity. Equation (16.2.7) is thus replaced by h 0=  Sh 1/vextendsingle/vextendsingle/vextendsingle/vextendsingle∆0 ∆1/vextendsingle/vextendsingle/vextendsingle/vextendsingle0.20 ∆0≥∆1 Sh 1/vextendsingle/vextendsingle/vextendsingle/vextendsingle∆ 0 ∆1/vextendsingle/vextendsingle/vextendsingle/vextendsingle0.25 ∆0<∆1(16.2.10 ) We have found this prescription to be a reliable one in practice. Here, then, is a stepper program that takes one “quality-controlled ”Runge- Kutta step. SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs) INTEGER n,NMAX REAL eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n) EXTERNAL derivsPARAMETER (NMAX=50) Maximum number of equations. C USES derivs,rkck Fifth-order Runge-Kutta step with monitoring of local truncation error to ensure accuracy and adjust stepsize. Input are the dependent variable vector y(1:n) and its derivative dydx(1:n) at the starting value of the independent variable x. Also input are the stepsize to be attempted htry , the required accuracy eps, and the vector yscal(1:n) against which the error is scaled. On output, yandxare replaced by their new values, hdid is the stepsize that was actually accomplished, and hnext is the estimated next stepsize. derivs is the user-supplied subroutine that computes the right-hand side derivatives. INTEGER iREAL errmax,h,htemp,xnew,yerr(NMAX),ytemp(NMAX),SAFETY,PGROW, * PSHRNK,ERRCON PARAMETER (SAFETY=0.9,PGROW=-.2,PSHRNK=-.25,ERRCON=1.89e-4) The value ERRCON equals (5/SAFETY)**(1/PGROW) , see use below. h=htry Set stepsize to the initial trial value. 1 call rkck(y,dydx,n,x,h,ytemp,yerr,derivs) Take a step. 16.2AdaptiveStepsizeControlforRunge-Kutta 713Sample 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).errmax=0. Evaluate accuracy. do11i=1,n errmax=max(errmax,abs(yerr(i)/yscal(i))) enddo 11 errmax=errmax/eps Scale relative to required tolerance. if(errmax.gt.1.)then Truncation error too large, reduce stepsize. htemp=SAFETY*h*(errmax**PSHRNK)h=sign(max(abs(htemp),0.1*abs(h)),h) No more than a factor of 10. xnew=x+h if(xnew.eq.x)pause ’stepsize underflow in rkqs’ goto 1 For another try. else Step succeeded. Compute size of next step. if(errmax.gt.ERRCON)then hnext=SAFETY*h*(errmax**PGROW) else No more than a factor of 5 increase. hnext=5.*h endif hdid=hx=x+hdo 12i=1,n y(i)=ytemp(i) enddo 12 return endif END Theroutine rkqscallstheroutine rkcktotakeaCash-KarpRunge-Kuttastep: SUBROUTINE rkck(y,dydx,n,x,h,yout,yerr,derivs) INTEGER n,NMAXREAL h,x,dydx(n),y(n),yerr(n),yout(n)EXTERNAL derivs PARAMETER (NMAX=50) Set to the maximum number of functions. C USES derivs Given values for nvariables yand their derivatives dydx known at x, use the fifth-order Cash-Karp Runge-Kutta method to advance the solution over an interval hand return the incremented variables as yout . Also return an estimate of the local truncation er- ror in yout using the embedded fourth-order method. The user supplies the subroutine derivs(x,y,dydx) , which returns derivatives dydx atx. INTEGER i REAL ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX),ak6(NMAX), * ytemp(NMAX),A2,A3,A4,A5,A6,B21,B31,B32,B41,B42,B43,B51,* B52,B53,B54,B61,B62,B63,B64,B65,C1,C3,C4,C6,DC1,DC3, * DC4,DC5,DC6 PARAMETER (A2=.2,A3=.3,A4=.6,A5=1.,A6=.875,B21=.2,B31=3./40., * B32=9./40.,B41=.3,B42=-.9,B43=1.2,B51=-11./54.,B52=2.5, * B53=-70./27.,B54=35./27.,B61=1631./55296.,B62=175./512., * B63=575./13824.,B64=44275./110592.,B65=253./4096.,* C1=37./378.,C3=250./621.,C4=125./594.,C6=512./1771.,* DC1=C1-2825./27648.,DC3=C3-18575./48384., * DC4=C4-13525./55296.,DC5=-277./14336.,DC6=C6-.25) do 11i=1,n First step. ytemp(i)=y(i)+B21*h*dydx(i) enddo 11 call derivs(x+A2*h,ytemp,ak2) Second step. do12i=1,n ytemp(i)=y(i)+h*(B31*dydx(i)+B32*ak2(i)) enddo 12 call derivs(x+A3*h,ytemp,ak3) Third step. do13i=1,n ytemp(i)=y(i)+h*(B41*dydx(i)+B42*ak2(i)+B43*ak3(i)) 714 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).enddo 13 call derivs(x+A4*h,ytemp,ak4) Fourth step. do14i=1,n ytemp(i)=y(i)+h*(B51*dydx(i)+B52*ak2(i)+B53*ak3(i)+ * B54*ak4(i)) enddo 14 call derivs(x+A5*h,ytemp,ak5) Fifth step. do15i=1,n ytemp(i)=y(i)+h*(B61*dydx(i)+B62*ak2(i)+B63*ak3(i)+ * B64*ak4(i)+B65*ak5(i)) enddo 15 call derivs(x+A6*h,ytemp,ak6) Sixth step. do16i=1,n Accumulate increments with proper weights. yout(i)=y(i)+h*(C1*dydx(i)+C3*ak3(i)+C4*ak4(i)+ * C6*ak6(i)) enddo 16 do17i=1,n Estimate error as difference between fourth and fifth order methods. yerr(i)=h*(DC1*dydx(i)+DC3*ak3(i)+DC4*ak4(i)+DC5*ak5(i) * +DC6*ak6(i)) enddo 17 return END Notingthattheaboveroutinesareallinsingleprecision,don ’tbetoogreedyin specifying eps. Thepunishmentforexcessivegreedinessisinterestingandworthyof GilbertandSullivan ’sMikado: Theroutinecanalwaysachieveanapparent zeroerror bymakingthestepsizesosmallthatquantitiesoforder hy/primeaddtoquantitiesoforder yas if they were zero. Then the routine chugs happily along taking in finitely many infinitesimal steps and neverchangingthe dependentvariablesoneiota. (Youguard against this catastrophic loss of your computer budget by signaling on abnormallysmall stepsizes or on the dependent variable vector remaining unchangedfrom step to step. On a personal workstation you guard against it by not taking too long a lunch hour while your program is running.) Here is a full- fledged“driver”for Runge-Kutta with adaptive stepsize control. Wewarmlyrecommendthisroutine,oronelikeit,foravarietyofproblems,notablyincluding garden-varietyODEs or sets of ODEs, and de finite integrals (augmenting the methods of Chapter 4). For storage of intermediate results (if you desire to inspectthem)we assumea commonblock path, whichcan holdupto KMAXXsteps. Because steps occur at unequal intervals results are stored only at intervals greater than dxsav. Also in the block is kmax, indicating the number of steps that can be stored. If kmax=0thereisnointermediatestorage,andtherestofthecommonblock need not exist. Otherwise you should set kmax =KMAXX. Storage of steps stops ifkmaxis exceeded, except that the ending values are always stored. Again, these controls are merely indicative of what you might need. The routine odeintshould be customized to the problem at hand. SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs) INTEGER nbad,nok,nvar,KMAXX,MAXSTP,NMAX REAL eps,h1,hmin,x1,x2,ystart(nvar),TINYEXTERNAL derivs,rkqsPARAMETER (MAXSTP=10000,NMAX=50,KMAXX=200,TINY=1.e-30) Runge-Kutta driver with adaptive stepsize control. Integrate the starting values ystart(1:nvar) from x1tox2with accuracy eps, storing intermediate results in the common block /path/ . h1should be set as a guessed first stepsize, hmin as the minimum allowed stepsize (can be zero). On output nok andnbad are the number of good and bad (but retried and 16.2AdaptiveStepsizeControlforRunge-Kutta 715Sample 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).fixed) steps taken, and ystart is replaced by values at the end of the integration interval. derivs is the user-supplied subroutine for calculating the right-hand side derivative, while rkqs is the name of the stepper routine to be used. /path/ contains its own information about how often an intermediate value is to be stored. INTEGER i,kmax,kount,nstp REAL dxsav,h,hdid,hnext,x,xsav,dydx(NMAX),xp(KMAXX),y(NMAX), * yp(NMAX,KMAXX),yscal(NMAX) COMMON /path/ kmax,kount,dxsav,xp,yp User storage for intermediate results. Preset dxsav andkmax . x=x1 h=sign(h1,x2-x1)nok=0nbad=0 kount=0 do 11i=1,nvar y(i)=ystart(i) enddo 11 if (kmax.gt.0) xsav=x-2.*dxsav Assures storage of first step. do16nstp=1,MAXSTP Take at most MAXSTP steps. call derivs(x,y,dydx) do12i=1,nvar Scaling used to monitor accuracy. This general-purpose choice can be modified if needbe. yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY enddo 12 if(kmax.gt.0)then if(abs(x-xsav).gt.abs(dxsav)) then Store intermediate results. if(kount.lt.kmax-1)then kount=kount+1xp(kount)=x do 13i=1,nvar yp(i,kount)=y(i) enddo 13 xsav=x endif endif endifif((x+h-x2)*(x+h-x1).gt.0.) h=x2-x If stepsize can overshoot, decrease. call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs) if(hdid.eq.h)then nok=nok+1 else nbad=nbad+1 endifif((x-x2)*(x2-x1).ge.0.)then Are we done? do 14i=1,nvar ystart(i)=y(i) enddo 14 if(kmax.ne.0)then kount=kount+1 Save final step. xp(kount)=xdo 15i=1,nvar yp(i,kount)=y(i) enddo 15 endif return Normal exit. endifif(abs(hnext).lt.hmin) pause ’stepsize smaller than minimum in odeint’h=hnext enddo 16 pause ’too many steps in odeint’ returnEND 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,onecouldusethemodi fiedmidpointmethodinitsownrightasanODE 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 modi fied 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 thefinal 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 quali fier“modified.” Themodi fiedmidpointmethodisasecond-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 modi fiedmidpointmethod by itself will be dominated by the embedded Runge-Kutta method with adaptive stepsize control, as implemented in the preceding section. Theusefulnessofthemodi fiedmidpointmethodtotheBulirsch-Stoertechnique (§16.4)derivesfroma “deep”result aboutequations(16.3.2),dueto Gragg. Itturns