Simulating gravitational dynamics—such as the famous three-body problem or solar system orbits—requires robust numerical integration. The Runge-Kutta-Fehlberg method (RK45 or RKF45) is an adaptive step-size integrator that balances accuracy and performance by estimating local truncation errors.

However, when moving from a single particle to an $N$-body system, many developers run into issues where simulations get stuck in infinite loops, step sizes collapse to min_step, or bodies veer off into deep space. Here is why this happens and how to build a working, robust RK45 $N$-body integrator from scratch.

The Critical Flaws in the Decoupled Approach

In the question's original code, three primary bugs prevent the simulation from advancing properly:

  • Coupled System Evaluated Independently: The gravitational acceleration of body $A$ depends directly on the current position of body $B$. In intermediate Runge-Kutta stages ($k_1, k_2, ext{...}, k_6$), all bodies must evaluate their positions simultaneously. Integrating one body while reading the static, un-stepped positions of the other bodies completely breaks the physics and creates massive artificial truncation errors.
  • Premature Return Inside the Step Loop: Inside rk45_step, the statement return out_arr, dt is indented inside the for loop. As a result, only body 0 gets updated, while bodies 1 and 2 are left untouched before the function exits.
  • Python List Mutability (Pass-by-Reference): Assigning out_arr = inp_arr creates a shallow alias. Mutating elements in out_arr directly mutates inp_arr during intermediate failed step trials, corrupting the initial conditions.

The Solution: Treating the System as a Single State Vector

An $N$-body system is not $N$ separate ordinary differential equations (ODEs); it is a single coupled system of $4N$ (or $6N$ in 3D) first-order differential equations. The entire state vector $\mathbf{Y}$ should be integrated simultaneously:

# Y is structured as an (N, 4) array: [x, y, vx, vy] for each body
state = np.array([
    [-1.0, 0.0,  0.151837333844,  0.139577137808],
    [ 0.0, 0.0, -0.495911893796, -0.455869192302],
    [ 1.0, 0.0,  0.151837333844,  0.139577137808]
])

Complete Working Implementation

Here is an end-to-end, functional implementation of the RKF45 integrator configured for the periodic "Golden Pancake Drama" three-body orbit:

import numpy as np
import matplotlib.pyplot as plt

# System parameters
big_G = 1.0
softening = 5e-5
masses = np.array([1.0, 0.6123560888279949, 1.0])
duration_t = 8.8457131110
err_tolerance = 1e-5
max_step = 1e-2
min_step = 1e-9

# Initial state: shape (N, 4) -> [x, y, vx, vy]
state = np.array([
    [-1.000000000000, 0.000000000000,  0.151837333844,  0.139577137808],
    [ 0.000000000000, 0.000000000000, -0.495911893796, -0.455869192302],
    [ 1.000000000000, 0.000000000000,  0.151837333844,  0.139577137808]
], dtype=np.float64)

def compute_derivatives(y, masses):
    """
    Computes dY/dt for all N bodies simultaneously.
    Returns an array of shape (N, 4): [vx, vy, ax, ay]
    """
    n_bodies = y.shape[0]
    derivs = np.zeros_like(y)
    
    # Velocity is the time derivative of position
    derivs[:, 0:2] = y[:, 2:4]
    
    # Pairwise gravitational accelerations
    pos = y[:, 0:2]
    for i in range(n_bodies):
        diff = pos - pos[i]  # Vectors pointing from body i to all other bodies
        r2 = np.sum(diff**2, axis=1) + softening**2
        inv_r3 = 1.0 / (r2 ** 1.5)
        
        # Exclude self-interaction
        inv_r3[i] = 0.0
        
        # Acceleration: a_i = sum_j (G * m_j * r_ij / |r_ij|^3)
        acc = big_G * np.sum(diff * (masses * inv_r3)[:, np.newaxis], axis=0)
        derivs[i, 2:4] = acc
        
    return derivs

def rkf45_step(y, dt, tol, masses):
    """
    Performs a single adaptive Runge-Kutta-Fehlberg step.
    """
    # Fehlberg 7(6) Butcher coefficients for 4th/5th order pairs
    k1 = dt * compute_derivatives(y, masses)
    k2 = dt * compute_derivatives(y + (1/4)*k1, masses)
    k3 = dt * compute_derivatives(y + (3/32)*k1 + (9/32)*k2, masses)
    k4 = dt * compute_derivatives(y + (1932/2197)*k1 - (7200/2197)*k2 + (7296/2197)*k3, masses)
    k5 = dt * compute_derivatives(y + (439/216)*k1 - 8*k2 + (3680/513)*k3 - (845/4104)*k4, masses)
    k6 = dt * compute_derivatives(y - (8/27)*k1 + 2*k2 - (3544/2565)*k3 + (1859/4104)*k4 - (11/40)*k5, masses)

    # 4th-order estimate
    y4 = y + (25/216)*k1 + (1408/2565)*k3 + (2197/4104)*k4 - (1/5)*k5
    
    # Truncation error between 4th and 5th order
    error = np.abs((1/360)*k1 - (128/4275)*k3 - (2197/4104)*k4 + (1/50)*k5 + (2/55)*k6)
    max_error = np.max(error)

    # Safety scaling factor
    safety = 0.84
    if max_error > 0:
        scale = safety * (tol / max_error) ** 0.25
    else:
        scale = 4.0

    scale = np.clip(scale, 0.1, 4.0)

    if max_error <= tol:
        # Step accepted
        return True, y4, scale * dt
    else:
        # Step rejected, try again with smaller dt
        return False, y, scale * dt

# Simulation loop
time = 0.0
dt = 1e-4
trajectory = [state.copy()]
timestamps = [time]

while time < duration_t:
    dt = np.clip(dt, min_step, max_step)
    accepted, new_state, next_dt = rkf45_step(state, dt, err_tolerance, masses)
    
    if accepted:
        time += dt
        state = new_state
        trajectory.append(state.copy())
        timestamps.append(time)
        
    dt = next_dt

trajectory = np.array(trajectory)

# Visualization
plt.figure(figsize=(7, 7))
colors = ['blue', 'crimson', 'green']
for i in range(3):
    plt.plot(trajectory[:, i, 0], trajectory[:, i, 1], label=f'Body {i+1}', color=colors[i], lw=1)

plt.title('Golden Pancake Drama (RKF45 Integration)')
plt.xlabel('X')
plt.ylabel('Y')
plt.axis('equal')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.show()

Key Takeaways for Orbital Simulations

  • Always vectorize the entire system: Never step bodies individually in high-order coupled ODE integrators. Intermediate stage points must sample all bodies at corresponding intermediate moments.
  • Step acceptance vs. rejection: If a step fails the error tolerance (max_error > tol), the simulation clock must not advance; only the time step size $dt$ should be adjusted.
  • Energy Conservation Notice: While RK45 is excellent for trajectories with changing timescales (e.g., hyperbolic flybys or space probes), standard Runge-Kutta schemes are not symplectic. For very long integrations (thousands of orbits), look into geometric integrators like the Velocity Verlet or Yoshida 4th-order Symplectic methods to avoid artificial orbital decay.