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:

Calculate the net force on the particle
The net force on particle j (with mass m_j) 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 F_i,x and y component F_i,y. Calculate this value for all particles.
Calculate Acceleration on Each Particle
Again, this calculation is done independently in the x and y direction. The equations are a_i,x = F_i,x / m_i and a_i,y = F_i,x / m_i. Again, these scalar values are calculated for every particle.
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 v_x' = v_x + a_i,xh and v_y' = v_y + a_i,yh .
Update Positions
We update the position using the new velocities, namely x' = x + v_x'h and y' = y + v_y'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:

Direct Comparison of Data
I used the RK4 method with a very small timestep to create a data set that I then treated as a standard of comparison. This data set was then compared to data generated by the RK4 and Euler methods using a timestep 20 times larger. This analysis showed the expected results -- the RK4 method was indeed much more accurate. As time passed, the Euler method began to diverge dramatically from the standard data set, while the RK4 method stayed fairly consistent. The reason that I am not giving exact numbers here is that I discovered a flaw in the data recording facility of my program. There is no inaccuracy, but because of the use of threads it is impossible to control exactly when data recording begins and ends. This made it impossible to exactly align the data sets with the standard set. Thus, some of the error was do to an offset between the two data sets, not to any real numerical error. However, the general trend of increasing divergence can be considered accurate.
Comparison of Behavior on Known System
A tested several systems where the behavior could be predicted. The first system, and perhaps the most interesting, is an 8 body system under a spring-like force (imagine a spring between each pair of bodies). The system is laid out in a symmetrical fashion, and this symmetry should be maintained indefinitely in a perfectly accurate simulation. However, after a certain amount of time the system was observed to break dramatically and suddenly from this symmetry. The time it takes to reach this disruption should serve as an indicator of the accuracy of the simulation. The disruption occurred after approximately 13,000 timesteps (using the default timestep for the 8bodyspring.txt file) using Euler's method, and after about 14,000 timesteps using the RK4 method. Even more interestingly, using half the default timestep, the breakup occurred after 24,000 iterations using the RK4 method. That is, in less simulation time than with the larger timestep. This counterintuitive result suggests that the disruption has to do with the number of calculations and hence the buildup of roundoff errors, not so much with the accuracy of the numerical methods.

There are a number of special cases for which an analytic solution can be found to a 3 body problem. One of these is a collinear system. That is, given the correct initial conditions it can be shown that the 3 bodies should always orbit in a collinear fashion. The simulation file collinear.txt sets up such a simulation. The RK4 method preserves the collinear property approximately 4800 iterations (using the default timestep), and the Euler method preserves it approximately 4100 iterations. Thus, this simulation does seem to demonstrate the superior accuracy of the RK4 method.
Behavior on Analytically Solvable Systems
An effective means of analysis would be to compare data generated by the two methods to data from an analytical solution to a two-body problem. Unfortunately, I did not have time to produce a analytical 2D data set to compare with the output from my program.

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