hamiltonian-mechanics

Hamiltonian mechanics — symplectic integrators (leapfrog, Yoshida), Hamilton's equations, Poisson brackets, canonical transformations, action-angle variables, and KAM theory analysis. Use for energy-conserving long-time integration of conservative systems.

Jump to install

Repository: synthetic-sciences/openscience
Last source activity: July 4, 2026 at 06:25
Detected SKILL.md language: English
Stars: 3,362
Forks: 454

Install options

The review-first prompt is selected by default. You can switch to a direct command or download a local copy.

Review the source files

Read SKILL.md and any companion files shown by SkillsMP before deciding whether to install.

Back to search Back to installation

name	hamiltonian-mechanics
description	Hamiltonian mechanics — symplectic integrators (leapfrog, Yoshida), Hamilton's equations, Poisson brackets, canonical transformations, action-angle variables, and KAM theory analysis. Use for energy-conserving long-time integration of conservative systems.
category	physics
version	1.0.0
author	Synthetic Sciences
license	MIT
tags	["Hamiltonian","Symplectic","Leapfrog","Classical Mechanics","Energy Conservation","Physics"]
dependencies	["scipy>=1.11.0","numpy>=1.24.0","matplotlib>=3.7.0"]

Hamiltonian Mechanics

Overview

Solve Hamilton's equations using symplectic integrators that exactly conserve the symplectic structure and approximately conserve energy for exponentially long times. Essential for N-body problems, celestial mechanics, molecular dynamics, and any conservative system needing long-time accuracy.

When to Use

Conservative (energy-preserving) systems
Long-time integration (thousands of orbits, molecular dynamics)
N-body gravitational or Coulomb problems
When energy drift from standard RK45 is unacceptable
Phase space structure analysis (KAM tori, resonances)

Do NOT Use When

System has dissipation (use ode-solver with RK45/Radau)
System is stiff (symplectic methods are explicit → CFL-limited)
You need adaptive time-stepping (symplectic methods use fixed dt)

Core Workflows

1. Leapfrog / Stormer-Verlet (2nd Order Symplectic)

import numpy as np
import matplotlib.pyplot as plt

def leapfrog(dH_dq, dH_dp, q0, p0, dt, n_steps):
    """
    Leapfrog (Stormer-Verlet) symplectic integrator.

    H(q, p) is the Hamiltonian.
    dH_dq = ∂H/∂q (returns force: dp/dt = -∂H/∂q)
    dH_dp = ∂H/∂p (returns velocity: dq/dt = ∂H/∂p)
    """
    n = len(q0)
    q = np.zeros((n_steps + 1, n))
    p = np.zeros((n_steps + 1, n))
    q[0] = q0
    p[0] = p0

    for i in range(n_steps):
        # Half-step momentum
        p_half = p[i] - 0.5 * dt * dH_dq(q[i], p[i])
        # Full-step position
        q[i+1] = q[i] + dt * dH_dp(q[i], p_half)
        # Half-step momentum
        p[i+1] = p_half - 0.5 * dt * dH_dq(q[i+1], p_half)

    return q, p



 ():
    r = np.linalg.norm(q)
     q / r**  

 ():
     p  


q0 = np.array([, ])
p0 = np.array([, ])  
dt = 
n_steps =   

q, p = leapfrog(dH_dq, dH_dp, q0, p0, dt, n_steps)


H =  * np.(p**, axis=) -  / np.linalg.norm(q, axis=)
()

fig, axes = plt.subplots(, , figsize=(, ))
axes[].plot(q[:, ], q[:, ], , linewidth=)
axes[].plot(, , , markersize=, label=)
axes[].set_xlabel()
axes[].set_ylabel()
axes[].set_title()
axes[].set_aspect()
axes[].legend()
axes[].grid(, alpha=)

t = np.arange(n_steps + ) * dt
axes[].plot(t, (H - H[])/(H[]), , linewidth=)
axes[].set_xlabel()
axes[].set_ylabel()
axes[].set_title()
axes[].grid(, alpha=)

plt.tight_layout()
plt.savefig(, dpi=, bbox_inches=)

