Why Does RK45 Deform and Shift N-Body Orbits?

When implementing Runge-Kutta adaptive methods from scratch for orbital mechanics, it is common to encounter trajectories that look roughly correct in shape but appear stretched, shifted, or fail to close periodic loops (such as the Golden Pancake Drama or Broucke A2 orbits). If your fixed-step RK4 works fine but your adaptive RK45 solver distorts the orbits, two critical bugs are almost certainly at play:

  • Decoupled RK Stages: You are likely stepping each body through intermediate Runge-Kutta stages individually while holding the other bodies fixed at their initial positions. In an N-body system, the gravitational acceleration on body $A$ at intermediate stage $k_n$ depends on the position of body $B$ at that exact same intermediate stage.
  • Step Size ($ ext{d}t$) Mutation Bug: The step size is scaled before accepting the step or returning it to the main loop, meaning the accepted state increment (computed with the old $ ext{d}t$) is matched with an incorrect time step increment.

1. The Decoupled Intermediate Stages Bug

In the original implementation, the function rk45_i calculates $k_1$ through $k_6$ for a single body by querying the positions of the other bodies in state_arr. While calculating stage $k_2$ for body 1, body 2 is still at stage $k_0$.

A system of $N$ interacting bodies is a single coupled system of ordinary differential equations (ODEs). The state vector $oldsymbol{Y}$ must represent all bodies simultaneously:

# Flattened or structured unified state
Y = [x1, y1, vx1, vy1, x2, y2, vx2, vy2, x3, y3, vx3, vy3]

When computing $k_2 = f(t + c_2 ext{d}t, oldsymbol{Y} + a_{21} k_1)$, every body's position must be displaced by its corresponding component of $k_1$ before evaluating accelerations.

2. The Adaptive Step Size Logic Bug

In standard Runge-Kutta-Fehlberg (RKF45), the local truncation error per unit step is calculated to determine whether a step is accepted or rejected:

delta = 0.84 * (err_tol / r_val) ** 0.25  # Safety factor included

If the step succeeds (r_val <= err_tol), the state must be updated using the current $ ext{d}t$ that was used to compute those $k$-values. Only then should the step size for the next step be scaled by $ ext{delta}$. Modifying $ ext{d}t$ beforehand causes time accumulation to drift out of sync with position and velocity.

Corrected Full Implementation

Below is a refactored, vectorised implementation of the RKF45 solver from scratch. It solves the coupled system together and accurately traces periodic three-body solutions:

import numpy as np
import matplotlib.pyplot as plt

# System constants
BIG_G = 1.0
SOFTENING = 1e-9  # Keep small to avoid artificial damping

# Initial conditions for Golden Pancake Drama (GPD)
masses = np.array([1.0, 0.6123560888279949, 1.0])
# State array: [x, y, vx, vy] for each body
initial_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)

duration_t = 8.8457131110
err_tolerance = 1e-10
max_step = 1e-2
min_step = 1e-12
initial_dt = 1e-4

def compute_derivatives(state, m):
    """
    Evaluates the derivatives d/dt [x, y, vx, vy] for all bodies coupled together.
    """
    n_bodies = len(m)
    derivatives = np.zeros_like(state)
    
    # Positions and velocities
    pos = state[:, 0:2]
    vel = state[:, 2:4]
    
    # dx/dt = vx, dy/dt = vy
    derivatives[:, 0:2] = vel
    
    # Compute coupled gravitational acceleration
    acc = np.zeros_like(pos)
    for i in range(n_bodies):
        for j in range(n_bodies):
            if i != j:
                diff = pos[j] - pos[i]
                dist_sq = np.sum(diff**2) + SOFTENING**2
                inv_dist_cube = 1.0 / (dist_sq * np.sqrt(dist_sq))
                acc[i] += BIG_G * m[j] * diff * inv_dist_cube
                
    derivatives[:, 2:4] = acc
    return derivatives

