Energy Drift Playground — Simple Pendulum

The Foundation of Computational Physics

In the realm of computational physics, one fundamental challenge towers above all others: how do we preserve the essential properties of physical systems when translating continuous mathematics into discrete computations? This question lies at the heart of every physics simulation, from climate models predicting global warming to spacecraft navigation systems guiding interplanetary missions.

Energy drift represents perhaps the most critical manifestation of numerical error in conservative systems. When simulating a frictionless pendulum, total mechanical energy must remain constant for all time—this is not merely a numerical preference, but a fundamental law of physics. Yet many integration schemes systematically violate this conservation, leading to completely unphysical behavior that can invalidate years of computational work.

What if we could visualize this phenomenon in real-time? What if we could feel how different mathematical approaches either preserve or destroy the fundamental structure of physical reality?

The Mathematical Foundation of Conservative Systems

Hamiltonian Mechanics and Phase Space

A pendulum with length \(L\) and mass \(m\) under gravity \(g\) is governed by the Hamiltonian:

where \(q = \theta\) is the angular position and \(p = mL^2\dot{\theta}\) is the canonical momentum.

Hamilton’s equations provide the time evolution:

This formulation reveals the symplectic structure underlying all Hamiltonian systems—a geometric property that governs how phase space volumes evolve over time.

Liouville’s Theorem and Phase Space Conservation

Liouville’s theorem states that Hamiltonian flow preserves phase space volume:

This fundamental result has profound implications for numerical integration: any integration scheme that preserves phase space volume will maintain bounded energy errors over arbitrarily long time scales.

The Symplectic Condition

A transformation \((q_n, p_n) \rightarrow (q_{n+1}, p_{n+1})\) is symplectic if it preserves the symplectic form:

Equivalently, if we write the transformation as \(\mathbf{z}_{n+1} = \mathbf{M}\mathbf{z}_n\) where \(\mathbf{z} = (q, p)^T\), then \(\mathbf{M}\) must satisfy:

where \(\mathbf{J} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}\) is the symplectic matrix.

Visualizing Phase Space: The Geometric Perspective

Before exploring individual integration methods, we must understand how they behave in phase space—the natural coordinate system for Hamiltonian dynamics.

Phase Space Trajectories and Conservation

In phase space, each point represents a complete state $(\theta, \omega)$ of the pendulum. As time evolves, the system traces out a trajectory in this space. For the conservative pendulum, these trajectories have remarkable properties:

Energy Conservation as Geometric Constraint: Each trajectory lies on a curve of constant energy: \(H(\theta, \omega) = \frac{1}{2}\omega^2 + \frac{g}{L}(1 - \cos\theta) = E_0\)

Three Regimes of Motion:

• Small oscillations (\(E_0 < 2g/L\)): Closed elliptical orbits around \((\theta, \omega) = (0, 0)\)
• Large oscillations (\(E_0 > 2g/L\)): Open trajectories representing continuous rotation
• Separatrix (\(E_0 = 2g/L\)): The critical boundary between oscillation and rotation

The Symplectic Advantage: Geometric Structure Preservation

The key insight is that symplectic methods preserve the geometric structure of phase space, while non-symplectic methods systematically distort it:

Explicit Euler: Non-symplectic transformation causes trajectories to spiral outward, violating energy conservation.

Symplectic Euler: Despite first-order accuracy, preserves the topological structure of phase space, keeping trajectories on correct energy surfaces.

Velocity Verlet: Combines symplectic structure preservation with second-order accuracy, producing nearly perfect phase space trajectories.

RK4: High local accuracy but lacks symplectic structure, leading to slow drift off energy surfaces.

The Critical Challenge: Long-Term Stability vs. Local Accuracy

The naive approach to numerical integration prioritizes local truncation error—how well each individual step approximates the true solution. However, for conservative systems, this focus can be catastrophically misguided.

Consider two fundamental questions:

1. Should we choose a method with fourth-order local accuracy that allows energy to drift systematically?
2. Or should we prefer a first-order method that maintains perfect energy conservation structure?

For short simulations, local accuracy dominates. For long-term evolution—climate modeling, orbital mechanics, molecular dynamics—geometric structure preservation becomes paramount.

