c16-6 (779598), страница 4
Текст из файла (страница 4)
. ,(16.6.35)where each member differs from its predecessor by the smallest multiple of 4 thatmakes the ratio of successive terms be ≤ 57 . The parameter KMAXX is taken to be 7.• The work per unit step now includes the cost of Jacobian evaluations as wellas function evaluations. We count one Jacobian evaluation as equivalent to Nfunction evaluations, where N is the number of equations.• Once again the user-supplied routine derivs is a dummy argument and so can haveany name. However, to maintain “plug-compatibility” with rkqs, bsstep andstiff, the routine jacobn is not an argument and must have exactly this name. Itis called once per step to return f 0 (dfdy) and ∂f/∂x (dfdx) as functions of x and y.Here is the routine, with comments pointing out only the differences from bsstep:#include <math.h>#include "nrutil.h"#define KMAXX 7#define IMAXX (KMAXX+1)#define SAFE1 0.25#define SAFE2 0.7#define REDMAX 1.0e-5#define REDMIN 0.7#define TINY 1.0e-30#define SCALMX 0.1float **d,*x;void stifbs(float y[], float dydx[], int nv, float *xx, float htry, float eps,float yscal[], float *hdid, float *hnext,void (*derivs)(float, float [], float []))Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).++a[i][i];}ludcmp(a,n,indx,&d);LU decomposition of the matrix.for (i=1;i<=n;i++)Set up right-hand side for first step.
Use youtyout[i]=h*(dydx[i]+h*dfdx[i]);for temporary storage.lubksb(a,n,indx,yout);for (i=1;i<=n;i++)First step.ytemp[i]=y[i]+(del[i]=yout[i]);x=xs+h;(*derivs)(x,ytemp,yout);Use yout for temporary storage of derivatives.for (nn=2;nn<=nstep;nn++) {General step.for (i=1;i<=n;i++)Set up right-hand side for general step.yout[i]=h*yout[i]-del[i];lubksb(a,n,indx,yout);for (i=1;i<=n;i++)ytemp[i] += (del[i] += 2.0*yout[i]);x += h;(*derivs)(x,ytemp,yout);}for (i=1;i<=n;i++)Set up right-hand side for last step.yout[i]=h*yout[i]-del[i];lubksb(a,n,indx,yout);for (i=1;i<=n;i++)Take last step.yout[i] += ytemp[i];free_vector(ytemp,1,n);free_vector(del,1,n);free_matrix(a,1,n,1,n);free_ivector(indx,1,n);16.6 Stiff Sets of Equations745d=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 workfor (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++) {Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Semi-implicit extrapolation step for integrating stiff o.d.e.’s, with monitoring of local truncationerror to adjust stepsize.
Input are the dependent variable vector y[1..n] and its derivativedydx[1..n] at the starting value of the independent variable x. Also input are the stepsize tobe attempted htry, the required accuracy eps, and the vector yscal[1..n] against whichthe error is scaled. On output, y and x are replaced by their new values, hdid is the stepsizethat was actually accomplished, and hnext is the estimated next stepsize. derivs is a usersupplied routine that computes the derivatives of the right-hand side with respect to x, whilejacobn (a fixed name) is a user-supplied routine that computes the Jacobi matrix of derivativesof the right-hand side with respect to the components of y. Be sure to set htry on successivesteps to the value of hnext returned from the previous step, as is the case if the routine iscalled 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 fromint reduct,exitflag=0;bsstep.746Chapter 16.Integration of Ordinary Differential Equations}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);Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (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 topzextr(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;}}16.7 Multistep, Multivalue, and Predictor-Corrector Methods747free_vector(ysav,1,nv);free_vector(yerr,1,nv);free_vector(x,1,KMAXX);free_vector(err,1,KMAXX);free_matrix(dfdy,1,nv,1,nv);free_vector(dfdx,1,nv);free_matrix(d,1,nv,1,KMAXX);The routine stifbs is an excellent routine for all stiff problems, competitive withthe best Gear-type routines.
stiff is comparable in execution time for moderate N and−4−8<∼ 10 . By the time ∼ 10 , stifbs is roughly an order of magnitude faster. Thereare further improvements that could be applied to stifbs to make it even more robust. Forexample, very occasionally ludcmp in simpr will encounter a singular matrix. You couldarrange for the stepsize to be reduced, say by a factor of the current nseq[k]. There arealso certain stability restrictions on the stepsize that come into play on some problems. Fora discussion of how to implement these automatically, see [6].CITED REFERENCES AND FURTHER READING:Gear, C.W.
1971, Numerical Initial Value Problems in Ordinary Differential Equations (EnglewoodCliffs, NJ: Prentice-Hall). [1]Kaps, P., and Rentrop, P. 1979, Numerische Mathematik, vol. 33, pp. 55–68. [2]Shampine, L.F. 1982, ACM Transactions on Mathematical Software, vol. 8, pp. 93–113. [3]Enright, W.H., and Pryce, J.D. 1987, ACM Transactions on Mathematical Software, vol. 13,pp. 1–27. [4]Bader, G., and Deuflhard, P. 1983, Numerische Mathematik, vol.
41, pp. 373–398. [5]Deuflhard, P. 1983, Numerische Mathematik, vol. 41, pp. 399–422.Deuflhard, P. 1985, SIAM Review, vol. 27, pp. 505–535.Deuflhard, P. 1987, “Uniqueness Theorems for Stiff ODE Initial Value Problems,” Preprint SC87-3 (Berlin: Konrad Zuse Zentrum für Informationstechnik). [6]Enright, W.H., Hull, T.E., and Lindberg, B. 1975, BIT, vol. 15, pp.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