def rkf45_step(state, dt, m, tol, min_h, max_h):
    """
    Performs a single adaptive Runge-Kutta-Fehlberg 4(5) step across all bodies.
    """
    # Coefficients for Fehlberg RK45
    c20, c21 = 1/4, 1/4
    c30, c31, c32 = 3/8, 3/32, 9/32
    c40, c41, c42, c43 = 12/13, 1932/2197, -7200/2197, 7296/2197
    c50, c51, c52, c53, c54 = 1.0, 439/216, -8.0, 3680/513, -845/4104
    c60, c61, c62, c63, c64, c65 = 1/2, -8/27, 2.0, -3544/2565, 1859/4104, -11/40

    # 4th and 5th order weighting coefficients
    b4 = np.array([25/216, 0.0, 1408/2565, 2197/4104, -1/5, 0.0])
    b5 = np.array([16/135, 0.0, 6656/12825, 28561/56430, -9/50, 2/55])
    
    current_dt = dt

    while True:
        k1 = current_dt * compute_derivatives(state, m)
        k2 = current_dt * compute_derivatives(state + c21*k1, m)
        k3 = current_dt * compute_derivatives(state + c31*k1 + c32*k2, m)
        k4 = current_dt * compute_derivatives(state + c41*k1 + c42*k2 + c43*k3, m)
        k5 = current_dt * compute_derivatives(state + c51*k1 + c52*k2 + c53*k3 + c54*k4, m)
        k6 = current_dt * compute_derivatives(state + c61*k1 + c62*k2 + c63*k3 + c64*k4 + c65*k5, m)

        # Difference between 4th and 5th order approximations: local error
        error_vector = (b5[0]-b4[0])*k1 + (b5[2]-b4[2])*k3 + (b5[3]-b4[3])*k4 + (b5[4]-b4[4])*k5 + (b5[5]-b4[5])*k6
        error = np.max(np.abs(error_vector))

        # Compute scale factor with a 0.84 safety margin
        if error == 0:
            factor = 2.0
        else:
            factor = 0.84 * (tol * current_dt / error) ** 0.25

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

        if error <= tol * current_dt or current_dt <= min_h:
            # Accept the step using the 5th order solution
            new_state = state + (b5[0]*k1 + b5[2]*k3 + b5[3]*k4 + b5[4]*k5 + b5[5]*k6)
            next_dt = np.clip(current_dt * factor, min_h, max_h)
            return new_state, current_dt, next_dt
        else:
            # Reject step and shrink dt
            current_dt = np.clip(current_dt * factor, min_h, max_h)

# Main simulation loop
time = 0.0
dt = initial_dt
state = initial_state.copy()

traj = [[], [], []]

while time < duration_t:
    # Prevent stepping past total duration
    if time + dt > duration_t:
        dt = duration_t - time

    state, step_used, dt = rkf45_step(state, dt, masses, err_tolerance, min_step, max_step)
    time += step_used

    for i in range(3):
        traj[i].append((state[i, 0], state[i, 1]))

# Visualization
plt.figure(figsize=(6, 6))
colors = ['blue', 'red', 'green']
for i in range(3):
    t = np.array(traj[i])
    plt.plot(t[:, 0], t[:, 1], color=colors[i], label=f'm{i+1}', linewidth=0.8)

plt.title('GOLDEN PANCAKE DRAMA (RKF45 Corrected)')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.axis('equal')
plt.show()

Key Takeaways for Orbital Simulations

  • Always couple full phase spaces: Avoid evaluating Runge-Kutta stages on isolated bodies. Derivative evaluations require simultaneous positions of all interacting celestial bodies.
  • Step acceptance before step updates: Use the step size that successfully generated the coefficients for advancing simulation time, then compute the next suggested step size.
  • Energy conservation in long orbits: While RKF45 is accurate for non-stiff initial value problems, Runge-Kutta methods are not symplectic. For simulating chaotic or multi-period orbits over large time horizons, consider using symplectic integrators such as the Velocity Verlet or Yoshida 4th-order integrator to prevent orbital energy drift.