16.6 Stiff Sets of Equations 743 Convert this equation into semi-implicit form by linearizing the right-hand side about f(y). The result is the semi-implicit midpoint rule: -+2h.)-】 af of yn+1=1+h (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 y is replaced by 了n三(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. For practical implementation,it is better to rewrite the equations using=-y 8 With h=H/m,start by calculating 4=1-hr] oy .hf(yo) (16.6.32) y1=yo+△0 Then for 1,...,m-1,set 室2 6 △k=△k-1+21-h .[hfy)-△k-] (16.6.33) yk+1=yk+△k 8”, America Press. Finally compute △m=1- f- 9 ·[hfym）-△m-i] Programs (16.6.34) 了m=ym+△m It is easy to incorporate the replacement(16.6.19)in the above formulas.The additional terms in the Jacobian that come from of/o all 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),the term hf becomes hf+hof/Or.The remaining equations are all unchanged. 豪 SCIENTIFIC COMPUTING(ISBN This algorithm is implemented in the routine simpr: 19891892 #include "nrutil.h" void simpr(float y[],float dydx],float dfdx],float **dfdy,int n, float xs,float htot,int nstep,float yout☐， void (*derivs)(float,float []float []) Numerical Recipes 10-621 43108 Performs one step of semi-implicit midpoint rule.Input are the dependent variable y [1..n],its derivative dydx[1..n],the derivative of the right-hand side with respect to x,dfdx[1..n] and the Jacobian dfdy [1..n][1..n]at xs.Also input are htot,the total step to be taken. (outside and nstep,the number of substeps to be used.The output is returned as yout [1..n] Software. derivs is the user-supplied routine that calculates dydx. void lubksb(float **a,int n,int *indx,float b[]); Ame void ludcmp(float **a,int n,int *indx,float *d); int i,i,nn,*indx; float d,h,x,**a,*del,*ytemp; indx=ivector(1,n); a=matrix(1,n,1,n); del=vector(1,n); ytemp=vector(1,n) h=htot/nstep; Stepsize this trip. for(1=1:1<n;it+){ Set up the matrix 1-hf. for (j=1;j<=n;j++)a[i][j]=-h*dfdy[i][j];16.6 Stiff Sets of Equations 743 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). Convert this equation into semi-implicit form by linearizing the right-hand side about f(yn). The result is the semi-implicit midpoint rule:  1 − h ∂f ∂y  · yn+1 =  1 + h ∂f ∂y  · yn−1 + 2h  f(yn) − ∂f ∂y · yn  (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 yn is 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. For practical implementation, it is better to rewrite the equations using ∆k ≡ yk+1 −yk. With h = H/m, start by calculating ∆0 =  1 − h ∂f ∂y −1 · hf(y0) y1 = y0 + ∆0 (16.6.32) Then for k = 1,...,m − 1, set ∆k = ∆k−1 + 2  1 − h ∂f ∂y −1 · [hf(yk) − ∆k−1] yk+1 = yk + ∆k (16.6.33) Finally compute ∆m =  1 − h ∂f ∂y −1 · [hf(ym) − ∆m−1] ym = ym + ∆m (16.6.34) It is easy to incorporate the replacement (16.6.19) in the above formulas. The additional terms in the Jacobian that come from ∂f/∂x all 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), the term hf becomes hf + h2∂f/∂x. The remaining equations are all unchanged. This algorithm is implemented in the routine simpr: #include "nrutil.h" void simpr(float y[], float dydx[], float dfdx[], float **dfdy, int n, float xs, float htot, int nstep, float yout[], void (*derivs)(float, float [], float [])) Performs one step of semi-implicit midpoint rule. Input are the dependent variable y[1..n], its derivative dydx[1..n], the derivative of the right-hand side with respect to x, dfdx[1..n], and the Jacobian dfdy[1..n][1..n] at xs. 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]. derivs is the user-supplied routine that calculates dydx. { void lubksb(float **a, int n, int *indx, float b[]); void ludcmp(float **a, int n, int *indx, float *d); int i,j,nn,*indx; float d,h,x,**a,*del,*ytemp; indx=ivector(1,n); a=matrix(1,n,1,n); del=vector(1,n); ytemp=vector(1,n); h=htot/nstep; Stepsize this trip. for (i=1;i<=n;i++) { Set up the matrix 1 − hf . for (j=1;j<=n;j++) a[i][j] = -h*dfdy[i][j];