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 machine￾readable 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); }