c16-2 (779594), страница 3
Текст из файла (страница 3)
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).static float a2=0.2,a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2,b31=3.0/40.0,b32=9.0/40.0,b41=0.3,b42 = -0.9,b43=1.2,b51 = -11.0/54.0, b52=2.5,b53 = -70.0/27.0,b54=35.0/27.0,b61=1631.0/55296.0,b62=175.0/512.0,b63=575.0/13824.0,b64=44275.0/110592.0,b65=253.0/4096.0,c1=37.0/378.0,c3=250.0/621.0,c4=125.0/594.0,c6=512.0/1771.0,dc5 = -277.00/14336.0;float dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0,dc4=c4-13525.0/55296.0,dc6=c6-0.25;float *ak2,*ak3,*ak4,*ak5,*ak6,*ytemp;16.2 Adaptive Stepsize Control for Runge-Kutta721#include <math.h>#include "nrutil.h"#define MAXSTP 10000#define TINY 1.0e-30extern int kmax,kount;extern float *xp,**yp,dxsav;User storage for intermediate results.
Preset kmax and dxsav in the calling program. If kmax 6=0 results are stored at approximate intervals dxsav in the arrays xp[1..kount], yp[1..nvar][1..kount], where kount is output by odeint. Defining declarations for these variables, withmemory allocations xp[1..kmax] and yp[1..nvar][1..kmax] for the arrays, should be inthe calling program.void odeint(float ystart[], int nvar, float x1, float x2, float eps, float h1,float hmin, int *nok, int *nbad,void (*derivs)(float, float [], float []),void (*rkqs)(float [], float [], int, float *, float, float, float [],float *, float *, void (*)(float, float [], float [])))Runge-Kutta driver with adaptive stepsize control.
Integrate starting values ystart[1..nvar]from x1 to x2 with accuracy eps, storing intermediate results in global variables. h1 shouldbe set as a guessed first stepsize, hmin as the minimum allowed stepsize (can be zero). Onoutput nok and nbad are the number of good and bad (but retried and fixed) steps taken, andystart is replaced by values at the end of the integration interval.
derivs is the user-suppliedroutine for calculating the right-hand side derivative, while rkqs is the name of the stepperroutine to be used.{int nstp,i;float xsav,x,hnext,hdid,h;float *yscal,*y,*dydx;yscal=vector(1,nvar);y=vector(1,nvar);dydx=vector(1,nvar);x=x1;h=SIGN(h1,x2-x1);*nok = (*nbad) = kount = 0;for (i=1;i<=nvar;i++) y[i]=ystart[i];if (kmax > 0) xsav=x-dxsav*2.0;Assures storage of first step.for (nstp=1;nstp<=MAXSTP;nstp++) {Take at most MAXSTP steps.(*derivs)(x,y,dydx);for (i=1;i<=nvar;i++)Scaling used to monitor accuracy. This general-purpose choice can be modifiedif need be.yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;if (kmax > 0 && kount < kmax-1 && fabs(x-xsav) > fabs(dxsav)) {xp[++kount]=x;Store intermediate results.for (i=1;i<=nvar;i++) yp[i][kount]=y[i];xsav=x;}if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;If stepsize can overshoot, decrease.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).including garden-variety ODEs or sets of ODEs, and definite integrals (augmentingthe methods of Chapter 4).
For storage of intermediate results (if you desire toinspect them) we assume that the top-level pointer references *xp and **yp have beenvalidly initialized (e.g., by the utilities vector() and matrix()). Because stepsoccur at unequal intervals results are only stored at intervals greater than dxsav. Thetop-level variable kmax indicates the maximum number of steps that can be stored.If kmax=0 there is no intermediate storage, and the pointers *xp and **yp need notpoint to valid memory.
Storage of steps stops if kmax is exceeded, except that theending values are always stored. Again, these controls are merely indicative of whatyou might need. The routine odeint should be customized to the problem at hand.722Chapter 16.Integration of Ordinary Differential Equations}nrerror("Too many steps in routine odeint");}CITED REFERENCES AND FURTHER READING:Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (EnglewoodCliffs, NJ: Prentice-Hall). [1]Cash, J.R., and Karp, A.H.
1990, ACM Transactions on Mathematical Software, vol. 16, pp. 201–222. [2]Shampine, L.F., and Watts, H.A. 1977, in Mathematical Software III, J.R. Rice, ed. (New York:Academic Press), pp. 257–275; 1979, Applied Mathematics and Computation, vol. 5,pp. 93–121.Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for MathematicalComputations (Englewood Cliffs, NJ: Prentice-Hall).16.3 Modified Midpoint MethodThis section discusses the modified midpoint method, which advances a vectorof dependent variables y(x) from a point x to a point x + H by a sequence of nsubsteps each of size h,h = H/n(16.3.1)In principle, one could use the modified midpoint method in its own right as an ODEintegrator. In practice, the method finds its most important application as a part ofthe more powerful Bulirsch-Stoer technique, treated in §16.4.
You can thereforeconsider this section as a preamble to §16.4.The number of right-hand side evaluations required by the modified midpointmethod is n + 1. The formulas for the method arez0 ≡ y(x)z1 = z0 + hf(x, z0 )zm+1 = zm−1 + 2hf(x + mh, zm )y(x + H) ≈ yn ≡1[zn + zn−1 + hf(x + H, zn )]2for m = 1, 2, . . . , n − 1(16.3.2)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).(*rkqs)(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,derivs);if (hdid == h) ++(*nok); else ++(*nbad);if ((x-x2)*(x2-x1) >= 0.0) {Are we done?for (i=1;i<=nvar;i++) ystart[i]=y[i];if (kmax) {xp[++kount]=x;Save final step.for (i=1;i<=nvar;i++) yp[i][kount]=y[i];}free_vector(dydx,1,nvar);free_vector(y,1,nvar);free_vector(yscal,1,nvar);return;Normal exit.}if (fabs(hnext) <= hmin) nrerror("Step size too small in odeint");h=hnext;.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
