Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 13
Текст из файла (страница 13)
Suppose we would like tokeep the single-step error in the neighborhood of 10−8 . We might then adopt, say 5 × 10−8as the upper limit of tolerable error and, for instance, 10−9 as the lower limit.Why should we have a lower limit? If the calculation is being done with more precisionthan necessary, the step size will be smaller than needed, and we will be wasting computertime as well as possibly building up roundoff error.58The Numerical Solution of Differential EquationsNow that we have fixed these limits we should be under no delusion that our computedsolution will depart from the true solution by no more than 5 × 10−8 , or whatever.
Whatwe are controlling is the one-step truncation error, a lot bettern than nothing, but certainlynot the same as the total accumulated error.With the upper and lower tolerances set, we embark on the computation. First, sincethe midpoint method needs two backward values before it can be used, something specialwill have to be done to get the procedure moving at the start.
This is typical of numericalsolution of differential equations. Special procedures are needed at the beginning to buildup enough computed values so that the predictor-corrector formulas can be used thereafter,or at least until the step size is changed.In the present case, since we’re going to use the trapezoidal corrector anyway, we mightas well use the trapezoidal rule, unassisted by midpoint, with complete iteration to convergence, to get the value of y at the first point x0 + h beyond the initial point x0 .Now we have two consecutive values of y, and the generic calculation can begin. Fromany two consecutive values, the midpoint rule is used to predict the next value of y.
Thispredicted value is also saved for future use in error estimation. The predicted value is thenrefined by the trapezoidal rule.With the trapezoidal value in hand, the local error is then estimated by calculatingone-fifth of the difference between that value and the midpoint guess.If the absolute value of the local error lies between the preset limits 10−9 and 5 × 10−8 ,then we just go on to the next step. This means that we augment x by h, and move back thenewer values of the unknown function to the locations that hold older values (we remember,at any moment, just two past values).Otherwise, suppose the local error was too large.
Then we must reduce the step sizeh, say by cutting it in half. When all this is done, some message should be printed outthat announces the change, and then we should restart the procedure, with the new valueof h, from the “farthest backward” value of x for which we still have the corresponding yin memory.
One reason for this is that we may find out right at the outset that our veryfirst choice of h is too large, and perhaps it may need to be halved, say, three times beforethe errors are tolerable. Then we would like to restart each time from the same originallygiven data point x0 , rather than let the computation creep forward a little bit with stepsizes that are too large for comfort.Finally, suppose the local error was too small. Then we double the step size, print amessage to that effect, and restart, again from the smallest possible value of x.Now let’s apply the philosophy of structured programming to see how the whole thingshould be organized. We ask first for the major logical blocks into which the computationis divided. In this case we see(i) a procedure midpt.
Input to this procedure will be x, h, yn−1 , yn . Output from itwill be yn+1 computed from the midpoint formula. No arrays are involved. The threevalues of y in question occupy just three memory locations. The leading statement inthis routine might be2.9 Controlling the step size59midpt:=proc(x,h,y0,ym1,n);and its return statement would be return(yp1);. One might think that it is scarcelynecessary to have a separate subroutine for such a simple calculation. The spirit ofstructured programming dictates otherwise.
Someday one might want to change fromthe midpoint predictor to some other predictor. If organized as a subrouting, then it’squite easy to disentangle it from the program and replace it. This is the “modular”arpproach.(ii) a procedure trapez. This routine will be called from two or three different places in themain routine: when starting, when restarting with a new value of h, and in a genericstep, wher eit is the corrector formula. Operation of thie routine is different, too,depending on the circumstances. When starting or restarting, there is no externallysupplied guess to help it.
It must find its own way to convergence. On a generic stepof the integration, however, we want it to use the prediction supplied by midpt, andthen just do a single correction.One way to handle this is to use a logical variable start. If trapez is called with start= TRUE, then the subroutine would supply a final converged value without looking for anyinput guess. Suppose the first line of trapez istrapez:=proc(x,h,yin,yguess,n,start,eps);and its return statement is return(yout);. When called with start = TRUE, then theroutine might use yin as its first guess to yout, and iterate to convergence from there, usingeps as its convergence criterion.
If start = FALSE, it will take yguess as an estimate ofyout, then use the trapezoidal rule just once, to refine the value as yout.The combination of these two modules plus a small main program that calls them asneeded, constitutes the whole package.
Each of the subroutines and the main routine shouldbe heavily documented in a self-contained way. That is, the descriptions of trapez, and ofmidpt, should precisely explain their operation as separate modules, with their own inputsand outputs, quite independently of the way they happen to be used by the main programin this problem. It should be possible to unplug a module from this application and use itwithout change in another.
The documentation of a procedure should never make referenceto how it is used by a certain calling program, but should describe only how it transformsits own inputs into its own outputs.In the next section we are going to take an intermission from the study of integrationrules in order to discuss an actual physical problem, the flight of a spacecraft to the moon.This problem will need the best methods that we can find for a successful, accurate solution.Then in section 2.12, we’ll return to the modules discussed above, and will display completecomputer programs that carry them out. In the meantine, the reader might enjoy trying towrite such a complete program now, and comparing it with the one in that section.60The Numerical Solution of Differential Equationsx=0x=Rx=Dx(t)lx=D−rs -EarthRocketfMoonFigure 2.2: 1D Moon Rocket2.10Case study: Rocket to the moonNow we have a reasonably powerful apparatus for integration of initial-value problems andsystems, including the automatic regulation of step size and built-in error estimation.
Inorder to try out this software on a problem that will use all of its capability, in this sectionwe are going to derive the differential equations that govern the flight of a rocket to themoon. We do this first in a one-dimensional model, and then in two dimensions. It willbe very useful to have these equations available for testing various proposed integrationmethods. Great accuracy will be needed, and the ability to chenge the step size, bothto increase it and the devrease it, will be essential, or else the computation will becomeintolerably long.
The variety of solutions that can be obtained is quite remarkable.First, in the one-dimensional simplified model, we place the center of the earth at theorigin of the x-axis, and let R denote the earth’s radius. At the point x = D we place themoon, and we let its radius be r. Finally, at a position x = x(t) is our rocket, making itsway towards the moon.We will use Newton’s law of gravitation to find the net gravitational force on therocket,and equate it to the mass of the rocket times its acceleration (Newton’s second law ofmotion). According to Newton’s law of gravitation, the gravitational force exerted by onebody on another is proportional to the product of their masses and inversely proportionalto the square of the distance between them.
If we use K for the constant of proportionality,then the force on the rocket due to the earth is−KME m,x2(2.10.1)whereas the force on the rocket due to the moon’s gravity isKMM m(D − x)2(2.10.2)where ME , MM and m are, respectively, the masses of the earth, the moon and the rocket.The acceleration of the rocket is of course x00 (t), and so the assertion that the net forceis equal to mass times acceleration takes the form:mx00 = −KME mMM m+K.x2(D − x)2(2.10.3)This is a (nasty) differential equation of second order in the unknown function x(t), theposition of the rocket at time t. Note the nonlinear way in which this unknown functionappears on the right-hand side.2.10 Case study: Rocket to the moon61A second-order differential equations deserves two initial values, and we will oblige.First.
let’s agree that at time t = 0 the rocket was on the surface of the earth, and second,that the rocket was fired at the moon with a certain initial velocity V . Hence, the initialconditions that go with (2.10.3) arex0 (0) = V.x(0) = R ;(2.10.4)Now, just a quick glance at (2.10.3) shows that m cancels out, so let’s remove it, butnot before pointing out the immense significance of that fact.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
