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)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)).

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

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.

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.

**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.

- 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.