f16-4

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 16.4RichardsonExtrapolationandtheBulirsch-StoerMethod 719Sample 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).6 steps2 steps4 steps ⊗ extrapolation to ∞ steps x x + Hy Figure 16.4.1. Richardson extrapolation as used in the Bulirsch-Stoer method. A large interval His spanned by different sequences of finer andfiner substeps. Their results are extrapolated to an answer that is supposed to correspond to in finitelyfine substeps. In the Bulirsch-Stoer method, the integrations are done by the modi fied midpoint method, and the extrapolation technique is rational function or polynomial extrapolation. plane. Rational function fits can remain good approximations to analytic functions even after the various terms in powers of hall have comparable magnitudes. In other words, hcan be so large as to make the whole notion of the “order”of the methodmeaningless —and the methodcanstill worksuperbly. Nevertheless,more recent experience suggests that for smooth problems straightforward polynomial extrapolationis slightly more ef ficient than rational functionextrapolation. We will accordingly adopt polynomial extrapolation as the default, but the routine bsstep below allows easy substitution of one kind of extrapolation for the other. You mightwishatthispointtoreview §3.1–§3.2,wherepolynomialandrationalfunction extrapolation were already discussed. The third idea was discussed in the section before this one, namely to use a method whose error function is strictly even, allowing the rational function or polynomialapproximationto be in terms of the variable h2instead of just h. Put these ideas together and you have the Bulirsch-Stoer method [1]. A single Bulirsch-Stoersteptakesusfrom xtox+H,where Hissupposedtobequitealarge —not at all in finitesimal —distance. That single step is a grand leap consisting of many (e.g., dozens to hundreds) substeps of modi fied midpoint method, which are then extrapolated to zero stepsize. Thesequenceofseparateattemptstocrosstheinterval Hismadewithincreasing values of n, the number of substeps. Bulirsch and Stoer originally proposed the sequence n=2 ,4,6,8,12,16,24,32,48,64,96,..., [nj=2 nj−2],... (16.4.1 ) More recent work by Deu flhard[2,3]suggests that the sequence n=2 ,4,6,8,10,12,14,..., [nj=2 j],... (16.4.2 ) is usually more ef ficient. For each step, we do not know in advance how far up this sequencewe will go. After each successive nis tried, a polynomialextrapolationis 720 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).attempted. That extrapolation returns both extrapolated values and error estimates. If the errors are not satisfactory, we go higher in n. If they are satisfactory, we go on to the next step and begin anew with n=2. Ofcoursetheremustbesomeupperlimit,beyondwhichweconcludethatthere is someobstacle in ourpathin the interval H, so that we must reduce Hratherthan just subdivide it more finely. In the implementations below, the maximum number ofn’s to be tried is called KMAXX. For reasons describedbelow we usually take this equal to 8; the 8th value of the sequence (16.4.2) is 16, so this is the maximum number of subdivisions of Hthat we allow. Weenforceerrorcontrol,asintheRunge-Kuttamethod,bymonitoringinternal consistency,andadaptingstepsizetomatchaprescribedboundonthelocaltruncation error. Eachnewresultfromthesequenceofmodi fiedmidpointintegrationsallowsa tableaulikethatin §3.1tobeextendedbyoneadditionalsetofdiagonals. Thesizeof the new correction added at each stage is taken as the (conservative)error estimate. How shouldwe use this errorestimate to adjust the stepsize? The best strategy now known is due to Deu flhard[2,3]. For completeness we describe it here: Supposetheabsolutevalueoftheerrorestimatereturnedfromthe kthcolumn(andhence thek+1st row) of the extrapolation tableau is /epsilon1k+1 ,k. Errorcontrol isenforced by requiring /epsilon1k+1 ,k</epsilon1 (16.4.3 ) as the criterion for accepting the current step, where /epsilon1is the required tolerance. For the even sequence (16.4.2) the order of the method is 2k+1: /epsilon1k+1 ,k∼H2k+1(16.4.4 ) Thusasimpleestimateofanewstepsize Hktoobtainconvergenceina fixedcolumn kwouldbe Hk=H/parenleftbigg/epsilon1 /epsilon1k+1 ,k/parenrightbigg1/(2k+1) (16.4.5 ) Which column kshould we aim to achieve convergence in? Let ’s compare the work required for different k. Suppose Akisthe work to obtain row kof the extrapolation tableau, soAk+1is the work to obtain column k. We will assume the work is dominated by the cost of evaluating the functions de fining the right-hand sides of the differential equations. For nk subdivisions in H, the number of function evaluations can be found from the recurrence A1=n1+1 Ak+1=Ak+nk+1(16.4.6 ) The work per unit step to get column kisAk+1/H k, which we nondimensionalize with a factor of Hand write as W k=Ak+1 HkH (16.4.7 ) =Ak+1/parenleftBig/epsilon1k+1 ,k /epsilon1/parenrightBig1/(2k+1) (16.4.8 ) The quantities W kcan be calculated during the integration. The optimal column index q is then de fined by W q=m i n k=1 ,...,k fW k (16.4.9 ) where kfis thefinal column, in which the error criterion (16.4.3) was satis fied. The q determined from (16.4.9) de fines the stepsize Hqto be used as the next basic stepsize, so that we can expect to get convergence in the optimal column q. Two important re finements have to be made to the strategy outlined so far: 16.4RichardsonExtrapolationandtheBulirsch-StoerMethod 721Sample 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).•If the current His“too small, ”thenkfwill be“too small, ”and so qremains “too small. ”It may be desirable to increase Hand aim for convergence in a column q>k f. •Ifthecurrent His“toobig,”wemaynotconvergeatallonthecurrentstepandwe willhave todecrease H. Wewould liketo detectthis by monitoring the quantities /epsilon1k+1 ,kfor each kso we can stop the current step as soon as possible. Deuflhard’s prescription fordealing withthesetwo problems uses ideas fromcommuni- cation theory to determine the “average expected convergence behavior ”of the extrapolation. His model produces certain correction factors α(k,q)by which Hkis to be multiplied to try to get convergence in column q. The factors α(k,q)depend only on /epsilon1and the sequence {ni} and so can be computed once during initialization: α(k,q)=/epsilon1Ak+1−Aq+1 (2k+1)( Aq+1−A1+1)fork<q (16.4.10 ) withα(q,q)=1. Now to handle the first problem, suppose convergence occurs in column q=kf. Then ratherthantaking Hqforthenextstep,wemightaimtoincreasethestepsizetogetconvergence incolumn q+1. Sincewedon ’t have Hq+1availablefromthecomputation, weestimateitas Hq+1=Hqα(q,q+1 ) ( 16.4.11 ) By equation (16.4.7) this replacement is ef ficient, i.e., reduces the work per unit step, if Aq+1 Hq>Aq+2 Hq+1(16.4.12 ) or Aq+1α(q,q+1 )>A q+2 (16.4.13 ) During initialization, this inequality can be checked for q=1,2,...to determine kmax, the largest allowed column. Then when (16.4.12) is satis fied it will always be ef ficient to use Hq+1. (In practice we limit kmaxto 8 even when /epsilon1is very small as there is very littlefurther gain in ef ficiency whereas roundoff can become a problem.) The problem of stepsize reduction is handled by computing stepsize estimates ¯Hk≡Hkα(k, q),k =1,...,q −1( 16.4.14 ) duringthecurrentstep. The ¯H’sareestimatesofthestepsizetogetconvergenceintheoptimal column q.I fa n y ¯Hkis“too small, ”we abandon the current step and restart using ¯Hk. The criterion of being “too small ”is taken to be Hkα(k, q+1 )<H (16.4.15 ) Theα’s satisfy α(k,q+1 ) >α(k, q). During the first step, when we have no information about the solution, the stepsize reduction check is made for all k. Afterwards, we test for convergence and for possible stepsize reduction only in an “order window ” max(1 ,q−1)≤k≤min(kmax,q+1 ) ( 16.4.16 ) The rationale for the order window is that if convergence appears to occur for k<q −1it is often spurious, resulting from some fortuitously small error estimate in the extrapolation.On the other hand, if you need to go beyond k=q+1to obtain convergence, your local model of the convergence behavior is obviously not very good and you need to cut thestepsize and reestablish it. In the routine bsstep, these various tests are actually carried out using quantities /epsilon1(k)≡H Hk=/parenleftBig/epsilon1k+1 ,k /epsilon1/parenrightBig1/(2k+1) (16.4.17 ) called err(k)in the code. As usual, we include a “safety factor ”in the stepsize selection. This is implemented by replacing /epsilon1by0.25/epsilon1. Other safety factors are explained in the program comments. Note that while the optimal convergence column is restricted to increase by at most one on each step, a sudden drop in order is allowed by equation (16.4.9). This gives the methoda degree of robustness for problems with discontinuities. 722 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).Let us remind you once again that scalingof the variables is often crucial for successful integration of differential equations. The scaling “trick”suggested in the discussion following equation (16.2.8) is a good general purpose choice, but not foolproof. Scaling by the maximum values of the variables is more robust, but requires you to have some prior information. The following implementation of a Bulirsch-Stoer step has exactly the same calling sequence as the quality-controlled Runge-Kutta stepper rkqs. This means that the driver odeintin§16.2 can be used for Bulirsch-Stoer as well as Runge- Kutta: Just substitute bsstepforrkqsinodeint’s argument list. The routine bsstepcalls mmidtotakethemodi fiedmidpointsequences,andcalls pzextr,given below, to do the polynomial extrapolation. SUBROUTINE bsstep(y,dydx,nv,x,htry,eps,yscal,hdid,hnext,derivs) INTEGER nv,NMAX,KMAXX,IMAX REAL eps,hdid,hnext,htry,x,dydx(nv),y(nv),yscal(nv),SAFE1, * SAFE2,REDMAX,REDMIN,TINY,SCALMX PARAMETER (NMAX=50,KMAXX=8,IMAX=KMAXX+1,SAFE1=.25,SAFE2=.7, * REDMAX=1.e-5,REDMIN=.7,TINY=1.e-30,SCALMX=.1) C USES derivs,mmid,pzextr Bulirsch-Stoer step with monitoring of local truncation error to ensure accuracy and adjuststepsize. Input are the dependent variable vector y(1:nv) and its derivative dydx(1:nv) at the starting value of the independent variable x. Also input are the stepsize to be at- tempted htry, the required accuracy eps, and the vector yscal(1:nv) against which the error is scaled. On output, yandxare replaced by their new values, hdidis the stepsize that was actually accomplished, and hnextis the estimated next stepsize. derivsis the user-supplied subroutine that computes the right-hand side derivatives. Be sure to set htry on successive steps to the value of hnextreturned from the previous step, as is the case if the routine is called by odeint. Parameters: NMAXis the maximum value of nv;KMAXXis the maximum row number used in the extrapolation; IMAXis the next row number; SAFE1andSAFE2are safety factors; REDMAXis the maximum factor used when a stepsize is reduced, REDMINthe minimum; TINYprevents division by zero; 1/ SCALMXis the maximum factor by which a stepsize can be increased. INTEGER i,iq,k,kk,km,kmax,kopt,nseq(IMAX) REAL eps1,epsold,errmax,fact,h,red,scale,work,wrkmin,xest, * xnew,a(IMAX),alf(KMAXX,KMAXX),err(KMAXX),yerr(NMAX),* ysav(NMAX),yseq(NMAX) LOGICAL first,reduct SAVE a,alf,epsold,first,kmax,kopt,nseq,xnew EXTERNAL derivsDATA first/.true./,epsold/-1./DATA nseq /2,4,6,8,10,12,14,16,18/ if(eps.ne.epsold)then A new tolerance, so reinitialize. hnext=-1.e29 “Impossible” values. xnew=-1.e29 eps1=SAFE1*eps a(1)=nseq(1)+1 Compute work coefficients A k. do11k=1,KMAXX a(k+1)=a(k)+nseq(k+1) enddo 11 do13iq=2,KMAXX Compute α(k, q ). do12k=1,iq-1 alf(k,iq)=eps1**((a(k+1)-a(iq+1))/ * ((a(iq+1)-a(1)+1.)*(2*k+1))) enddo 12 enddo 13 epsold=epsdo 14kopt=2,KMAXX-1 Determine optimal row number for conver- gence. if(a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt))goto 1 enddo 14 16.4RichardsonExtrapolationandtheBulirsch-StoerMethod 723Sample 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 kmax=kopt endif h=htry do15i=1,nv Save the starting values. ysav(i)=y(i) enddo 15 if(h.ne.hnext.or.x.ne.xnew)then Anewstepsizeoranewintegration: re-establish the order window. first=.true. kopt=kmax endif reduct=.false. 2d o 17k=1,kmax Evaluate the sequence of modified midpoint integrations. xnew=x+h if(xnew.eq.x)pause ’step size underflow in bsstep’ call mmid(ysav,dydx,nv,x,h,nseq(k),yseq,derivs)xest=(h/nseq(k))**2 Squared, since error series is even. call pzextr(k,xest,yseq,y,yerr,nv) Perform extrapolation. if(k.ne.1)then Compute normalized error estimate /epsilon1(k). errmax=TINYdo 16i=1,nv errmax=max(errmax,abs(yerr(i)/yscal(i))) enddo 16 errmax=errmax/eps Scale error relative to tolerance. km=k-1 err(km)=(errmax/SAFE1)**(1./(2*km+1)) endifif(k.ne.1.and.(k.ge.kopt-1.or.first))then In order window. if(errmax.lt.1.)goto 4 Converged. if(k.eq.kmax.or.k.eq.kopt+1)then Checkforpossiblestepsizereduction. red=SAFE2/err(km) goto 3 else if(k.eq.kopt)then if(alf(kopt-1,kopt).lt.err(km))then red=1./err(km) goto 3 endif else if(kopt.eq.kmax)then if(alf(km,kmax-1).lt.err(km))then red=alf(km,kmax-1)* * SAFE2/err(km) goto 3 endif else if(alf(km,kopt).lt.err(km))then red=alf(km,kopt-1)/err(km)goto 3 endif endif enddo 17 3 red=min(red,REDMIN) Reduce stepsize by at least REDMINand at most REDMAX. red=max(red,REDMAX) h=h*redreduct=.true. goto 2 Try again. 4 x=xnew Successful step taken. hdid=h first=.false. wrkmin=1.e35 Compute optimal row for convergence and corresponding stepsize. do 18kk=1,km fact=max(err(kk),SCALMX) work=fact*a(kk+1) if(work.lt.wrkmin)then scale=factwrkmin=work kopt=kk+1 724 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).endif enddo 18 hnext=h/scaleif(kopt.ge.k.and.kopt.ne.kmax.and..not.reduct)then Check for possible order in- crease, but not if step- size was just reduced.fact=max(scale/alf(kopt-1,kopt),SCALMX) if(a(kopt+1)*fact.le.wrkmin)then hnext=h/factkopt=kopt+1 endif endif returnEND Thepolynomialextrapolationroutineisbasedonthesamealgorithmas polint §3.1. Itissimplerinthatitisalwaysextrapolatingtozero,ratherthantoanarbitrary value. However,it is more complicatedin that it must individuallyextrapolateeach component of a vector of quantities. SUBROUTINE pzextr(iest,xest,yest,yz,dy,nv) INTEGER iest,nv,IMAX,NMAX REAL xest,dy(nv),yest(nv),yz(nv)PARAMETER (IMAX=13,NMAX=50) Use polynomial extrapolation to evaluate nvfunctions at x=0by fitting a polynomial to a sequence of estimates with progressively smaller values x=xest, and corresponding func- tion vectors yest(1:nv) . This call is number iestin the sequence of calls. Extrapolated function values are output as yz(1:nv) , and their estimated error is output as dy(1:nv) . Parameters: Maximum expected value of iestisIMAX;o fnvisNMAX. INTEGER j,k1REAL delta,f1,f2,q,d(NMAX),qcol(NMAX,IMAX),x(IMAX) SAVE qcol,x x(iest)=xest Save current independent variable. do 11j=1,nv dy(j)=yest(j) yz(j)=yest(j) enddo 11 if(iest.eq.1) then S t o r efi r s te s t i m a t ei nfi r s tc o l u m n . do12j=1,nv qcol(j,1)=yest(j) enddo 12 else do13j=1,nv d(j)=yest(j) enddo 13 do15k1=1,iest-1 delta=1./(x(iest-k1)-xest) f1=xest*deltaf2=x(iest-k1)*delta do 14j=1,nv Propagate tableau 1 diagonal more. q=qcol(j,k1)qcol(j,k1)=dy(j)delta=d(j)-q dy(j)=f1*delta d(j)=f2*deltayz(j)=yz(j)+dy(j) enddo 14 enddo 15 do16j=1,nv qcol(j,iest)=dy(j) enddo 16 endif return END 16.4RichardsonExtrapolationandtheBulirsch-StoerMethod 725Sample 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).Current wisdom favors polynomialextrapolationover rational functionextrap- olationintheBulirsch-Stoermethod. However,ourfeelingisthatthisviewisguidedmore by the kinds of problems used for tests than by one method being actually“better.”Accordingly,we providethe optionalroutine rzextrfor rational function extrapolation, an exact substitution for pzextrabove. SUBROUTINE rzextr(iest,xest,yest,yz,dy,nv) INTEGER iest,nv,IMAX,NMAX REAL xest,dy(nv),yest(nv),yz(nv) PARAMETER (IMAX=13,NMAX=50) Exact substitute for pzextr, but uses diagonal rational function extrapolation instead of polynomial extrapolation. INTEGER j,k REAL b,b1,c,ddy,v,yy,d(NMAX,IMAX),fx(IMAX),x(IMAX)SAVE d,x x(iest)=xest Save current independent variable. if(iest.eq.1) then do 11j=1,nv yz(j)=yest(j) d(j,1)=yest(j) dy(j)=yest(j) enddo 11 else do12k=1,iest-1 fx(k+1)=x(iest-k)/xest enddo 12 do14j=1,nv Evaluate next diagonal in tableau. yy=yest(j)v=d(j,1) c=yy d(j,1)=yydo 13k=2,iest b1=fx(k)*v b=b1-c if(b.ne.0.) then b=(c-v)/bddy=c*b c=b1*b else Care needed to avoid division by 0. ddy=v endif if (k.ne.iest) v=d(j,k)d(j,k)=ddyyy=yy+ddy enddo 13 dy(j)=ddy yz(j)=yy enddo 14 endif returnEND CITED REFERENCES AND FURTHER READING: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §7.2.14. [1] Gear,C.W.1971, NumericalInitialValueProblemsinOrdinaryDifferentialEquations (Englewood Cliffs, NJ: Prentice-Hall), §6.2. Deuflhard, P. 1983, Numerische Mathematik , vol. 41, pp. 399–422. [2] Deuflhard, P. 1985, SIAM Review , vol. 27, pp. 505–535. [3] 726 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).16.5 Second-Order Conservative Equations Usually when you have a system of high-order differential equations to solve it is best to reformulate them as a system of first-order equations, as discussed in §16.0. There is a particular class of equations that occurs quite frequently in practice where you can gainabout a factor of two in ef ficiency by differencing the equations directly. The equations are second-order systems where the derivative does not appear on the right-hand side: y/prime/prime=f(x, y),y (x0)=y0,y/prime(x0)=z0 (16.5.1 ) As usual, ycan denote a vector of values. Stoermer’s rule, dating back to 1907, has been a popular method for discretizing such systems. With h=H/mwe have y1=y0+h[z0+1 2hf(x0,y0)] yk+1−2yk+yk−1=h2f(x0+kh,y k),k =1,...,m −1 zm=(ym−ym−1)/h+1 2hf(x0+H,y m)(16.5.2 ) Herezmisy/prime(x0+H). Henricishowed howtorewriteequations (16.5.2) toreduce roundoff error by using the quantities ∆k≡yk+1−yk. Start with ∆0=h[z0+1 2hf(x0,y0)] y1=y0+∆ 0(16.5.3 ) Then for k=1,...,m −1, set ∆k=∆ k−1+h2f(x0+kh,y k) yk+1=yk+∆ k(16.5.4 ) Finally compute the derivative from zm=∆ m−1/h+1 2hf(x0+H,y m)( 16.5.5 ) Gragg again showed that the error series for equations (16.5.3) –(16.5.5) contains only evenpowersof h,andsothemethodisalogicalcandidateforextrapolation `alaBulirsch-Stoer. We replace mmidby the following routine stoerm: SUBROUTINE stoerm(y,d2y,nv,xs,htot,nstep,yout,derivs) INTEGER nstep,nv,NMAXREAL htot,xs,d2y(nv),y(nv),yout(nv) EXTERNAL derivs PARAMETER (NMAX=50) Maximum value of nv. C USES derivs Stoermer’s rule for integrating y/prime/prime=f(x, y )for a system of n=nv/2equations. On input y(1:nv) contains yin its first nelements and y/primein its second nelements, all evaluated atxs.d2y(1:nv) contains the right-hand side function f(also evaluated at xs)i ni t s firstnelements. Its second nelements are not referenced. Also input is htot, the total step to be taken, and nstep, the number of substeps to be used. The output is returned asyout(1:nv) , with the same storage arrangement as y.derivsis the user-supplied subroutine that calculates f. INTEGER i,n,neqns,nn REAL h,h2,halfh,x,ytemp(NMAX)h=htot/nstep Stepsize this trip. halfh=0.5*h neqns=nv/2 Number of equations. do 11i=1,neqns First step. n=neqns+i ytemp(n)=h*(y(n)+halfh*d2y(i))