Homework #5: N-body problems: it's too crowded in here

Posted 4/2; Due 4/16

Background

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:

• individual molecules (e.g., the protein folding problem),
• many molecules (e.g., protein clusters, cell organelles, viruses, cell membrane channels),
• crack propagation at the nano- to meso-scale,
• contact angle at and moving point of contact between two fluids and a solid, and
• rare gas dynamics.

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 acceleration is:

a = -Gx/||x||³

See also the discussion in the introduction to Hairer, Lubich, and Wanner's "Geometric Numerical Integration: Structure-preserving algorithms for ordinary differential equations".

Some preliminary analytics

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:

• Heun's method, the trapezoidal rule modified to make it an explicit integrator
• the Euler-Cromer method where positions are calculated explicity and velocities implicitly (but because accelerations depend only on positions --- which are calculated explicitly --- they can be calculated without solving an equation)
• full Taylor (lowest order Nystrom method)

• leapfrog integration
• (semi-implicit) trapezoidal rule integration
• Greenspan's method, an energy conserving implicit method

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?

• translation invariance
• rotation invariance
• Galilean invariance
• center of mass doesn't accelerate
• conservation of linear momentum
• conservation of angular momentum
• conservation of energy
• time reversibility
• equal areas swept out in equal time

Some numerics

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.

• The exact period is "known". Take this and divide by an integer N to get a Δt.
• Track the mean orbital position orbit-by-orbit, and plot it. For the continuum problem, the mean orbital position is the center of the ellipse (for a closed, elliptical orbit); it is also the average position over one orbit. For your numerical approximation you can use the second idea: average all the positions the earth takes on over one "orbit". (I put "orbit" in quotes because it won't be closed: the earth won't return to its starting position. But you can average over the calculated "period" of the orbit.)
• The mean position should follow a circular orbit. The lack of angular momentum conservation (the non-zero torque) can be estimated by the turning rate.

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.

• Track energy gain/loss over a whole orbit. (Again, the orbital period and energy can be calculated exactly.)

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. I'll provide 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 look here 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).

Other problems within your reach (but that are perhaps time-consuming)

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.