Numerical Methods Cheatsheet

🎯 Why Numerical Methods Matter

Every physics simulation — a bouncing ball, a cloth sheet, an orbiting planet — is governed by differential equations that almost never have a closed-form solution once collisions, non-linear forces, or many interacting bodies are involved. Numerical methods turn continuous equations into discrete update rules a computer can run every frame, trading exactness for speed and controllable error.

Choosing the right method is a trade-off between three axes: accuracy (how close to the true trajectory), stability (does error blow up over time), and cost (function evaluations per step). This page collects the methods actually used across the simulations on this site, with the practical guidance that textbooks often bury.

🌀 ODE Integrators

Simulation state (position, velocity) evolves according to a first-order ODE $\dot{\mathbf{y}} = f(t, \mathbf{y})$. An integrator advances $\mathbf{y}$ from $t$ to $t+h$ using evaluations of $f$.

Explicit (Forward) Euler

The simplest integrator: one function evaluation per step, trivial to implement, but it systematically gains energy in oscillatory systems (a pendulum simulated with Euler slowly spirals outward). Local error is $O(h^2)$, global error $O(h)$.

function eulerStep(state, force, h) {
  // state = { x, v }, force(x, v) returns acceleration
  const a = force(state.x, state.v);
  state.x += h * state.v;
  state.v += h * a;
  return state;
}

Semi-implicit (Symplectic) Euler

Swap the order — update velocity first, then use the new velocity to update position. Same cost as forward Euler, still only first-order accurate, but it is symplectic: energy oscillates around the true value instead of drifting away. This is the default integrator in almost every real-time game physics engine (Box2D, Cannon-es, Rapier).

function symplecticEulerStep(state, force, h) {
  const a = force(state.x, state.v);
  state.v += h * a;        // velocity first
  state.x += h * state.v;  // then position, using new v
  return state;
}

Velocity Verlet

The workhorse of particle systems, molecular dynamics, and cloth/rope solvers. Second-order accurate, symplectic (bounded energy error), and time-reversible. Costs two acceleration evaluations per step when acceleration depends on velocity; cheaper when it doesn't.

function verletStep(state, accel, h) {
  const a0 = accel(state.x);
  state.x += h * state.v + 0.5 * h * h * a0;
  const a1 = accel(state.x);
  state.v += 0.5 * h * (a0 + a1);
  return state;
}

Runge-Kutta 4 (RK4)

Four function evaluations per step buy fourth-order accuracy — global error $O(h^4)$ versus $O(h)$ for Euler. Excellent for smooth, non-collisional dynamics (orbital mechanics, camera paths) where you want high accuracy at a larger step size. RK4 is not symplectic: over very long integrations (millions of orbits) energy can still drift, just far more slowly than with Euler.

function rk4Step(y, t, h, f) {
  const k1 = f(t, y);
  const k2 = f(t + h / 2, add(y, scale(k1, h / 2)));
  const k3 = f(t + h / 2, add(y, scale(k2, h / 2)));
  const k4 = f(t + h, add(y, scale(k3, h)));
  return add(y, scale(
    add(add(k1, scale(k2, 2)), add(scale(k3, 2), k4)),
    h / 6
  ));
}

⚖ Stability & Step Size

A method is stable for a given step size $h$ if small errors do not grow unboundedly over many steps. For the test equation $\dot{y} = \lambda y$ (with $\lambda$ complex, $\text{Re}(\lambda) \le 0$), each explicit method has a bounded region of absolute stability in the complex plane $h\lambda$.

For a damped spring-mass system with stiffness $k$ and mass $m$, this translates to a practical rule of thumb:

Exceed this and a stiff spring (high k) simulated with an explicit integrator will oscillate with growing amplitude and eventually explode — the classic "jittery exploding cloth" bug. Three fixes, in order of preference:

• Fixed, small timestep — decouple physics update rate from render framerate (e.g. 60 Hz physics via an accumulator, even at 30 fps render).
• Substepping — run N smaller physics steps per render frame instead of one large one.
• Implicit integration — backward Euler or position-based dynamics (PBD) are unconditionally stable for linear stiff terms, at the cost of solving a linear (or iterative) system per step.

// Fixed-timestep accumulator loop (decouples physics from render rate)
let accumulator = 0;
const FIXED_DT = 1 / 60;

function frame(now) {
  accumulator += (now - lastTime) / 1000;
  lastTime = now;
  while (accumulator >= FIXED_DT) {
    physicsStep(FIXED_DT);
    accumulator -= FIXED_DT;
  }
  render(accumulator / FIXED_DT); // interpolation alpha
  requestAnimationFrame(frame);
}

🔍 Root Finding

Used for collision time-of-impact queries, inverse kinematics, and solving implicit constraint equations $g(\mathbf{x}) = 0$.

Bisection method

Requires $f(a)$ and $f(b)$ to have opposite signs. Always converges, halving the bracket each iteration — reliable but slow ($O(\log_2(\epsilon^{-1}))$ iterations for tolerance $\epsilon$).

Newton-Raphson method

