f16-5

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 efficiency by differencing the equations directly. The equations aresecond-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 ) Here zmisy/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 first nelements. 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.derivs is 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)) 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) Useyout for temporary 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=−999 u−1999 v(16.6.1 ) with boundary conditions u(0) = 1 v(0) = 0 ( 16.6.2 )