Interactive Laboratory: Four Integration Paradigms

The simple pendulum provides an ideal testbed because:

• Conservative system: Energy conservation is mathematically guaranteed
• Nonlinear dynamics: Small errors can compound exponentially
• Well-understood physics: We know the ground truth behavior
• Rich phase space: Exhibits libration, circulation, and separatrix dynamics

Method 1: Explicit Euler - The Unstable Baseline

The Explicit Euler method represents the most straightforward approach to numerical integration:

Algorithm: \(\omega_{n+1} = \omega_n + h \cdot f(\theta_n, \omega_n)\) \(\theta_{n+1} = \theta_n + h \cdot \omega_n\)

where \(f(\theta, \omega) = -\frac{g}{L}\sin(\theta)\) is the angular acceleration.

Mathematical Analysis:

• Order of accuracy: \(O(h)\)
• Stability: Conditionally stable with severe restrictions on \(h\)
• Energy behavior: Systematic energy growth due to the use of “stale” velocity information
• Phase space: Non-symplectic transformation that violates Liouville’s theorem

The Fundamental Flaw: Explicit Euler uses the velocity at time \(t_n\) to update position to time \(t_{n+1}\), while using the acceleration at time \(t_n\) to update velocity to time \(t_{n+1}\). This temporal mismatch systematically overestimates kinetic energy, leading to the characteristic exponential energy growth.

Method 2: Symplectic Euler - The Geometric Insight

Symplectic Euler makes a subtle but revolutionary change: update momentum first, then use the updated momentum to advance position.

Algorithm: \(p_{n+1} = p_n + h \cdot F(q_n)\) \(q_{n+1} = q_n + h \cdot \frac{p_{n+1}}{m}\)

For the pendulum: \(\omega_{n+1} = \omega_n + h \cdot \left(-\frac{g}{L}\sin(\theta_n)\right)\), then \(\theta_{n+1} = \theta_n + h \cdot \omega_{n+1}\)

Symplectic Structure Preservation: For the linearized harmonic oscillator, the transformation matrix for symplectic Euler is: \(\mathbf{M} = \begin{pmatrix} 1 - h^2 \omega_0^2 & h \\ - h \omega_0^2 & 1 \end{pmatrix}\)

One can verify that \(\det(\mathbf{M}) = 1\), ensuring the method is symplectic (area-preserving in phase space).

Key Properties:

• Order of accuracy: Still \(O(h)\) locally
• Global behavior: Energy oscillates with bounded amplitude
• Reversibility: Time-reversible integration
• Long-term stability: No secular drift in energy

The Geometric Miracle: Despite having the same local accuracy as Explicit Euler, Symplectic Euler captures the essential geometric structure of Hamiltonian flow. This leads to fundamentally different long-term behavior—energy remains bounded rather than growing exponentially.

Method 3: Velocity Verlet - The Optimal Balance

Velocity Verlet (also known as Leapfrog) combines symplectic structure with higher-order accuracy through a three-stage process:

Algorithm: \(v_{n+1/2} = v_n + \frac{h}{2} a(x_n)\) \(x_{n+1} = x_n + h \cdot v_{n+1/2}\) \(v_{n+1} = v_{n+1/2} + \frac{h}{2} a(x_{n+1})\)

For the pendulum: \(\omega_{n+1/2} = \omega_n + \frac{h}{2} \left(-\frac{g}{L}\sin(\theta_n)\right)\) \(\theta_{n+1} = \theta_n + h \cdot \omega_{n+1/2}\) \(\omega_{n+1} = \omega_{n+1/2} + \frac{h}{2} \left(-\frac{g}{L}\sin(\theta_{n+1})\right)\)

Superior Properties:

• Order of accuracy: \(O(h^2)\) - second-order accurate
• Symplectic: Preserves phase space structure exactly
• Time-reversible: \(\mathbf{M}^{-1} = \mathbf{M}^T\)
• Stability: Excellent long-term energy conservation

Computational Cost Analysis:

• Function evaluations: 2 per timestep (vs. 1 for Euler methods)
• Memory overhead: Minimal (stores half-step velocity)
• Efficiency: Optimal balance of accuracy and computational cost

