Nash - Scientific Computing with PCs (523165), страница 48
Текст из файла (страница 48)
In general, the vector quantities at and ac are given as follows:(16.2.3)at= 2 Ω × (dr/dt)where Ω (Greek capital omega) is the angular velocity of rotation of the coordinate system. Note thevector cross-product. Further, using a vector dot-product,(16.2.4)ac= - Ω2 r + (dΩ/dt) × r + Ω (Ω · r)By considering only trajectories in the plane of the earth-moon rotation, we greatly simplify thesecorrections. Consider the coordinates defined in Figure 16.1.1.
Note that the earth has its center at (-d, 0),the moon at (D-d, 0), the spacecraft at (x,y). Newton’s gravitation law gives the total gravitational forceon the capsule as(16.2.5)F = - G w {(M / r12) i + (m / r22) j}where i and j are the unit vectors along the earth-spacecraft and spacecraft-moon lines, M ≈ 5.977E+24kg. and m ≈ 7.349E+22 kg. are the earth and moon masses, and G ≈ 6.6733E-11 m3 kg-1 s-2 is thegravitational constant. The distance center-to-center from the earth to the moon is D ≈ 3.844E+8 m. Thecenter of mass of the earth-moon system, which we use for our coordinate origin, is therefore at a distanced ≈ 4.669E+4 m from the center of the earth, and is actually inside the earth.
The angular velocity of theearth-moon rotation Ω ≈ 9.583E-3 hr-1 = 2.661E-6 s-1. This gives a period T = 2 π / Ω ≈ 655.66 hrs ≈ 27.32days (Tennent, 1971).Combining and simplifying the equations yields the linked second order ordinary differential equations140Copyright © 1984, 1994 J C & M M NashNash Information Services Inc., 1975 Bel Air Drive, Ottawa, ON K2C 0X1 Canada(16.2.6a)d 2 x/dt 2SCIENTIFIC COMPUTING WITH PCsCopy for:Dr. Dobb’s Journal= - G M (x + a)/r1 3 - G m (x - (D-a))/r2 3+ Ω 2 x + 2 Ω dy/dt= g1 (x, x’, y, y’, t)(16.2.6b)d 2y/dt 2= - G M y/r1 3 - G m y/r2 3 + Ω2 y - 2 Ω dx/dt= g2 (x, x’, y, y’, t)where we use the functionsg1 and g2 as a shorthand for later equations.16.3 Methods and AlgorithmsSolution of equations (16.2.6) from an initial position with a known initial velocity (i.e., x, y, dx/dt, dy/dtknown at t=0) is not always straightforward because programs and packages such as RKF45 (Forsythe,Malcolm and Moler, 1977) are designed to solve a set of first order differential equations.
We cantransform our pair of equations into four first-order ordinary differential equations (ODEs) by using newvariables u and v(16.3.1)u = x’ and v = y’(16.3.2a)dx/dt = u(16.3.2b)du/dt = g1 (x, u, y, v, t)(16.3.2c)du/dt = v(16.3.2d)du/dt = g2 (x, u, y, v, t)We therefore need methods that handle either linked first or second order differential equations. Thereis a large literature on such problems (e.g., Kahaner, Moler, and Nash S G, 1989, Chapter 8).16.4 Programs or PackagesThe software tools for the solution of such a problem fall into three categories:•Existing library subroutines or programs to which we supply our problem by means of customprogram code;•Problem solving environments of a general nature in which our particular problem is expressed insome algebraic language;•Special-purpose software for the solution of differential equation problems.The first approach was the only option as recently as 1980.
For example, in 1983 we used a BASICtranslation (by Kris Stewart) of the FORTRAN program RKF45 given in Forsythe, Malcolm and Moler(1977). The user has to write driver code and write or link to graphical routines.The last approach, using special-purpose software, is clearly the least work for the user if the software iswell-designed. Since spacecraft trajectories are needed by many agencies, such software likely exists, butaccess to it may be less easy to obtain.
For the more general task of solving sets of differential equationswe have examples like PLOD (PLotted solutions of Ordinary Differential equations, Kahaner et al., 1989).This package builds and executes FORTRAN programs from a description of the first order differentialequations (16.3.2). The equations are entered almost as they appear, though we must also supply16: SPACECRAFT TRAJECTORIES141expressions for the distances r1 and r2 as well as physical quantities that appear. PLOD lets thisinformation be placed in a file for reuse.
Figure 16.4.1 shows the file APOLC.MOD with the starting values(16.4.1)x = 461.28; y = 0.0; u = 0.0; v = -3.865412Figure 16.4.1PLOD input file for Apollo spacecraft trajectoriesAPOLC.TMe = 5169.26643mm = 63.5585om = 9.5827E-03tmu = 379.73r1 = sqrt( (x + mu)**2 + y**2 )r2 = sqrt( (x - tmu)**2 + y**2 )x’=uy’=vu’ = 2*om*v + x*om**2 - Me*(x+mu)/r1**3 - mm*(x-tmu)/r2**3v’ = -2*om*u + y*om**2 - Me*y/r1**3 - mm*y/r2**316.5 Some solutionsFigure 16.5.1 shows a plot of the trajectory computed by PLOD for initial conditions (16.4.1).
The PLODinput and output illustrate the "minor" but very annoying difficulties that accompany any real-worldscientific problem. First, we must ensure all the constants in the equations are in consistent units.Attaching the units to all quantities and performing careful algebra with these and the conversion factors,we obtained the input in Figure 16.4.1. A second annoyance arose here — PLOD forced us to use distancesand velocities expressed in terms of megameters (106 meters) because it would not accept the largenumbers present when kilometers were used. A third annoyance is that PLOD does not appear to allowus to save high quality image information.
Nor can we draw the earth and moon on the output or putlabels on. Finally, while it is helpful to plot the position of the spacecraft at regular time intervals (to giveus an idea of how fast the spacecraft is moving at different points along its trajectory), PLOD (and mostother differential equation integrators) vary the integration time step for efficiency and error control. Thuswe have a set of data (t, x, y) with irregular time steps.Avoiding scale and numerical errors in input quantities is important.
A spreadsheet or symbolic algebrapackage lets us perform unit conversion calculations. By making a program display the data that has beenentered we may catch data storage "bugs". This is how we found that PLOD was "truncating" the largenumbers for distances in kilometers.Equispaced time points could be obtained by outputting the (t, x, y) data and using a separateinterpolation and graphing package. For each new time point new_t, with two adjacent original timepoints such that(16.5.1)ti-1 <= new_t < tia linear interpolation gives a corresponding approximate(16.5.2)x valuenew_x = xi-1 + (xi - xi-1) * (new_t - ti-1)/(ti - ti-1)with similar expressions for any other output variables.
In a program, check that the inequality 16.5.2 istruly satisfied, and avoid asking for time points outside the available data range. The linear interpolationdoes result in small but perceptible "jaggies" in the plot, but alternatives require more sophisticatedprogramming. We do not show results here.A different difficulty is illustrated if we choose to initialize the trajectory at the same place and directionbut with a velocity of 5000 kilometers per hour. That is, the initial conditions are(15.6.4)x = 461.28; y = 0.0; u = 0.0; v = -5000Figure 16.5.2 shows the output from PLOD. Our second software choice, a problem solving environment,142Copyright © 1984, 1994 J C & M M NashNash Information Services Inc., 1975 Bel Air Drive, Ottawa, ON K2C 0X1 CanadaFigure 16.5.1SCIENTIFIC COMPUTING WITH PCsCopy for:Dr.
Dobb’s JournalPLOD screen plot for Apollo spacecraft trajectory from the starting point given asEquation 16.4.1.lets us lets us understand the situation. Using the Student Edition of MATLAB, the setup and solution ofour problem is only slightly more work than with PLOD. The MATLAB integrator (ode23) stopped at aboutT = 154 hours with the message "SUSPECTED SINGULARITY". Drawing the earth and moon on theMATLAB plot we see why in Figure 16.5.3.16.6 AssessmentIn this case study, we have not used special-purpose programs that directly carry out the trajectorycalculation.
Clearly, such an approach would be preferred if we had many calculations to perform, as longas the software were known to be reliable. We would not, for instance, want the program to compute atrajectory through the earth as PLOD did. That MATLAB found a singularity was likely accidental. Whatis more important is that the MATLAB plot of the trajectory includes the earth and moon.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
