746 Chapter 16. Integration of Ordinary Differential Equations xnew=(*xx)+h; if (xnew =(*xx))nrerror("step size underflow in stifbs"); simpr(ysav,dydx,dfdx,dfdy,nv,*xx,h,nseq[k],yseq,derivs); Semi-implicit midpoint rule. xest=SQR(h/nseq [k]); The rest of the routine is identical to pzextr(k,xest,yseq,y,yerr,nv); bsstep. 1f(k1=1)[ errmax=TINY; for (i=1;i<=nv;i++)errmax=FMAX(errmax,fabs(yerr[i]/yscal[i])); errmax /eps; kmsk-1; err [km]=pow(errmax/SAFE1,1.0/(2*km+1)); http://www.nr. readable files Permission is opyright(C Sample page (k !1&(k >kopt-1 II first)){ if (errmax 1.0){ exitflag=1; break: 3 if (k =kmax Il k =kopt+1){ .com or call 1-800-872- (including this one) granted fori internet red=SAFE2/err [km]; break; else if (k =kopt&alf [kopt-1][kopt]err[km]){ red=1.0/err[km]; break; 2 else if (kopt =kmax&&alf [km][kmax-1]err[km]){ -7423 (North America red-alf [km][kmax-1]*SAFE2/err [km]; break; 2 t users to make one paper copy for their else if (alf [km][kopt]err[km]){ red=alf [km][kopt-1]/err[km]; to any server computer,is strictly prohibited. break: if (exitflag)break: red=FMIN(red,REDMIN): red=FMAX(red,REDMAX) h *=red; 1988-1992 by Cambridge University Press.Programs Copyright (C)1988-1992 by Numerical Recipes from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) reduct=1; *xx=xnew; *hdid=h; first=0; wrkmin=1.0e35; only),orsend email to directcustserv@cambridge.org(outside North America). personal use.Further reproduction,or for (kk=1;kk<=km;kk++){ fact=FMAX(err [kk],SCALMX); work=fact*a[kk+1]; if (work wrkmin){ scale=fact; Software. wrkmin=work; kopt=kk+1; a如 ying of machine *hnext=h/scale; if (kopt >=k&&kopt !kmax&&!reduct){ fact=FMAX(scale/alf [kopt-1][kopt]SCALMX); if (a[kopt+1]*fact <wrkmin){ *hnext=h/fact; kopt++; free_vector(yseq,1,nv);746 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). xnew=(*xx)+h; if (xnew == (*xx)) nrerror("step size underflow in stifbs"); simpr(ysav,dydx,dfdx,dfdy,nv,*xx,h,nseq[k],yseq,derivs); Semi-implicit midpoint rule. xest=SQR(h/nseq[k]); The rest of the routine is identical to pzextr(k,xest,yseq,y,yerr,nv); bsstep. if (k != 1) { errmax=TINY; for (i=1;i<=nv;i++) errmax=FMAX(errmax,fabs(yerr[i]/yscal[i])); errmax /= eps; km=k-1; err[km]=pow(errmax/SAFE1,1.0/(2*km+1)); } if (k != 1 && (k >= kopt-1 || first)) { if (errmax < 1.0) { exitflag=1; break; } if (k == kmax || k == kopt+1) { red=SAFE2/err[km]; break; } else if (k == kopt && alf[kopt-1][kopt] < err[km]) { red=1.0/err[km]; break; } else if (kopt == kmax && alf[km][kmax-1] < err[km]) { red=alf[km][kmax-1]*SAFE2/err[km]; break; } else if (alf[km][kopt] < err[km]) { red=alf[km][kopt-1]/err[km]; break; } } } if (exitflag) break; red=FMIN(red,REDMIN); red=FMAX(red,REDMAX); h *= red; reduct=1; } *xx=xnew; *hdid=h; first=0; wrkmin=1.0e35; for (kk=1;kk<=km;kk++) { fact=FMAX(err[kk],SCALMX); work=fact*a[kk+1]; if (work < wrkmin) { scale=fact; wrkmin=work; kopt=kk+1; } } *hnext=h/scale; if (kopt >= k && kopt != kmax && !reduct) { fact=FMAX(scale/alf[kopt-1][kopt],SCALMX); if (a[kopt+1]*fact <= wrkmin) { *hnext=h/fact; kopt++; } } free_vector(yseq,1,nv);