c16-2 (779594)
Текст из файла
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);16.2 Adaptive Stepsize Control for Runge-Kutta715y(x + 2h) = y1 + (2h)5 φ + O(h6 ) + .
. .y(x + 2h) = y2 + 2(h5 )φ + O(h6 ) + . . .(16.2.1)where, to order h5 , the value φ remains constant over the step. [Taylor seriesexpansion tells us the φ is a number whose order of magnitude is y(5) (x)/5!.] Thefirst expression in (16.2.1) involves (2h)5 since the stepsize is 2h, while the secondexpression involves 2(h5 ) since the error on each step is h5 φ. The differencebetween the two numerical estimates is a convenient indicator of truncation error∆ ≡ y2 − y1(16.2.2)It is this difference that we shall endeavor to keep to a desired degree of accuracy,neither too large nor too small.
We do this by adjusting h.It might also occur to you that, ignoring terms of order h6 and higher, we cansolve the two equations in (16.2.1) to improve our numerical estimate of the truesolution y(x + 2h), namely,y(x + 2h) = y2 +∆+ O(h6 )15(16.2.3)This estimate is accurate to fifth order, one order higher than the original Runge-Kuttasteps. However, we can’t have our cake and eat it: (16.2.3) may be fifth-orderaccurate, but we have no way of monitoring its truncation error.
Higher order isnot always higher accuracy! Use of (16.2.3) rarely does harm, but we have noway of directly knowing whether it is doing any good. Therefore we should use∆ as the error estimate and take as “gravy” any additional accuracy gain derivedfrom (16.2.3). In the technical literature, use of a procedure like (16.2.3) is called“local extrapolation.”An alternative stepsize adjustment algorithm is based on the embedded RungeKutta formulas, originally invented by Fehlberg.
An interesting fact about RungeKutta formulas is that for orders M higher than four, more than M functionevaluations (though never more than M + 2) are required. This accounts for theSample 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).the calculation of this information will add to the computational overhead, but theinvestment will generally be repaid handsomely.With fourth-order Runge-Kutta, the most straightforward technique by far isstep doubling (see, e.g., [1]).
We take each step twice, once as a full step, then,independently, as two half steps (see Figure 16.2.1). How much overhead is this,say in terms of the number of evaluations of the right-hand sides? Each of thethree separate Runge-Kutta steps in the procedure requires 4 evaluations, but thesingle and double sequences share a starting point, so the total is 11.
This is to becompared not to 4, but to 8 (the two half-steps), since — stepsize control aside —we are achieving the accuracy of the smaller (half) stepsize. The overhead cost istherefore a factor 1.375. What does it buy us?Let us denote the exact solution for an advance from x to x + 2h by y(x + 2h)and the two approximate solutions by y1 (one step 2h) and y2 (2 steps each of sizeh). Since the basic method is fourth order, the true solution and the two numericalapproximations are related by716Chapter 16.Integration of Ordinary Differential Equationsbig steptwo small stepsFigure 16.2.1. Step-doubling as a means for adaptive stepsize control in fourth-order Runge-Kutta.Points where the derivative is evaluated are shown as filled circles. The open circle represents the samederivatives as the filled circle immediately above it, so the total number of evaluations is 11 per twosteps.
Comparing the accuracy of the big step with the two small steps gives a criterion for adjusting thestepsize on the next step, or for rejecting the current step as inaccurate.popularity of the classical fourth-order method: It seems to give the most bangfor the buck. However, Fehlberg discovered a fifth-order method with six functionevaluations where another combination of the six functions gives a fourth-ordermethod. The difference between the two estimates of y(x + h) can then be used asan estimate of the truncation error to adjust the stepsize. Since Fehlberg’s originalformula, several other embedded Runge-Kutta formulas have been found.Many practitioners were at one time wary of the robustness of Runge-KuttaFehlberg methods.
The feeling was that using the same evaluation points to advancethe function and to estimate the error was riskier than step-doubling, where theerror estimate is based on independent function evaluations. However, experiencehas shown that this concern is not a problem in practice. Accordingly, embeddedRunge-Kutta formulas, which are roughly a factor of two more efficient, havesuperseded algorithms based on step-doubling.The general form of a fifth-order Runge-Kutta formula isk1 = hf(xn , yn )k2 = hf(xn + a2 h, yn + b21 k1 )···(16.2.4)k6 = hf(xn + a6 h, yn + b61 k1 + · · · + b65 k5 )yn+1 = yn + c1 k1 + c2 k2 + c3 k3 + c4 k4 + c5 k5 + c6 k6 + O(h6 )The embedded fourth-order formula is∗yn+1= yn + c∗1 k1 + c∗2 k2 + c∗3 k3 + c∗4 k4 + c∗5 k5 + c∗6 k6 + O(h5 )(16.2.5)and so the error estimate is∗=∆ ≡ yn+1 − yn+16X(ci − c∗i )ki(16.2.6)i=1The particular values of the various constants that we favor are those found by Cashand Karp [2], and given in the accompanying table.
These give a more efficientmethod than Fehlberg’s original values, with somewhat better error properties.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).x71716.2 Adaptive Stepsize Control for Runge-KuttaCash-Karp Parameters for Embedded Runga-Kutta Methodiaibij1151533103409404353109− 106551− 115452− 702735276781631552961755125751382444275110592253409612345j=c∗i373782825276480025062118575483841255941352555296027714336512177114Now that we know, at least approximately, what our error is, we need toconsider how to keep it within desired bounds.
What is the relation between ∆and h? According to (16.2.4) – (16.2.5), ∆ scales as h5 . If we take a step h1and produce an error ∆1 , therefore, the step h0 that would have given some othervalue ∆0 is readily estimated as ∆0 0.2h0 = h1 ∆1 (16.2.7)Henceforth we will let ∆0 denote the desired accuracy. Then equation (16.2.7) isused in two ways: If ∆1 is larger than ∆0 in magnitude, the equation tells howmuch to decrease the stepsize when we retry the present (failed) step. If ∆1 issmaller than ∆0 , on the other hand, then the equation tells how much we can safelyincrease the stepsize for the next step.
Local extrapolation consists in acceptingthe fifth order value yn+1 , even though the error estimate actually applies to the∗fourth order value yn+1.Our notation hides the fact that ∆0 is actually a vector of desired accuracies,one for each equation in the set of ODEs. In general, our accuracy requirement willbe that all equations are within their respective allowed errors. In other words, wewill rescale the stepsize according to the needs of the “worst-offender” equation.How is ∆0 , the desired accuracy, related to some looser prescription like “get asolution good to one part in 106 ”? That can be a subtle question, and it depends onexactly what your application is! You may be dealing with a set of equations whosedependent variables differ enormously in magnitude.
In that case, you probablywant to use fractional errors, ∆0 = y, where is the number like 10−6 or whatever.On the other hand, you may have oscillatory functions that pass through zero butare bounded by some maximum values. In that case you probably want to set ∆0equal to times those maximum values.A convenient way to fold these considerations into a generally useful stepperroutine is this: One of the arguments of the routine will of course be the vector ofdependent variables at the beginning of a proposed step.
Call that y[1..n]. Letus require the user to specify for each step another, corresponding, vector argumentyscal[1..n], and also an overall tolerance level eps. Then the desired accuracySample 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.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
