2. N-Body Simulation

Goals

  • To learn about reading input using the StdIn library and printing formatted output using the StdOut library.

  • To learn about graphics and animation using the StdDraw library.

  • To learn about using the command-line to redirect standard input to read from a file.

  • To gain more experience using arrays and loops.

  • To learn how to decompose a large program into small, manageable steps - key to becoming a power programmer!

Getting Started

  • Download the project zip file for this assignment from TigerFile , which contains the files you will need for this assignment.

  • Expand the zip file, which will create the folder nbody. We recommend this folder be placed in your COS126-Assignments folder.

  • To test your program, you will need planets.txt and the accompanying image and sound files. The project folder contains these files along with a readme.txt template and additional test universes.

  • Get familiar with the command-line. Use the javac-introcs and java-introcs commands to access the standard libraries.

  • Make sure you understand the precept exercises before proceeding, including the ones using the command-line with redirection.

Background

In this assignment, you will write a program to simulate the motion of \(n\) particles (bodies) in the plane, 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.

In 1687, Isaac Newton formulated the principles governing the motion of two particles under the influence of their mutual gravitational attraction in his famous Principia. However, Newton was unable to solve the problem for three particles. Indeed, in general, solutions to systems of three or more particles must be approximated via numerical simulations.

Approach / Design

Simulating the universe: 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 physics is a bit rusty; all of the necessary formulas are included below. We’ll assume for now that the position \( (p_x,p_y) \) and velocity \( (v_x,v_y) \) 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.

$$ G = 6.67 \times 10^{-11} N \cdot m^2 \cdot kg^{-2} $$

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 \((F_x,F_y)\) as illustrated below.

Forces between bodies

  • 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. Calculating the force (between all pairs of bodies at time \(t\)) involves adding all the pairwise forces.

    For example, the forces on earth (at time \(t\)) are:

$$ \vec{F}_{earth} = \vec{F}_{mars \to earth} + \vec{F}_{mercury \to earth} + \vec{F}_{sun \to earth} + \vec{F}_{venus \to earth} $$

  • Acceleration. Newton’s second law of motion postulates that the accelerations in the \(x\)- and \(y\)-directions are given by: $$ a_x = \frac{F_x}{m} , a_y = \frac{F_y}{m} $$

Simulating the universe: 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\) (measured in seconds). We maintain the position \( (p_x,p_y) \) and velocity \( (v_x,v_y) \) of each particle at each time step. The steps below illustrate how to evolve the positions and velocities of the particles.

  • Step A: calculate the forces. For each particle, calculate the net force \((F_x,F_y) \) at the current time \(t\) acting on that particle using Newton’s law of gravitation and the principle of superposition. Note that force is a vector (i.e., it has direction). In particular, ∆x and ∆y are signed (positive or negative). In the diagram above, when you compute the force the sun (yellow-orange) exerts on the earth (dark blue), the sun is pulling the earth up (∆y positive) and to the right (∆x positive). Likewise the earth is pulling the sun down and to the left (negative ∆x and negative ∆y), albeit with a much smaller effect due to the sun’s substantially greater mass.

  • Step B: update the velocities and positions. For each particle:

    1. Calculate its acceleration \((a_x, a_y)\) at time \(t\) using the net force computed in Step A and Newton’s second law of motion: \( a_x = \frac{F_x}{m} \) and \( a_y = \frac{F_y}{m} \)
    2. Calculate its new velocity \( (v_x,v_y) \) at the next time step (time \(t\) + \(∆t\)) by using the acceleration computed in (1) and the velocity from the old time step: Assuming the acceleration remains constant in this interval, the new velocity is \( (v_x + a_x∆t, v_y + a_y∆t). \)
    3. Calculate its new position \( (p_x,p_y) \) at time \(t + ∆t\) by using the velocity computed in (2) and its old position at time \(t\): Assuming the velocity remains constant in this interval, the new position is \( (p_x+v_x∆t, p_y+v_y∆t) \).
  • Step C: draw the universe. Draw each particle, using the position computed in Step B.

Do not interleave steps A and B; otherwise, you will be computing the forces at time \(t\) using the positions of some of the particles at time \(t\) and others at time \(t + ∆t\). Also, the simulation is more accurate when \(∆t\) is very small, but this comes at the price of more computation.

