Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Modeling the Universe
ASTR 596 - Fall 2025

Project 2: N-body Dynamics with Monte Carlo Sampling

San Diego State University

ASTR 596 Fall 2025
Instructor: Anna Rosen
Due: Wednesday September 24, 2025 by 11:59 pm

Pair Programming Assignments for Project 2

Project 2 Pairings:

Individual Implementation Requirements: While you work together during Friday sessions:

Review the Pair Programming Guidelines for best practices.

N-body Dynamics with Monte Carlo Initial Conditions

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

Your code must be modular and object-oriented, allowing users to integrate the orbits of NN 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).

Relevant Course Materials:

Numerical Problem Description

The acceleration of particle ii in an NN-body gravitational system, including a softening parameter ϵ\epsilon, is given by:

a⃗i=∑j=1j≠iNG mjr⃗j−r⃗i(∣r⃗j−r⃗i∣2+ϵ2)3/2\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}}
(1)

where mim_i and mjm_j are the particle masses, and the distance rijr_{ij} between particles ii and jj is defined as:

rij2=(xj−xi)2+(yj−yi)2+(zj−zi)2r_{ij}^2 = (x_j - x_i)^2 + (y_j - y_i)^2 + (z_j - z_i)^2
(2)

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

ϵ≈0.01×Rcluster/N1/3,\epsilon \approx 0.01 \times R_{\text{cluster}}/N^{1/3},
(3)

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

Choosing the Softening Parameter ε\varepsilon:

Required units:

Mass:M⊙∣Length:AU∣Time:yr\begin{aligned} \text{Mass}: &\quad M_{\odot} \quad | \quad \text{Length}: \quad \text{AU} \quad | \quad \text{Time}: &\quad \text{yr} \end{aligned}
(4)

In these units, the gravitational constant becomes:

G=4π2 AU3 M⊙−1 yr−2G = 4\pi^2~\text{AU}^3 \text{ M}_{\odot}^{-1} \text{ yr}^{-2}
(5)

Energy and Virial Diagnostics: See Section 7 of the Kroupa IMF & Plummer Sphere Guide for complete formulae and physical interpretation.

You must monitor and plot the evolution of:

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<ji<j in the sum).

ODE Integration Methods Reference

For N-body dynamics, you’ll implement four integration schemes. Full derivations in course materials:

Key Algorithmic Structures for the system y⃗˙=f⃗(t,y⃗)\dot{\vec{y}} = \vec{f}(t, \vec{y}) where y⃗=[r⃗,v⃗]\vec{y} = [\vec{r}, \vec{v}] and f⃗\vec{f} returns [v⃗,a⃗][\vec{v}, \vec{a}]:

Euler’s Method (1st order, simplest):

y⃗n+1=y⃗n+Δt⋅f⃗(tn,y⃗n)\vec{y}_{n+1} = \vec{y}_n + \Delta t \cdot \vec{f}(t_n, \vec{y}_n)
(6)

RK2 (Midpoint Method; 2nd order):

k⃗1=f⃗(tn,y⃗n)k⃗2=f⃗(tn+Δt2,y⃗n+Δt2k⃗1)y⃗n+1=y⃗n+Δt⋅k⃗2\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{y}_{n+1} &= \vec{y}_n + \Delta t \cdot \vec{k}_2 \end{aligned}
(7)

RK4 (4th order, high accuracy):

k⃗1=f⃗(tn,y⃗n)k⃗2=f⃗(tn+Δt2,y⃗n+Δt2k⃗1)k⃗3=f⃗(tn+Δt2,y⃗n+Δt2k⃗2)k⃗4=f⃗(tn+Δt,y⃗n+Δt⋅k⃗3)y⃗n+1=y⃗n+Δt6(k⃗1+2k⃗2+2k⃗3+k⃗4)\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}
(8)

Leapfrog (2nd order, symplectic - conserves energy):

v⃗n+1/2=v⃗n+Δt2⋅a⃗(r⃗n)(half-step velocity)r⃗n+1=r⃗n+Δt⋅v⃗n+1/2(full-step position)v⃗n+1=v⃗n+1/2+Δt2⋅a⃗(r⃗n+1)(half-step velocity)\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}
(9)

Critical Implementation Notes:

Steps to Complete the Problem

Why This Progression?

We build complexity incrementally to isolate and debug different aspects:

Why Implement All Four Methods?

You’re implementing methods of increasing sophistication to understand the trade-offs:

Building from simple to complex helps you:

Each step builds on the previous one - don’t skip ahead until the current step works correctly!

