N-Body Particle Dynamics


Dynamics is the study of particle motion. In classical, or Newtonian mechanics the dynamics of one and two particle systems is completely understood. For these simple cases we can derive equations that will give an exact description of the system at any time once we are given the positions and velocities of the particles at some reference time. The case of n-body dynamics, in which there is an arbitrary number of particles, is different. We cannot derive any analogous equations. -- The problem is not merely difficult, it is mathematically impossible. Hence we must turn to computation to solve our equations. We will study a system of three twelve identical particles moving in space. All have the same mass and feel no forces other than the gravitational attraction of the other particles. We begin by choosing an initial configuration. A configuration is simply the positions and velocities of each particle. In classical mechanics, a configuration uniquely identifies a system.

To determine the motion of the particles over time we must integrate Newton's equations. You may be familiar with:

F = ma

Where F is the force on a particle, m its mass, and a the acceleration of the particle at that instant. Note that force and acceleration, as well as velocity and position are vectors. That is, they have a direction as well as a magnitude. We represent vectors as a triple of numbers, which are the x, y and z components of the vector. Because the force is due to other, moving particles, the acceleration vector changes both magnitude and direction at every instant of time.Acceleration is the derivative of velocity with respect to time. Once we know the acceleration we can make a small step in time, to find the new velocity. The simplest step we can make is

Vnew = Vold + a * dt

where dt is a small increment of time. We will use a more sophisticated approach that attempts to compensate for the fact that even during this small time step, dt, the acceleration has changed. The method we will use is Simplectic Verlet-type Integrator with Adaptive Step Size. Similarly, velocity is the derivative of position, and we compute new positions from the velocity in the same manner as computing new velocities from accelerations. Although we cannot find equations to solve our problem, we can make some general claims about any system of particles. The most important is that energy is conserved. In our system energy is of two forms: kinetic, which is due to the motion of the particles; and potential, due to the attractive forces between particles. Energy is a scalar quantity, it has no direction, but can be positive or negative. By convention, kinetic energy (denoted by the symbol T) is a positive quantity, while potential energy (denoted by V) is negative if the force is attractive, and positive if it is repulsive. If we determine the energy of our initial configuration, we can monitor it during the course of the simulation. The individual kinetic and potential terms may change in value, but if our program is correct, their sum must always be the same. Another conserved quantity is momentum, which again comes in two forms: linear momentum is due to motion along a line; angular momentum is due to rotational motion. Momentum is a vector. We will choose (or force, if necessary!) our initial configurations to have 0 linear momentum - this is convenient when it comes time to graphically illustrating the dynamics. If the linear momentum were not 0, our particles would all eventually drift off the screen!

Below is an outline of the program flow:

  • Choose an initial configuration.
    This is a random choice within some limits that make the computation and graphics easier.
  • Take a small step in time
    We determine how large a time step to take by considering the velocities of the particles. To obtain smooth motion in the graphics, each time step should be small enough that a particle moving with average velocity moves approximately one pixel.
  • You will also write routines that compute the total energy and momenta before and after the time step.
  • In some cases, such as when particles are close together, the time step determined above is too large.
  • So you will also write a recursive routine that controls the Verlet integrator by taking smaller time steps until an accurate answer is obtained.
    The Verlet routine needs to know what the forces are on each particle, so you will write a routine that computes this. You will be given the necessary equations. Generate graphics output
  • Repeat taking time steps and updating the graphics 180 iterations. You may find that some initial configurations result in collisions of particles, after which the simulation cannot be continued.

Dave Ennis is the OSC coordinator for the Particle Dynamics project. Dave's office is in 420-4. Please contact Dave to set up appointment(s) for consultation.

For assistance, write si-contact@osc.edu or call 614-292-0890.