Creating an animation. Draw each particle at its current position, and repeat this process at each time step until the designated stopping time. By displaying this sequence of snapshots (or frames) in rapid succession, you will create the illusion of movement.

Assignment Requirements

  • Implement a single class:
    • NBody.java
  • You MUST follow the assignment specifications for reading input from command-line arguments and standard input, and writing output to standard drawing and standard output.
  • There are two optional challenges:
    • deluxe-universe.txt
    • DeluxeNBody.java
  • Submit a readme.txt.
  • Complete the acknowledgments.txt file.
  • Read the COS 126 Style Guide to ensure your code follows our conventions. Style is an important component of writing code, and not following guidelines will result in deductions.

Program specification

Your program should:

  • Take two double command-line arguments: the duration of the simulation, \(Τ\) and the simulation time increment, \(∆t\).
  • Read in information about the universe from standard input using StdIn, using several parallel arrays to store the data.
  • Initialize the standard drawing scale and coordinate system; play the music (i.e., 2001.wav) on standard audio.
  • Simulate the universe, starting at time \(t = 0.0\), and continuing in \(Δt\) increments as long as \(t < Τ\), using the leapfrog scheme described above.
  • Animate the results using StdDraw. That is, for each time increment, draws each planet’s image.
  • After the simulation completes (reaches time \(T\)), print the state of the universe at the end of the simulation (in the same format as the input file) to standard output using StdOut.

Input format

The input format is a text file that contains the information for a particular universe.

  • The first value is an integer n that represents the number of particles.
  • The second value is a real number radius that represents the radius of the universe; it is used to determine the scaling of the drawing window.
  • Next, there are n lines, one per particle.
    • The first two values are the \(x\)- and \(y\)-coordinates of the initial position.
    • The next two values are the \(x\)- and \(y\)-components of the initial velocity.
    • The fifth value is the mass.
    • The sixth value is a String that is the name of an image file used to display the particle.
  • The remainder of the file optionally contains a description of the universe, which your program must ignore.

As an example, planets.txt contains real data from part of our Solar System.

Input file format

Compiling and executing your program

To compile your program from the command line, type the following in your embedded terminal:

javac-introcs NBody.java

To execute your program from the command line, redirecting the file planets.txt to standard input, type:

java-introcs NBody 157788000.0 25000.0 < planets.txt

After the animation stops, your program must output the final state of the universe in the same format as the input.

Possible Progress Steps

Click to show possible progress steps

The key to composing a large program is decomposing it into smaller steps that can be implemented and tested incrementally. Here is the decomposition we recommend for this assignment.

  1. Start with comments
  2. Parse command-line arguments
  3. Read universe from standard input
  4. Initialize standard drawing
  5. Play music on standard audio
  6. Simulate the universe (main time loop)
  • A Calculate net forces
  • B. Update velocities and positions
  • C. Draw universe to standard drawing
  1. Print universe to standard output

These are the order in which the components will appear in your code. Although your final code will appear in order 1–6, we recommend implementing these steps in the order 0, 1, 2, 6, 3, 4, 5, 5B, 5C, 5A. Why? Because this order will facilitate testing and debugging along the way.

Step 0: Start with comments

This outline should help you organize your program.

// Step 1. Parse command-line arguments.

// Step 2. Read the universe from standard input.

// Step 3. Initialize standard drawing.

// Step 4. Play music on standard audio.

// Step 5. Simulate the universe (main time loop).

  // Step 5A. Calculate net forces.
  // Step 5B. Update velocities and positions.
  // Step 5C. Draw universe to standard drawing.

// Step 6. Print universe to standard output.

Step 1: Parse command-line arguments

  • Parse the two command-line arguments Τ and Δt and store their values. Name the first command-line argument tau or stoppingTime since you should not begin a Java variable name with an uppercase letter.
  • Print the variables to make sure you parsed them correctly (and in the specified order). For example:
 java-introcs NBody 10 1					java-introcs NBody 157788000.0 25000.0 < planets.txt
T = 10.0							T = 1.57788E8
dt = 1.0							dt = 25000.0
  • Remove the print statements after you confirm you are parsing the command-line arguments correctly.

