15-859B: Assignment 3: Boundary Value Problem: Put a Spaceship in Orbit

Due: Thu, 16 Nov, 5pm at my office (Newell-Simon 4205), or in class earlier that day.

Summary

Simulate 2-D gravitational attraction on three bodies: the earth, the moon, and a spaceship. Solve a boundary value problem to determine the launch velocity of a spaceship to put it into orbit around the moon.

Detail

I recommend you do the simulations in 2-D, since it's much easier to visualize, and use the following formulas for pairwise potential and force between particles i and j:

potential energy = phi = G*mi*mj*log(rij)
force = -gradient(phi) = -G*mi*mj*(Xi-Xj)/rij^2

where mi and Xi are the mass and 2-D position of particle i, rij=||Xi-Xj|| is the distance between particles i and j, and G is the gravitational constant. Note that the use of the (non-unit) vector Xi-Xj in the numerator makes the force law a 1/r law. This formula is discussed in the paper: Greengard-Rokhlin, J. Comput. Phys. 73, 325-348, 1987, if you're curious. Note that this is not Newton's gravitation law, for which the potential goes as 1/r and the force goes as 1/r^2. This 2-D law tends to give more interesting behavior for 2-D simulations, and it obeys Laplace's equation in 2-D. Whatever potential you choose, make sure that your potential and force are consistent (force = -gradient(phi)).

The differential equation says that the mass of particle i (mi) times its acceleration (Xi'') equals the sum of the forces on it, namely the sum of the gravitational forces given by the formula above, due to each of the other particles. Each force is a 2-D vector.

I don't care particularly about units. We can use G=1 and fictitious masses and distances. (Feel free to use more physically accurate quantities if you like, but preserve the qualitative nature of the problem.)

What we have here is a system of ODE's: for n particles it's a system of 2n 2nd order ODE's, or a system of 4n 1st order ODE's. Implement a boundary value solution method -- you can choose from one of the methods discussed in Heath's book. I believe the shooting method will be easiest, since this is a nonlinear differential equation.

Part 1

First, put the moon in approximately circular orbit around the earth in a 2-particle simulation. At the beginning of time (t=0), start with the Earth, a mass 1000 particle, at rest, at the origin. Start the moon, a mass 10 particle, at (1,0) at t=0. The other boundary condition is that the moon should pass through (0,1) at t=tq=pi*r/(2*sqrt(G*me)) where r is the radius of the desired orbit (I'm using r=1), and me is the mass of the earth. Solve for the initial velocity (a 2-vector) of the moon, while simulating the mutual attraction of earth and moon due to gravity, to satisfy this boundary value problem. Although 2-body gravitational motion is solvable analytically, use a numerical method to solve this; this is warmup for the more difficult part, to come.

Don't proceed until you have the above, and you're convinced (by graphing trajectories) that your moon is in a stable, approximately circular orbit of the chosen radius around the earth, and that the earth is moving only slowly. Note that the coordinate system I've recommended is not fixed relative to the center of mass of the two bodies, so after a long time, both will drift away from the origin.

Part 2

Now launch a spaceship of total mass .1, with two astronauts aboard, from a point in space (.8,0) at t=0, and solve for the initial velocity that will put it in orbit around the moon at time t=tq. Do a 3-body simulation; don't neglect the force of the spaceship on the earth and moon. Note: the spaceship motion might force you to use a small step size. Two of the possible strategies for this part are:

a) Turn the moon motion into an IVP using the initial velocity from part 1. Then the only BVP is the spaceship. It may require some trial and error to determine the appropriate position of the spaceship at t=tq to put it into orbit. That you can do manually. But solving for the initial velocity of the spaceship, once you've fixed the exact target position of the spaceship at t=tq, should be automated. Note that the moon probably won't pass through the same point at t=tq, because of the moon's gravitational attraction to the spaceship.

b) Keep the moon motion a BVP and add the spaceship BVP, so you'd be solving for both the (new, spaceship-perturbed) moon initial velocity and the spaceship initial velocity.

If you're unable to put your spaceship into orbit around the moon, you can crash it into the moon at t=tq, but there will be a grading penalty, and it will be your job to notify the families of the deceased :-) .

The assignment can be implemented in Matlab or any other language.

Turn in:

Discussion of what you did, what worked and didn't work, of 2 pages or less.
Graphs from the first part, of the paths of the moon and earth for approximately one orbit. These graphs should be in the xy plane, not tx and ty. Visualization tip: I suggest you use a small step size h to do the simulation and draw the curve, and separately put x's or o's on the curve with a time increment of h2, so it is clear how fast the two are moving. Use h2>h, so that there are more than 10 but fewer than 100 markers around the orbit.
Graphs from the second part, of the paths of the earth, moon, and spaceship, also all on one plot, in the xy plane. Use of different colors, line patterns, or widths might help to avoid confusion between curves here. Tick marks on the curves at equal time increments h2 is critical. It should be possible to figure out by looking carefully at your plots, and counting tick marks (better yet, label them with small numbers!) the positions of each particle at a given time. You can separately graph x(t) and y(t) of the spaceship, if you want.
Additional graphs if you feel they're necessary to show your solution approach (e.g. graphs of some of the trajectories tried by your shooting algorithm, should you attack the problem that way).
Code listing of your BVP solution method.

15-859B, Introduction to Scientific Computing
Paul Heckbert
Revision 1: 3 Nov. 2000.
Revision 2: 3 Nov. 2000. Clarified what the diff.eq. is.
Revision 3: 5 Nov. 2000. Clarified how to build part 2 from results of part 1.