744 Chapter 16. Integration of Ordinary Differential Equations +a[1][1]; ludcmp(a,n,indx,&d); LU decomposition of the matrix. for(i=1;1<=n;i+) Set up right-hand side for first step.Use yout yout [i]=h*(dydx[i]+h*dfdx[i]); for temporary storage. lubksb(a,n,indx,yout); for (i=1;i<=n;i++) First step. ytemp[i]=y [i]+(del[i]=yout [i]); x=xs+h; (*derivs)(x,ytemp,yout); Use yout for temporary storage of derivatives. for (nn=2;nn<=nstep;nn++) General step. for (i=1;i<=n;i++) Set up right-hand side for general step yout [i]=h*yout [i]-del[i]; lubksb(a,n,indx,yout); for(1=1;1<=n;i+) ytemp[i]+(del[i]+2.0*yout [i]); granted for X+=h: (*derivs)(x,ytemp,yout); 1600 (including this one) for(i=1;1<=n;i++) Set up right-hand side for last step. 872 yout [i]=h*yout [i]-del[i]; lubksb(a,n,indx,yout); to any Cambridge from NUMERICAL RECIPES IN 18881992 for(i=1:1<n;i++) Take last step. yout [i]+ytemp[i]; free_vector(ytemp,1,n); free_vector(del,1,n); free_matrix(a,1,n,1,n); server computer, (North America to make free_ivector(indx,1,n); one paper UnN电.t THE ART The routine simpr is intended to be used in a routine stifbs that is almost exactly the Programs same as bsstep.The only differences are: .The stepsize sequence is n=2,6,10,14,22,34,50, (16.6.35) to dir where each member differs from its predecessor by the smallest multiple of 4 that makes the ratio of successive terms be <The parameter KMAXX is taken to be 7. The work per unit step now includes the cost of Jacobian evaluations as well 1788-1982 OF SCIENTIFIC COMPUTING(ISBN as function evaluations.We count one Jacobian evaluation as equivalent to N function evaluations,where N is the number of equations. v@cam Once again the user-supplied routine derivs is a dummy argument and so can have any name.However,to maintain "plug-compatibility"with rkqs,bsstep and .Further reproduction, 10-621 stiff,the routine jacobn is not an argument and must have exactly this name.It is called once per step to return f (dfdy)and of/or (dfdx)as functions ofx and y. Numerical Recipes 43108 Here is the routine,with comments pointing out only the differences from bsstep: #include <math.h> (outside #include "nrutil.h" Software. #define KMAXX 7 #define IMAXX (KMAXX+1) #define SAFE1 0.25 #define SAFE2 0.7 #define REDMAX 1.0e-5 #define REDMIN 0.7 #define TINY 1.0e-30 #define SCALMX 0.1 float **d,*x; void stifbs(float y],float dydx[],int nv,float *xx,float htry,float eps, f1 oat yscal口，f1 oat *hdid,f1 oat *hnext, void(*derivs)(float,f1oat☐，f1oat))744 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 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). ++a[i][i]; } ludcmp(a,n,indx,&d); LU decomposition of the matrix. for (i=1;i<=n;i++) Set up right-hand side for first step. Use yout yout[i]=h*(dydx[i]+h*dfdx[i]); for temporary storage. lubksb(a,n,indx,yout); for (i=1;i<=n;i++) First step. ytemp[i]=y[i]+(del[i]=yout[i]); x=xs+h; (*derivs)(x,ytemp,yout); Use yout for temporary storage of derivatives. for (nn=2;nn<=nstep;nn++) { General step. for (i=1;i<=n;i++) Set up right-hand side for general step. yout[i]=h*yout[i]-del[i]; lubksb(a,n,indx,yout); for (i=1;i<=n;i++) ytemp[i] += (del[i] += 2.0*yout[i]); x += h; (*derivs)(x,ytemp,yout); } for (i=1;i<=n;i++) Set up right-hand side for last step. yout[i]=h*yout[i]-del[i]; lubksb(a,n,indx,yout); for (i=1;i<=n;i++) Take last step. yout[i] += ytemp[i]; free_vector(ytemp,1,n); free_vector(del,1,n); free_matrix(a,1,n,1,n); free_ivector(indx,1,n); } The routine simpr is intended to be used in a routine stifbs that 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 the ratio of successive terms be ≤ 5 7 . The parameter KMAXX is taken to be 7. • 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 N is the number of equations. • Once again the user-supplied routine derivs is a dummy argument and so can have any name. However, to maintain “plug-compatibility” with rkqs, bsstep and stiff, the routine jacobn is not an argument and must have exactly this name. It is called once per step to return f (dfdy) and ∂f/∂x (dfdx) as functions of x and y. Here is the routine, with comments pointing out only the differences from bsstep: #include <math.h> #include "nrutil.h" #define KMAXX 7 #define IMAXX (KMAXX+1) #define SAFE1 0.25 #define SAFE2 0.7 #define REDMAX 1.0e-5 #define REDMIN 0.7 #define TINY 1.0e-30 #define SCALMX 0.1 float **d,*x; void stifbs(float y[], float dydx[], int nv, float *xx, float htry, float eps, float yscal[], float *hdid, float *hnext, void (*derivs)(float, float [], float []))