16.4 Richardson Extrapolation and the Bulirsch-Stoer Method 729 to the value of hnext returned from the previous step,as is the case if the routine is called by odeint. void mmid(float y[],float dydx[],int nvar,float xs,float htot, int nstep,float yout[],void (*derivs)(float,float []float[])); vo1 d pzextr(int iest,f1 oat xest,f1 oat yest[☐，f1 oat yz☐，f1 loat dy[], int nv); int i,iq,k,kk,km; static int first=1,kmax,kopt; static float epsold =-1.0,xnew; float eps1,errmax,fact,h,red,scale,work,wrkmin,xest; float *err,*yerr,*ysav,*yseq; static float a[IMAXX+1]; static float alf [KMAXX+1][KMAXX+1]; stat1c1 nt nseg[IMAXX+1]=[0,2,4,6,8,10,12,14,16,18}: int reduct,exitflag=0; d=matrix(1,nv,1,KMAXX) 1600 (including this one) granted for 18881992 err=vector(1,KMAXX); x=vector(1,KMAXX); 872 yerr=vector(1,nv); Cambridge ysav-vector(1,nv); n NUMERICAL RECIPES IN yseg=vector(1,nv); if (eps !epsold){ A new tolerance,so reinitialize. to make *hnext xnew =-1.0e29; "Impossible"values. (Nor eps1=SAFE1*eps; a[1]=nseg[1]+1; Compute work coefficients A America server computer, one paper UnN电.t THE for (k=1;k<=KMAXX;k++)a[k+1]=a[k]+nseq[k+1]; ART for (iq=2;iq<=KMAXX;iq++){ Compute a(k,g) for (k=1;k<ig;k++) alf [k][iq]=pow(eps1,(a[k+1]-a[iq+1])/ Programs (a[iq+1]-a[1]+1.0)*(2*k+1)); epsold=eps: for (kopt=2;kopt<KMAXX;kopt++) Determine optimal row number for if (a[kopt+1]a[kopt]*alf [kopt-1][kopt])break; convergence. to dir kmax-kopt; h=htry; 19881982 OF SCIENTIFIC COMPUTING(ISBN for (i=1;i<=nv;i++)ysav [i]-y[i]; Save the starting values if (*xx !xnew Il h !(*hnext)){ A new stepsize or a new integration: re-establish the order window. kopt=kmax; 10-621 reduct=0; for (;;) -43108 for (k=1;k<=kmax;k++){ Evaluate the sequence of modified Numerical Recipes xnew=(*xx)+h; midpoint integrations. if (xnew ==(*xx))nrerror("step size underflow in bsstep"); (outside mmid(ysav,dydx,nv,*xx,h,nseq[k],yseq,derivs); xest=SQR(h/nseq[k]); Squared,since error series is even North Software. pzextr(k,xest,yseq,y,yerr,nv); Perform extrapolation. 1f(k1=1) Compute normalized error estimate Ame errmax=TINY; c() for (i=1;i<=nv;i++)errmax=FMAX(errmax,fabs(yerr[i]/yscal[i])); errmax /eps; Scale error relative to tolerance. km=k-1; err [km]=pow(errmax/SAFE1,1.0/(2*km+1)); if (k !=1&&(k >kopt-1 II first)) In order window. if (errmax 1.0){ Converged. exitflag=1; break;16.4 Richardson Extrapolation and the Bulirsch-Stoer Method 729 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). to the value of hnext returned from the previous step, as is the case if the routine is called by odeint. { void mmid(float y[], float dydx[], int nvar, 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; static float epsold = -1.0,xnew; float eps1,errmax,fact,h,red,scale,work,wrkmin,xest; float *err,*yerr,*ysav,*yseq; static float a[IMAXX+1]; static float alf[KMAXX+1][KMAXX+1]; static int nseq[IMAXX+1]={0,2,4,6,8,10,12,14,16,18}; int reduct,exitflag=0; d=matrix(1,nv,1,KMAXX); err=vector(1,KMAXX); x=vector(1,KMAXX); yerr=vector(1,nv); ysav=vector(1,nv); yseq=vector(1,nv); if (eps != epsold) { A new tolerance, so reinitialize. *hnext = xnew = -1.0e29; “Impossible” values. eps1=SAFE1*eps; a[1]=nseq[1]+1; Compute work coefficients Ak. for (k=1;k<=KMAXX;k++) a[k+1]=a[k]+nseq[k+1]; for (iq=2;iq<=KMAXX;iq++) { Compute α(k, q). 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; for (kopt=2;kopt<KMAXX;kopt++) Determine optimal row number for if (a[kopt+1] > a[kopt]*alf[kopt-1][kopt]) break; convergence. kmax=kopt; } h=htry; for (i=1;i<=nv;i++) ysav[i]=y[i]; Save the starting values. if (*xx != xnew || h != (*hnext)) { A new stepsize or a new integration: first=1; re-establish the order window. kopt=kmax; } reduct=0; for (;;) { for (k=1;k<=kmax;k++) { Evaluate the sequence of modified xnew=(*xx)+h; midpoint integrations. if (xnew == (*xx)) nrerror("step size underflow in bsstep"); mmid(ysav,dydx,nv,*xx,h,nseq[k],yseq,derivs); xest=SQR(h/nseq[k]); Squared, since error series is even. pzextr(k,xest,yseq,y,yerr,nv); Perform extrapolation. if (k != 1) { Compute normalized error estimate errmax=TINY; (k). for (i=1;i<=nv;i++) errmax=FMAX(errmax,fabs(yerr[i]/yscal[i])); errmax /= eps; Scale error relative to tolerance. km=k-1; err[km]=pow(errmax/SAFE1,1.0/(2*km+1)); } if (k != 1 && (k >= kopt-1 || first)) { In order window. if (errmax < 1.0) { Converged. exitflag=1; break; }