🔬 Numerical Stability: Why Simulations Explode and How to Fix Them

Floating Point: Not Real Numbers

64-bit floating point numbers (IEEE 754 double precision) represent numbers with about 15–17 significant decimal digits of precision. They are not uniformly spaced: the gap between representable numbers near 1.0 is about 2.2 × 10⁻¹⁶ (machine epsilon), while near 10⁶ the gap is about 2.2 × 10⁻¹⁰. This is why order of operations matters numerically even for "exact" arithmetic.

The classic failure is catastrophic cancellation: subtracting two nearly-equal numbers loses significant digits. Consider computing (1 + x) - 1 for x = 1e-16. In exact arithmetic, this equals x. In floating point, 1 + 1e-16 rounds to exactly 1.0 (below machine epsilon), so the result is 0.0 — a 100% error. This bites in many physics formulas:

// Bad: catastrophic cancellation for small x
const bad = (1 - Math.cos(x)) / (x * x);

// Good: use the trigonometric identity
// (1 - cos x) = 2 sin²(x/2)
const good = 2 * Math.sin(x/2)**2 / (x * x);

Other floating point hazards: summing many numbers in a naive loop accumulates rounding errors O(n · ε). Kahan summation corrects this by maintaining a running compensation term, reducing accumulated error to O(ε) regardless of n. For physics simulations summing forces from thousands of particles, Kahan summation can make the difference between stable and drifting energy.

Condition Numbers: Ill-Posed Problems

The condition number of a problem measures how much the output changes relative to small changes in input. A condition number of 100 means a relative input error of 10⁻¹⁵ can produce a relative output error of 10⁻¹³ — still fine. A condition number of 10¹⁵ means floating-point precision alone makes the output meaningless.

Solving a linear system Ax = b has condition number κ(A) = ||A|| · ||A⁻¹||. Nearly-singular matrices (where two rows are almost parallel) have enormous condition numbers. This is why inverting a matrix directly is ill-advised for numerical computation — use LU decomposition, QR factorization, or iterative methods that exploit the problem structure instead.

For simulations, ill-conditioning often arises from extreme parameter ratios. A fluid simulation coupling fast acoustic waves (propagating at 340 m/s) with slow flow (0.1 m/s) has a Mach number of 0.0003 and the acoustic-to-flow timescale ratio of 3400. Explicit time-stepping must resolve the fast acoustic waves even when we only care about the slow flow — enormously wasteful. The condition number of the coupled system is proportional to this ratio.

Stiffness: When Explicit Methods Fail

A stiff ODE contains processes with vastly different timescales. Chemical kinetics is the classic example: a fast reaction with timescale 10⁻⁹ s and a slow reaction with timescale 10³ s are coupled. An explicit integrator must take timesteps < 10⁻⁹ s to remain stable — billions of steps to simulate a millisecond, even though the fast process reaches quasi-equilibrium almost immediately and contributes no useful dynamics at the slow timescale.

The stability region of an explicit method is a bounded region in the complex plane. For Euler's method, stability requires |1 + Δt·λ| ≤ 1 for all eigenvalues λ of the Jacobian. For a stiff system with λ = −10⁶, stability requires Δt ≤ 2/10⁶ = 2 μs. Explicit RK4 has a larger stability region but the same fundamental constraint.

Implicit methods include the future state in the derivative evaluation, making stability conditions much less restrictive. The backward Euler method y(t+Δt) = y(t) + Δt · f(y(t+Δt)) requires solving a nonlinear system at each step (using Newton's method) but is A-stable: stable for any Δt for problems with negative-real eigenvalues. This allows huge timesteps for stiff problems, often making implicit methods 10–1000× faster despite the higher cost per step.

Diagnosing an Unstable Simulation

When a simulation explodes, the root cause is almost always one of three things:

• Timestep too large: Check Δt · max(|eigenvalues|). If this exceeds the stability region of your integrator, reduce Δt or switch to an implicit method.
• Floating point catastrophe: Log intermediate values and look for sudden jumps or NaN propagation. Common culprits: dividing by a near-zero quantity (particle collision distances), square roots of negative numbers (floating-point noise making a should-be-zero quantity slightly negative), or taking logs of numbers that should be positive.
• Energy drift from non-symplectic integration: Monitor total energy over time. If it increases monotonically (not oscillating), you have systematic energy injection from your integrator. Switch to a symplectic method.

// Defensive simulation loop
function step(dt) {
  const newState = integrate(state, dt);
  if (!isFinite(newState.energy) || newState.energy > MAX_ENERGY) {
    console.warn('Instability detected, halving timestep');
    dt /= 2;
    return step(dt);  // retry with smaller step
  }
  state = newState;
}

Practical Stability Techniques

For fluid simulations, velocity field advection using semi-Lagrangian methods (tracing particle paths backward in time) is unconditionally stable for any timestep — at the cost of some numerical diffusion. The Courant-Friedrichs-Lewy (CFL) condition bounds the timestep for explicit advection: Δt ≤ CFL · Δx / max(|u|). With CFL = 0.9, a fluid cell crossed in one timestep remains stable; CFL > 1 means information propagates more than one cell per step and instability follows.

For spring-mass systems (cloth, soft bodies), the critical timestep scales as Δt_crit ~ √(m/k): stiffer springs (higher k) require smaller timesteps. Position-based dynamics (PBD) avoids this by projecting constraints after each step rather than integrating forces directly — trading physical accuracy for unconditional stability at any timestep, which is why it's widely used in game physics.