Fixing Three-Body Leapfrog Integration: Half-Step Scheme and Hamiltonian Energy Tracking in Python

Learn how to implement a time‑centered leapfrog for the three‑body problem in Python, with forces and unique-pair potential for reliable Hamiltonian tracking.

Tracking the Hamiltonian in a three-body setup is a straightforward way to sanity‑check a gravitational integrator. If the time integrator is implemented consistently, the energy time series should evolve in a controlled manner. Below is a compact walkthrough of a three-body leapfrog, what in the initial variant makes it hard to reason about, and how to align it with the canonical half‑step scheme to make the Hamiltonian readout more reliable.

Original implementation and the pain point

The goal is to evolve three bodies with unit masses, record the Hamiltonian T + V at each step, and stop early if two bodies get too close. The positions and velocities are 2D NumPy arrays, organized as lists of length three.

from numpy import sum as npsum
from numpy.linalg import norm
from copy import deepcopy

def energy_fn(q, u, w):
    T = npsum([w[i] * norm(u[i])**2 / 2 for i in range(3)])
    V = -npsum([w[0 - i] * w[1 - i] / norm(q[0 - i] - q[1 - i]) for i in range(3)])
    return T + V

def accel_single(q, idx, w):
    return w[idx - 1] * (q[idx - 1] - q[idx]) / (norm(q[idx - 1] - q[idx])**3) \
         + w[idx - 2] * (q[idx - 2] - q[idx]) / (norm(q[idx - 2] - q[idx])**3)

def has_collision(q):
    for k in range(3):
        if norm(q[0 - k] - q[1 - k]) < 0.1:
            return True

# leapfrog flavor
def leapfrog_run(Q, U, h, iters, w):
    pos, vel = deepcopy(Q), deepcopy(U)
    energy_log = [energy_fn(pos, vel, w)]

    for _ in range(iters):
        acc_now = [accel_single(pos, i, w) for i in range(3)]

        for i in range(3):
            pos[i] = pos[i] + vel[i] * h + 0.5 * acc_now[i] * h**2

        for i in range(3):
            vel[i] = vel[i] + 0.5 * (acc_now[i] + accel_single(pos, i, w)) * h

        if has_collision(pos):
            return energy_log
        energy_log.append(energy_fn(pos, vel, w))

    return energy_log

This variant advances positions using acceleration at the current time and then updates velocities with the average of old and new accelerations. The potential energy is computed over the three unique pairs via negative indexing, and the early‑stop condition is a proximity threshold.

What’s actually off

The confusion stems from the update order. The code mixes full‑step positions with a velocity update that implicitly references two different times for the acceleration. This pattern is easy to get wrong when you want the clean, time‑centered behavior that a “textbook” leapfrog is known for. Another source of cognitive overhead is pair handling via negative indices: while it hits each unique pair once for three bodies, it is opaque and easy to misread. When the objective is to document Hamiltonian variation, it helps to make both the force assembly and the half‑step structure explicit.

Corrected, time‑centered leapfrog

The adjusted routine follows the half‑step velocity update, full‑step position update, and then the second half‑step for velocities. Forces are computed for all bodies at once, and the Hamiltonian uses a clear unique‑pair sum.

import numpy as np
from numpy.linalg import norm
from copy import deepcopy

def energy_total(r, v, mass):
    T = sum([mass[i] * norm(v[i])**2 / 2 for i in range(3)])
    V = -sum([mass[i] * mass[j] / norm(r[i] - r[j])
              for i in range(3) for j in range(i + 1, 3)])
    return T + V

def accel_all(pos, mass):
    acc = [np.zeros(2) for _ in range(3)]
    for i in range(3):
        for j in range(3):
            if i != j:
                diff = pos[i] - pos[j]
                acc[i] -= mass[j] * diff / (norm(diff)**3)
    return acc

def touch(pos):
    for i in range(3):
        for j in range(i + 1, 3):
            if norm(pos[i] - pos[j]) < 0.1:
                return True
    return False

def leapfrog_solver(R0, V0, dt, steps, mass):
    r, v = deepcopy(R0), deepcopy(V0)
    H_series = [energy_total(r, v, mass)]

    acc_prev = accel_all(r, mass)

    for _ in range(steps):
        v_half = [v[i] + 0.5 * acc_prev[i] * dt for i in range(3)]
        r = [r[i] + v_half[i] * dt for i in range(3)]
        acc_new = accel_all(r, mass)
        v = [v_half[i] + 0.5 * acc_new[i] * dt for i in range(3)]

        if touch(r):
            return H_series
        H_series.append(energy_total(r, v, mass))
        acc_prev = acc_new

    return H_series

Why the adjusted version is better for energy tracking

In this structure, velocities live on half steps and positions on full steps, which makes the ordering unambiguous. Acceleration is computed once per step for all bodies, avoiding timing mismatches. The potential energy is a clear sum over unique pairs, which directly matches the intent when reviewing energy terms. Together, these choices make the Hamiltonian series easier to interpret over time.

Wrap‑up

If your target is a trustworthy Hamiltonian time series for the three‑body problem, stick to an explicit half‑step leapfrog. Compute forces for all bodies in one place, sum potential over unique pairs, and evaluate energy after each completed step. With these pieces in place, the energy log you produce is much more straightforward to analyze.