Step 2: Read universe from standard input

  • Read the number of bodies n as an int from standard input.
  • Read the radius of the universe radius as a double from standard input.
  • Create six (6) parallel arrays, each of length n, to store the six (6) pieces of information characterizing a body. Each array px[i], py[i], vx[i], vy[i], and mass[i] is a double arrays and usused to store the current position (\(x\)- and \(y\)-coordinates), velocity (\(x\)- and \(y\)-components), and mass of particle i, respectively; let image[i] be a String that represents the filename of the image used to display particle i.
  • Read data associated with each body and store it in parallel arrays.

Hint!
Recall the Students.java exercise from your precept.

  • Run your program with planets.txt.
java-introcs NBody 157788000.0 25000.0 < planets.txt

Step 6: Print universe to standard output

  • Print the universe in the specified format.
  • Our input files were created with the following StdOut.printf() statements.
StdOut.printf("%d\n", n);
StdOut.printf("%.2e\n", radius);
for (int i = 0; i < n; i++) {
    StdOut.printf("%11.4e %11.4e %11.4e %11.4e %11.4e %12s\n",
                  px[i], py[i], vx[i], vy[i], mass[i], image[i]);
}

Per the course collaboration policy, copying or adapting code that is not yours is permitted only if it comes from the course materials. If you do so, you must cite any code that you copy or adapt (with the exception of code that is included with the assignment). So, you are free to copy or adapt this code fragment, with or without citation.

  • Run your program with planets.txt. You should see exactly the output below.
java-introcs NBody 0.0 25000.0 < planets.txt
5
2.50e+11
 1.4960e+11  0.0000e+00  0.0000e+00  2.9800e+04  5.9740e+24    earth.gif
 2.2790e+11  0.0000e+00  0.0000e+00  2.4100e+04  6.4190e+23     mars.gif
 5.7900e+10  0.0000e+00  0.0000e+00  4.7900e+04  3.3020e+23  mercury.gif
 0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.9890e+30      sun.gif
 1.0820e+11  0.0000e+00  0.0000e+00  3.5000e+04  4.8690e+24    venus.gif
  • Test your program on another input file.
java-introcs NBody 0.0 1.0 < 3body-zero-gravity.txt
3
5.12e+02
 0.0000e+00  0.0000e+00  1.0000e+00  1.0000e+00  1.0000e-30    earth.gif
 1.2800e+02  0.0000e+00  2.0000e+00  1.0000e+00  1.0000e-40    venus.gif
 0.0000e+00  1.2800e+02  1.0000e+00  2.0000e+00  1.0000e-50     mars.gif
  • Submit your code to TigerFile and click the Check Submitted Files button. If everything is correct so far (reading and printing the universe), your code should pass Tests 1–7. Do not continue to the next step until it does.

Step 3: Initialize standard drawing

Initialize standard drawing

  • The default \(x\)- and \(y\)-scale supports coordinates between 0 and 1. Set the scale of the standard drawing window so that (0, 0) is at the center, (−radius, −radius) is at the lower-left corner, and (+radius, +radius) is at the upper-right corner.

    Hint!
    Recall StdDraw.setXscale() and StdDraw.setYscale()

  • Call StdDraw.enableDoubleBuffering() to enable double buffering and support animation.

Step 4: Play music on standard audio

  • Call StdAudio.playInBackground("2001.wav") to play the background music.

Step 5: Simulate the universe (the big time loop)

  • Create the time simulation loop. Each iteration will perform a single time step of the simulation, starting at \(t = 0.0\), and continuing in \(Δt\) increments as long as \(t < Τ\). This code goes between the code where you read the input and the code where you print the output.
  • Test your program by printing the time \(t\) at each time step.
java-introcs NBody 23.0 2.5 < planets.txt		java-introcs NBody 25.0 2.5 < planets.txt
t = 0.0							t = 0.0
t = 2.5							t = 2.5
t = 5.0							t = 5.0
t = 7.5							t = 7.5
t = 10.0						t = 10.0
t = 12.5						t = 12.5
t = 15.0						t = 15.0
t = 17.5						t = 17.5
t = 20.0						t = 20.0
t = 22.5						t = 22.5
  • Once you are confident that your time loop is correct, remove the code that prints the time \(t\) at each time step.

