c16-1 (779593), страница 2
Текст из файла (страница 2)
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).316.1 Runge-Kutta Method713}The Runge-Kutta method treats every step in a sequence of steps in identicalmanner. Prior behavior of a solution is not used in its propagation.
This ismathematically proper, since any point along the trajectory of an ordinary differentialequation can serve as an initial point. The fact that all steps are treated identically alsomakes it easy to incorporate Runge-Kutta into relatively simple “driver” schemes.We consider adaptive stepsize control, discussed in the next section, an essentialfor serious computing. Occasionally, however, you just want to tabulate a function atequally spaced intervals, and without particularly high accuracy. In the most commoncase, you want to produce a graph of the function. Then all you need may be asimple driver program that goes from an initial xs to a final xf in a specified numberof steps. To check accuracy, double the number of steps, repeat the integration, andcompare results. This approach surely does not minimize computer time, and it canfail for problems whose nature requires a variable stepsize, but it may well minimizeuser effort.
On small problems, this may be the paramount consideration.Here is such a driver, self-explanatory, which tabulates the integrated functionsin the global arrays *x and **y; be sure to allocate memory for them with theroutines vector() and matrix(), respectively.#include "nrutil.h"float **y,*xx;For communication back to main.void rkdumb(float vstart[], int nvar, float x1, float x2, int nstep,void (*derivs)(float, float [], float []))Starting from initial values vstart[1..nvar] known at x1 use fourth-order Runge-Kuttato advance nstep equal increments to x2.
The user-supplied routine derivs(x,v,dvdx)evaluates derivatives. Results are stored in the global variables y[1..nvar][1..nstep+1]and xx[1..nstep+1].{void rk4(float y[], float dydx[], int n, float x, float h, float yout[],void (*derivs)(float, float [], float []));int i,k;float x,h;float *v,*vout,*dv;v=vector(1,nvar);vout=vector(1,nvar);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).yt=vector(1,n);hh=h*0.5;h6=h/6.0;xh=x+hh;for (i=1;i<=n;i++) yt[i]=y[i]+hh*dydx[i];First step.(*derivs)(xh,yt,dyt);Second step.for (i=1;i<=n;i++) yt[i]=y[i]+hh*dyt[i];(*derivs)(xh,yt,dym);Third step.for (i=1;i<=n;i++) {yt[i]=y[i]+h*dym[i];dym[i] += dyt[i];}(*derivs)(x+h,yt,dyt);Fourth step.for (i=1;i<=n;i++)Accumulate increments with properyout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);weights.free_vector(yt,1,n);free_vector(dyt,1,n);free_vector(dym,1,n);714Chapter 16.Integration of Ordinary Differential Equations}CITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A.
1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.5. [1]Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (EnglewoodCliffs, NJ: Prentice-Hall), Chapter 2. [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. [3]Rice, J.R. 1983, Numerical Methods, Software, and Analysis (New York: McGraw-Hill), §9.2.16.2 Adaptive Stepsize Control for Runge-KuttaA good ODE integrator should exert some adaptive control over its own progress,making frequent changes in its stepsize.
Usually the purpose of this adaptive stepsizecontrol is to achieve some predetermined accuracy in the solution with minimumcomputational effort. Many small steps should tiptoe through treacherous terrain,while a few great strides should speed through smooth uninteresting countryside.The resulting gains in efficiency are not mere tens of percents or factors of two;they can sometimes be factors of ten, a hundred, or more. Sometimes accuracymay be demanded not directly in the solution itself, but in some related conservedquantity that can be monitored.Implementation of adaptive stepsize control requires that the stepping algorithmsignal information about its performance, most important, an estimate of its truncationerror. In this section we will learn how such information can be obtained.
Obviously,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).dv=vector(1,nvar);for (i=1;i<=nvar;i++) {Load starting values.v[i]=vstart[i];y[i][1]=v[i];}xx[1]=x1;x=x1;h=(x2-x1)/nstep;for (k=1;k<=nstep;k++) {Take nstep steps.(*derivs)(x,v,dv);rk4(v,dv,nvar,x,h,vout,derivs);if ((float)(x+h) == x) nrerror("Step size too small in routine rkdumb");x += h;xx[k+1]=x;Store intermediate steps.for (i=1;i<=nvar;i++) {v[i]=vout[i];y[i][k+1]=v[i];}}free_vector(dv,1,nvar);free_vector(vout,1,nvar);free_vector(v,1,nvar);.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
