Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 46
Текст из файла (страница 46)
It corresponds to the adaptivequadrature schemes of Chapter 5.232CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSEngland’s formulas are as follows:One drawback of conventional algorithms for solving initial value problems is thatthey produce a table of approximate values while the mathematical solution y(x) is acontinuous function. It is possible to approximate the solution for all x by interpolation. At the beginning of a step we have Yn and formSimilarly, at the start of the second half-step we haveandThe code Rke does local extrapolation, hence reports the fifth order solutionas itsapproximate solution at xn+1. It will be called Y n+1 on the next step.
By programming6.5 THE ALGORITHMS233the procedure so thatis evaluated in the current step and later used asthe K0 of the next step, we have an approximation to Y´ (xn+1) In this way we haveapproximations to Y(x) and Y´(x) at xn,xn + h/2, and xn+1 in the course of taking astep. It is natural to interpolate these data by quintic (degree 5) Hermite interpolation.It turns out that in a certain sense this interpolant is as accurate an approximation toY (x) on (xn,xn+1) as Y n+l is to Y(xn+1).
Notice that only information generated inthe current step is needed. Furthermore, the quintic polynomial on [xn,xn+l] agrees invalue and slope at xn+l with the quintic polynomial on [xn+l,xn+2]. Thus the piecewisepolynomial function resulting from this scheme is continuous and has a continuousderivative on all of [a,b]. An interpolation capability can greatly increase the efficiencyof an ordinary differential equation solver because the step size can be selected solelyfor reliability and efficiency.6.5 THE ALGORITHMSIt is easy enough to write a code based on a pair of Runge-Kutta formulas like England’s, but it is not easy to write a code of production quality.
A subdiscipline ofnumerical analysis called mathematical software has developed that concerns itselfwith such tasks. References [l], [2], [4], [6], [8], [9], and [10] discuss the issues forordinary differential equations at length. The code Rke that we provide is significantlymore complex than the codes in other chapters, so its description is far more detailed.In this section we consider briefly portions of the code with the aim of explainingsome of the care needed in converting a pair of formulas into a production code. It ishard to come to a full understanding of such complex codes because decisions aboutone algorithmic question usually depend on those made in connection with severalothers.
Perhaps a good place to start is with the way error is to be measured.For simplicity, only the scalar case is described; the vector case is handled similarly. At each step the code attempts to keep the local error less than a tolerancespecified by the user:|local error| < τ.How this error is measured is important. A reasonable error tolerance will depend onthe size of the solution. Because this size is usually not known in advance, a good wayto proceed is to measure the error relative to the size computed by the code:It is not so clear what we should take here as the “size” of y. Besides needing areasonable definition of size when a solution component vanishes, we need to avoidthe technical difficulty of dividing by zero. We have a value yn at the beginning of thestep, an approximate solution of order 4, yn+½, at half the step, and two approximatesolutions of orders 4 and 5,at the end of the step.
A reasonable way todefine the size of y over the step is to average these magnitudes, taking account of thefact that two of the values approximate y(x) at the same point:234CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSWith this definition, it is unlikely that a zero value for wt would arise unless the solution underflows identically to zero. The local error is approximated byIfthe error is estimated to be zero, and there is no need to compute theweighted error.
Ifthe definition of wt implies that wt > 0. Becausethere is then no difficulty in performing the testProceeding in this way, we have a good measure of “size” and we avoid numericaldifficulties. Nonetheless, a pure relative error control may not be appropriate when thesolution vanishes or becomes “small.” What constitutes “small” necessarily dependson the scale of the problem and must be supplied from the insight of the problem solveror a preliminary computation. The code Rke asks the user to specify a threshold valueand measures the error relative to max(wt,threshold). This tells the code that whenthe magnitude of the solution drops below the threshold, the user is interested only inan absolute error.Some attention must be paid to the arithmetic of the computer being used. Theerror control provided in Rke stresses pure relative error.
It makes no sense to askfor a numerical solution more accurate than the correctly rounded true solution. Toavoid difficulties the code insists that τ not be smaller than 10 units of roundoff andh not be too small for the precision available. These elementary precautions are veryhelpful. They are the kinds of things that distinguish mathematical software fromresearch codes.Suppose we have just tried a step of length h from xn and formed the local errorestimateA Taylor series expansion of the local error leads to an expression of the formlocal error = h5 Φ (xn,yn) + O(h6).Earlier we wrote out Φ explicitly for some low order formulas. If h is small enough,If the step is a success, meaning thatwe would like to estimate a suitable step size H for the next step.
The largest step sizepossible hasThe function Φ is (usually) smooth so that6.5 THE ALGORITHMS235Writing H = αh we then findso that the “optimal” H isIt is much too bold to try this H because if it is even slightly too large, the step willfail and this is expensive. Besides, after making all those approximations, we shouldnot take H too seriously. In practice a fraction of H is tried, or equivalently a fractionof τ is used in computing H, and an efficient fraction determined by experiment. InRke the tolerance aimed at is 0.6τ.
This is equivalent to using about nine-tenths of the“optimal” H.The same argument is used to obtain the step size for trying again after a failedstep. In either case we must program the process with some care. For example, wemust deal with the possibility that the estimated local error is zero. This technicaldifficulty highlights the fact that large increases or decreases cannot be justified by thearguments made. For this reason changes are limited to an order of magnitude. If alarge change is truly called for, this action allows a large change to be accumulatedover a few steps without the disastrous possibilities opened up by a large change in asingle step.
In the case of a failed step we must be especially cautious about taking theestimated error at face value. In Rke we try the predicted value once, but if it fails, wesimply halve the step size until we get success.For the numerical solution of ordinary differential equations there are two difficultranges of tolerances. It is to be expected that tolerances near limiting precision aredifficult, but it turns out that nominal tolerances are also difficult. Often users thinkthat because engineering accuracy of 10% will suffice in their use of the results, theycan keep their costs down by specifying such a large tolerance.
This may result incheap results that are not reliable because the local error test may not keep the stepsize small enough to justify the approximations used throughout the code. Even if theresults are reliable, they can be far from what is desired because at crude tolerances thetrue, or global, errors can be much larger than the local errors controlled by the code.It is prudent to ask for accuracy of at least a couple of digits, so that the error controlof Rke emphasizes relative error and it is required that the relative tolerance τ < 0.01.EXERCISES6.15 Implement Euler’s method and a local error estimatorbased on Heun’s method.
Apply it to the problemy´ = 10(y-x), y(0) = 1/10and compare the estimated local error to the true localerror. Also compare the global error at several pointsto the general size of the local errors made in the computations up to this point.236CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONS6.6 THE CODE RKEThe routine Rke solves the initial value problem for a system of first order ordinarydifferential equations of the formA typical call for Rke in FORTRAN isCALL RKE (F, NEQ, X, Y, H, FIRST, TOL, THRES, FLAG, STEP,YCOEFF, SCR, NDIM)In C it isRke (f, neq, &x, y, &h, &first, tol, threshold, &flag, &step, ycoeff);and it isRke (f, neq, x, y, h, first, tol, threshold, flag, step, ycoeff);in C++.Input parameters to Rke are F, the name of the routine that defines the differentialequations [i.e., F(x, Y)]; NEQ, the number of first order differential equations to be integrated; X, the initial value of the independent variable; Y, an array of dimension NEQcontaining the values of the solution components at X; H, step size for the current step(its sign determines the direction of integration); FIRST, a variable indicating whetherthis is the first or a subsequent step; a scalar TOL and a vector THRES (or thresholdin C and C++) are tolerances for the local error control; and NDIM > 6× NEQ, thedimension of the output vector YCOEFF and in FORTRAN of the auxiliary storagevector SCR.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
