Project 2: N-body Dynamics

COMP 536 | Short Projects

Author

Dr. Anna Rosen

Published

April 22, 2026

Due: Tuesday March 3, 2026 by 11:59 pm

Lab Check-ins:

Check-in 1: Wednesday February 18 (in lab)
Check-in 2: Wednesday February 25 (in lab)

Collaboration Guidelines

Individual Implementation Requirements:

Core algorithms should be independently written
Shared: Overall structure, test cases, debugging strategies
Independent: Actual implementation of integrators, force calculation
OK to share: Plotting code, file I/O, validation tests

Review the Project Submission Guide for collaboration best practices.

N-body Dynamics with Plummer Initial Conditions

Write a Python software package to simulate the dynamics of a star cluster using N-body algorithms. You will implement and compare three ODE integration schemes: Euler, RK4, and Leapfrog.

Your code must be modular and object-oriented, allowing users to integrate the orbits of \(N\) bodies under their mutual gravitational interactions, assuming the system is isolated. Your code should allow users to select any integration method (e.g., via a class ODEIntegrator or a method parameter).

One simulator, multiple verification tests

You’re building one N-body simulator that works for any \(N\). The 2-body and 3-body tests are verification steps — you run the same code with \(N=2\) or \(N=3\) to validate against known solutions before trusting it for \(N=100\).

Relevant Course Materials:

Numerical Problem Description

Gravitational Acceleration

The acceleration of particle \(i\) in an \(N\)-body gravitational system, including a softening parameter \(\epsilon\), is given by:

\[ \vec{a}_i = \sum_{\substack{j=1 \\ j \ne i}}^{N} G\,m_j \frac{\vec{r}_j - \vec{r}_i}{\left(\lvert\vec{r}_j - \vec{r}_i\rvert^2 + \epsilon^2\right)^{3/2}} \]

where \(m_i\) and \(m_j\) are the particle masses, and the distance \(r_{ij}\) between particles \(i\) and \(j\) is defined as:

\[ r_{ij}^2 = (x_j - x_i)^2 + (y_j - y_i)^2 + (z_j - z_i)^2 \]

Softening Parameter

The softening parameter \(\epsilon\) prevents numerical divergences when two particles approach closely, effectively replacing the point-mass potential at small scales with a smoothed potential core. Use:

\[\epsilon \approx 0.01 \times R_{\text{cluster}}/N^{1/3},\]

where \(R_{\text{cluster}}/N^{1/3}\) roughly corresponds to the mean inter-particle spacing.

Choosing \(\epsilon\) for your simulations:

For testing (\(a = 100\) AU, \(N \leq 50\)): \(\epsilon \approx 1\) AU
For production (\(a = 1000\) AU, \(N = 100\)–\(200\)): \(\epsilon \approx 5\) AU
Too small: numerical errors when particles approach closely
Too large: artificially smooths dynamics

Worked example: For \(N = 50\), \(R = 100\) AU: mean spacing \(\approx 100 / 50^{1/3} \approx 27\) AU, so \(\epsilon \approx 0.01 \times 27 \approx 0.27\) AU. We round up to \(\epsilon = 1\) AU for safety.

Required Units

\[ \begin{aligned} \text{Mass}: &\quad M_{\odot} \quad | \quad \text{Length}: \quad \text{AU} \quad | \quad \text{Time}: &\quad \text{yr} \end{aligned} \]

In these units, the gravitational constant becomes: \[ G = 4\pi^2~\text{AU}^3 \text{ M}_{\odot}^{-1} \text{ yr}^{-2} \]

Energy and Virial Diagnostics

You must monitor and plot the evolution of:

Total kinetic energy \(E_K\)
Total gravitational potential energy \(W\) (count each pair only once!)
Total energy \(E_{\rm tot} = E_K + W\)
Virial ratio \(Q\) (should be approximately 1 at equilibrium)
Virial residual \(\Delta\) (should be less than 0.01 near equilibrium)

Use:

\[ E_K = \frac{1}{2}\sum_{i=1}^{N} m_i v_i^2, \qquad W = -\sum_{i<j} \frac{G m_i m_j}{r_{ij}}, \qquad E_{\rm tot} = E_K + W \]

