ODE Solvers: Euler, RK4, and Leapfrog Compared

1. What Is an ODE?

An Ordinary Differential Equation (ODE) expresses a state's rate of change as a function of the state itself. Newton's second law is a classic example:

Newton's second law as an ODE system dv/dt = F(x, v, t) / m (velocity changes due to force) dx/dt = v (position changes due to velocity)

Given initial conditions x(0) and v(0), an ODE solver advances the system in small time steps Δt, approximating the exact solution. The error per step and the long-term behaviour depend heavily on the algorithm.

2. Euler's Method

The simplest possible scheme: take the derivative at the current point and step in that direction by Δt.

Forward (Explicit) Euler x(t + Δt) = x(t) + Δt \cdot ẋ(t) v(t + Δt) = v(t) + Δt \cdot a(t)

Order: 1st order (error ∝ Δt per step, ∝ Δt² global).
Cost: 1 force evaluation per step — very cheap.
Problem: The method artificially adds energy to conservative systems. Orbits spiral outward; springs oscillate with growing amplitude. The error accumulates even with tiny Δt.

Avoid Euler for long simulations: A planet orbiting with Euler integration will slowly spiral away from the Sun — not because physics is wrong, but because the integrator is non-conservative.

3. Midpoint & Heun Methods

A simple improvement: evaluate the derivative at an intermediate point. Heun's method (predictor–corrector):

Heun's method (RK2) k1 = f(t, y) k2 = f(t + Δt, y + Δt\cdotk1) (predicted endpoint) y(t + Δt) = y + (Δt/2)\cdot(k1 + k2)

Both are 2nd order — much more accurate than Euler with only 2× the cost. Still non-symplectic; energy still drifts over very long runs.

4. Runge–Kutta 4th Order (RK4)

The workhorse of scientific computing. RK4 evaluates four derivative estimates per step and takes their weighted average:

RK4 update k1 = f(t, y) k2 = f(t + Δt/2, y + (Δt/2)\cdotk1) k3 = f(t + Δt/2, y + (Δt/2)\cdotk2) k4 = f(t + Δt, y + Δt\cdotk3) y(t + Δt) = y + (Δt/6)\cdot(k1 + 2k2 + 2k3 + k4)

Order: 4th — global error ∝ Δt⁴. Extremely accurate for smooth problems.
Cost: 4 force evaluations per step — 4× more than Euler.
Problem: Also non-symplectic. Over millions of steps, energy drifts — slower than Euler, but it still drifts.

Rule of thumb: RK4 is ideal for trajectory calculations with a known end time (rocket manoeuvres, ODE-based animations). For indefinitely long simulations (molecular dynamics, long-term orbital mechanics), use Leapfrog or other symplectic solvers instead.

5. Leapfrog / Verlet Integration

The Leapfrog (Störmer–Verlet) integrator staggers position and velocity by half a time step — positions and velocities "leapfrog" past each other:

Leapfrog (velocity Verlet form) v(t + Δt/2) = v(t) + (Δt/2)\cdota(t) x(t + Δt) = x(t) + Δt \cdot v(t + Δt/2) a(t + Δt) = F(x(t + Δt)) / m (re-evaluate force) v(t + Δt) = v(t + Δt/2) + (Δt/2)\cdota(t + Δt)

Order: 2nd order — same as Heun, but with half the cost (1 force evaluation per step).
Symplectic: Yes — the total energy oscillates around the correct value but does not drift. Orbits stay closed forever at any fixed Δt.
Used in: All major molecular dynamics packages (LAMMPS, GROMACS), n-body codes, celestial mechanics.

Why is leapfrog symplectic? It preserves the volume of phase space (Liouville's theorem). This is the discrete analogue of Hamilton's equations, which also conserve phase-space volume.

6. Symplectic vs. Non-Symplectic

A symplectic integrator exactly conserves a slightly perturbed Hamiltonian H' (close to the true H). This prevents systematic energy growth or decay, making the simulation stable for arbitrarily long times.

Non-symplectic methods (Euler, RK4) do not satisfy this property. Their energy errors grow without bound, eventually ruining long simulations even if the per-step error is tiny.

Leapfrog/Verlet — 2nd order symplectic
Ruth's 3rd order — 3rd order symplectic
Yoshida 4th order — 4th order symplectic (combines 3 leapfrog sub-steps)
PEFRL — modern 4th order, often preferred in MD simulations

7. Comparison Table

Method	Order	Evals/step	Symplectic	Energy long-term	Best for
Euler (forward)	1st	1	No	Grows	Demos, learning only
Symplectic Euler	1st	1	Yes	Stable	Fast games, many particles
Midpoint / Heun	2nd	2	No	Slow drift	Short animations
RK4	4th	4	No	Very slow drift	Precise trajectories
Leapfrog / Verlet	2nd	1	Yes	Oscillates ≈ const	Orbital, molecular dynamics
Yoshida 4th	4th	3	Yes	Oscillates ≈ const	High-accuracy long runs

8. JavaScript Implementation

All three integrators for a simple harmonic oscillator a = −ω²·x (spring without damping):

const omega = 1.0;                  // angular frequency
const dt    = 0.05;                  // time step (seconds)

function accel(x) { return -omega * omega * x; }

// ── Forward Euler ───────────────────────────────────
function stepEuler(x, v) {
  const a = accel(x);
  return { x: x + dt * v,  v: v + dt * a };
}

// ── RK4 ────────────────────────────────────────────
function stepRK4(x, v) {
  const k1x = v,                     k1v = accel(x);
  const k2x = v + 0.5*dt*k1v,       k2v = accel(x + 0.5*dt*k1x);
  const k3x = v + 0.5*dt*k2v,       k3v = accel(x + 0.5*dt*k2x);
  const k4x = v + dt*k3v,            k4v = accel(x + dt*k3x);
  return {
    x: x + (dt/6) * (k1x + 2*k2x + 2*k3x + k4x),
    v: v + (dt/6) * (k1v + 2*k2v + 2*k3v + k4v)
  };
}

// ── Leapfrog (Velocity Verlet) ─────────────────────
function stepLeapfrog(x, v) {
  const vHalf = v   + 0.5 * dt * accel(x);   // half-step velocity
  const xNew  = x   + dt   * vHalf;             // full-step position
  const vNew  = vHalf + 0.5 * dt * accel(xNew); // half-step velocity again
  return { x: xNew, v: vNew };
}

// Energy of a harmonic oscillator (should be conserved)
function energy(x, v) { return 0.5*v*v + 0.5*omega*omega*x*x; }

Tip: Plug these into a simple render loop and plot energy(x, v) over time. You will see Euler's energy grow linearly, RK4's grow very slowly, and Leapfrog's oscillate tightly around the initial value.

9. Which Method Should You Use?

Game / real-time physics: Symplectic Euler (update velocity first, then position). One evaluation, symplectic, cheap.
N-body / molecular dynamics: Leapfrog/Verlet. Standard in the field for 50 years.
Precise short trajectory (rocket, projectile): RK4. Accuracy over long-term stability.
Long-run high-accuracy astronomy: Yoshida 4th order or higher-order symplectic.
Stiff equations (chemical kinetics, circuits): Implicit methods (backward Euler, SDIRK) — these are a separate topic.

Cannon-es & Rapier (the physics engines powering this site) use variations of symplectic Euler and impulse-based solvers. They sacrifice some accuracy in exchange for extreme stability and rigid-body collision resolution.

Explore the n-body and orbital simulations on this site — they use Leapfrog integration so that planets orbit stably for millions of simulated years without spiralling away.

🪐 Open Orbital Simulation →