Why Verlet Dominates Molecular Dynamics: The method’s symplectic nature combined with second-order accuracy makes it the gold standard for molecular dynamics simulations, where energy conservation over millions of timesteps is crucial.

Method 4: Runge-Kutta 4th Order - The Accuracy Champion

RK4 achieves exceptional local accuracy through four intermediate evaluations per timestep:

Algorithm: \(k_1 = f(t_n, y_n)\) \(k_2 = f(t_n + h/2, y_n + hk_1/2)\) \(k_3 = f(t_n + h/2, y_n + hk_2/2)\) \(k_4 = f(t_n + h, y_n + hk_3)\) \(y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)\)

Theoretical Foundation: RK4 achieves fourth-order accuracy by matching the first four terms of the Taylor series expansion. The method can be derived through the theory of trees (Butcher trees) which provide a systematic approach to constructing high-order Runge-Kutta methods.

The Symplectic Trade-off: While RK4 excels in local accuracy, it is not symplectic. For dissipative systems (with damping), this is often acceptable or even desirable. For conservative systems, the lack of geometric structure preservation can lead to subtle but systematic energy drift over long time scales.

Performance Characteristics:

• Order of accuracy: \(O(h^4)\) - fourth-order accurate
• Function evaluations: 4 per timestep
• Memory overhead: Minimal temporary storage
• Computational cost: 4× that of simple methods per timestep

Real-Time Method Comparison: Side-by-Side Analysis

Experience all four integration methods simultaneously to see their fundamental differences in action.

Comparative Analysis and Stability Theory

Computational Complexity Analysis

The Fundamental Trade-off: Accuracy vs. Structure

For conservative Hamiltonian systems, we face a fundamental choice:

Local Accuracy Priority (RK4, higher-order methods):

• Superior accuracy for individual timesteps
• Excellent for short-term, high-precision calculations
• May exhibit energy drift over extended simulations
• Ideal for dissipative systems where energy conservation is not required

Geometric Structure Priority (Symplectic methods):

• Preserve essential physical invariants
• Maintain bounded energy errors indefinitely
• Provide qualitatively correct long-term behavior
• Essential for conservative system simulation

Stability Analysis Through Linear Theory

For the linearized pendulum (\(\sin\theta \approx \theta\)), the equation becomes:

where \(\omega_0 = \sqrt{g/L}\) is the natural frequency.

Explicit Euler Stability: The amplification matrix eigenvalues are: \(\lambda_{\pm} = 1 \pm i h \omega_0\)

Their magnitudes are: \(|\lambda_{\pm}| = \sqrt{1 + (h\omega_0)^2} > 1 \quad \text{for any } h>0\)

Symplectic Euler Stability: The eigenvalues lie exactly on the unit circle: \(|\lambda| = 1\) for all \(h\). This guarantees bounded solutions for arbitrary timestep sizes.

Interactive Comparison: Complete Behavioral Analysis

Real-World Applications and Case Studies

Case Study 1: Climate Modeling - The ECMWF Experience

The European Centre for Medium-Range Weather Forecasts (ECMWF) faced a critical challenge: their atmospheric models were exhibiting systematic energy drift over long-term climate runs.

The Problem:

• Simulation duration: 100+ year climate projections
• Energy drift: 0.1% per decade using standard methods
• Consequence: Artificial warming trends corrupting climate predictions

The Solution: Implementation of energy-conserving finite difference schemes based on symplectic principles:

! Symplectic time stepping for atmospheric dynamics
do k = 1, nlevels
    ! Update momentum first (symplectic)
    u_new(k) = u_old(k) + dt * forcing_u(u_old, v_old, k)
    v_new(k) = v_old(k) + dt * forcing_v(u_old, v_old, k)
    
    ! Update positions using new momentum
    x(k) = x(k) + dt * u_new(k)
    y(k) = y(k) + dt * v_new(k)
end do

Results:

• Energy drift reduced to < 0.01% per century
• Improved long-term climate stability
• More reliable temperature trend predictions

Case Study 2: Molecular Dynamics - The GROMACS Success Story

GROMACS, one of the most widely used molecular dynamics packages, owes its success largely to the adoption of Velocity Verlet integration.

System Characteristics:

• Particle count: 10⁶ - 10⁹ atoms
• Simulation time: nanoseconds to microseconds
• Timestep: 1-2 femtoseconds
• Total steps: 10⁹ - 10¹² integration steps

Why Symplectic Methods Are Essential:

// Velocity Verlet implementation in GROMACS
void update_positions_and_velocities(int n, real dt,
                                   rvec x[], rvec v[], rvec f[], real mass[]) {
    for (int i = 0; i < n; i++) {
        // Half-step velocity update
        for (int d = 0; d < DIM; d++) {
            v[i][d] += (dt * 0.5) * f[i][d] / mass[i];
        }
        
        // Full-step position update
        for (int d = 0; d < DIM; d++) {
            x[i][d] += dt * v[i][d];
        }
    }
    
    // Compute new forces at updated positions
    compute_forces(x, f);
    
    // Complete velocity update
    for (int i = 0; i < n; i++) {
        for (int d = 0; d < DIM; d++) {
            v[i][d] += (dt * 0.5) * f[i][d] / mass[i];
        }
    }
}

Critical Performance Metrics:

• Energy conservation: < 0.01% drift over microsecond simulations
• Temperature stability: Canonical ensemble maintained without thermostats
• Computational efficiency: Optimal balance of accuracy and speed

Case Study 3: Spacecraft Navigation - The New Horizons Mission

NASA’s New Horizons mission to Pluto required extraordinary precision in trajectory calculation over a 9-year journey.

Navigation Challenges:

• Journey duration: 9 years (2.8 × 10⁸ seconds)
• Precision requirement: < 3000 km accuracy at Pluto encounter
• Perturbations: Gravitational effects from multiple bodies, solar radiation pressure

Integration Strategy: NASA used a combination of symplectic methods for long-term stability with adaptive high-order methods for critical maneuvers:

def spacecraft_propagator(state, dt, adaptive=False):
    if adaptive or critical_maneuver:
        # High-precision RK4 with adaptive timestep
        return rk4_adaptive(state, dt, tolerance=1e-12)
    else:
        # Symplectic integration for cruise phase
        return symplectic_euler(state, dt)

Results:

• Pluto encounter within 3000 km of predicted position
• Total trajectory error < 0.01% after 9-year journey
• Successful flyby of Arrokoth (2014 MU69) in 2019

Advanced Topics in Numerical Integration

Higher-Order Symplectic Methods

While Symplectic Euler and Velocity Verlet are first and second-order accurate respectively, higher-order symplectic methods can be constructed using composition methods.

Forest-Ruth Algorithm (4th Order Symplectic): \(\mathbf{S}_4 = \mathbf{S}_h^{\theta_1} \circ \mathbf{S}_h^{\theta_2} \circ \mathbf{S}_h^{\theta_3} \circ \mathbf{S}_h^{\theta_2} \circ \mathbf{S}_h^{\theta_1}\)

where \(\mathbf{S}_h\) is the basic symplectic map and:

• \[\theta_1 = \theta_3 = \frac{1}{2 - 2^{1/3}}\]
• \[\theta_2 = -\frac{2^{1/3}}{2 - 2^{1/3}}\]

Trade-offs:

• Pros: Fourth-order accuracy with exact energy conservation
• Cons: 5× computational cost per timestep vs. Velocity Verlet

Adaptive Time Stepping Strategies

For systems with multiple timescales, adaptive time stepping can dramatically improve efficiency:

Richardson Extrapolation:

def adaptive_symplectic_step(state, h_initial, tolerance):
    # Take one step of size h
    state_h = symplectic_step(state, h_initial)
    
    # Take two steps of size h/2
    state_temp = symplectic_step(state, h_initial/2)
    state_h2 = symplectic_step(state_temp, h_initial/2)
    
    # Estimate error
    error = norm(state_h - state_h2)
    
    if error < tolerance:
        # Accept step with Richardson extrapolation
        return (4*state_h2 - state_h) / 3, h_initial
    else:
        # Reduce timestep and retry
        return adaptive_symplectic_step(state, h_initial/2, tolerance)

Multiple Time Scale Methods

Systems like planetary motion (fast orbital motion + slow precession) benefit from multiple timestep methods:

Impulse-Velocity-Impulse (IVI) Decomposition:

def multiple_timestep_integrator(state, dt_fast, dt_slow, n_substeps):
    # Slow force impulse (half step)
    state = apply_slow_forces(state, dt_slow/2)
    
    # Fast evolution with multiple substeps
    dt_sub = dt_slow / n_substeps
    for i in range(n_substeps):
        state = symplectic_step_fast_forces(state, dt_sub)
    
    # Slow force impulse (half step)
    state = apply_slow_forces(state, dt_slow/2)
    
    return state

Implementation Best Practices

1. Timestep Selection Guidelines

Conservative Systems (Hamiltonian):

• Symplectic methods: \(h \leq \frac{T}{20}\) where \(T\) is the shortest characteristic period
• RK4: \(h \leq \frac{T}{100}\) for comparable accuracy to symplectic methods

Dissipative Systems:

• RK4 preferred: Higher local accuracy compensates for energy dissipation
• Adaptive methods: Essential for stiff systems with multiple timescales

2. Energy Conservation Monitoring

class EnergyMonitor:
    def __init__(self, initial_energy, tolerance=1e-6):
        self.E0 = initial_energy
        self.tolerance = tolerance
        self.max_error = 0
        
    def check_energy(self, current_energy, time):
        relative_error = abs(current_energy - self.E0) / abs(self.E0)
        self.max_error = max(self.max_error, relative_error)
        
        if relative_error > self.tolerance:
            warnings.warn(f"Energy error {relative_error:.2e} exceeds tolerance at t={time}")
            
        return relative_error

3. Numerical Precision Considerations

Double vs. Extended Precision:

• Standard double precision (64-bit): ~15 decimal digits
• For long-term simulations: Consider extended precision (80-bit) or arbitrary precision arithmetic
• GPU implementations: Be aware of reduced precision in single-precision calculations

4. Symplectic Integrator Verification

def verify_symplectic_property(integrator, state0, dt, n_steps=1000):
    """Verify that phase space volume is conserved"""
    # Create small perturbation in phase space
    epsilon = 1e-8
    states = [
        state0,
        state0 + [epsilon, 0],
        state0 + [0, epsilon],
        state0 + [epsilon, epsilon]
    ]
    
    # Evolve all perturbed states
    final_states = []
    for state in states:
        current = state.copy()
        for _ in range(n_steps):
            current = integrator(current, dt)
        final_states.append(current)
    
    # Compute initial and final phase space volumes
    initial_volume = compute_phase_space_volume(states)
    final_volume = compute_phase_space_volume(final_states)
    
    volume_change = abs(final_volume - initial_volume) / initial_volume
    
    print(f"Phase space volume change: {volume_change:.2e}")
    return volume_change < 1e-10  # Should be essentially zero for symplectic methods

Troubleshooting Common Issues

Problem: Energy Growth in Symplectic Methods

Symptoms: Energy increases systematically despite using symplectic integrator

Potential Causes:

1. Implementation error: Check that momentum is updated before position
2. Floating-point precision: Accumulation of round-off errors over millions of steps
3. Timestep too large: Linear stability analysis assumes small timesteps

Solutions:

# Correct symplectic Euler implementation
def symplectic_euler_correct(theta, omega, dt, g, L):
    # CORRECT: Update omega first
    omega_new = omega - (dt * g / L) * sin(theta)
    theta_new = theta + dt * omega_new  # Use NEW omega
    return theta_new, omega_new

# INCORRECT implementation (produces energy growth)
def symplectic_euler_wrong(theta, omega, dt, g, L):
    # WRONG: Updates done simultaneously
    omega_new = omega - (dt * g / L) * sin(theta)
    theta_new = theta + dt * omega  # Uses OLD omega
    return theta_new, omega_new

Problem: Stiff System Instability

Symptoms: Simulation becomes unstable for small timesteps

Cause: Stiff differential equations require implicit methods

Solution: Use implicit-explicit (IMEX) methods:

def imex_euler(q, p, dt, force_explicit, force_implicit):
    # Explicit update for non-stiff terms
    p_temp = p + dt * force_explicit(q, p)
    
    # Implicit update for stiff terms (requires solver)
    p_new = solve_implicit(p_temp, dt, force_implicit, q)
    q_new = q + dt * p_new
    
    return q_new, p_new