\[ Q = \frac{2E_K}{|W|}, \qquad \Delta = \frac{|2E_K + W|}{|W|} = |Q - 1| \]

Monitoring all components separately helps debug which type of energy causes conservation issues.

Note: These quantities represent the system values (i.e., summed over all particles). When computing the potential energy, count each pair only once (i.e., use \(i<j\) in the sum).

ODE Integration Methods Reference

Full derivations are in the course materials (linked above). For N-body dynamics, you’ll implement three integration schemes for the system \(\dot{\vec{y}} = \vec{f}(t, \vec{y})\) where \(\vec{y} = [\vec{r}, \vec{v}]\) and \(\vec{f}\) returns \([\vec{v}, \vec{a}]\):

Euler’s Method (1st order, simplest): \[\vec{y}_{n+1} = \vec{y}_n + \Delta t \cdot \vec{f}(t_n, \vec{y}_n)\]

RK4 (4th order, high accuracy): \[\begin{aligned} \vec{k}_1 &= \vec{f}(t_n, \vec{y}_n) \\ \vec{k}_2 &= \vec{f}\left(t_n + \frac{\Delta t}{2}, \vec{y}_n + \frac{\Delta t}{2} \vec{k}_1\right) \\ \vec{k}_3 &= \vec{f}\left(t_n + \frac{\Delta t}{2}, \vec{y}_n + \frac{\Delta t}{2} \vec{k}_2\right) \\ \vec{k}_4 &= \vec{f}(t_n + \Delta t, \vec{y}_n + \Delta t \cdot \vec{k}_3) \\ \vec{y}_{n+1} &= \vec{y}_n + \frac{\Delta t}{6}\left(\vec{k}_1 + 2\vec{k}_2 + 2\vec{k}_3 + \vec{k}_4\right) \end{aligned}\]

Leapfrog (2nd order, symplectic — preserves phase-space structure): \[\begin{aligned} \vec{v}_{n+1/2} &= \vec{v}_n + \frac{\Delta t}{2} \cdot \vec{a}(\vec{r}_n) \quad \text{(half-step velocity)} \\ \vec{r}_{n+1} &= \vec{r}_n + \Delta t \cdot \vec{v}_{n+1/2} \quad \text{(full-step position)} \\ \vec{v}_{n+1} &= \vec{v}_{n+1/2} + \frac{\Delta t}{2} \cdot \vec{a}(\vec{r}_{n+1}) \quad \text{(half-step velocity)} \end{aligned}\]

Critical Implementation Notes:

For N-body: \(\vec{f}\) computes accelerations from positions
Leapfrog updates positions and velocities separately (different structure from Euler/RK4!)
Store intermediate \(\vec{k}\) values efficiently — they’re expensive to compute
Test with \(N=2\) first where an analytical solution exists

How Leapfrog conserves energy

Leapfrog is symplectic — it conserves a shadow Hamiltonian \(\tilde{H} = H + O(\Delta t^2)\) that is close to the true Hamiltonian. This means:

The true energy oscillates with bounded amplitude \(\propto \Delta t^2\)
There is no secular (long-term) drift — the error does not grow over time
This is fundamentally different from RK4, where the energy error is smaller at any given instant but drifts monotonically

The advantage of Leapfrog is long-term stability, not instantaneous accuracy. This is why it is the standard choice for long N-body simulations.

Project Phases

The Big Picture

You’re building one N-body simulator and running it with different values of \(N\):

\(N=2\) (Sun-Earth): Verification test against a known analytical solution. Any errors here are purely from the integration method, not force calculation.
\(N=3\) (Sun-Earth-Jupiter): Verification test for multi-body force summation. You can visually verify Jupiter’s 12-year period.
\(N=10\) to \(N=1000+\): The full star cluster simulation where both integration accuracy and computational efficiency matter.

The pipeline-first approach

Phase 1 builds the complete pipeline — force calculation, integration, energy diagnostics, plotting — using the simplest integrator (Euler). Once your pipeline works end-to-end, Phase 2 adds RK4 and Leapfrog as drop-in replacements. Phase 3 adds realistic initial conditions. Phase 4 optimizes for performance.

Each phase adds one new idea on top of a validated foundation.