Steps 5B, 5C, and 5A (below) will go inside this time loop.

FAQ
In which order should I implement these 3 sub-steps?
5B, 5C, 5A because calculating forces is hardest.

FAQ
Can I interleave steps 5A, 5B, and 5C?
No. Not only is it bad design, but it ruins the physics: you need the position of all bodies at time \(t\) to calculate the forces between each pair of bodies, and updating some positions to time \(t + ∆t\) would mess up the forces for time \(t\).

Step 5B: Update the velocities and positions

  • Since we have not implemented Step 5A yet, the forces \(F_x\) and \(F_y\) are 0, so we can initially assume that the accelerations \(a_x\) and \(a_y\) are zero and the velocities are constant.
  • Write a loop to calculate the new velocity and position for each particle.
    • Update the velocity of each body: \(v_x = v_x + a_x∆t\) and \(v_y = v_y + a_y∆t\) (note for now \(a_x\) and \(a_y\) are 0).
    • Update the position of each body: \(p_x = p_x + v_x∆t\) and \(p_y = p_y + v_y∆t\).
  • Test on an artificial universe, where it is easy to check the expected results by hand.
java-introcs NBody 1 1 < 3body-zero-gravity.txt
3
5.12e+02
 1.0000e+00  1.0000e+00  1.0000e+00  1.0000e+00  1.0000e-30    earth.gif
 1.3000e+02  1.0000e+00  2.0000e+00  1.0000e+00  1.0000e-40    venus.gif
 1.0000e+00  1.3000e+02  1.0000e+00  2.0000e+00  1.0000e-50     mars.gif
java-introcs NBody 192 1 < 3body-zero-gravity.txt
3
5.12e+02
 1.9200e+02  1.9200e+02  1.0000e+00  1.0000e+00  1.0000e-30    earth.gif
 5.1200e+02  1.9200e+02  2.0000e+00  1.0000e+00  1.0000e-40    venus.gif
 1.9200e+02  5.1200e+02  1.0000e+00  2.0000e+00  1.0000e-50     mars.gif
  • Submit your code via TigerFile and click the Check Submitted Files button. If everything is correct so far (ignoring the gravitational forces and drawing), your code should pass Tests 1–10. Do not continue to the next step until it does.

Step 5C: Draw universe to standard drawing

  • Draw each particle at its current position to standard drawing, and repeat this process at each time step until the designated stopping time. By displaying this sequence of snapshots (or frames) in rapid succession, you will create the illusion of movement.
  • After each time step, draw the background image starfield.jpg; redraw all the particles in their new positions; and control the animation speed.
    • Use StdDraw.picture(x, y, filename) to draw the background image starfield.jpg, centered at (0, 0).
    • Then, write a loop to display the n particles.
    • Call StdDraw.show() to display results on screen.
    • Call StdDraw.pause(20) to control animation speed, which produces about 50 frames per second.
  • Run the program with planets.txt. You should now see the four planets moving off the screen in a straight line, with constant velocity.
java-introcs NBody 157788000.0 25000.0 < planets.txt
  • If the planets are flickering, be sure you are using double buffering, following the template in BouncingBallsDeluxe.java. In particular, call StdDraw.enableDoubleBuffering() once at the beginning of the program, and call StdDraw.show() and StdDraw.pause(20) at the end of each time step (frame), not after each call to StdDraw.picture().
  • Run your program on another input file, such as kaleidoscope.txt, which produces this animation:
java-introcs NBody 157788000.0 25000.0 < kaleidoscope.txt
  • Submit your code via TigerFile and click the Check Submitted Files button. If everything is correct so far (ignoring the gravitational forces), your code should pass Tests 1–17. Do not continue to the next step until it does.

Step 5A: Calculate the forces

  • This step is the most difficult one to implement. Let’s review:
    • Pairwise force. Recall this diagram from above:

Forces between bodies

  • In Java, you can represent the gravitational constant using scientific notation: 6.67e-11 (or 6.67E-11)

  • 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. Calculating the force (between all pairs of bodies at time \(t\)) involves adding all the pairwise forces. For example, the forces on earth (at time \(t\)) are: $$ \vec{F}_{earth} = \vec{F}_{mars \to earth} + \vec{F}_{mercury \to earth} + \vec{F}_{sun \to earth} + \vec{F}_{venus \to earth} $$

  • To calculate the net force acting on each body, you will need two additional arrays fx[i] and fy[i] to store the net force acting on body i.

  • The net force (arrays fx[i] and fy[i]) must be initialized to 0 at each time step.

  • Use a nested loop pattern to calculate the net force exerted by body j on body i. Add these values to fx[i] and fy[i], but skip the case when i equals j (you don’t want to calculate the force of body i on itself).

    Hint!
    Before calculating the net forces, can you implement a nested loop pattern that enumerates all pairs of bodies? For example: what pairs are not printed?

n == 4                  n == 5
0-1 0-2 0-3             0-1 0-2 0-3 0-4
1-0 1-2 1-3             1-0 1-2 1-3 1-4
2-0 2-1 2-3             2-0 2-1 2-3 2-4
3-0 3-1 3-2             3-0 3-1 3-2 3-4
                        4-0 4-1 4-2 4-3
  • Once you have these values computed, go back to Step 5B and use them to compute the acceleration (instead of assuming it is zero).
    • Calculate the acceleration \( (a_x, a_y) \) at time \(t\) using the net force computed in Step 5A and Newton’s second law of motion: \( a_x = \frac{F_x}{m} \) and \( a_y = \frac{F_y}{m} \).
    • The code from Step 5B should then compute the new velocities and positions.
  • Test your program on several data files.
  • Submit your code via TigerFile and click the Check Submitted Files button. Your code should now pass all of the tests!

Debugging

Checkstyle complains about a variable because it begins with an uppercase letter. Which name should I use instead? By convention, a variable in Java should begin with a lowercase letter and use camelCase, such as isLeapYear. A constant variable (a variable whose value does not change during the execution of a program or from one execution of the program to the next) should begin with an uppercase letter and use underscores to separate any word boundaries, such as GRAVITATIONAL_CONSTANT. So, a Java variable name might not always perfectly align with the corresponding domain-specific mathematical notation. In this case, tau or stoppingTime would be fine choices.
Why do I get an exception when I read data? InputMismatchException means that your program is attempting to read data of one type but the next piece of data on standard input is of an incompatible type. For example, if you use StdIn.readInt() and the next token on standard input is 3.14, you’ll get this error. NoSuchElementException means that your program is trying to read data from standard input when there is no more data available from standard input. For example, if you use StdIn.readInt() to read two integers but standard input contains only a single integer. (This can also happen if you wipe out an input file, such as planets.txt, by errantly typing > instead of < during a redirection command.)
My animation looks fine, but the numbers are off a little bit when I follow the tests below. What could be causing this? Check that you followed the assignment instructions exactly. In particular, you should have separate loops for Steps A, B, and C and you should update the velocities before the positions in Step B. Check that your gravitational constant is exactly 6.67e-11 as given.
I draw the bodies, but they do not appear on the screen. Why? Did you set the scale of standard drawing (StdDraw.setXscale() and StdDraw.setYscale()) to change the \(x\)- and \(y\)-coordinate systems to use the physics coordinates instead of the unit box? Since you want it centered on the origin with a square radius of radius, the minimum of each axis should be –radius and the maximum should be +radius.
My bodies repel each other. Why don’t they attract each other?

Make sure that you get the sign right when you apply Newton’s law of gravitation: (x[j] - x[i]) vs. (x[i] - x[j]). Note that \(Δx\) and \(Δy\) can be positive or negative, so do not use Math.abs(). Do not consider changing the universal gravitational constant \(G\) to patch your code!

Why do my bodies disappear?

Make sure you are using the correct index variables. It’s very easy to confuse an i with a j, especially when copying and pasting Java statements.

Testing Output

Here are our results for a few sample inputs.