Problem: Resonance-Induced Instability

Symptoms: Energy oscillates with growing amplitude over time

Cause: Timestep resonates with system’s natural frequency

Solution: Use irrational timestep ratios or variable timestep methods:

# Golden ratio timestep to avoid resonances
dt_base = 0.01
golden_ratio = (1 + sqrt(5)) / 2
dt_adaptive = dt_base / golden_ratio

Future Directions and Research Frontiers

Machine Learning Enhanced Integration

Recent developments combine traditional numerical methods with machine learning:

Physics-Informed Neural Networks (PINNs):

import torch
import torch.nn as nn

class HamiltonianNet(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(2, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 1)
        )
    
    def forward(self, q, p):
        return self.net(torch.cat([q, p], dim=1))
    
    def symplectic_gradient(self, q, p):
        H = self.forward(q, p)
        dH_dq = torch.autograd.grad(H.sum(), q, create_graph=True)[0]
        dH_dp = torch.autograd.grad(H.sum(), p, create_graph=True)[0]
        return -dH_dq, dH_dp  # Hamilton's equations

Quantum-Classical Hybrid Systems

Emerging applications require integration of quantum and classical dynamics:

Mean Field Theory Integration: For systems where quantum effects influence classical motion through expectation values.

Geometric Deep Learning

Using the geometric structure of phase space in neural network architectures:

Symplectic Neural Networks: Networks constrained to preserve symplectic structure exactly while learning complex dynamics.

Conclusion: The Art and Science of Numerical Integration

The exploration of energy drift in the simple pendulum reveals fundamental principles that extend far beyond this seemingly elementary system. We have witnessed how the choice of integration method can mean the difference between simulation success and catastrophic failure—between climate models that accurately predict future temperatures and those that exhibit spurious warming trends.

Key Insights from Our Journey

1. Geometric Structure Trumps Local Accuracy: Symplectic methods with first-order accuracy often outperform fourth-order methods that ignore geometric structure.

2. Conservation Laws Are Not Optional: For conservative systems, preserving invariants like energy, momentum, and phase space volume is essential for physical realism.

3. The Mathematics Has Real Consequences: The abstract concepts of symplectic geometry and Hamiltonian mechanics translate directly into practical simulation reliability.

4. Understanding Enables Optimization: Knowledge of stability theory, convergence analysis, and computational complexity allows informed selection of methods for specific applications.

The Broader Impact

These principles extend throughout computational science:

• Climate Science: Energy-conserving atmospheric models
• Astronomy: Long-term planetary motion calculations
• Materials Science: Molecular dynamics simulations
• Engineering: Structural dynamics and vibration analysis
• Biology: Protein folding and cellular dynamics
• Finance: Long-term economic modeling

Looking Forward

As computational power continues to grow exponentially, the temptation exists to simply use smaller timesteps with traditional methods. However, this brute-force approach fails for systems requiring integration over geological timescales or astronomical distances. The future belongs to methods that respect the underlying mathematical structure of physical systems.

The simple pendulum, in its elegant simplicity, embodies the rich mathematical structure that governs all of physics. By understanding how to preserve this structure numerically, we gain the power to simulate reality with unprecedented fidelity and extend our understanding of complex systems across all scales of nature.

The mathematics of numerical integration isn’t merely computational technique—it’s the bridge between theoretical understanding and practical prediction, between the elegance of continuous mathematics and the reality of discrete computation.

Further Reading and References

Essential Textbooks

• Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Springer.
• Leimkuhler, B., & Reich, S. (2004). Simulating Hamiltonian Dynamics. Cambridge University Press.

Research Papers

• Forest, E., & Ruth, R. D. (1990). Fourth-order symplectic integration. Physica D, 43(1), 105-117.
• Yoshida, H. (1990). Construction of higher order symplectic integrators. Physics Letters A, 150(5), 262-268.

Computational Tools

• GROMACS: gromacs.org - Molecular dynamics with Velocity Verlet
• NAMD: namd.org - Large-scale parallel MD simulations
• yt-project: yt-project.org - Analysis tools for astrophysical simulations