Phase 1: Build the Pipeline (Euler + Validation)

Goal: A complete, working N-body simulator with Euler integration, validated against known physics.

Implementation:

Force calculation (general \(N\), with softening)
Euler integrator
Energy and virial diagnostics
Plotting infrastructure

Validation — N=2 (Sun-Earth):

Your N-body code with \(N=2\) particles
Initial conditions: Sun at origin (\(m = 1\,M_\odot\)), Earth at \(r = 1\) AU with circular velocity \(v_c = \sqrt{GM_\odot/r}\)
Simulate for 10 orbital periods (10 years)
Use \(\Delta t = 0.01\) years
Observe: Euler produces a spiral trajectory and systematic energy drift — this is expected!

Validation — N=3 (Sun-Earth-Jupiter):

Same code with \(N=3\) particles
Sun at origin (\(m = 1\,M_{\odot}\))
Earth at 1 AU (\(m = 3 \times 10^{-6}\,M_{\odot}\)), circular velocity
Jupiter at 5.2 AU (\(m = 9.55 \times 10^{-4}\,M_{\odot}\)), circular velocity
Simulate for 24 years (captures 2 Jupiter orbits)
Check: No particles escape or crash. Earth completes \(\sim 24\) orbits, Jupiter completes \(\sim 2\).

After this phase: You have a working, validated N-body simulator. The force calculation, energy diagnostics, and output infrastructure are done. Everything that follows builds on this foundation.

Phase 2: Upgrade the Integrator

Goal: Implement RK4 and Leapfrog as drop-in replacements. Compare all three methods to understand their tradeoffs.

Just two new functions

Adding RK4 and Leapfrog means implementing two new integration functions with the same interface as your Euler step. Everything else — force calculation, energy diagnostics, plotting — stays the same. This is the payoff of modular software design.

Implementation:

Implement RK4 step function (same interface as Euler)
Implement Leapfrog step function (note: different structure — position/velocity split)
Re-run \(N=2\) (Sun-Earth) with all three methods using the same \(\Delta t = 0.01\) yr
Plot energy error comparison

What you should observe:

Euler: Systematic energy drift; orbit spirals outward
RK4: Extremely small energy error initially, but it drifts monotonically over time (grows slowly but never stops)
Leapfrog: Energy error oscillates with bounded amplitude — no long-term drift, even after many orbits

After this phase: You understand why Leapfrog is the standard choice for long-term gravitational simulations: not because it is the most accurate at any instant, but because its error is bounded rather than accumulating. Choose Leapfrog for your cluster simulations in Phase 3.

Phase 3: Star Cluster with Plummer Initial Conditions

Goal: Implement realistic star cluster initial conditions using a Plummer sphere and simulate cluster dynamics.

Phase 3A — Simple sanity check (\(N=10\), loops):

Uniform masses (\(1\,M_\odot\) each)
Random positions in a sphere of radius \(R = 100\) AU
Zero initial velocities (cold collapse)
Run with Leapfrog, \(\Delta t = 0.001\) yr
Purpose: Verify your code works for \(N > 3\) before adding Plummer sampling

Phase 3B — Plummer initial conditions (\(N = 10\)–\(100\), loops):

Plummer sphere position sampler (see Science Background for derivations):
- Sample radii from the Plummer CDF using inverse-transform sampling
- Sample angles uniformly on a sphere
- Convert to Cartesian coordinates
Virial equilibrium velocities: \(v_i \sim \mathcal{N}(0, \sigma_{1D}(r_i))\) where \(\sigma_{1D}^2(r) = \frac{GM}{6a}\left(1 + \frac{r^2}{a^2}\right)^{-1/2}\)
Uniform masses (\(1\,M_\odot\) each)
Center-of-mass corrections: subtract \(\vec{r}_{\rm cm}\) and \(\vec{v}_{\rm cm}\) from all particles
Validate your initial conditions (see validation checklist below)
Run with Leapfrog for \(N = [10, 50, 100]\)
Do NOT attempt \(N > 100\) with loops

Validation checklist for Plummer IC:

Radial density profile matches Plummer model (log-log plot)
\(|r_{\rm cm}| < 0.001 \times a\) (center of mass position)
\(|v_{\rm cm}| < 0.001 \times \sigma_{3D}\) (center of mass velocity)
Initial virial ratio \(|Q - 1| < 0.01\)

