16.6 Stiff Sets of Equations 745 Semi-implicit extrapolation 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..nv]and its derivative dydx [1..nv]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..nv]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.Be sure to set htry on successive steps to the value of hnext returned from the previous step,as is the case if the routine is called by odeint void jacobn(float x,float y],float dfdx[],float **dfdy,int n); void simpr(float y[],float dydx[],float dfdx[],float **dfdy, int n,float xs,float htot,int nstep,float yout void(*derivs)(f1oat,f1oat☐，float[☐))； granted for 18881892 void pzextr(int1est,f1 oat xest,f1 oat yest[☐，f1 oat yz☐，f1 loat dy[], int nv); 1.800 int i,iq,k,kk,km; static int first=1,kmax,kopt,nvold =-1; static float epsold =-1.0,xnew; float eps1,errmax,fact,h,red,scale,work,wrkmin,xest; from NUMERICAL RECIPES IN float *dfdx,**dfdy,*err,*yerr,*ysav,*yseqi static float a[IMAXX+1]; (North static float alf [KMAXX+1][KMAXX+1]; stat1c1 nt nseq[IMAXX+1]=[0,2,6,10,14,22,34,50,70]; Sequence is different from America to any server computer. one paper Cambridge University Press. THE int reduct,exitflag=0; bsstep. ART d=matrix(1,nv,1,KMAXX) dfdx=vector(1,nv); only), dfdy=matrix(1,nv,1,nv); err=vector(1,KMAXX); x=vector(1,KMAXX) yerr=vector(1,nv); st st Programs ysav-vector(1,nv); to dir Copyright (C) yseq=vector(1,nv); if(eps !epsold II nv !nvold){ Reinitialize also if nv has changed. 米nnext=Xnew= -1.0e29; eps1=SAFE1*eps; a[1]=nseq[1]+1: ectcustser To order OF SCIENTIFIC COMPUTING(ISBN for (k=1;k<=KMAXX;k++)a[k+1]=a[k]+nseq[k+1]; for (iq=2;iq<=KMAXX;iq++) for (k=1;k<iq;k++) 1988-1992 by Numerical Recipes 10-621 alf [k][ig]=pow(eps1,((a[k+1]-a[iq+1])/ ((a[ig+1]-a[1]+1.0)*(2*k+1)))); -43108 epsold=eps; nvold=nv; Save nv. a[1]+=nv; Add cost of Jacobian evaluations to work (outside for (k=1;k<=KMAXX;k++)a[k+1]=a[k]+nseq[k+1]; coefficients. for (kopt=2;kopt<KMAXX;kopt++) North Software. if (a[kopt+1]a[kopt]*alf [kopt-1][kopt])break; kmax=kopt; h=htry; for (i=1;i<-nv;i++)ysav[i]sy[i]; jacobn(*xx,y,dfdx,dfdy,nv); Evaluate Jacobian if (*xx !xnew II h !(*hnext)){ first=1; kopt=kmax; reduct=0; for (;;) for (k=1;k<=kmax;k++){16.6 Stiff Sets of Equations 745 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). Semi-implicit extrapolation step for integrating stiffo.d.e.’s, with monitoring of local truncation error to adjust stepsize. Input are the dependent variable vector y[1..nv] and its derivative dydx[1..nv] 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..nv] 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. Be sure to set htry on successive steps to the value of hnext returned from the previous step, as is the case if the routine is called by odeint. { void jacobn(float x, float y[], float dfdx[], float **dfdy, int n); 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 [])); void pzextr(int iest, float xest, float yest[], float yz[], float dy[], int nv); int i,iq,k,kk,km; static int first=1,kmax,kopt,nvold = -1; static float epsold = -1.0,xnew; float eps1,errmax,fact,h,red,scale,work,wrkmin,xest; float *dfdx,**dfdy,*err,*yerr,*ysav,*yseq; static float a[IMAXX+1]; static float alf[KMAXX+1][KMAXX+1]; static int nseq[IMAXX+1]={0,2,6,10,14,22,34,50,70}; Sequence is different from int reduct,exitflag=0; bsstep. d=matrix(1,nv,1,KMAXX); dfdx=vector(1,nv); dfdy=matrix(1,nv,1,nv); err=vector(1,KMAXX); x=vector(1,KMAXX); yerr=vector(1,nv); ysav=vector(1,nv); yseq=vector(1,nv); if(eps != epsold || nv != nvold) { Reinitialize also if nv has changed. *hnext = xnew = -1.0e29; eps1=SAFE1*eps; a[1]=nseq[1]+1; for (k=1;k<=KMAXX;k++) a[k+1]=a[k]+nseq[k+1]; for (iq=2;iq<=KMAXX;iq++) { for (k=1;k<iq;k++) alf[k][iq]=pow(eps1,((a[k+1]-a[iq+1])/ ((a[iq+1]-a[1]+1.0)*(2*k+1)))); } epsold=eps; nvold=nv; Save nv. a[1] += nv; Add cost of Jacobian evaluations to work for (k=1;k<=KMAXX;k++) a[k+1]=a[k]+nseq[k+1]; coefficients. for (kopt=2;kopt<KMAXX;kopt++) if (a[kopt+1] > a[kopt]*alf[kopt-1][kopt]) break; kmax=kopt; } h=htry; for (i=1;i<=nv;i++) ysav[i]=y[i]; jacobn(*xx,y,dfdx,dfdy,nv); Evaluate Jacobian. if (*xx != xnew || h != (*hnext)) { first=1; kopt=kmax; } reduct=0; for (;;) { for (k=1;k<=kmax;k++) {