Duke CS Logo CompSci 100e: Program Design & Analysis II
(Fall 2009)
Home   Assignments   N-Body   How To  

N-Body Problem

Taken from N-Body Simulation Assignment by Robert Sedgewick & Kevin Wayne.

Please read this assignment and the associated How-To before you start programming.

Background

In 1687 Sir Isaac Newton formulated the principles governing the the motion of two particles under the influence of their mutual gravitational attraction in his famous Principia[1]. However, Newton was unable to solve the problem for three particles. Indeed, in general, systems of three or more particles can only be solved numerically. The N-Body problem is a reoccuring problem in many disciplines - from Biology and Biochemistry to Physics and Astronomy to Hydrology and Aerodynamics[2].

In this assignment, you will write a program to simulate the motion of N particles, mutually affected by gravitational forces, and animate the results. Such methods are widely used in cosmology, semiconductors, and fluid dynamics to study complex physical systems. Scientists also apply the same techniques to other pairwise interactions including Coulombic, Biot-Savart, and van der Waals.

The physics

We review the equations governing the motion of the particles according to Newton's laws of motion and gravitation. Don't worry if your grounding in physics is a bit rusty or even nonexistent; all of the necessary formulas are included below. We'll assume for now that the position (pxpy) and velocity (vx, vy) of each particle is known. In order to model the dynamics of the system, we must know the net force exerted on each particle.

  • Pairwise force

    Newton's law of universal gravitation asserts that the strength of the gravitational force between two particles is given by the product of their masses divided by the square of the distance between them, scaled by the gravitational constant G (6.67 × 10-11 N m2 / kg2). The pull of one particle towards another acts on the line between them. Since we are using Cartesian coordinates to represent the position of a particle, it is convenient to break up the force into its x and y components (Fx, Fy) as illustrated below.

    force diagram

  • Net force

    The principle of superposition says that the net force acting on a particle in the x or y direction is the sum of the pairwise forces acting on the particle in that direction.

  • Acceleration. Newton's second law of motion postulates that the accelerations in the x and y directions are given by: ax = Fx / m, ay = Fy / m.

The numerics.  We use the leapfrog finite difference approximation scheme to numerically integrate the above equations: this is the basis for most astrophysical simulations of gravitational systems. In the leapfrog scheme, we discretize time, and update the time variable t in increments of the time quantum Δt. We maintain the position (px, py) and velocity (vx, vy) of each particle at each time step. The steps below illustrate how to evolve the positions and velocities of the particles.

  1. For each particle: Calculate the net force (Fx, Fy) at the current time t acting on that particle using Newton's law of gravitation and the principle of superposition.

  2. For each particle:

    1. Calculate its acceleration (ax, ay) at time t using the net force computed in Step 1 and Newton's second law of motion: ax = Fx / m, ay = Fy / m.

    2. Calculate its new velocity (vx, vy) at the next time step by using the acceleration computed in Step 2a and the velocity from the old time step: Assuming the acceleration remains constant in this interval, the new velocity is (vx + Δt ax, vy + Δt ay).

    3. Calculate its new position (px, py) at time t + Δt by using the velocity computed in Step 2b and its old position at time t: Assuming the velocity remains constant in this interval, the new position is (px + Δt vx, py + Δt vy).

  3. For each particle: Draw it using the position computed in Step 2.
The simulation is more accurate when Δt is very small, but this comes at the price of more computation. Use Δt = 25,000 in this assignment for a reasonable tradeoff.

Program Specifications

Input format

The input file is a text file that contains the information for a particular universe. The first value is an integer N which represents the number of particles. The second value is a real number R which represents the radius of the universe: assume all particles will have x- and y-coordinates that remain between -R and R. Finally, there are N rows, and each row contains 6 values. The first two values are the x- and y-coordinates of the initial position; the second two values are the x- and y-coordinates of the initial velocity; the third value is the mass; the last value is a String that is the name of an image file used to display the particle. As an example, the input file planets.txt contains data for our solar system (in SI units). 

