730 Chapter 16. Integration of Ordinary Differential Equations if (k =kmax Il k =kopt+1) Check for possible stepsize red-SAFE2/err [km]; reduction. 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; http://www.nr Permission is read able files else if (alf [km][kopt]err[km]){ red=alf [km][kopt-1]/err[km]; break; 839 granted for 19881992 11-600 if (exitflag)break: red=FMIN(red,REDMIN): Reduce stepsize by at least REDMIN red=FMAX(red,REDMAX) and at most REDMAX. /Cambridge h *red; from NUMERICAL RECIPES IN reduct=l: Try again. (Nor server computer, to make *xx=xnew; Successful step taken. *hdid=h; America one paper UnN电.t THE first=0; wrkmin=1.0e35; Compute optimal row for convergence ART for (kk=1;kk<=km;kk++){ and corresponding stepsize. fact=FMAX(err [kk],SCALMX); work=fact*a[kk+1]; for their Programs if (work wrkmin){ scale=fact; wrkmin-work; 1CIYP.ic kopt=kk+1; *hnext-h/scale; if (kopt >k&&kopt !kmax&&!reduct){ OF SCIENTIFIC COMPUTING(ISBN Check for possible order increase,but not if stepsize was just reduced. fact=FMAX(scale/alf [kopt-1][kopt],SCALMX); if (a[kopt+1]*fact <wrkmin){ *hnext=h/fact; kopt++; @cambridge.org 1988-1992 by Numerical Recipes 10-:6211 43106 free_vector(yseq,1,nv); free_vector(ysav,1,nv); free_vector(yerr,1,nv); free_vector(x,1,KMAXX) free_vector(err,1,KMAXX); (outside North Amer Software. free_matrix(d,1,nv,1,KMAXX); visit website The polynomial extrapolation routine is based on the same algorithm as polint 83.1.It is simpler in that it is always extrapolating to zero,rather than to an arbitrary value.However,it is more complicated in that it must individually extrapolate each component of a vector of quantities.730 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). if (k == kmax || k == kopt+1) { Check for possible stepsize red=SAFE2/err[km]; reduction. 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); Reduce stepsize by at least REDMIN red=FMAX(red,REDMAX); and at most REDMAX. h *= red; reduct=1; } Try again. *xx=xnew; Successful step taken. *hdid=h; first=0; wrkmin=1.0e35; Compute optimal row for convergence for (kk=1;kk<=km;kk++) { and corresponding stepsize. 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) { Check for possible order increase, but not if stepsize was just reduced. fact=FMAX(scale/alf[kopt-1][kopt],SCALMX); if (a[kopt+1]*fact <= wrkmin) { *hnext=h/fact; kopt++; } } free_vector(yseq,1,nv); free_vector(ysav,1,nv); free_vector(yerr,1,nv); free_vector(x,1,KMAXX); free_vector(err,1,KMAXX); free_matrix(d,1,nv,1,KMAXX); } The polynomial extrapolation routine is based on the same algorithm as polint §3.1. It is simpler in that it is always extrapolating to zero, rather than to an arbitrary value. However, it is more complicated in that it must individually extrapolate each component of a vector of quantities