Converges much faster than bisection when it converges, but can diverge or cycle if the initial guess is far from the root or $f'(x_n) \approx 0$. Used heavily in continuous collision detection to refine the exact time of impact.

function newtonRaphson(f, fPrime, x0, tol = 1e-6, maxIter = 50) {
  let x = x0;
  for (let i = 0; i < maxIter; i++) {
    const fx = f(x);
    if (Math.abs(fx) < tol) return x;
    const dfx = fPrime(x);
    if (Math.abs(dfx) < 1e-12) break; // derivative too flat, bail out
    x -= fx / dfx;
  }
  return x; // best effort — caller should validate
}

Bisection vs Newton — when to use which

🔢 Linear Systems

Solving $A\mathbf{x} = \mathbf{b}$ appears throughout simulation: constraint solvers (rigid body contacts, ropes), implicit integration, and inverse kinematics Jacobians.

Gauss-Seidel iteration

Physics engines (Box2D, Bullet, Cannon-es) solve contact and joint constraints with a handful of Gauss-Seidel (or projected Gauss-Seidel, for inequality constraints like non-penetration) iterations per frame rather than an exact solve — a few iterations of an approximate answer every frame beats one exact solve that misses the deadline.

// N Gauss-Seidel sweeps over a sparse constraint system
function solveConstraints(constraints, iterations = 10) {
  for (let iter = 0; iter < iterations; iter++) {
    for (const c of constraints) {
      c.solve(); // updates velocities in place, uses latest values
    }
  }
}

LU decomposition (dense, direct)

For small dense systems (e.g. a 3×3 or 4×4 inertia tensor solve), a direct LU decomposition is cheap and exact: factor $A = LU$ once, then solve two triangular systems per right-hand side via forward/back substitution — $O(n^3)$ to factor, $O(n^2)$ per solve.

📈 Interpolation

Linear interpolation (lerp)

Cubic Hermite / Catmull-Rom spline

Catmull-Rom passes exactly through every control point (unlike Bézier) and is the standard choice for camera paths and procedural roads/tracks in the simulations here.

Spherical linear interpolation (slerp)

Unlike component-wise lerp of quaternions, slerp interpolates along the great-circle arc so rotation speed stays constant — essential for smooth camera and character orientation blending.

∫ Numerical Differentiation & Integration

Central difference (derivative)

Used to numerically approximate a Jacobian or normal vector (e.g. gradient of a signed-distance field) when no analytic derivative is available. The central difference cancels the first-order error term that a naive forward difference $\frac{f(x+h)-f(x)}{h}$ carries.

Trapezoidal rule (integration)

Simpson's rule (integration)

Simpson's rule fits a parabola through every three points and is the cheapest big accuracy win over the trapezoidal rule for smooth integrands — useful when precomputing arc-length tables for splines (needed for constant-speed path following).

🌊 FFT Basics

The Fast Fourier Transform converts a signal from the time/space domain to the frequency domain in $O(n\log n)$ instead of the $O(n^2)$ of a naive Discrete Fourier Transform (DFT) — the backbone of ocean-wave simulation (Tessendorf/FFT ocean), audio-reactive visuals, and GPU-based fluid solvers.

The Cooley-Tukey FFT algorithm recursively splits the sum into even- and odd-indexed terms, halving the problem size each recursion level — hence $N$ must typically be a power of two for the classic radix-2 implementation used in most WebGL ocean shaders.

In practice, ocean and water simulations here run the FFT on a GPU compute shader (or a WebGL ping-pong render target pass) to transform a heightfield spectrum into displacement — a 256×256 grid costs roughly 16 passes of radix-2 butterflies per dimension.

📊 Error Order Quick Table

Halving the step size $h$ for a method of order $p$ reduces the error by roughly $2^p$. Use this to sanity-check a method's convergence empirically (Richardson extrapolation): halve $h$, re-run, and confirm the error drops by the expected factor.

⚠ Common Pitfalls

• Variable timestep with Euler — running forward Euler at whatever deltaTime the browser gives you makes energy drift depend on framerate, so the same simulation behaves differently on different machines.
• Catastrophic cancellation — subtracting two nearly-equal floating point numbers (e.g. b*b - 4*a*c in the quadratic formula when discriminant ≈ 0) amplifies rounding error; use the numerically stable form of the quadratic formula when precision matters.
• Ignoring stiffness — a single very stiff spring or constraint in an otherwise soft system forces the entire simulation to use a tiny step size (or go unstable); isolate stiff constraints and solve them implicitly.
• Wrong integrator for the job — using RK4 for colliding rigid bodies wastes cycles (it isn't symplectic and collisions violate its smoothness assumption anyway); using Euler for accurate orbital mechanics accumulates visible drift.
• Aliasing without a low-pass filter — sampling an FFT-based heightfield at too coarse a grid, or not windowing a signal, introduces high-frequency artifacts.
• Comparing floats for equality — always use an epsilon tolerance (Math.abs(a - b) < 1e-6) instead of a === b for any quantity derived from iterative numerical methods.