740 Chapter 16. Integration of Ordinary Differential Equations #define C1X (1.0/2.0) #def1neC2x(-3.0/2.0) #define C3X (121.0/50.0) #define C4X (29.0/250.0) #define A2X 1.0 #define A3X (3.0/5.0) void stiff(f1oaty▣，f1 oat dydx，intn,f1oat*x,float htry,f1 oat eps, f1 oat yscal☐，f1oat*hdid,f1oat*hnext, void(*derivs)(f1oat,f1oat[，f1oat[])) Fourth-order Rosenbrock step for integrating stiff o.d.e.'s,with monitoring of local truncation error to adjust stepsize.Input are the dependent variable vector y[1..n]and its derivative dydx [1..n]at the starting value of the independent variable x.Also input are the stepsize to be attempted htry.the required accuracy eps,and the vector yscal[1..n]against which the error is scaled.On output,y and x are replaced by their new values,hdid is the stepsize 83g that was actually accomplished,and hnext is the estimated next stepsize.derivs is a user- granted for 18881992 supplied routine that computes the derivatives of the right-hand side with respect to x,while jacobn (a fixed name)is a user-supplied routine that computes the Jacobi matrix of derivatives 1-600 of the right-hand side with respect to the components of y void jacobn(float x,float y[],float dfdx[],float **dfdy,int n); void lubksb(float **a,int n,int *indx,float b[]); from NUMERICAL RECIPES IN void ludcmp(float **a,int n,int *indx,float *d); int i,j,jtry,*indx; float d,errmax,h,xsav,**a,*dfdx,**dfdy,*dysav,*err; float *g1,*g2,*g3,*g4,*ysav; (North America to any server computer,is sto make one paper THE indx=ivector(1,n); /Cambridge University Press.Programs ART a=matrix(1,n,1,n); dfdx-vector(1,n); dfdy=matrix(1,n,1,n); dysav=vector(1,n); copy for their err=vector(1,n); gi=vector(1,n); g2-vector(1,n); g3=vector(1,n); g4=vector(1,n); ysav=vector(1,n); xsav=(*x); Save initial values to directcustser for(1=1;1<=n;i++){ ysav[i]=y[i]; dysav[i]=dydx[i]; v@cambr jacobn(xsav,ysav,dfdx,dfdy,n): The user must supply this routine to return the n-by-n matrix dfdy and the vector dfdx h=htry; Set stepsize to the initial trial value. G7/Repe气maR 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING (ISBN 0-521 -43108 for (jtry=1;jtry<=MAXTRY;jtry++) for(i=1;i<=n;1++){ Set up the matrix 1-hf. for (j=1;j<=n;j++)a[i][j]=-dfdy[i][j]; (outside a[1][i]+=1.0/(GAM*h); North Software. ludcmp(a,n,indx,&d); LU decomposition of the matrix. for (i=1;i<=n;i++) Set up right-hand side for g1. Amer g1[i]=dysav[i]+h*C1X*dfdx[i]; lubksb(a,n,indx,g1); Solve for g visit website for (i=1;i<=n;i++) Compute intermediate values of y and x. 3 y[i]=ysav[i]+A21*g1[1]; *x=xsav+A2X*h: (*derivs)(*x,y,dydx)i Compute dydx at the intermediate values. for (i=1;i<=n;i++) Set up right-hand side for g2. g2[i]=dydx [i]+h*C2X*dfdx [i]+C21*g1[i]/h; lubksb(a,n,indx,g2); Solve for g2- for (i=1;i<=n;i++) Compute intermediate values of y and x. y[1]=ysav[1]+A31*g1[1]+A32*g2[1]:740 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). #define C1X (1.0/2.0) #define C2X (-3.0/2.0) #define C3X (121.0/50.0) #define C4X (29.0/250.0) #define A2X 1.0 #define A3X (3.0/5.0) void stiff(float y[], float dydx[], int n, float *x, float htry, float eps, float yscal[], float *hdid, float *hnext, void (*derivs)(float, float [], float [])) Fourth-order Rosenbrock step for integrating stiffo.d.e.’s, with monitoring of local truncation error to adjust stepsize. Input are the dependent variable vector y[1..n] and its derivative dydx[1..n] at the starting value of the independent variable x. Also input are the stepsize to be attempted htry, the required accuracy eps, and the vector yscal[1..n] against which the error is scaled. On output, y and x are replaced by their new values, hdid is the stepsize that was actually accomplished, and hnext is the estimated next stepsize. derivs is a user￾supplied routine that computes the derivatives of the right-hand side with respect to x, while jacobn (a fixed name) is a user-supplied routine that computes the Jacobi matrix of derivatives of the right-hand side with respect to the components of y. { void jacobn(float x, float y[], float dfdx[], float **dfdy, int n); void lubksb(float **a, int n, int *indx, float b[]); void ludcmp(float **a, int n, int *indx, float *d); int i,j,jtry,*indx; float d,errmax,h,xsav,**a,*dfdx,**dfdy,*dysav,*err; float *g1,*g2,*g3,*g4,*ysav; indx=ivector(1,n); a=matrix(1,n,1,n); dfdx=vector(1,n); dfdy=matrix(1,n,1,n); dysav=vector(1,n); err=vector(1,n); g1=vector(1,n); g2=vector(1,n); g3=vector(1,n); g4=vector(1,n); ysav=vector(1,n); xsav=(*x); Save initial values. for (i=1;i<=n;i++) { ysav[i]=y[i]; dysav[i]=dydx[i]; } jacobn(xsav,ysav,dfdx,dfdy,n); The user must supply this routine to return the n-by-n matrix dfdy and the vector dfdx. h=htry; Set stepsize to the initial trial value. for (jtry=1;jtry<=MAXTRY;jtry++) { for (i=1;i<=n;i++) { Set up the matrix 1 − γhf . for (j=1;j<=n;j++) a[i][j] = -dfdy[i][j]; a[i][i] += 1.0/(GAM*h); } ludcmp(a,n,indx,&d); LU decomposition of the matrix. for (i=1;i<=n;i++) Set up right-hand side for g1. g1[i]=dysav[i]+h*C1X*dfdx[i]; lubksb(a,n,indx,g1); Solve for g1. for (i=1;i<=n;i++) Compute intermediate values of y and x. y[i]=ysav[i]+A21*g1[i]; *x=xsav+A2X*h; (*derivs)(*x,y,dydx); Compute dydx at the intermediate values. for (i=1;i<=n;i++) Set up right-hand side for g2. g2[i]=dydx[i]+h*C2X*dfdx[i]+C21*g1[i]/h; lubksb(a,n,indx,g2); Solve for g2. for (i=1;i<=n;i++) Compute intermediate values of y and x. y[i]=ysav[i]+A31*g1[i]+A32*g2[i];