Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 49
Текст из файла (страница 49)
References [4], [8], and [9] provide usefulorientation, especially in regard to the software available. State-of-the-art codes areavailable from netlib: RKSUITE [l] makes available explicit Runge-Kutta formulasof three different orders. ODE/STEP/INTRP [10] is a variable step, variable orderAdams-Bashforth-Moulton code. The methods of RKSUITE and ODE/STEP/INTRPare not appropriate for stiff problems. Both suites of codes diagnose stiffness when itis responsible for unsatisfactory performance.
VODE [2] is a variable step, variableorder code that makes available two kinds of methods, Adams-Moulton formulas anda variant of the BDFs. The computing environments MATLAB and MATHCAD providecodes based on a variety of methods, including some not mentioned here, but the codethat is to be tried first (assuming that the problem is not stiff) is an explicit RungeKutta code. Mathematics provides a single code that, like VODE, makes availableboth Adams-Moulton methods and the BDFs. It is unusual in that the code attemptsto recognize stiffness and select an appropriate method automatically.6.8 CASE STUDY 6The restricted three-body problem is obtained from Newton’s equations of motionfor the gravitational attraction of three bodies when one has a mass infinitesimal incomparison to the other two.
For example, the position (x,y) of a spaceship or satellitemoving under the influence of the earth and the moon in a coordinate system rotatingso as to keep the positions of the earth and moon fixed changes according toHereand µ = 1/82.45, µ* = 1 - µ. More insight is possible when the general equationsof motion are reduced to those of the restricted three-body problem, but it is still notpossible to determine orbits analytically. A search using high precision computationidentified several periodic orbits. One has initial conditions x(0) = 1.2, x´(0) = 0,y(0) = 0, and y´(0) = -1.04935750983031990726.
The period of this orbit is aboutT = 6.19216933131963970674. Integration of this problem with Rke is straightforward after the equations are written as a first order system by introducing the vectorof unknowns y(t) = (x(t),y(t),x´(t),y´(t))T. The orbit displayed in Figure 6.3 wascomputed using 10-6 for TOL and all components of THRESHOLD. Although the6.8 CASE STUDY 6245Figure 6.3 Rke solution of the satellite problem.analytical solution is not known, we do know y(T) = y(0) = Y 0 because the orbithas period T. Using this known value, we measured the global error of the approximation YN to y(T) and tested whether the computed orbit is periodic by computing||Y N - Y0||.
The discrepancy turned out to be about 6.1 × 10-5, which is about whatwe would expect for the local error tolerances given the code. The figure shows the“natural” output, the values at each step, connected by straight lines in the manner typical of plotting packages. At this tolerance the natural output is sufficiently dense thatthe curve appears to be smooth. However, at less stringent tolerances it was found thatin portions of the integration, the step size is so large that the curve is visibly composedof straight lines.
This unsatisfactory situation can be remedied easily by computing inexpensively with Yvalue the additional output values needed for a smooth graph. Aless efficient alternative that is acceptable in this particular instance is to limit the stepsize so as to force Rke to take more steps.Conservation of energy has a special form for the restricted three-body problem.The Jacobi integral isA little calculation shows that the derivative dJ/dt is zero when it is evaluated at arguments x(t), y(t) that satisfy the differential equations. This leads to the conservationlawexpressing the fact that the Jacobi integral is constant along a solution of the restrictedthree-body problem.
We monitored the integration by computing G(t) at each step.For the tolerances specified, it was never larger in magnitude than about 1.8 × 10-5.Many differential equations have conservation laws that arise naturally from physical principles, but others satisfy laws that are not so readily interpreted. Recall, for246CHAPTER 6 ORDINARY DIFFERENTIAL EQUATIONSexample, the conservation law we used in Case Study 4 to solve the Lotka-Volterraequations in the phase plane. Conservation laws are consequences of the form of thedifferential equations. It does not follow that numerical approximations will satisfythe laws, and generally speaking they do not.
However, they must satisfy the lawsapproximately. To see this, suppose that G(t) = G(t,y (t)) = 0 for any solution y(t) ofa differential equation y´ = F(t, y ). If y ny(tn), then linearization tells us thatEvidently the residual of the numerical solution in the conservation law is of the sameorder as the global error of the numerical solution, y(tn) - y n. This observation helpsus understand the size of the residual we found in integrating the periodic earth-moonorbit.
It is worth remarking that solutions of a system of differential equations mightsatisfy several conservation laws.The conservation laws mentioned so far are nonlinear, but others arising fromconservation of mass, charge balance, and the like are linear. A linear conservationlaw for the equation y´ = F( t, y) arises mathematically when there is a constant vectorv such that vT F( t,u) = 0 for all arguments (t,u). If y(t) is a solution of the equation,thenimplying that G(t) = v T y (t) - v T y (0) = 0. A simple example is provided by a systemof equations that describes a certain radioactive decay chain:The right-hand sides here sum to zero, hence the system satisfies a linear conservationlaw with v T = (1, 1,. . .
, l) T. Figure 6.4 shows the solution of this system with initialTcondition y(0) = (1, 0,. . . ,0) computed using Rke with TOL equal to 10-3 and allcomponents of THRESHOLD equal to 10-10. Despite the modest relative accuracyrequested, the error in the conservation law was at the roundoff level, specifically amaximum of 4.4 × 10-l6 in the MATLAB environment on the workstation we used.This illustrates an interesting fact: all standard numerical methods produce approximations that satisfy linear conservation laws exactly. This is easy to show for explicitRunge-Kutta formulas.
In advancing from xn to xn + h, such a formula has the formEach stage K i has the form F( x*,Y*) for arguments (x*,Y*) that are defined in termsof xn, h, and the previous stages. The details do not matter, for all we need here is thatv T K i = 0 because vT F( t,u) = 0 for all arguments (t,u). It then follows immediately6.8 CASE STUDY 6247Figure 6.4 Rke solution of the radioactive decay problem.that v T Y n + l = v T Y n for all n. We start the integration with the given initial values sothat Y0 = y(0). This implies that for all n,Put differently, G( xn ,Y n ) = 0 for all n.
Accordingly, it was no accident that the numerical solution computed by Rke satisfied the conservation law to roundoff error, that iswhat we should expect of a linear conservation law.Returning now to the satellite problem, suppose we would like to know when thesatellite is nearest and farthest from earth. The distance to the satellite isso we look for the extrema of this function by finding the times t for which d’(t) =0.
Let us avoid square roots by working with the square of the distance, D(t). Thederivative of this function isin the original variables and those of the first order system, respectively. Notice that forthe orbit we study, the initial distance d(0) = 1.2 is an extremum because D´(0) = 0.This will afford a check on our computations because the orbit is periodic and thesame must be true of d(T).
We want to compute the roots of F(t) = D´ (t) = 0. Ateach step we have an approximation to the solution that allows us to evaluate F(t),so when we reach x after a step of size step, we ask if F(x - step) F(x) < 0. If so,we have found just the kind of bracket for a root that we need for applying Zero tolocate a root precisely. Evaluation of F(t) is easy enough; we just invoke Yvalue toget an approximation to y at t, and then use these approximations in the expressionfor D´(t).
There is one snag, which is a very common one when combining items248CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSof mathematical software. We invoke Zero with the name of a function F of justone argument t. However, to evaluate F(t) we must invoke Yvalue, and it requiresthree other arguments, namely x, step, and ycoeff, the array of coefficients definingthe interpolating polynomial over [x - step,x], returned by Rke. Somehow we mustcommunicate this output from Rke in the main program to the function for evaluatingF. There can be more than one way to do this, depending on the language, but it isalways possible to do it by means of global variables.
As a specific example, we codedthe function in MATLAB asfunction Ft = F(t)global x step ycoeffyt = Yvalue(t,x,step,ycoeff);Ft = 2*(yt(1:2)´*yt(3:4));The quantities listed in the second line of this code are computed by Rke in the mainprogram. By duplicating the line there, the quantities are made accessible from thefunction F.
Proceeding in the manner described with the same tolerances in both Rkeand Zero, we foundtime1.458573.096064.733546.19210distance0.033381.262440.033381.19998The last extremum occurs at a time that agrees to about five figures with the period andagrees with the initial distance to about the same accuracy, which is quite reasonablegiven the tolerances. In Figure 6.3 the points nearest the earth and farthest away thatwe found in this way are indicated by circles.REFERENCES1. R. Brankin, I. Gladwell, and L. Shampine, “RKSUITE: a suite of Runge-Kutta codes for theinitial value problem for ODES,” Softreport 92-S1, Math.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
