The N-Body Problem, Numerical Methods, and Simulation Details

The n-body problem is this: for a system of n particles with known interactions, find the evolution of the system from a given initial state. The Newtonian n-body problem, for example, is the problem where the interactions are described by Newton's Universal Law of Gravitation. The problem can be reduced to finding the solutions of a system of coupled nonlinear differential equations coming from Newton's laws. The n-body problem is of importance for the understanding of various real applications and astronomical phenomena, and much effort has been made to characterize it. Unfortunately, a general solution does not appear to exist for systems of more than two bodies.

Because it is not possible to find a general solution for most interesting systems, numerical techniques are often used to model n body systems. My program, JavaPlanets, models the behavior of n-bodies in 2-dimensional space when they are acted upon by a known force. The program allows the user to give particles an initial mass, position, and velocity, and then watch the system evolve over time.

Different Approaches to the n-body Problem

The JavaPlanets program approaches the n-body problem with a "Particle-Particle" method. That is, the exact force between every pair of particles is calculated. This approach is fine for the type of simulations we are interested in. The particle-particle method basically reduces to the numerical solution of differential equations. More details of how this is done are provided in the next section. First, it is useful to consider some of the limitations of this approach and alternative techniques. The main disadvantage of particle-particle methods are that they become exceptionally expensive computationally as the number of particles grows. In a system of n particles, there are n(n-1)/2 different pairs of particles. Thus, the number of force calculations increases as the number of particles squared. That is, the method is O(n^2). For systems with large number of particles, many other approaches are possible. Each has its strengths and weakness. Some of the major methods are: Particle-Particle, Particle-Mesh, Particle-Particle/Particle-Mesh, Particle Multiple-Mesh, Nested Grid Particle-Mesh, Tree-Code (such as Barnes-Hut), Fast-Multipole-Method, Tree-Code Particle Mesh, Self-Consistent Field, and the Symplectic Method. A good general introduction to these methods can be found at http://www.amara.com/papers/nbody.html.

Many of these approaches are quite complex, but all try to make some kind of approximation to reduce the number of calculations required, while maintaining accuracy. All of these methods are directed almost entirely at modeling 1/r^2 forces. These forces are "nice" because the force between two particles is negligible when sufficient distance separates the particles.

Simulation Details

As indicated above, one way to numerical solve the n-body problem is to solve the differential equations that arise from considering all particle-particle interactions. A wide variety of techniques are available to tackle this problem. The original version of the program used perhaps the simplest -- Euler's method. This method is less accurate than more sophisticated methods but it is easy to implement correctly, and a well chosen timestep provides enough accuracy for the kinds of simulations we wanted to run. The current version of the program also implements 4th-order Runge-Kutta integration, a much more accurate technique. Another important way to gain accuracy is to use a variable timestep. However, if implemented in the traditional way, this would make it very difficult to interpret the movement of the system on the screen because the passage of simulation time would be a nonlinear function of real time. Thus, this technique was originally ruled out. A discussion of one possible way to implement a variable timestep system is described at the bottom of this page. The steps in the algorithm for calculating the movement of a particular particle using Euler's method are as follows:
  1. Calculate the net force on the particle
    The net force on particle j (with mass mj) will be the sum of the forces on the particle from every other particle in the system, and any frictional force and position dependent forces (though no position dependent forces were investigated in this paper). Let the net force on the ith particle be Fi, with x component Fi,x and y component Fi,y. Calculate this value for all particles.
  2. Calculate Acceleration on Each Particle
    Again, this calculation is done independently in the x and y direction. The equations are ai,x = Fi,x / mi and ai,y = Fi,x / mi. Again, these scalar values are calculated for every particle.
  3. Update Velocities
    Let the time step chosen be h. Then, the new velocity is the current velocity plus the change in velocity, namely acceleration multiplied by the timestep. We have vx' = vx + ai,xh and vy' = vy + ai,yh .
  4. Update Positions
    We update the position using the new velocities, namely x' = x + vx'h and y' = y + vy'h.

This routine is repeated approximately 20-30 times per second when the simulation is running in its graphical mode. However, most of this time is taken up by drawing the screen, not numerical calculations.

Error Analysis

Quantifying error in the various systems modeled and comparing the accuracies of the different integration methods proved more difficult than I anticipated, and my lack of background in numerical analysis made any rigorous investigation impossible. One primary problem was my inability to distinguish between round-off error and error inherent in the numerical integration method used. The following is a brief summary of techniques tried:

The other thing to note about this error analysis is that it does not take into account computer time. My RK4 method is much slower than my Euler implementation -- much more so than necessary, probably. Hence, even though the RK4 method is more accurate for a given timestep, it is possible that a more accurate simulation could be produced in the same amount of time by simply using Euler's method with a smaller timestep. I did not investigate this idea because I had no good way to measure the time spent on computation as opposed to painting the screen or the thread sleeping.

An Ideal Implementation: A Buffered Discrete-Event Simulation

My original plan was to incorporate the concept of a DES into my implementation, but I abandoned the idea because I did not see a good way to maintain the real-time display property that is important to the program. Dr. Conery suggested an extension of the basic DES idea that would have made this implementation feasible. The idea is to treat drawing events as events on the event queue, just like updates for the planets. At these fixed intervals, a "snapshot" of the system would be taken. This does not solve the real time issue. But, for systems with small number of bodies (the only type this simulation is real designed to handle) very little information is needed to render the screen. Thus, rather than drawing at each draw event, the current state can simply be pushed on a queue. The drawing thread can then take from this queue at regular intervals. The framerate and simulation parameters could then be adjusted to insure a steady framerate and realistic simulation, while still taking advantage of dynamically chosen time steps.




Page Maintained By:  Brendan McMahan
Email to:  mcmahahb@whitman.edu
Last Modified:  June 4th, 1999