java-introcs NBody 0.0 25000.0 < planets.txt           ## zero steps
5
2.50e+11
 1.4960e+11  0.0000e+00  0.0000e+00  2.9800e+04  5.9740e+24    earth.gif
 2.2790e+11  0.0000e+00  0.0000e+00  2.4100e+04  6.4190e+23     mars.gif
 5.7900e+10  0.0000e+00  0.0000e+00  4.7900e+04  3.3020e+23  mercury.gif
 0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  1.9890e+30      sun.gif
 1.0820e+11  0.0000e+00  0.0000e+00  3.5000e+04  4.8690e+24    venus.gif
java-introcs NBody 25000.0 25000.0 < planets.txt       ## one step
5
2.50e+11
 1.4960e+11  7.4500e+08 -1.4820e+02  2.9800e+04  5.9740e+24    earth.gif
 2.2790e+11  6.0250e+08 -6.3860e+01  2.4100e+04  6.4190e+23     mars.gif
 5.7875e+10  1.1975e+09 -9.8933e+02  4.7900e+04  3.3020e+23  mercury.gif
 3.3087e+01  0.0000e+00  1.3235e-03  0.0000e+00  1.9890e+30      sun.gif
 1.0819e+11  8.7500e+08 -2.8329e+02  3.5000e+04  4.8690e+24    venus.gif
java-introcs NBody 50000.0 25000.0 < planets.txt       ## two steps
5
2.50e+11
 1.4959e+11  1.4900e+09 -2.9640e+02  2.9799e+04  5.9740e+24    earth.gif
 2.2790e+11  1.2050e+09 -1.2772e+02  2.4100e+04  6.4190e+23     mars.gif
 5.7826e+10  2.3945e+09 -1.9789e+03  4.7880e+04  3.3020e+23  mercury.gif
 9.9262e+01  2.8198e-01  2.6470e-03  1.1279e-05  1.9890e+30      sun.gif
 1.0818e+11  1.7499e+09 -5.6660e+02  3.4998e+04  4.8690e+24    venus.gif
java-introcs NBody 60000.0 25000.0 < planets.txt       ## three steps
5
2.50e+11
 1.4958e+11  2.2349e+09 -4.4460e+02  2.9798e+04  5.9740e+24    earth.gif
 2.2789e+11  1.8075e+09 -1.9158e+02  2.4099e+04  6.4190e+23     mars.gif
 5.7752e+10  3.5905e+09 -2.9682e+03  4.7839e+04  3.3020e+23  mercury.gif
 1.9852e+02  1.1280e+00  3.9705e-03  3.3841e-05  1.9890e+30      sun.gif
 1.0816e+11  2.6248e+09 -8.4989e+02  3.4993e+04  4.8690e+24    venus.gif
java-introcs NBody 31557600.0 25000.0 < planets.txt    ## one year
5
2.50e+11
 1.4959e+11 -1.6531e+09  3.2949e+02  2.9798e+04  5.9740e+24    earth.gif
-2.2153e+11 -4.9263e+10  5.1805e+03 -2.3640e+04  6.4190e+23     mars.gif
 3.4771e+10  4.5752e+10 -3.8269e+04  2.9415e+04  3.3020e+23  mercury.gif
 5.9426e+05  6.2357e+06 -5.8569e-02  1.6285e-01  1.9890e+30      sun.gif
-7.3731e+10 -7.9391e+10  2.5433e+04 -2.3973e+04  4.8690e+24    venus.gif
# this test should take only a few seconds
# 4.294E9 is bigger than the largest int
java-introcs NBody 4.294E9 2.147E9 < 3body.txt
3
1.25e+11
 2.1470e+12 -7.8082e-03  5.0000e+02 -3.6368e-12  5.9740e+24    earth.gif
 1.2882e+14 -1.5100e+17  3.0000e+04 -3.5165e+07  1.9890e+30      sun.gif
-1.2882e+14  1.5100e+17 -3.0000e+04  3.5165e+07  1.9890e+30      sun.gif

Optional Challenges