Step 1: Two-Body System (Sun-Earth)

Goal: Validate your four integration methods against known physics.

Implementation:

What You Should Observe:

After This Step: Select your two best methods for the remaining work. You should find that Leapfrog (for energy conservation) and RK4 (for accuracy) perform best. If your results don’t show this, debug before proceeding!

Step 2: Three-Body System (Sun-Earth-Jupiter)

Goal: Verify correct force summation between multiple bodies.

Implementation:

Success Criteria:

Step 3: N-Body Star Cluster

Goal: Implement realistic star cluster with Monte Carlo initial conditions, first with loops then with vectorization.

Implementation Phase A - Simple Test with Loops (N=10N=10):

Implementation Phase B - Full Physics with Loops:

Implementation Phase C - Vectorized Implementation: After your loop version works correctly:

Success Criteria:

Required Plots and Analysis

Standard Energy Diagnostics (Apply to ALL Steps)

For every simulation phase, create a 3-panel figure showing:

  1. Energy Components: Plot EKE_K, WW, and EtotE_{tot} vs time on same axes

  2. Relative Energy Error: Plot (Etot(t)−Etot,0)/∣Etot,0∣(E_{tot}(t) - E_{tot,0})/|E_{tot,0}| vs time (log scale)

  3. Virial Ratio: Plot Q=∣2EK+W∣/∣W∣Q = |2E_K + W|/|W| vs time

Expected values:

Step-Specific Additional Plots

Step 1 (Two-Body):

Step 2 (Three-Body):

Step 3A (Simple N-Body):

Step 3B (Full Physics):

Step 3C (Vectorized):

Key Implementation Steps

  1. Implement Kroupa mass sampling using analytical normalization (Section 2 of guide)

  2. Implement Plummer position sampling (Section 3 of guide)

  3. Assign virial equilibrium velocities (Section 3.6 of guide): draw each component from N(0,σ1D(r))\mathcal{N}(0, \sigma_{1D}(r))

  4. Apply center-of-mass corrections: subtract r⃗cm\vec{r}_{\text{cm}} and v⃗cm\vec{v}_\text{cm} from all particles

Validation Checklist

Before production runs, Validate your initial conditions by verifying:

Simulation Parameters

Timestep Selection:

Timestep Stability Criterion:

Δt<ηmin⁡(ϵ∣a⃗∣max⁡,0.01×min orbital period2π)\Delta t < \eta \min\left(\sqrt{\frac{\epsilon}{|\vec{a}|_{\max}}}, \frac{0.01 \times \text{min orbital period}}{2\pi}\right)
(10)

where η≈0.01−0.1\eta \approx 0.01-0.1 is a safety factor.

Simulation Duration:

Output Management:

Implementation Tips

Data Structure Convention:

Testing Strategy

Common Bug Patterns

Watch out for these frequent mistakes:

  1. 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

  2. Double Counting Bug: Counting each pair twice in potential energy

    • Symptom: Potential energy exactly 2× what it should be

    • Fix: Use i<j not i≠j in summation

  3. Sign Error Bug: Wrong sign in force or acceleration

    • Symptom: Particles repel instead of attract

    • Fix: Check force points from i to j

  4. Unit Vector Bug: Not normalizing direction vectors

    • Symptom: Forces scale incorrectly with distance

    • Fix: Divide by |r_ij| not |r_ij|²

  5. Timestep Bug: Using different dt in different parts of integrator

    • Symptom: Energy not conserved even with symplectic integrator

    • Fix: Use consistent dt throughout each method

  6. Array Shape Bug: Mixing (N,3) and (3,N) arrays

    • Symptom: Broadcasting errors or wrong results

    • Fix: Be consistent with array shapes throughout

  7. Accumulation Bug: Not resetting force array each timestep

    • Symptom: Forces grow without bound

    • Fix: Zero force array before computing new forces

Suggested Timeline & Milestones

To ensure steady progress and avoid last-minute rushing, aim for these checkpoints:

By End of Day 3 (Wednesday)

✓ Two-body system working with at least Euler and one other integrator
✓ Energy conservation plots showing method differences
✓ Repository set up with initial commits

By End of Week 1 (Sunday)

✓ All four integrators implemented and tested on two-body problem
✓ Three-body (Sun-Earth-Jupiter) system working correctly
✓ Basic N-body working with loops (N=10, uniform masses)
✓ At least 3-4 meaningful Git commits

By End of Day 10 (Wednesday of Week 2)

