Posted 4/2; Due 4/16
We've been talking about the N-body problem in class. It has applications in astrophysics --- simulations of planetary systems, galaxies, clusters, and so on --- but the mathematics is equally applicable to simulations of:
We'll be sticking to a relatively simple problem for now: a two-body
gravitational problem in 2-D. We'll make it even simpler by making one of the
masses heavy (say, m=1) and the other mass light (with m=0).
That way, the heavy mass doesn't move (the light mass exerts no pull on it),
plopping the heavy mass down as our center for our coordinate system makes
sense, and we only have to track the movement of the light particle. Its
See also the discussion in the introduction to Hairer, Lubich, and Wanner's "Geometric Numerical Integration: Structure-preserving algorithms for ordinary differential equations".
Calculate using pen and paper the solution to the above initial value problem. What is the orbital period? What is the eccentricity of the orbit? What energy does the system have? What is the escape velocity of the massless particle?
You'll be asked to implement Euler's method and one other. Some easier choices are:
In class we showed that Euler's method had some mimetic (feature-preserving) properties. We also showed that there were a few physical properties of the continuum problem that Euler's method did not preserve. For your chosen alternate to Euler's method, how many of the following mimetic properties can you prove and/or computationally demonstrate?
Our goal will be to verify some of the mimetics for Euler's method and your chosen alternate. Specifically, we'll verify the lack of angular momentum and energy conservation for the two-body problem with a heavy sun and a zero-weight earth.
Exercise #1: show the elliptical orbit precesses. This will demonstrate the lack of angular momentum conservation: in the continuum problem there is no torque on the system.
Exercise #2: show the orbital energy increases (show that the sum of kinetic and potential energy increases). This will demonstrate the lack of energy conservation: in the continuum problem there is no external forcing, no energy input. The energy of the system is constant.
Repeat these exercises with orbits of varying eccentricity. Try at least eccentricities of 0.1, 0.9, and 0.99. For exercise #1, how many orbital periods must pass before the orbit precesses by 30 degrees? How does this change with varying N and varying eccentricity? For exercise #2, what is the change in eccentricity per orbit? (And does this change with changing initial eccentricity? Changing N?) What is the escape time? (Is there one?)
So ... how do you pick a time step? (How do you choose N, the number of time steps per orbit?) Do you want to vary N as the eccentricity changes? Why or why not?
You have a reference solution in your analytics.
I'll also provide a
high-order adaptive method (for a sanity check on your analytics) in both C
and Java. Matlab's ODE45 now uses DOPRI5 and produces dense output.
Given Matlab's plotting routines, it might be easiest to just work there. A
call such as
[ts,ys] = ode45(odefun,0.0:h:T,y0) would work; "h"
is your step size, "T" is the orbital period (or several times it), "odefun"
is a function that returns a vector whose first two components are the
velocity and whose second two components are the acceleration, and "y0" is a
vector whose first two components are the initial position (1,0) and whose
second two components are the initial velocity (0,v0). If you don't trust
your analytics and don't want to use Matlab, ask me for a reference code in
another language. See the plots (Figure 2.2, for example) from Hairer,
Lubich, and Wanner above for ideas on visualizing the orbits.
code for visualization in Java only (though you can do all these exercises
without "looking" at the solutions). I'll try my best to get you a little
applet to play with. In the mean time, you can
and here for some ideas; the
source code for both should be available on their respective pages.
"Do you care?" (Do you care about the non-mimetic nature of Euler's method?) Well, maybe not ... depends on what questions you're asking of the system. And euler's method is simple to implement. And understand (the devil you know may be a more acceptable choice than the devil you don't).
Make plots showing the accuracy of your method(s). Using your analytic solutions (or the reference code I provide), calculate the error in the position (at a certain fixed time) for a variety of time steps. Plot the log of the norm of the error versus the log of the time step. Do the same for the velocities. Verify the expected order of convergence.
Make a plot of the time your method takes to solve a given problem. Either build a timing routine into your code, or count the number of times you have to evaluate an acceleration. Plot the log of the time taken versus the log of the time step size. This should be a straight line.
More interesting, though, is a work/accuracy plot. Plot the log of the time taken versus the log of the norm of the error. This plot tells you (for a given method) how much time it takes to achieve a given accuracy.
If you're interested in trying a many-body problem, I can suggest some ways to initialize the problem and some interesting statistics to gather.
You can also check out Chapter 31 of of the GPU Gems 3 book. They also have some sample code you can pick apart. There's a little bit of information in the Silly-pedia article on CUDA, too.