732 Chapter 16. Integration of Ordinary Differential Equations dy[j]=yest[j]; else for (k=1;k<iest;k++) fx[k+1]-x[iest-k]/xest; for (j=1;j<=nv;j++) Evaluate next diagonal in tableau. vd[j][1]: d[j][1]=yy=c=yest [j]; for (k=2;k<=iest;k++) bi=fx[k]*v; b=b1-c; 1f(b)[ b=(c-v)/b; ddy=c*b; c=b1*b; 83 else Care needed to avoid division by 0. granted for 19881992 ddy=v; if (k !iest)v=d[i][k]; (including this one) d[j][k]-ddy; yy +=ddy; 111800.872 Cambridge n NUMERICAL RECIPES dy[j]=ddy; yz[j]=yy; free_vector(fx,1,iest); (North America 州bMe se to make one paper University Press. THE ART CITED REFERENCES AND FURTHER READING: ictly pror Programs Stoer.J..and Bulirsch.R.1980.Introduction to Numerical Analysis(New York:Springer-Verlag). 57.2.14.[1] Gear,C.W.1971,Numerical Initial Value Problems in Ordinary Differential Equations(Englewood 6 Cliffs,NJ:Prentice-Hall),86.2. Deufhard,P.1983,Numerische Mathematik,vol.41,pp.399-422.[2] Deuflhard,P.1985,SIAM Review,vol.27,pp.505-535.[3] OF SCIENTIFIC COMPUTING(ISBN 1920 16.5 Second-Order Conservative Equations Numerical 1021 Recipes 43108 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 816.0.There is (outside a particular class of equations that occurs quite frequently in practice where you can gain Software. about a factor of two in efficiency by differencing the equations directly.The equations are North second-order systems where the derivative does not appear on the right-hand side: y”=(x,y),y(xo)=0,y(xo)=20 (16.5.1) As usual,y can denote a vector of values. Stoermer's rule,dating back to 1907,has been a popular method for discretizing such systems.With h=H/m we have 1=0+h[20+hf(xo,o] 张+1-2k+k-1=2f(xo+kh,),k=1,,m-1 (16.5.2) 2m =(9m -ym-1)/h+hf(o+H,ym)
732 Chapter 16. Integration of Ordinary Differential Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). dy[j]=yest[j]; } else { for (k=1;k<iest;k++) fx[k+1]=x[iest-k]/xest; for (j=1;j<=nv;j++) { Evaluate next diagonal in tableau. v=d[j][1]; d[j][1]=yy=c=yest[j]; for (k=2;k<=iest;k++) { b1=fx[k]*v; b=b1-c; if (b) { b=(c-v)/b; ddy=c*b; c=b1*b; } else Care needed to avoid division by 0. ddy=v; if (k != iest) v=d[j][k]; d[j][k]=ddy; yy += ddy; } dy[j]=ddy; yz[j]=yy; } } free_vector(fx,1,iest); } CITED REFERENCES AND FURTHER READING: Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), §7.2.14. [1] Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (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] 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 gain about a factor of two in efficiency by differencing the equations directly. The equations are second-order systems where the derivative does not appear on the right-hand side: y = f(x, y), y(x0) = y0, y (x0) = z0 (16.5.1) As usual, y can denote a vector of values. Stoermer’s rule, dating back to 1907, has been a popular method for discretizing such systems. With h = H/m we have y1 = y0 + h[z0 + 1 2 hf(x0, y0)] yk+1 − 2yk + yk−1 = h2 f(x0 + kh, yk), k = 1,...,m − 1 zm = (ym − ym−1)/h + 1 2 hf(x0 + H, ym) (16.5.2)
16.5 Second-Order Conservative Equations 733 Here zm is y(zo+H).Henrici showed how to rewrite equations(16.5.2)to reduce roundoff error by using the quantities Ak=y-yk.Start with △o=hz0+5hf(ro,yo] (16.5.3) 1=0+△0 Then for 1,...,m-1,set △k=△k-1+h2f(xo+kh,yk) (16.5.4) k+1=k十△k 83 Finally compute the derivative from zm=△m-1/h+5hf(xo+H,ym） (16.5.5) Gragg again showed that the error series for equations (16.5.3)-(16.5.5)contains only even powers of h,and so the method is a logical candidate for extrapolationa la Bulirsch-Stoer. gg2侯 We replace mmid by the following routine stoerm: (Nort server 6 #include "nrutil.h" America computer one paper University Press. THE void stoerm(float y[],float d2y[],int nv,float xs,float htot,int nstep, float yout[],void (*derivs)(float,float []float []) ART Stoermer's rule for integrating y"=f(,y)for a system of n=nv/2 equations.On input 是 y[1..nv]contains y in its first n elements and y'in its second n elements,all evaluated at Programs xs.d2y [1..nv]contains the right-hand side function f (also evaluated at xs)in its first n elements.Its second n elements 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 as yout [1..nv], with the same storage arrangement as y.derivs is the user-supplied routine that calculates f. { int i,n,negns,nn; float h,h2,halfh,x,*ytemp; ectcustser OF SCIENTIFIC COMPUTING(ISBN ytemp=vector(1,nv); h=htot/nstep; Stepsize this trip. halfh=0.5*h; negns=nv/2; Number of equations. for (i=1;i<=negns;i++){ First step. 10-621 n=negns+i; ytemp[i]=y [i]+(ytemp[n]=h*(y [n]+halfh*d2y [i])); 1988-1992 by Numerical Recipes 43108 x=xs+h; (*derivs)(x,ytemp,yout); Use yout for temporary storage of derivatives. h2-h*h; (outside for (nn=2;nn<=nstep;nn++){ General step. Software. for (i=1;i<=neqns;i++) ytemp[i]+(ytemp[(n=neqns+i)]+h2*yout[i]); x+=h; (*derivs)(x,ytemp,yout); for (i=1;i<=negns;i++){ Last step n=negns+i; yout [n]=ytemp[n]/h+halfh*yout [i]; yout [i]=ytemp[i]; free_vector(ytemp,1,nv);
16.5 Second-Order Conservative Equations 733 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Here zm is y (x0 + H). Henrici showed how to rewrite equations (16.5.2) to reduce roundoff error by using the quantities ∆k ≡ yk+1 − yk. Start with ∆0 = h[z0 + 1 2 hf(x0, y0)] y1 = y0 + ∆0 (16.5.3) Then for k = 1,...,m − 1, set ∆k = ∆k−1 + h2 f(x0 + kh, yk) yk+1 = yk + ∆k (16.5.4) Finally compute the derivative from zm = ∆m−1/h + 1 2 hf(x0 + H, ym) (16.5.5) Gragg again showed that the error series for equations (16.5.3)–(16.5.5) contains only even powers of h, and so the method is a logical candidate for extrapolationa la Bulirsch-Stoer. ` We replace mmid by the following routine stoerm: #include "nrutil.h" void stoerm(float y[], float d2y[], int nv, float xs, float htot, int nstep, float yout[], void (*derivs)(float, float [], float [])) Stoermer’s rule for integrating y = f(x, y) for a system of n = nv/2 equations. On input y[1..nv] contains y in its first n elements and y in its second n elements, all evaluated at xs. d2y[1..nv] contains the right-hand side function f (also evaluated at xs) in its first n elements. Its second n elements 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 as yout[1..nv], with the same storage arrangement as y. derivs is the user-supplied routine that calculates f. { int i,n,neqns,nn; float h,h2,halfh,x,*ytemp; ytemp=vector(1,nv); h=htot/nstep; Stepsize this trip. halfh=0.5*h; neqns=nv/2; Number of equations. for (i=1;i<=neqns;i++) { First step. n=neqns+i; ytemp[i]=y[i]+(ytemp[n]=h*(y[n]+halfh*d2y[i])); } x=xs+h; (*derivs)(x,ytemp,yout); Use yout for temporary storage of derivatives. h2=h*h; for (nn=2;nn<=nstep;nn++) { General step. for (i=1;i<=neqns;i++) ytemp[i] += (ytemp[(n=neqns+i)] += h2*yout[i]); x += h; (*derivs)(x,ytemp,yout); } for (i=1;i<=neqns;i++) { Last step. n=neqns+i; yout[n]=ytemp[n]/h+halfh*yout[i]; yout[i]=ytemp[i]; } free_vector(ytemp,1,nv); }
734 Chapter 16.Integration of Ordinary Differential Equations Note that for compatibility with bsstep the arrays y and d2y are of length 2n for a system of n second-order equations.The values of y are stored in the first n elements of y, while the first derivatives are stored in the second n elements.The right-hand side f is stored in the first n elements of the array d2y,the second n elements are unused.With this storage arrangement you can use bsstep simply by replacing the call to mmid with one to stoerm using the same arguments;just be sure that the argument nv of bsstep is set to 2n.You should also use the more efficient sequence of stepsizes suggested by Deufhard: n=1,2,3,4,5,… (16.5.6) and set KMAXX 12 in bsstep. ://www.nr CITED REFERENCES AND FURTHER READING: Deuflhard,P.1985,SIAM Review,vol.27,pp.505-535. 八g2 1.800 NUMERICAL RECIPES 16.6 Stiff Sets of Equations (North America to make one paper server computer, University Press. THE As soon as one deals with more than one first-order differential equation,the ART possibility of a stiff set of equations arises.Stiffness occurs in a problem where there are two or more very different scales of the independent variable on which 9 Program the dependent variables are changing.For example,consider the following set of equations [1]: 6 、三 OF SCIENTIFIC( u'=998u+1998U (16.6.1) v'=-999u-1999u with boundary conditions 10621 u(0)=1 u(0)=0 (16.6.2) Numerica 431 By means of the transformation Recipes u=2y-z v=-y+z (16.6.3) North we find the solution u=2e-z-e-1000z v=-e-z+e-1000z (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 term would require a stepsize h1/1000 for the method to be stable (the reason for this is explained below).This is so even
734 Chapter 16. Integration of Ordinary Differential Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Note that for compatibility with bsstep the arrays y and d2y are of length 2n for a system of n second-order equations. The values of y are stored in the first n elements of y, while the first derivatives are stored in the second n elements. The right-hand side f is stored in the first n elements of the array d2y; the second n elements are unused. With this storage arrangement you can use bsstep simply by replacing the call to mmid with one to stoerm using the same arguments; just be sure that the argument nv of bsstep is set to 2n. You should also use the more efficient sequence of stepsizes suggested by Deuflhard: n = 1, 2, 3, 4, 5,... (16.5.6) and set KMAXX = 12 in bsstep. 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 stiff set of equations arises. Stiffness occurs in a problem where there are two or more very different scales of the independent variable on which the dependent variables are changing. For example, consider the following set of equations [1]: u = 998u + 1998v v = −999u − 1999v (16.6.1) with boundary conditions u(0) = 1 v(0) = 0 (16.6.2) By means of the transformation u = 2y − z v = −y + z (16.6.3) we find the solution u = 2e−x − e−1000x v = −e−x + e−1000x (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−1000x term would require a stepsize h 1/1000 for the method to be stable (the reason for this is explained below). This is so even