What to expect physically: For a virialized Plummer sphere, the cluster should remain roughly stable, with some particles exchanging energy. Some dynamical evolution is normal — the cluster is not static.

Phase 4: Vectorize and Scale Up

Goal: Rewrite the force calculation for performance. Scale to large \(N\).

After your loop version works correctly:

Rewrite force calculation using NumPy broadcasting
No explicit for i in range(N): for j in range(i+1, N): loops
Use array operations: r_ij shape (N, N, 3) for all pairwise vectors at once
First verify: matches loop results exactly (same random seed, same \(N\))
Then scale up: \(N = [200, 500, 1000]\) or more if performance allows
Collect timing data: loop vs. vectorized at several \(N\) values

Success Criteria:

Vectorized matches loop results exactly (same random seed)
Vectorization provides significant speedup (\(10\)–\(100\times\))
Energy conservation maintained across all implementations
Can simulate \(N \geq 1000\) particles with vectorized code

Required Plots and Analysis

Phase 2: Integrator Comparison (the key figure)

Figure: \(N=2\) method comparison — the most important figure in the project:

Panel A: Orbital trajectories for all three methods (different line styles or colors)
Panel B: Relative energy error \(\frac{E_{\rm tot}(t) - E_{\rm tot,0}}{|E_{\rm tot,0}|}\) vs. time for all three methods (log scale on \(y\)-axis)

Document your observations: which method drifts, which oscillates, and why this matters for long-term simulations.

Figure: \(N=3\) verification (can go in appendix):

Orbital trajectories showing Earth’s \(\sim 24\) orbits during Jupiter’s 2

Phase 3: Cluster Simulation

Figure: Plummer IC validation:

Radial density profile (log-log) compared to the analytical Plummer profile

Figure: Cluster energy diagnostics (3-panel figure):

Energy components (\(E_K\), \(W\), \(E_{\rm tot}\)) vs. time
Relative energy error vs. time (log scale)
Virial ratio \(Q(t)\) with reference line at \(Q = 1\)

Figure: Cluster evolution snapshots:

3 panels showing cluster at \(t/t_{\rm final} = [0,\, 0.5,\, 1.0]\)
Use consistent markers and colors
Fixed axis limits across panels for comparison
Include time labels

Phase 4: Performance and Showcase

Figure: Performance comparison:

Wall-clock time vs. \(N\) for loop and vectorized implementations (log-log scale)

Figure: Large-\(N\) showcase (your largest stable \(N\)):

3 evolution snapshots at \(t/t_{\rm final} = [0,\, 0.5,\, 1.0]\)
Include energy diagnostics showing conservation at large \(N\)

Simulation Parameters

Timestep Selection:

Start with \(\Delta t = 0.01\) years for \(N=2\) (Sun-Earth)
Use \(\Delta t = 0.001\) years for \(a = 100\) AU clusters
Use \(\Delta t = 0.01\) years for \(a = 1000\) AU clusters
Rule of thumb: \(\Delta t \approx 0.01 \times\) shortest orbital period
If energy conservation poor (\(> 1\%\) drift with Leapfrog), reduce \(\Delta t\) by factor of 2

Simulation Duration:

\(N=2\): 10 orbital periods (10 years)
\(N=3\): 24 years (2 Jupiter orbits)
Clusters with \(a = 100\) AU: at least 100 years
Clusters with \(a = 1000\) AU: at least 1000 years
Stop early if any particle reaches \(r > 10a\) (escaper — this can be physical or numerical; note which)

Output Management:

Save snapshots every 10–100 timesteps (balance file size vs. temporal resolution)
Store at minimum: positions, velocities, energies (total, kinetic, potential), virial ratio
For large \(N\) and long integration times, consider saving only every 100th timestep

Implementation Tips

Data Structure Convention

Positions: shape (N, 3) — each row is [x, y, z] for one particle
Velocities: shape (N, 3) — each row is [vx, vy, vz] for one particle
Accelerations: shape (N, 3) — each row is [ax, ay, az]
Masses: shape (N,) — 1D array of masses

