Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 14
Текст из файла (страница 14)
It implies that the motionof the rocket is independent of its mass. For performing a now-legendary experiment withrocks of different sizes dropping from the Tower of Pisa, Galileo demonstrated that fact toan incredulous world.At any rate, (2.10.3) now reads asx00 = −KMEKMM.+x(D − x)2(2.10.5)We can make this equation a good bit prettier by changing the units of distance and timefrom miles and seconds (or whatever) to a set of more natural units for the problem.For our unit of distance we choose R, the radius of the earth.
If we divide (2.10.5)through by R, we can write the result asxR00KMEKMMR33= − R 2 + XRDR−xR2 .(2.10.6)Now instead of the unknown function x(t), we define y(t) = x(t)/R. Then y(t) is theposition of the rocket, expressed in earth radii, at time t. Further, the ratio D/R thatoccurs in (2.10.6) is a dimensionless quantity, whose numerical value is about 60. Hence(2.10.6) has now been transformed to00y =KME3− R2y+KMMR3(60 − y)2.(2.10.7)Next we tackle the new time units.
Since y is now dimensionless, the dimension of theleft side of the equation is the reciprocal of the square of a time. If we look next at thefirst term on the right, which of course has the same dimension, we see that the quantityR3 /KME is the square of a time, sosT0 =R3KME(2.10.8)is a time. Its numerical value is easier to calculate if we change the formula first, as follows.Consider a body of mass m on the surface of the earth. Its weight is the magnitude ofthe force exerted on it by the earth’s gravity, namely KME m/R2 .
Its weight is also equalto m times the acceleration of the body, namely the acceleration due to gravity, usuallydenoted by g, and having the value 32.2 feet/sec2 .62The Numerical Solution of Differential EquationsIt follows thatKMe m= mg,R2and if we substitute into (2.10.8) we find that our time unit issT0 =R.g(2.10.9)(2.10.10)We take R = 4000 miles, and find T0 is about 13 minutes and 30 seconds. We propose tomeasure time in units of T0 . To that end, we multiply through equation (2.10.7) by T0 andgetMM1MET02 y 00 = − 2 +.(2.10.11)y(60 − y)2The ratio of the mass MM of the moon to the mass ME of the earth is about 0.012.Furthermore, we will now introduce a new independent variable τ and a new dependentvariable u = u(τ ) by the relationsu(τ ) = y(τ T0 ) ;t = τ T0 .(2.10.12)Thus, u(τ ) represents the position of the rocket, measured in units of the radius of theearth, at a time τ that is measured in units of T0 , i.e., in units of 13.5 minutes.The substitution of (2.10.12) into (2.10.11) yields the differential equation for the scaleddistance u(τ ) as a function of the scaled time τ in the formu00 = −10.012+.u2 (60 − u)2(2.10.13)Finally we must translate the initial conditions (2.10.4) into conditions on the newvariables.
The first condition is easy: u(0) = 1. Next, if we differentiate (2.10.12) and setτ = 0 we getT0 VVu0 (0) =.(2.10.14)=RR/T0This is a ratio of two velocities. In the numerator is the velocity with which the rocket islaunched. What is the significance of the velocity R/T0 ?We claim that it is, aside from a numerical factor, the escape velocity from the earth, ifthere were no moon.
Perhaps the quickest way to see this is to go back to equation (2.10.11)and drop the second term on the right-hand side (the one that comes from the moon). Thenwe will be looking at the differential equation that would govern the motion if the moonwere absent. This equation can be solved. Multiply bth sides by 2y 0 , and it becomesT02 (y 0 )20 0=2y,(2.10.15)and integration yieldsT02 (y 0 )2 =2+ C.y(2.10.16)2.10 Case study: Rocket to the moon63Now let t = 0 and find that C = T02 V 2 /R2 − 2, soT02 (y 0 )2!T02 V 2−2 .R22= −y(2.10.17)Suppose the rocket is launched with sufficient initial velocity to escape from the earth.
Thenthe function y(t) will grow without bound. Hence let y → ∞ on the right side of (2.10.17).For all values of y, the left side is a square, and therefore a nonnegative quantity. Hencethe right side, which approaches the constant C, must also be nonnegative. Thus C ≥ 0 or,equivalently√ RV ≥ 2 .(2.10.18)T0Thus, if the rocket√ escapes, then (2.10.18) is true, and the converse is easy to show also.Hence the quantity 2 R/T0 is the escape velocity from the earth. We shall denote it byVesc . Its numerical value is approximately 25, 145 miles per hour.0 (t) into initialNow we can return to (2.10.12) to translate the initial conditions on x√conditions on u0 (τ ).
In terms of the escape velocity, it becomes u0 (0) = 2 V /Vesc . Wemight say that if we choose to measure distance in units of earth radii, and time in units√ ofT0 , then velocities turn out to be measured in units of escape velocity, aside from the 2.In summary, the differential equation and the initial conditions have the final formu00 = −10.012+2u(60 − u)2u(0) = 1√ Vu0 (0) = 2Vesc(2.10.19)Since that was all so easy, let’s try the two-dimensional case next. Here, the earth iscentered at the origin of the xy-plane, and the moon is moving. Let the coordinates of themoon at time t be (xm (t), ym (t)).
For example, if we take the orbit of the moon to be acircle of radius D, then we would have xm = D cos(ωt) and ym (t) = D sin(ωt).If we put the rocket at a generic position (x(t), y(t)) on the way to the moon, then wehave the configuration shown in figure 1.16.2.Consider the net force on the rocket in the x direction. It is given byFx = −KME m cos θKMM m cos ψ+,x2 + y 2(x − xm )2 + (y − ym )2(2.10.20)where the angles θ and ψ are shown in figure 1.16.2.
From that figure, we see thatqandcos θ = x x2 + y 2(2.10.21)xm − xcos ψ = p.(xm − x)2 + (ym − y)2(2.10.22)64The Numerical Solution of Differential EquationsvMoonyM (t)y(t)s ψ rocket θzEarthx(t)xM (t)Figure 2.3: The 2D Moon RocketNow we substitute into (2.10.20), and equate the force in the x direction to mx00 (t), toobtain the differential equationmx00 (t) = −KME mxKMM m(xm − x)+.223/2(x + y )((xm − x)2 + (ym − y)2 )3/2(2.10.23)If we carry out a similar analysis for the y-component of the force on the rocket, we getmy 00 (t) = −KME myKMM m(ym − y)+.(x2 + y 2 )3/2 ((xm − x)2 +)ym − y)2 )3/2(2.10.24)We are now looking at two (even nastier) simultaneous differential equations of thesecond order in the two unknown functions x(t), y(t) that describe the position of therocket.
To go with these equations, we need four initial conditions. We will suppose that attime t = 0, the rocket is on the earth’s surface, at the point (R, 0). Further, at time t = 0,it will be fired with an initial speed of V , in a direction that makes an angle α with thepositive x-axis. Thus, our initial conditions are(x(0) = R ;y(0) = 0.x0 (0) = V cos α ; y 0 (0) = V sin α(2.10.25)The problem has now been completely defined. It remains to change the units into the samenatural dimensions of distance and time that were used in the one-dimensional problem.This time we leave the details to the reader, and give only the results.
If u(τ ) and v(τ )denote the x and y coordinates of the rocket, measured in units of earth radii, at a time τmeasured in units of T0 (see (2.10.10)), then it turns out the u and v satisfy the differentialequationsu0.012(um − u)u00 = − 2+(u + v 2 )3/2((um − u)2 + (vm − v)2 )3/2(2.10.26)v0.012(vm − v)00v =− 2+.(u + v 2 )3/2((um − u)2 + (vm − v)2 )3/2Furthermore, the initial data (2.10.25) take the form u(0) = 1 ;v(0) = 0√ V cos α√ V sin α u0 (0) = 2; v 0 (0) = 2.VescVesc(2.10.27)2.11 Maple programs for the trapezoidal rule65In these equations, the functions um (τ ) and vm (τ ) are the x and y coordinates of the moon,in units of R, at the time τ .
Just to be specific, let’s decree that the moon is in a circularorbit of radius 60R, and that it completes a revolution every twenty eight days. Then, aftera brief session with a hand calculator or a computer, we discover that the equationsum (τ ) = 60 cos (0.002103745τ )vm (τ ) = 60 sin (0.002103745τ )(2.10.28)represent the position of the moon.2.11Maple programs for the trapezoidal ruleIn this section we will first display a complete Maple program that can numerically solve asystem of ordinary differential equations of the first order together with given initial values.After discussing those programs, we will illustrate their operation by doing the numericalsolution of the one dimensional moon rocket problem.We will employ Euler’s method to predict the values of the unknowns at the next pointx + h from their values at x, and then we will apply the trapezoidal rule to correct thesepredicted values until sufficient convergence has occurred.First, here is the program that does the Euler method prediction.>>>>>>>>>>>>eulermethod:=proc(yin,x,h,f)local yout,ll,i:# Given the array yin of unknowns at x, uses Euler method to return# the array of values of the unknowns at x+h.
The function f(x,y) is# the array-valued right hand side of the given system of ODE’s.ll:=nops(yin):yout:=[]:for i from 1 to ll doyout:=[op(yout),yin[i]+h*f(x,yin,i)];od:RETURN(yout):end:Next, here is the program that takes as input an array of guessed values of the unknownsat x + h and refines the guess to convergence using the trapezoidal rule.>>>>>>>>>traprule:=proc(yin,x,h,eps,f)local ynew,yfirst,ll,toofar,yguess,i,allnear,dist;# Input is the array yin of values of the unknowns at x.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
