Numerical Integration for Physics — Euler, Verlet & RK4

Physics simulations evolve state over time by approximating continuous equations with discrete time steps. The choice of integrator determines stability, energy conservation, and accuracy. This guide compares four integrators on a spring-mass system, from the simplest forward Euler to the gold-standard RK4.

1The problem: integrating ODEs

A physics body has state: position x and velocity v. Newton's second law gives acceleration: a = F/m. Integration means advancing the state forward by a small time step dt:

          // State: position x, velocity v // Acceleration: a = F(x, v) / m //
          Spring-mass: F = -k * x - b * v (spring + damping) const k = 50; //
          spring constant const b = 0.5; // damping coefficient const m = 1; //
          mass function acceleration(x, v) { return (-k * x - b * v) / m; } //
          Initial conditions let x = 1.0; // stretched 1 unit from equilibrium
          let v = 0.0; // at rest const dt = 0.016; // 16ms / frame ≈ 60fps
        

2Forward (Explicit) Euler — simple but unstable

Uses current derivative to extrapolate state. Fast and simple, but adds energy to the system — springs eventually explode at large dt.

          function stepEuler(x, v, dt) { const a = acceleration(x, v); const
          xNew = x + v * dt; // position uses OLD velocity const vNew = v + a *
          dt; // velocity uses OLD acceleration return [xNew, vNew]; } // Energy
          drifts upward over time — amplitude increases // STABLE only when: dt
          < 2 / sqrt(k/m) ≈ 0.28s for k=50, m=1
        

Forward Euler is first-order accurate (error ∝ dt). For a spring with k=50 and m=1, the critical timestep is ~0.28s. A 60fps game loop (dt≈0.016s) is safely inside this limit, but large springs at low framerates will blow up.

3Symplectic Euler — energy conserving

One small change: update velocity first using current force, then update position using the new velocity. This makes the integrator symplectic — it conserves a modified Hamiltonian and oscillators don't gain energy.

          function stepSymplecticEuler(x, v, dt) { const a = acceleration(x, v);
          const vNew = v + a * dt; // velocity first (NEW) const xNew = x + vNew
          * dt; // position uses NEW velocity return [xNew, vNew]; } // Energy
          is bounded — the spring oscillates stably forever // This is what most
          game physics engines (Bullet, Box2D) use internally // Same cost as
          forward Euler — just two lines swapped!
        

4Verlet integration — position-based

Velocity Verlet stores the previous position rather than explicit velocity, giving second-order accuracy and good energy conservation. Popular in cloth, molecular dynamics (GROMACS), and position-based dynamics.

          // Velocity Verlet (leapfrog variant) function stepVerlet(x, v, dt) {
          const a = acceleration(x, v); const xNew = x + v * dt + 0.5 * a * dt *
          dt; const aNew = acceleration(xNew, v); // acceleration at new
          position const vNew = v + 0.5 * (a + aNew) * dt; // average
          acceleration return [xNew, vNew]; } // Classic Verlet (no explicit
          velocity — uses x_prev): let xPrev = x - v * dt; // initialise x_prev
          function stepClassicVerlet(x, xPrev, dt) { const a = acceleration(x,
          0); // velocity-independent forces only const xNew = 2 * x - xPrev + a
          * dt * dt; const vApprox = (xNew - xPrev) / (2 * dt); // central
          difference return [xNew, x, vApprox]; }
        

Classic Verlet cannot use velocity-dependent forces (like damping) directly — use Velocity Verlet or the leapfrog variant for damped systems.

5Runge-Kutta 4 — the gold standard

RK4 computes four slope estimates per step and takes a weighted average. Fourth-order accuracy means halving dt reduces error 16×. The cost is four force evaluations per step.

          function stepRK4(x, v, dt) { // Each k = [dx, dv] derivative pair
          const k1 = deriv(x, v ); const k2 = deriv(x + k1[0]*dt/2, v +
          k1[1]*dt/2); const k3 = deriv(x + k2[0]*dt/2, v + k2[1]*dt/2); const
          k4 = deriv(x + k3[0]*dt, v + k3[1]*dt ); const xNew = x + (dt/6) *
          (k1[0] + 2*k2[0] + 2*k3[0] + k4[0]); const vNew = v + (dt/6) * (k1[1]
          + 2*k2[1] + 2*k3[1] + k4[1]); return [xNew, vNew]; } // Derivative
          function [ẋ = v, v̇ = a(x,v)] function deriv(x, v) { return [v,
          acceleration(x, v)]; } // RK4 is 4× more expensive than Euler but
          allows much larger dt // while maintaining accuracy — ideal for
          orbital mechanics and stiff systems
        

6Choosing an integrator and fixed timestep

Integrator	Order	Energy	Eval/step	Use case
Forward Euler	1st	Gains	1	Teaching only
Symplectic Euler	1st	Conserves	1	Games, particles
Velocity Verlet	2nd	Conserves	2	Cloth, bodies
RK4	4th	~Conserves	4	Orbital, precision

Always use a fixed timestep with an accumulator to decouple rendering from physics:

          const FIXED_DT = 1 / 120; // 120 Hz physics let accumulator = 0; let
          prevTime = performance.now() / 1000; function loop(nowMs) { const now
          = nowMs / 1000; const frameTime = Math.min(now - prevTime, 0.1); //
          cap at 100ms (tab unfocus) prevTime = now; accumulator += frameTime;
          while (accumulator >= FIXED_DT) { [x, v] = stepSymplecticEuler(x, v,
          FIXED_DT); accumulator -= FIXED_DT; } // Interpolate for smooth
          rendering (alpha = leftover fraction) const alpha = accumulator /
          FIXED_DT; const xRender = x; // or lerp(xPrev, x, alpha) for sub-frame
          smoothness mesh.position.x = xRender; requestAnimationFrame(loop); }
        

The accumulator pattern is the industry standard (Gaffer On Games: "Fix Your Timestep!"). Without it, physics runs differently at 30fps vs 144fps and simulation outcomes change with display refresh rate.