We offer two challenges, for fame and glory only (not extra credit).

  • Challenge - deluxe-universe.txt

    There are limitless opportunities for additional excitement and discovery in this assignment. Create an alternate universe (using the same input format). Submit the deluxe-universe.txt file (using the same input format) as well as any additional image files. We provide some alternative universes written by former students (e.g., planetsparty.txt, twinbinaries.txt, chaosblossum.txt, galaxy.txt, among others). To try these universes, use your NBody.java with a corresponding input file named deluxe-universe.txt. For example:

    java-introcs NBody 157788000.0 25000.0 < deluxe-universe.txt
    

    Since this is a challenge, you may not ask for design/debugging help from COS 126 staff or the lab TAs. Submit deluxe-universe.txt (along with any image files) and briefly explain your approach in the readme.txt file.

  • Challenge - DeluxeNBody.java

    Design and implement a DeluxeNBody.java that adds other features, such as supporting elastic or inelastic collisions. Or, make the simulation three-dimensional by doing calculations for \(x\)-, \(y\)-, and \(z\)-coordinates, then using the \(z\)-coordinate to vary the sizes of the planets. Add a rocket ship that launches from one planet and has to land on another. Allow the rocket ship to exert force with the consumption of fuel.

    Since this is a challenge, you may not ask for design/debugging help from COS 126 staff or the lab TAs. Submit DeluxeNBody.java and any input file(s); include instructions in the readme.txt file. Also, record a two-minute video using loom.com to showcase your universe; include the link to the video along with an explanation of your approach in the readme.txt file.

Submission

Submit NBody.java, readme.txt and acknowledgments.txt files to TigerFile . If you attempted the challenges, submit any files that make up your alternative universe (deluxe-universe.txt, image files) and/or DeluxeNBody.java (along with any additional data files). Please document your work in the readme.txt.

Enrichment

Did Prof. Andrew Appel do his senior thesis based on the n-body problem? Yes, Appel graduated summa cum laude with an A.B. in physics from Princeton University in 1981 after completing a senior thesis, titled Investigation of galaxy clustering using an asymptotically fast N-body algorithm, under the supervision of Nobel laureate James Peebles (Wikipedia).
What is the music in 2001.wav? It’s the fanfare to Also sprach Zarathustra by Richard Strauss. It was popularized as the key musical motif in Stanley Kubrick’s 1968 film 2001: A Space Odyssey.
I’m a physicist. Why should I use the leapfrog method instead of the formula I derived in high school? In other words, why does the position update formula use the velocity at the updated time step rather than the previous one? Why not use the \( \frac{1}{2}at^2 \) formula?

The leapfrog method is more stable for integrating Hamiltonian systems than conventional numerical methods like Euler’s method or Runge-Kutta. 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!

Here’s a more complete explanation of how you should interpret the variables. The classic Euler method updates the position using 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, \( (p_x, p_y) \) represents the position at time \(t\) and \( (v_x, v_y) \) represents the velocity at time \(t − Δt / 2\). Interpreting the position and velocity in this way, the updated position \( (p_x + v_xΔt, p_y+v_yΔt) \) 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 velocity at time \(t = −Δt / 2 \) (instead of \(t = 0.0\) ), and you can assume that we have already done this for you. 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.

Are there analytic solutions known for n-body systems? Yes, Sundman and Wang developed global analytic solutions using convergent power series. However, the series converge so slowly that they are not useful in practice. See The Solution of the N-body Problem for a brief history of the problem.
Who uses n-body simulation?

N-body simulations play a crucial role in our understanding of the universe. Astrophysicists use it to study stellar dynamics at the galactic center, stellar dynamics in a globular cluster, colliding galaxies, and the formation of the structure of the Universe. The strongest evidence we have for the belief that there is a black hole in the center of the Milky Way comes from very accurate n-body simulations. Many of the problems that astrophysicists want to solve have millions or billions of particles. More sophisticated computational techniques are needed.

The same methods are also widely used in molecular dynamics, except that the heavenly bodies are replaced by atoms, gravity is replaced by some other force, and the leapfrog method is called Verlet’s method. With van der Waals forces, the interaction energy may decay as \( 1/R^6 \) instead of an inverse square law. Occasionally, 3-way interactions must be taken into account, e.g., to account for the covalent bonds in diamond crystals.

Kepler’s Orrery uses an n-body simulator to compose and play music.

What techniques are used to do n-body simulation in practice? Here’s a wealth of information in n-body simulation.

Copyright © 1999–2022, Robert Sedgewick and Kevin Wayne.