5
2.50e11
1.496e11 0.000e00 0.000e00 2.980e04 5.974e24 earth.gif
2.279e11 0.000e00 0.000e00 2.410e04 6.419e23 mars.gif
5.790e10 0.000e00 0.000e00 4.790e04 3.302e23 mercury.gif
0.000e00 0.000e00 0.000e00 0.000e00 1.989e30 sun.gif
1.082e11 0.000e00 0.000e00 3.500e04 4.869e24 venus.gif
The nbody project contains the planets.txt file, images of the planets, and many other sample universes.

Creating an animation

Draw each particle at its current position using standard draw (princeton.StdDraw), and repeat this process at each time step. By displaying this sequence of snapshots (or frames) in rapid succession, you will create the illusion of movement. After each time step:
  1. draw the background image starfield.jpg,
  2. redraw all the bodies in their new positions, and
  3. control the animation speed using StdDraw.show().

Your program

Write a program NBody.java that reads in the universe from standard input using StdIn, simulates its dynamics using the leapfrog scheme described above, and animates it using StdDraw. Maintain several arrays to store the data. To make the computer simulation, write an infinite loop that repeatedly updates the position and velocity of the particles. When plotting, use StdDraw.setXscale(-R, +R) and StdDraw.setYscale(-R, +R) to scale the physics coordinates to the screen coordinates.

Finishing touch. For a finishing touch, play the theme to 2001: A Space Odyssey using StdAudio and the file 2001.mid or my favorite superman.mid. It's a one-liner using the method StdAudio.play(). Running your program on planets.txt

N-body simulation of planets.txt

It produces this animation.

Exaple Animation

For the physicists. Are you surprised that the update formula for position uses the velocity at the updated time step rather than the previous one? Or did you expect the formula for position to be px = vx Δt + 1/2 ax (Δt)2? Here's an explanation.

The classic Euler method updates the position uses the velocity at time t instead of using the updated velocity at time t + Δt. A better idea is to use the velocity at the midpoint t + Δt / 2. The leapfrog method does this in a clever way. It maintains the position and velocity one-half time step out of phase: At the beginning of an iteration, (px, py) represents the position at time t and (vx, vy) represents the velocity at time t - Δt / 2. Interpreting the position and velocity in this way, the updated position (px + Δt vx, py + Δt vy). uses the velocity at time t + Δt / 2. Almost magically, the only special care needed to deal with the half time-steps is to initialize the system's the velocity at time t = -Δt / 2 (instead of t = 0). Note also that the acceleration is computed at time t so that when we update the velocity, we are using the acceleration at the midpoint of the interval under consideration.

The leapfrog method has several important advantages (over class numerical integration methods such as Euler or Runge-Kutta) for integrating Hamiltonian systems. (See Feynman Lectures on Physics, Vol. I, Sec. 9.6.) The leapfrog method is symplectic, which means it preserves properties specific to Hamiltonian systems (conservation of linear and angular momentum, time-reversibility, and conservation of energy of the discrete Hamiltonian). In contrast, ordinary numerical methods become dissipative and exhibit qualitatively different long-term behavior. For example, the earth would slowly spiral into (or away from) the sun. For these reasons, symplectic methods are extremely popular for N-body calculations in practice. You asked!

Extra credit. Submit a universe in our input format along with the necessary image files. If its behavior is sufficiently interesting, we'll award extra credit.

Challenges for the bored. There are limitless opportunities for additional excitement and discovery here. Try adding other features, such as:

References

[1] Newton, Isaac. Philosophiae naturalis principia mathematica. 1713.

[2] Lubomir Ivanov, "The N-Body {roblem throughout the Computer Science Curriculum,"

[3] J. Barnes, P. Hut, "A Hierarchical O(n log n) force calculation algorithm", Nature, vol.324, December 1986

Last updated Mon Aug 24 14:16:01 EST 2009