Testing Strategy

Unit tests: Test each component separately (force calculation, energy computation, single ODE integration step)
Conservation tests: Energy should be bounded for Leapfrog and slowly drifting for RK4
Known solutions: Two-body problem has an analytical solution
Scaling tests: Verify \(O(N^2)\) scaling for force calculations (use timeit)

Common Bug Patterns

Watch out for these frequent mistakes:

Self-force Bug: Computing force of particle on itself
- Symptom: NaN values or particles shooting to infinity
- Fix: Skip when i == j in force loop
Double Counting Bug: Counting each pair twice in potential energy
- Symptom: Potential energy exactly \(2\times\) what it should be
- Fix: Use i < j not i != j in summation
Sign Error Bug: Wrong sign in force or acceleration
- Symptom: Particles repel instead of attract
- Fix: Check force points from \(i\) toward \(j\)
Unit Vector Bug: Not normalizing direction vectors correctly
- Symptom: Forces scale incorrectly with distance
- Fix: The softened force already includes the correct \(r^{-2}\) scaling via the \(r^3\) in the denominator — don’t divide by \(|r_{ij}|\) again
Timestep Bug: Using different \(\Delta t\) in different parts of integrator
- Symptom: Energy not conserved even with symplectic integrator
- Fix: Use consistent \(\Delta t\) throughout each method
Array Shape Bug: Mixing (N, 3) and (3, N) arrays
- Symptom: Broadcasting errors or wrong results
- Fix: Be consistent with array shapes throughout
Accumulation Bug: Not resetting force array each timestep
- Symptom: Forces grow without bound
- Fix: Zero force array before computing new forces

Suggested Timeline

Week	Milestone
Week 1 (by Feb 18)	Phase 1 complete: pipeline working with Euler, \(N=2\) and \(N=3\) validated. Check-in 1.
Week 2 (by Feb 25)	Phase 2 complete: RK4 + Leapfrog working, comparison plots done. Phase 3 started (Plummer IC). Check-in 2.
Week 3 (by Mar 3)	Phases 3–4 complete: cluster simulations, vectorization, all plots, both memos. Submit by 11:59 PM.

Key Advice: If you’re behind at any checkpoint, come to office hours or post on Slack immediately. Don’t wait until the last weekend.

Deliverables

Following the Project Submission Guide, you must submit:

1. Code Repository

Well-organized modular code with clear separation of physics, numerics, and visualization
Both loop-based and vectorized implementations
Clear README with installation and usage instructions
At least 5–10 meaningful Git commits showing incremental development

2. Research Memo (2–3 pages)

Focus on your N-body cluster simulation results. Your memo should include:

Main Content:

Executive Summary: What you accomplished and key findings
Integrator Analysis: Comparison of Euler, RK4, and Leapfrog on the \(N=2\) problem — energy conservation behavior, why you chose Leapfrog for cluster simulations
Cluster Simulation Results: Key findings from Phase 3 with embedded plots:
- Plummer IC validation (radial density profile)
- Energy conservation diagnostics (\(E_K\), \(W\), \(E_{\rm tot}\) vs. time)
- Virial diagnostics: \(Q(t)\) and relative energy error
- Cluster evolution snapshots
Computational Performance: Loop vs. vectorized runtime analysis, speedup
Challenges & Solutions: What was difficult and how you solved it

Appendices (not counted toward page limit):

Appendix A: \(N=3\) verification — orbital plot confirming Jupiter’s 12-year period
Appendix B: Additional diagnostics (optional)

3. Growth Memo (1 page)

Reflect on your learning process:

Key learning moments and breakthroughs
Skill evolution from Project 1 to Project 2
How you used AI tools appropriately (within course guidelines)
Next learning goals

Submission Notes

Main memo focuses on cluster results; verification work goes in appendices
All plots must have labeled axes with units
Include figure captions explaining what each plot shows
Code should be executable with clear installation instructions

Appendix: Example Code Structure

Here’s a skeleton of how you might structure your N-body integrator. You will need to fill in the methods with the appropriate logic for force calculation, integration steps, and energy computation.

import numpy as np