def yoshida4(dH_dq, dH_dp, q0, p0, dt, n_steps):
    """
    4th-order symplectic integrator (Yoshida 1990).
    Composes three leapfrog steps with specific coefficients.
    """
    # Yoshida coefficients
    w1 = 1 / (2 - 2**(1/3))
    w0 = -2**(1/3) * w1
    c = [w1/2, (w0+w1)/2, (w0+w1)/2, w1/2]
    d = [w1, w0, w1, 0]  # last is unused

    n = len(q0)
    q = np.zeros((n_steps + 1, n))
    p = np.zeros((n_steps + 1, n))
    q[0] = q0
    p[0] = p0

    for i in range(n_steps):
        qi, pi = q[i].copy(), p[i].copy()
        for k in range(4):
            qi = qi + c[k] * dt * dH_dp(qi, pi)
            if k < 3:
                pi = pi - d[k] * dt * dH_dq(qi, pi)
        q[i+1] = qi
        p[i+1] = pi

    return q, p

# Compare leapfrog vs Yoshida on same problem
q_y, p_y = yoshida4(dH_dq, dH_dp, q0, p0, dt, n_steps)
H_y = 0.5 * np.sum(p_y**2, axis=1) - 1 / np.linalg.norm(q_y, axis=1)
print(f"Yoshida4 energy: max |ΔH/H₀| = {np.max(np.abs((H_y - H_y[0])/H_y[0])):.2e}")

def nbody_leapfrog(masses, q0, p0, dt, n_steps, G=1.0):
    """
    N-body gravitational simulation with leapfrog.
    q0: (N, 3) positions
    p0: (N, 3) momenta
    """
    N = len(masses)

    def acceleration(q):
        a = np.zeros_like(q)
        for i in range(N):
            for j in range(i+1, N):
                rij = q[j] - q[i]
                r = np.linalg.norm(rij)
                force = G * masses[i] * masses[j] / r**3 * rij
                a[i] += force / masses[i]
                a[j] -= force / masses[j]
        return a

    q = q0.copy()
    p = p0.copy()
    v = p / masses[:, None]

    trajectory = [q.copy()]

    a = acceleration(q)
    for step in range(n_steps):
        v += 0.5 * dt * a
        q += dt * v
        a = acceleration(q)
        v += 0.5 * dt * a

        if step % 100 == 0:
            trajectory.append(q.copy())

    return np.array(trajectory)

# 3-body problem (figure-8 solution)
masses = np.array([1.0, 1.0, 1.0])
# Chenciner-Montgomery initial conditions
q0 = np.array([
    [-0.97000436, 0.24308753, 0],
    [0.97000436, -0.24308753, 0],
    [0, 0, 0]
])
v0_base = np.array([0.4662036850, 0.4323657300, 0])
p0 = np.array([-v0_base/2, -v0_base/2, v0_base])

traj = nbody_leapfrog(masses, q0, p0, dt=0.001, n_steps=100000)

def poincare_section_hamiltonian(dH_dq, dH_dp, H_func, q0, p0, dt, n_steps,
                                  section_idx=0, section_val=0.0):
    """
    Compute Poincare section for a Hamiltonian system.
    Collects (q_j, p_j) when q[section_idx] crosses section_val.
    """
    q, p = leapfrog(dH_dq, dH_dp, q0, p0, dt, n_steps)

    crossings_q = []
    crossings_p = []

    for i in range(1, len(q)):
        q_prev = q[i-1, section_idx] - section_val
        q_curr = q[i, section_idx] - section_val

        if q_prev < 0 and q_curr >= 0:  # positive crossing
            # Linear interpolation
            frac = -q_prev / (q_curr - q_prev)
            q_cross = q[i-1] + frac * (q[i] - q[i-1])
            p_cross = p[i-1] + frac * (p[i] - p[i-1])

            # Record other phase space coordinates
            other_idx = [j for j in range(len(q0)) if j != section_idx]
            crossings_q.append(q_cross[other_idx])
            crossings_p.append(p_cross[other_idx])

    return np.array(crossings_q), np.array(crossings_p)

Method	Order	Energy Error (1000 orbits)	Cost/Step	Best For
Euler (non-symplectic)	1	O(1) — drifts!	1×	Never use for Hamiltonian
Leapfrog/Verlet	2	O(dt²) bounded	1×	Default choice
Yoshida 4th order	4	O(dt⁴) bounded	3×	High accuracy needed
Ruth 3rd order	3	O(dt³) bounded	3×	Rarely used
RK4 (non-symplectic)	4	O(dt⁴) but DRIFTS	4×	Short integrations only

Symptom	Fix
Energy drifts monotonically	You're using a non-symplectic method — switch to leapfrog
Close encounters cause blowup	Soften potential: 1/r → 1/√(r² + ε²)
Need adaptive timestep	Use time-transformed leapfrog (Mikkola & Aarseth)
Phase space structure smeared	Reduce dt (symplectic structure preserved at any dt, but accuracy improves)