✓ Vectorized implementation complete and tested
✓ Performance comparison data collected
✓ Kroupa IMF sampling with analytical normalization working

By End of Week 2 (Sunday)

✓ Plummer spatial distribution implemented
✓ Center-of-mass corrections applied
✓ Chosen extension implemented
✓ All plots generated
✓ Both memos written

Monday (Due Date)

✓ Final code review and cleanup
✓ Verify all deliverables complete
✓ Submit by 11:59 PM

Key Advice: If you’re behind at any checkpoint, come to Student Hacking Hours (Monday/Wednesday 1-2 pm) or post on Slack immediately. Don’t wait until the last weekend!

Example Extension Ideas

You must implement at least one extension beyond the base requirements. Extensions are graded on depth of exploration, not just implementation. A thoughtful analysis of one extension is better than superficial implementation of three. Let your curiosity drive you! Talk to me if you want to discuss other extensions of interest.

Performance Extensions

Physics Extensions

Analysis Extensions

Visualization Extensions

Deliverables

Following the Project Submission Guide, you must submit:

1. Code Repository

2. Research Memo (2-3 pages)

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

Main Content:

Appendices (not counted toward page limit):

3. Growth Memo (1 page)

Reflect on your learning process:

Submission Notes

Appendix: Example Class Structure

Appendix: Example Code Structure

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

import numpy as np
import ODEIntegrators as odeint  # Your module containing integration methods

class NBodyIntegrator:
    """Base class for N-body integration."""

    def __init__(self, masses, positions, velocities, epsilon=1.0, method='leapfrog'):
        # TODO: Implement
        pass

    def compute_accelerations(self, positions):
        # TODO: Implement
        pass

    def step(self, dt):
        # TODO: Implement
        pass

    def compute_energy(self):
        # TODO: Implement
        pass

Appendix: Visualization Tips for Large N

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

Plotting Strategy by Scale

Tips for N ≈ 10,000 Scatter Plots

  1. Sort by depth: Plot particles sorted by z-coordinate so nearer stars appear on top

  2. Size scaling: Try size ∝ m^(2/3) to make massive stars visible without overwhelming

  3. Transparency: Use alpha=0.5-0.7 with lw=0 (no edge lines)

  4. Performance: Set rasterized=True when saving PDFs to reduce file size

  5. Highlight interesting features: Overlay the top 100-200 most massive stars with slightly larger, opaque markers

When Scatter Plots Struggle

If your cluster core looks like an ink blot:

Consistent Visualization

Appendix: Essential NumPy Functions for N-body

FunctionPurposeExample Usage
np.zerosInitialize arrayspositions = np.zeros((N, 3))
np.random.uniformRandom samplingu = np.random.uniform(0, 1, size=N)
np.random.randomRandom [0,1) valuesrand_vals = np.random.random(N)
np.linalg.normVector magnitudesr = np.linalg.norm(pos[i] - pos[j])
np.crossCross product for perpendicular vectorsv_perp = np.cross(r_vec, z_hat)
np.sumSummation with axis controltotal_KE = np.sum(0.5 * m * v**2, axis=0)
np.sqrtVectorized square rootdistances = np.sqrt(dx**2 + dy**2 + dz**2)
np.newaxisAdd dimension for broadcastingr_ij = pos[:, np.newaxis, :] - pos[np.newaxis, :, :]
np.whereConditional operationsmassive = masses[np.where(masses > 10)]
np.arccosInverse cosine for anglestheta = np.arccos(1 - 2*u)
np.sin, np.cosTrig functionsx = r * np.sin(theta) * np.cos(phi)
np.arangeInteger sequencesindices = np.arange(N)
np.linspaceEvenly spaced valuestimes = np.linspace(0, t_end, n_steps)
np.logspaceLog-spaced values for mass binsmass_bins = np.logspace(-1, 2, 50)
np.histogramBinning data for analysiscounts, bins = np.histogram(masses, bins=mass_bins)
np.triu_indicesUpper triangular indicesi, j = np.triu_indices(N, k=1) for unique pairs
np.fill_diagonalZero self-forcesnp.fill_diagonal(force_matrix, 0)
np.clipClamp values to rangeu_safe = np.clip(u, 1e-12, 1-1e-12)
np.save/np.loadSave/load arraysnp.save('positions.npy', pos)
timeit.timeitPerformance testingtime = timeit.timeit(lambda: compute_forces(pos), number=100)

Vectorization Tips:

Project 1
Starter Scripts for Project 1
Project 2
Project 2: Science Background