class NBodySimulator:
    """N-body gravitational simulator with selectable integration method."""

    def __init__(self, masses, positions, velocities, epsilon=1.0,
                 method='leapfrog'):
        self.masses = masses          # shape (N,)
        self.positions = positions    # shape (N, 3)
        self.velocities = velocities  # shape (N, 3)
        self.epsilon = epsilon
        self.method = method
        self.G = 4 * np.pi**2        # AU^3 M_sun^-1 yr^-2

    def compute_accelerations(self, positions):
        """Compute gravitational accelerations for all particles."""
        # TODO: Implement — returns shape (N, 3)
        pass

    def euler_step(self, dt):
        """Single Euler integration step."""
        # TODO: Implement
        pass

    def rk4_step(self, dt):
        """Single RK4 integration step."""
        # TODO: Implement
        pass

    def leapfrog_step(self, dt):
        """Single Leapfrog integration step."""
        # TODO: Implement — note different structure!
        pass

    def step(self, dt):
        """Advance simulation by one timestep using selected method."""
        if self.method == 'euler':
            self.euler_step(dt)
        elif self.method == 'rk4':
            self.rk4_step(dt)
        elif self.method == 'leapfrog':
            self.leapfrog_step(dt)

    def compute_energies(self):
        """Compute kinetic, potential, and total energy."""
        # TODO: Implement — return E_K, W, E_tot
        pass

Note how the three integrators share the same interface — this is the modular design that makes Phase 2 straightforward.

Appendix: Visualization Tips for Large N

When visualizing your star clusters, choose your plotting strategy based on particle count:

Plotting Strategy by Scale

\(N \leq 10{,}000\): Use plt.scatter() with tuning
\(N > 20{,}000\): Consider density maps (hexbin/hist2d) with massive stars overlaid

Tips for Scatter Plots

Sort by depth: Plot particles sorted by \(z\)-coordinate so nearer stars appear on top
Transparency: Use alpha=0.5–0.7 with lw=0 (no edge lines)
Performance: Set rasterized=True when saving PDFs to reduce file size
Consistent visualization: Use ax.set_aspect('equal') for proper spatial representation
Fix axis limits across snapshots for easy comparison
Always include colorbars (if coloring by mass) and axis labels with units

Appendix: Essential NumPy Functions for N-body

Function	Purpose	Example Usage
`np.zeros`	Initialize arrays	`positions = np.zeros((N, 3))`
`np.random.uniform`	Random sampling	`u = np.random.uniform(0, 1, size=N)`
`np.linalg.norm`	Vector magnitudes	`r = np.linalg.norm(pos[i] - pos[j])`
`np.sum`	Summation with axis control	`total_KE = 0.5 * np.sum(m[:, np.newaxis] * v**2)`
`np.sqrt`	Vectorized square root	`distances = np.sqrt(dx2 + dy2 + dz**2)`
`np.newaxis`	Add dimension for broadcasting	`r_ij = pos[:, np.newaxis, :] - pos[np.newaxis, :, :]`
`np.arccos`	Inverse cosine for angles	`theta = np.arccos(1 - 2*u)`
`np.sin`, `np.cos`	Trig functions	`x = r * np.sin(theta) * np.cos(phi)`
`np.linspace`	Evenly spaced values	`times = np.linspace(0, t_end, n_steps)`
`np.logspace`	Log-spaced values	`mass_bins = np.logspace(-1, 2, 50)`
`np.triu_indices`	Upper triangular indices	`i, j = np.triu_indices(N, k=1)` for unique pairs
`np.fill_diagonal`	Zero self-forces	`np.fill_diagonal(force_matrix, 0)`
`np.clip`	Clamp values to range	`u_safe = np.clip(u, 1e-12, 1-1e-12)`
`np.save`/`np.load`	Save/load arrays	`np.save('positions.npy', pos)`
`timeit.timeit`	Performance testing	`time = timeit.timeit(lambda: compute_forces(pos), number=100)`

Vectorization Tips:

Use broadcasting to avoid loops: r_ij shape (N, N, 3) computes all pairwise differences at once
Mask diagonal for self-forces: np.fill_diagonal(force_matrix, 0)
Use axis parameter in reductions: np.sum(forces, axis=1) sums forces on each particle
For advanced users: np.einsum provides elegant notation for tensor operations