f16-6

16.6StiffSetsofEquations 727Sample 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).ytemp(i)=y(i)+ytemp(n) enddo 11 x=xs+hcall derivs(x,ytemp,yout) Useyoutfortemporary storage of derivatives. h2=h*h do 13nn=2,nstep General step. do12i=1,neqns n=neqns+iytemp(n)=ytemp(n)+h2*yout(i) ytemp(i)=ytemp(i)+ytemp(n) enddo 12 x=x+hcall derivs(x,ytemp,yout) enddo 13 do14i=1,neqns Last step. n=neqns+i yout(n)=ytemp(n)/h+halfh*yout(i) yout(i)=ytemp(i) enddo 14 return END Note that for compatibility with bsstepthe arrays yandd2yare of length 2nfor a system of nsecond-order equations. The values of yare stored in the first nelements of y, while the firstderivatives are stored in the second nelements. The right-hand side fis stored in the first nelements of the array d2y; the second nelements are unused. With this storage arrangement you can use bsstepsimply by replacing the call to mmidwith one to stoerm using the same arguments; just be sure that the argument nvofbsstepis set to 2n.Y o u should also use the more efficient sequence of stepsizes suggested by Deuflhard: n=1,2,3,4,5,... (16.5.6 ) and set KMAXX =1 2inbsstep. CITED REFERENCES AND FURTHER READING: Deuflhard, P. 1985, SIAM Review , vol. 27, pp. 505–535. 16.6 Stiff Sets of Equations As soon as one deals with more than one first-order differential equation, the possibility of a stiffset of equations arises. Stiffness occurs in a problem where there are two or more very different scales of the independent variable on whichthe dependent variables are changing. For example, consider the following set of equations [1]: u/prime= 998 u+ 1998 v v/prime=−999u−1999 v(16.6.1 ) with boundary conditions u(0) = 1 v(0) = 0 ( 16.6.2 ) 728 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).xy Figure 16.6.1. Example of an instability encountered in integrating a stiff equation (schematic). Here it is supposed that the equation has two solutions, shown as solid and dashed lines. Although the initial conditions are such as to give the solid solution, the stability of the integration (shown as the unstable dotted sequence of segments) is determined by the more rapidly varying dashed solution, even after thatsolution has effectively died away to zero. Implicit integration methods are the cure. By means of the transformation u=2y−zv =−y+z (16.6.3 ) we find the solution u=2e−x−e−1000 x v=−e−x+e−1000 x(16.6.4 ) If we integrated the system (16.6.1) with any of the methods given so far in this chapter, the presence of the e−1000 xterm would require a stepsize h/lessmuch1/1000for the method to be stable (the reason for this is explained below). This is so eventhoughthe e −1000 xterm is completelynegligiblein determiningthevaluesof uand vas soon as one is away from the origin (see Figure 16.6.1). This is the generic disease of stiff equations: we are required to follow the variation in the solution on the shortest length scale to maintain stability of the integration,even thoughaccuracyrequirementsallow a much larger stepsize. To see how we might cure this problem, consider the single equation y/prime=−cy (16.6.5 ) where c> 0is a constant. The explicit (or forward) Euler scheme for integrating this equation with stepsize his yn+1=yn+hy/prime n=( 1−ch)yn (16.6.6 ) The method is called explicit because the new value yn+1is given explicitly in terms of the old value yn. Clearly the method is unstable if h> 2/c, for then |yn|→∞asn→∞. 16.6StiffSetsofEquations 729Sample 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).Thesimplestcureistoresortto implicitdifferencing,wheretheright-handside is evaluatedat the newylocation. Inthis case, we getthe backwardEuler scheme: yn+1=yn+hy/prime n+1 (16.6.7 ) or yn+1=yn 1+ch(16.6.8 ) The method is absolutely stable: even as h→∞,yn+1→0, which is in fact the correct solution of the differential equation. If we think of xas representing time, thentheimplicitmethodconvergestothetrueequilibriumsolution(i.e.,thesolutionat late times) for large stepsizes. This nice feature of implicit methods holds only forlinearsystems,buteveninthegeneralcaseimplicitmethodsgivebetterstability. Of course, we give up accuracy in following the evolution towards equilibrium if we use large stepsizes, but we maintain stability. These considerations can easily be generalized to sets of linear equations with constant coefficients: y /prime=−C·y (16.6.9 ) whereCis a positive definite matrix. Explicit differencing gives yn+1=(1−Ch)·yn (16.6.10 ) Now a matrix Antends to zero as n→∞only if the largest eigenvalue of A has magnitude less than unity. Thus ynis bounded as n→∞only if the largest eigenvalue of 1−Chis less than 1, or in other words h<2 λmax(16.6.11 ) where λmaxis the largest eigenvalue of C. On the other hand, implicit differencing gives yn+1=yn+hy/prime n+1 (16.6.12 ) or yn+1=(1+Ch)−1·yn (16.6.13 ) If the eigenvalues of Careλ, then the eigenvalues of (1+Ch)−1are (1 + λh)−1, which has magnitude less than one for all h. (Recall that all the eigenvalues of a positivedefinite matrix are nonnegative.) Thus the methodis stable for all stepsizes h. The penalty we pay for this stability is that we are required to invert a matrix at each step. Not all equations are linear with constant coefficients, unfortunately! For the system y/prime=f(y)( 16.6.14 ) implicit differencing gives yn+1=yn+hf(yn+1)( 16.6.15 ) 730 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).Ingeneralthisissomenastysetofnonlinearequationsthathastobesolvediteratively at each step. Suppose we try linearizingthe equations,as in Newton’s method: yn+1=yn+h/bracketleftBigg f(yn)+∂f ∂y/vextendsingle/vextendsingle/vextendsingle/vextendsingley n·(yn+1−yn)/bracketrightBigg (16.6.16 ) Here ∂f/∂yisthematrixofthepartialderivativesoftheright-handside(theJacobian matrix). Rearrange equation (16.6.16) into the form yn+1=yn+h/bracketleftbigg 1−h∂f ∂y/bracketrightbigg−1 ·f(yn)( 16.6.17 ) Ifhis not too big, only one iteration of Newton’s method may be accurate enough to solve equation (16.6.15) using equation (16.6.17). In other words, at each stepwe have to invert the matrix 1−h∂f ∂y(16.6.18 ) to findyn+1. Solving implicit methods by linearization is called a “semi-implicit” method,soequation(16.6.17)isthe semi-implicitEulermethod . Itisnotguaranteed to be stable, but it usually is, because the behavior is locally similar to the case of a constant matrix Cdescribed above. So far we have dealt only with implicit methods that are first-order accurate. While these are veryrobust, most problemswill benefitfromhigher-ordermethods. There are three important classes of higher-ordermethods for stiff systems: •Generalizations of the Runge-Kutta method, of which the most useful are the Rosenbrock methods. The first practical implementation of these ideas was by Kaps and Rentrop, and so these methods are also calledKaps-Rentrop methods. •GeneralizationsoftheBulirsch-Stoermethod,inparticularasemi-implicit extrapolation method due to Bader and Deuflhard. •Predictor-corrector methods, most of which are descendants of Gear’s backward differentiation method. We shall give implementations of the first two methods. Note that systems where the right-hand side depends explicitly on x,f(y,x), can be handled by adding xto the list of dependent variables so that the system to be solved is /parenleftbigg y x/parenrightbigg /prime =/parenleftbigg f 1/parenrightbigg (16.6.19 ) In both the routines to be given in this section, we have explicitly carried out this replacement for you, so the routines can handle right-handsides of the form f(y,x) without any special effort on your part. We nowmentionanimportantpoint: Itis absolutelycrucialtoscaleyourvari- ables properly when integrating stiff problems with automatic stepsize adjustment. As in our nonstiff routines, you will be asked to supply a vector yscalwith which the error is to be scaled. For example, to get constant fractional errors, simply set yscal =|y|. You can get constant absolute errors relative to some maximum values 16.6StiffSetsofEquations 731Sample 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).by setting yscalequal to those maximum values. In stiff problems, there are often strongly decreasing pieces of the solution which you are not particularly interestedin following once they are small. You can control the relative error above some threshold Cand the absolute error below the threshold by setting y scal =m a x (C,|y|)( 16.6.20 ) Ifyouareusingappropriatenondimensionalunits,theneachcomponentof Cshould be of order unity. If you are not sure what values to take for C, simply try setting each component equal to unity. We strongly advocate the choice (16.6.20) for stiff problems. One final warning: Solving stiff problems can sometimes lead to catastrophic precision loss. Be alert for situations where double precision is necessary. Rosenbrock Methods These methods have the advantage of being relatively simple to understand and imple- ment. For moderate accuracies ( /epsilon1<∼10−4–10−5in the error criterion) and moderate-sized systems ( N<∼10), they are competitive with the more complicated algorithms. For more stringent parameters, Rosenbrock methods remain reliable; they merely become less efficientthan competitors like the semi-implicit extrapolation method (see below). A Rosenbrock method seeks a solution of the form y(x 0+h)=y0+s/summationdisplay i=1ciki (16.6.21 ) where the corrections kiare found by solving slinear equations that generalize the structure in (16.6.17): (1−γhf/prime)·ki=hf/parenleftBigg y0+i−1/summationdisplay j=1αijkj/parenrightBigg +hf/prime·i−1/summationdisplay j=1γijkj,i =1,...,s (16.6.22 ) Here we denote the Jacobian matrix by f/prime. The coefficients γ,ci,αij, and γijare fixed constants independent of the problem. If γ=γij=0, this is simply a Runge-Kutta scheme. Equations (16.6.22) can be solved successively for k1,k2,.... Crucial to the success of a stiff integration scheme is an automatic stepsize adjustment algorithm. Kaps and Rentrop [2]discovered an embedded or Runge-Kutta-Fehlberg method asdescribed in §16.2: Twoestimatesoftheform(16.6.21)arecomputed, the“real”one yand a lower-order estimate /hatwideywith different coefficients ˆci,i=1,..., ˆs, where ˆs<sbut theki arethesame. Thedifferencebetween yand/hatwideyleadstoanestimateofthelocaltruncationerror, whichcan thenbe used forstepsizecontrol. Kapsand Rentrop showed thatthesmallestvalueofsfor which embedding is possible is s=4,ˆs=3, leading to a fourth-order method. To minimize the matrix-vector multiplications on the right-hand side of (16.6.22), we rewrite the equations in terms of quantities g i=i−1/summationdisplay j=1γijkj+γki (16.6.23 ) The equations then take the form (1/γh−f/prime)·g1=f(y0) (1/γh−f/prime)·g2=f(y0+a21g1)+c21g1/h (1/γh−f/prime)·g3=f(y0+a31g1+a32g2)+(c31g1+c32g2)/h (1/γh−f/prime)·g4=f(y0+a41g1+a42g2+a43g3)+(c41g1+c42g2+c43g3)/h (16.6.24 ) 732 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).In our implementation stiffof the Kaps-Rentrop algorithm, we have carried out the replacement (16.6.19) explicitly in equations (16.6.24), so you need not concern yourselfabout it. Simply provide a subroutine (called derivsinstiff) that returns f(called dydx) as a function of xandy. Also supply a subroutine jacobnthat returns f /prime(dfdy) and ∂f/∂x (dfdx)asfunctionsof xandy.I fxdoesnotoccurexplicitlyontheright-handside,then dfdx willbezero. UsuallytheJacobianmatrixwillbeavailabletoyoubyanalyticdifferentiationoftheright-handside f. Ifnot,yoursubroutinewillhavetocomputeitbynumericaldifferencing with appropriate increments ∆y. Kaps and Rentrop gave two different sets of parameters, which have slightly different stability properties. Several other sets have been proposed. Our default choice is that ofShampine [3],butwealsogiveyouoneoftheKaps-Rentropsetsasanoption. Someproposed parameter sets require function evaluations outside the domain of integration; we prefer toavoid that complication. The calling sequence of stiffis exactly the same as the nonstiff routines given earlier in this chapter. It isthus “plug-compatible” with them in the general ODEintegrating routine odeint. This compatibility requires, unfortunately, one slight anomaly: While the user- supplied routine derivsis a dummy argument (which can therefore have any actual name), the other user-supplied routine is notan argument and must be named (exactly) jacobn. stiffbegins by saving the initial values, in case the step has to be repeated because the error tolerance is exceeded. The linear equations (16.6.24) are solved by first computingtheLUdecomposition of the matrix 1/γh−f /primeusing the routine ludcmp. Then the four giare found by back-substitution of the four different right-hand sides using lubksb. Note that each step of the integration requires one call to jacobnand three calls to derivs(one call to get dydxbefore calling stiff, and two calls inside stiff). The reason only three calls are needed and not four is that the parameters have been chosen so that the last twocalls in equation (16.6.24) are done with the same arguments. Counting the evaluation ofthe Jacobian matrix as roughly equivalent to Nevaluations of the right-hand side f, we see that the Kaps-Rentrop scheme involves about N+3function evaluations per step. Note that ifNis large and the Jacobian matrix is sparse, you should replace the LUdecomposition by a suitable sparse matrix procedure. Stepsize control depends on the fact that y exact =y+O(h5) yexact =/hatwidey+O(h4)(16.6.25 ) Thus |y−/hatwidey|=O(h4)( 16.6.26 ) Referring back to the steps leading from equation (16.2.4) to equation (16.2.10), we see that the new stepsize should be chosen as in equation (16.2.10) but with the exponents 1/4and 1/5 replaced by 1/3 and 1/4, respectively. Also, experience shows that it is wise toprevent too large a stepsize change in one step, otherwise we will probably have to undo the large change in the next step. We adopt 0.5 and 1.5 as the maximum allowed decrease and increase of hin one step. SUBROUTINE stiff(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs) INTEGER n,NMAX,MAXTRY REAL eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n),SAFETY,GROW, * PGROW,SHRNK,PSHRNK,ERRCON,GAM,A21,A31,A32,A2X,A3X,C21,* C31,C32,C41,C42,C43,B1,B2,B3,B4,E1,E2,E3,E4,C1X,C2X,C3X, * C4X EXTERNAL derivsPARAMETER (NMAX=50,SAFETY=0.9,GROW=1.5,PGROW=-.25, * SHRNK=0.5,PSHRNK=-1./3.,ERRCON=.1296,MAXTRY=40) PARAMETER (GAM=1./2.,A21=2.,A31=48./25.,A32=6./25.,C21=-8., * C31=372./25.,C32=12./5.,C41=-112./125.,C42=-54./125.,* C43=-2./5.,B1=19./9.,B2=1./2.,B3=25./108.,B4=125./108., * E1=17./54.,E2=7./36.,E3=0.,E4=125./108.,C1X=1./2., 16.6StiffSetsofEquations 733Sample 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).* C2X=-3./2.,C3X=121./50.,C4X=29./250.,A2X=1.,A3X=3./5.) C USES derivs,jacobn,lubksb,ludcmp Fourth-order Rosenbrock step for integrating stiff o.d.e.’s, with monitoring of local trun- cation error to 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 thestepsizetobeattempted htry,therequiredaccuracy eps,andthevector yscal(1:n) againstwhichtheerrorisscaled. Onoutput, yandxarereplacedbytheirnewvalues, hdid is the stepsize that was actually accomplished, and hnextis the estimated next stepsize. derivsis a user-supplied subroutine that computes the derivatives of the right-hand side withrespectto x,while jacobn(afixedname)isauser-suppliedsubroutinethatcomputes theJacobimatrixofderivativesoftheright-handsidewithrespecttothecomponentsof y. Parameters: NMAXisthemaximumvalueof n;GROWandSHRNKarethelargestandsmallest factorsbywhichstepsizecanchangeinonestep; ERRCON=(GROW/SAFETY)**(1/PGROW) and handles the case when errmax /similarequal0. INTEGER i,j,jtry,indx(NMAX)REAL d,errmax,h,xsav,a(NMAX,NMAX),dfdx(NMAX),dfdy(NMAX,NMAX), * dysav(NMAX),err(NMAX),g1(NMAX),g2(NMAX),g3(NMAX), * g4(NMAX),ysav(NMAX) xsav=x Save initial values. do 11i=1,n ysav(i)=y(i) dysav(i)=dydx(i) enddo 11 call jacobn(xsav,ysav,dfdx,dfdy,n,NMAX) Theusermustsupplythissubroutine toreturnthe n-by-nmatrix dfdyandthevector dfdx. h=htry Set stepsize to the initial trial value. do23jtry=1,MAXTRY do13i=1,n Set up the matrix 1−γhf/prime. do12j=1,n a(i,j)=-dfdy(i,j) enddo 12 a(i,i)=1./(GAM*h)+a(i,i) enddo 13 call ludcmp(a,n,NMAX,indx,d) LU decomposition of the matrix. do14i=1,n Set up right-hand side for g1. g1(i)=dysav(i)+h*C1X*dfdx(i) enddo 14 call lubksb(a,n,NMAX,indx,g1) Solve for g1. do15i=1,n Compute intermediate values of yandx. y(i)=ysav(i)+A21*g1(i) enddo 15 x=xsav+A2X*hcall derivs(x,y,dydx) Compute dydxat the intermediate values. do 16i=1,n Set up right-hand side for g2. g2(i)=dydx(i)+h*C2X*dfdx(i)+C21*g1(i)/h enddo 16 call lubksb(a,n,NMAX,indx,g2) Solve for g2. do17i=1,n Compute intermediate values of yandx. y(i)=ysav(i)+A31*g1(i)+A32*g2(i) enddo 17 x=xsav+A3X*hcall derivs(x,y,dydx) Compute dydxat the intermediate values. do 18i=1,n Set up right-hand side for g3. g3(i)=dydx(i)+h*C3X*dfdx(i)+(C31*g1(i)+ * C32*g2(i))/h enddo 18 call lubksb(a,n,NMAX,indx,g3) Solve for g3. do19i=1,n Set up right-hand side for g4. g4(i)=dydx(i)+h*C4X*dfdx(i)+(C41*g1(i)+ * C42*g2(i)+C43*g3(i))/h enddo 19 call lubksb(a,n,NMAX,indx,g4) Solve for g4. do21i=1,n Get fourth-order estimate of yand errorestimate. y(i)=ysav(i)+B1*g1(i)+B2*g2(i)+B3*g3(i)+B4*g4(i) 734 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).err(i)=E1*g1(i)+E2*g2(i)+E3*g3(i)+E4*g4(i) enddo 21 x=xsav+hif(x.eq.xsav)pause ’stepsize not significant in stiff’errmax=0. Evaluate accuracy. do 22i=1,n errmax=max(errmax,abs(err(i)/yscal(i))) enddo 22 errmax=errmax/eps Scale relative to required tolerance. if(errmax.le.1.)then Step succeeded. Compute size of next step and re- turn. hdid=h if(errmax.gt.ERRCON)then hnext=SAFETY*h*errmax**PGROW else hnext=GROW*h endif return else Truncation error too large, reduce stepsize. hnext=SAFETY*h*errmax**PSHRNKh=sign(max(abs(hnext),SHRNK*abs(h)),h) endif enddo 23 Go back and re-try step. pause ’exceeded MAXTRY in stiff’ END Here are the Kaps-Rentrop parameters, which can be substituted for those of Shampine simply by replacing the PARAMETER statement: PARAMETER (GAM=.231,A21=2.,A31=4.52470820736,A32=4.16352878860, * C21=-5.07167533877,C31=6.02015272865,C32=.159750684673, * C41=-1.856343618677,C42=-8.50538085819,C43= * -2.08407513602,B1=3.95750374663,B2=4.62489238836,B3=* .617477263873,B4=1.282612945268,E1=-2.30215540292, * E2=-3.07363448539,E3=.873280801802,E4=1.282612945268, * C1X=GAM,C2X=-.396296677520e-01,C3X=.550778939579,* C4X=-.553509845700e-01,A2X=.462,A3X=.880208333333) As an example of how stiffis used, one can solve the system y/prime 1=−.013y1−1000y1y3 y/prime 2=−2500y2y3 y/prime 3=−.013y1−1000y1y3−2500y2y3(16.6.27 ) with initial conditions y1(0) = 1 ,y 2(0) = 1 ,y 3(0) = 0 ( 16.6.28 ) (This istest problem D4 in [4].) Weintegrate the system up to x=5 0with an initialstepsize ofh=2.9×10−4using odeint. The components of Cin (16.6.20) are all set to unity. The routines derivsandjacobnfor this problem are given below. Even though the ratio of largest to smallest decay constants for this problem is around 106,stiffsucceeds in integrating this set in only 29 steps with /epsilon1=1 0−4. By contrast, the Runge-Kutta routine rkqsrequires 51,012 steps! SUBROUTINE jacobn(x,y,dfdx,dfdy,n,nmax) INTEGER n,nmax,i REAL x,y(*),dfdx(*),dfdy(nmax,nmax) do11i=1,3 dfdx(i)=0. enddo 11 16.6StiffSetsofEquations 735Sample 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).dfdy(1,1)=-.013-1000.*y(3) dfdy(1,2)=0. dfdy(1,3)=-1000.*y(1) dfdy(2,1)=0.dfdy(2,2)=-2500.*y(3) dfdy(2,3)=-2500.*y(2) dfdy(3,1)=-.013-1000.*y(3)dfdy(3,2)=-2500.*y(3)dfdy(3,3)=-1000.*y(1)-2500.*y(2) return END SUBROUTINE derivs(x,y,dydx) REAL x,y(*),dydx(*) dydx(1)=-.013*y(1)-1000.*y(1)*y(3)dydx(2)=-2500.*y(2)*y(3) dydx(3)=-.013*y(1)-1000.*y(1)*y(3)-2500.*y(2)*y(3) returnEND Semi-implicitExtrapolation Method TheBulirsch-Stoermethod,whichdiscretizesthedifferentialequationusingthemodified midpoint rule, does not work for stiff problems. Bader and Deuflhard [5]discovered a semi- implicit discretization that works very well and that lends itself to extrapolation exactly asin the original Bulirsch-Stoer method. The starting point is an implicit form of the midpoint rule: y n+1−yn−1=2hf/parenleftbiggyn+1+yn−1 2/parenrightbigg (16.6.29 ) Convert this equation into semi-implicit form by linearizing the right-hand side about f(yn). The result is the semi-implicit midpoint rule : /bracketleftbigg 1−h∂f ∂y/bracketrightbigg ·yn+1=/bracketleftbigg 1+h∂f ∂y/bracketrightbigg ·yn−1+2h/bracketleftbigg f(yn)−∂f ∂y·yn/bracketrightbigg (16.6.30 ) It is used with a special first step, the semi-implicit Euler step (16.6.17), and a special “smoothing” last step in which the last ynis replaced by yn≡1 2(yn+1+yn−1)( 16.6.31 ) Bader and Deuflhard showed that the error series for this method once again involves only even powers of h. Forpracticalimplementation,itisbettertorewritetheequationsusing ∆k≡yk+1−yk. With h=H/m, start by calculating ∆0=/bracketleftbigg 1−h∂f ∂y/bracketrightbigg−1 ·hf(y0) y1=y0+∆ 0(16.6.32 ) Then for k=1,...,m −1, set ∆k=∆ k−1+2/bracketleftbigg 1−h∂f ∂y/bracketrightbigg−1 ·[hf(yk)−∆k−1] yk+1=yk+∆ k(16.6.33 ) 736 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).Finally compute ∆m=/bracketleftbigg 1−h∂f ∂y/bracketrightbigg−1 ·[hf(ym)−∆m−1] ym=ym+∆ m(16.6.34 ) Itiseasy toincorporate thereplacement (16.6.19) intheabove formulas. Theadditional terms in the Jacobian that come from ∂f/∂xall cancel out of the semi-implicit midpoint rule (16.6.30). In the special first step (16.6.17), and in the corresponding equation (16.6.32), thetermhfbecomes hf+h 2∂f/∂x. The remaining equations are all unchanged. This algorithm is implemented in the routine simpr: SUBROUTINE simpr(y,dydx,dfdx,dfdy,nmax,n,xs,htot,nstep,yout, * derivs) INTEGER n,nmax,nstep,NMAXX REAL htot,xs,dfdx(n),dfdy(nmax,nmax),dydx(n),y(n),yout(n) EXTERNAL derivsPARAMETER (NMAXX=50) Maximum expected value of n. C USES derivs,lubksb,ludcmp Performsonestepofsemi-implicitmidpointrule. Inputarethedependentvariable y(1:n), itsderivative dydx(1:n),thederivativeoftheright-handsidewithrespectto x,dfdx(1:n), and the Jacobian dfdy(1:nmax,1:nmax) atxs. Also input are htot, the total step to be taken, and nstep, the number of substeps to be used. The output is returned as yout(1:n).derivsis the user-supplied subroutine that calculates dydx. INTEGER i,j,nn,indx(NMAXX)REAL d,h,x,a(NMAXX,NMAXX),del(NMAXX),ytemp(NMAXX) h=htot/nstep Stepsize this trip. do 12i=1,n Set up the matrix 1−hf/prime. do11j=1,n a(i,j)=-h*dfdy(i,j) enddo 11 a(i,i)=a(i,i)+1. enddo 12 call ludcmp(a,n,NMAXX,indx,d) LU decomposition of the matrix. do13i=1,n Set up right-hand side for first step. Use youtfor temporary storage. yout(i)=h*(dydx(i)+h*dfdx(i)) enddo 13 call lubksb(a,n,NMAXX,indx,yout) do14i=1,n First step. del(i)=yout(i) ytemp(i)=y(i)+del(i) enddo 14 x=xs+hcall derivs(x,ytemp,yout) Useyoutfortemporary storage of derivatives. do 17nn=2,nstep General step. do15i=1,n Set up right-hand side for general step. yout(i)=h*yout(i)-del(i) enddo 15 call lubksb(a,n,NMAXX,indx,yout) do16i=1,n del(i)=del(i)+2.*yout(i) ytemp(i)=ytemp(i)+del(i) enddo 16 x=x+h call derivs(x,ytemp,yout) enddo 17 do18i=1,n Set up right-hand side for last step. yout(i)=h*yout(i)-del(i) enddo 18 call lubksb(a,n,NMAXX,indx,yout) do19i=1,n Take last step. yout(i)=ytemp(i)+yout(i) 16.6StiffSetsofEquations 737Sample 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 19 return END The routine simpris intended to be used in a routine stifbsthat is almost exactly the same as bsstep. The only differences are: •The stepsize sequence is n=2,6,10,14,22,34,50,..., (16.6.35 ) where each member differs from its predecessor by the smallest multiple of 4 that makes theratioofsuccessivetermsbe ≤5 7. Theparameter KMAXXistakentobe7. •The work per unit step now includes the cost of Jacobian evaluations as well as function evaluations. We count one Jacobian evaluation as equivalent to N function evaluations, where Nis the number of equations. •Onceagaintheuser-suppliedroutine derivsisadummyargumentandsocanhave any name. However, to maintain “plug-compatibility” with rkqs,bsstepand stiff,theroutine jacobnisnot anargument and musthave exactly thisname. It iscalledoncepersteptoreturn f/prime(dfdy)and∂f/∂x(dfdx)asfunctionsof xandy. Here is the routine, with comments pointing out only the differences from bsstep: SUBROUTINE stifbs(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 EXTERNAL derivs PARAMETER (NMAX=50,KMAXX=7,IMAX=KMAXX+1,SAFE1=.25,SAFE2=.7, * REDMAX=1.e-5,REDMIN=.7,TINY=1.e-30,SCALMX=.1) C USES derivs,jacobn,simpr,pzextr Semi-implicit extrapolation step forintegrating stiff o.d.e.’s, with monitoring of local trun-cation error to adjust stepsize. Input are the dependent variable vector y(1:nv)and its derivative dydx(1:nv) atthe starting valueofthe independent variable x. Also input are thestepsizetobeattempted htry,therequiredaccuracy eps,andthevector yscal(1:nv) againstwhichtheerrorisscaled. Onoutput, yandxarereplacedbytheirnewvalues, hdid is the stepsize that was actually accomplished, and hnextis the estimated next stepsize. derivsis a user-supplied subroutine that computes the derivatives of the right-hand side withrespectto x,while jacobn(afixedname)isauser-suppliedsubroutinethatcomputes theJacobimatrixofderivativesoftheright-handsidewithrespecttothecomponentsof y. Be sure to set htryon successive steps to the valueof hnextreturned from the previous step, as is the case if the routine is called by odeint. INTEGER i,iq,k,kk,km,kmax,kopt,nvold,nseq(IMAX) REAL eps1,epsold,errmax,fact,h,red,scale,work,wrkmin,xest,xnew, * a(IMAX),alf(KMAXX,KMAXX),dfdx(NMAX),dfdy(NMAX,NMAX),* err(KMAXX),yerr(NMAX),ysav(NMAX),yseq(NMAX) LOGICAL first,reduct SAVE a,alf,epsold,first,kmax,kopt,nseq,nvold,xnewDATA first/.true./,epsold/-1./,nvold/-1/ DATA nseq /2,6,10,14,22,34,50,70/ Sequence is different from bsstep. if(eps.ne.epsold.or.nv.ne.nvold)then Reinitialize alsoif nvhas changed. hnext=-1.e29xnew=-1.e29 eps1=SAFE1*eps a(1)=nseq(1)+1do 11k=1,KMAXX a(k+1)=a(k)+nseq(k+1) enddo 11 do13iq=2,KMAXX 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 738 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).epsold=eps nvold=nv Save nv. a(1)=nv+a(1) AddcostofJacobianevaluationstoworkco- efficients. do14k=1,KMAXX a(k+1)=a(k)+nseq(k+1) enddo 14 do15kopt=2,KMAXX-1 if(a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt))goto 1 enddo 15 1 kmax=kopt endifh=htrydo 16i=1,nv ysav(i)=y(i) enddo 16 call jacobn(x,y,dfdx,dfdy,nv,nmax) Evaluate Jacobian. if(h.ne.hnext.or.x.ne.xnew)then first=.true.kopt=kmax endif reduct=.false. 2d o 18k=1,kmax xnew=x+h if(xnew.eq.x)pause ’stepsize underflow in stifbs’ call simpr(ysav,dydx,dfdx,dfdy,nmax,nv,x,h,nseq(k),yseq, * derivs) Semi-implicit midpoint rule. xest=(h/nseq(k))**2 Therestoftheroutineisidenticalto bsstep. call pzextr(k,xest,yseq,y,yerr,nv) if(k.ne.1)then errmax=TINY do17i=1,nv errmax=max(errmax,abs(yerr(i)/yscal(i))) enddo 17 errmax=errmax/eps km=k-1 err(km)=(errmax/SAFE1)**(1./(2*km+1)) endifif(k.ne.1.and.(k.ge.kopt-1.or.first))then if(errmax.lt.1.)goto 4 if(k.eq.kmax.or.k.eq.kopt+1)then 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 18 3 red=min(red,REDMIN) red=max(red,REDMAX)h=h*redreduct=.true. goto 2 16.6StiffSetsofEquations 739Sample 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).4 x=xnew hdid=h first=.false. wrkmin=1.e35do 19kk=1,km fact=max(err(kk),SCALMX) work=fact*a(kk+1)if(work.lt.wrkmin)then scale=fact wrkmin=work kopt=kk+1 endif enddo 19 hnext=h/scaleif(kopt.ge.k.and.kopt.ne.kmax.and..not.reduct)then fact=max(scale/alf(kopt-1,kopt),SCALMX) if(a(kopt+1)*fact.le.wrkmin)then hnext=h/factkopt=kopt+1 endif endif returnEND The routine stifbsis an excellent routine for all stiff problems, competitive with the best Gear-type routines. stiffis comparable in execution time for moderate Nand /epsilon1<∼10−4. By the time /epsilon1∼10−8,stifbsis roughly an order of magnitude faster. There are further improvements that could be applied to stifbsto make it even more robust. For example, very occasionally ludcmpinsimprwill encounter a singular matrix. You could arrange for the stepsize to be reduced, say by a factor of the current nseq(k). There are also certain stability restrictions on the stepsize that come into play on some problems. Fora discussion of how to implement these automatically, see [6]. CITED REFERENCES AND FURTHER READING: Gear,C.W.1971, NumericalInitialValueProblemsinOrdinaryDifferentialEquations (Englewood Cliffs, NJ: Prentice-Hall). [1] Kaps, P., and Rentrop, P. 1979, Numerische Mathematik , vol. 33, pp. 55–68. [2] Shampine, L.F. 1982, ACM Transactions on Mathematical Software , vol. 8, pp. 93–113. [3] Enright, W.H., and Pryce, J.D. 1987, ACM Transactions on Mathematical Software , vol. 13, pp. 1–27. [4] Bader, G., and Deuflhard, P. 1983, Numerische Mathematik , vol. 41, pp. 373–398. [5] Deuflhard, P. 1983, Numerische Mathematik , vol. 41, pp. 399–422. Deuflhard, P. 1985, SIAM Review , vol. 27, pp. 505–535. Deuflhard, P. 1987, “Uniqueness Theorems for Stiff ODE Initial Value Problems,” Preprint SC- 87-3(Berlin: Konrad Zuse Zentrum f¨ ur Informationstechnik). [6] Enright, W.H., Hull, T.E., and Lindberg, B. 1975, BIT, vol. 15, pp. 10–48. Wanner,G.1988,in NumericalAnalysis1987 ,PitmanResearchNotesinMathematics,vol.170, D.F. Griffiths and G.A. Watson, eds. (Harlow, Essex, U.K.: Longman Scientific and Tech-nical). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag). 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,y n), 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