Stability Limits of Explicit Integration in Dynamics

A simulation advances state by a time step, and the size of that step is not chosen for accuracy alone. Above a threshold the simulation stops being an approximation and becomes a different, oscillating, or exploding system, even when the differential equation it models is perfectly stable. The threshold is a property of the integrator and the stiffest force in the scene, and it is the reason a cloth simulation can run at sixty steps per second while a stiff constraint needs hundreds or a different method entirely. Understanding the bound is the difference between tuning a solver and being surprised by it.

The test problem that reveals the bound

Every stability argument about an integrator reduces to the harmonic oscillator. A spring with stiffness k and unit mass obeys $\ddot{x} = -\omega^2 x$ with $\omega^2 = k$, and the exact solution is a bounded sinusoid of frequency $\omega$. An integrator that produces unbounded growth from this input is unstable, and the step size at which it first does so is the stability limit. The harmonic oscillator is the right test because every force in a scene is locally a spring near an equilibrium, so the stiffest spring sets the global limit.

Forward Euler integrates $\dot{x} = f(x)$ as $x_{n+1} = x_n + h f(x_n)$. Applied to the oscillator it produces a recurrence whose amplification per step has magnitude $\sqrt{1 + (h\omega)^2}$, which is greater than one for any positive step. Forward Euler therefore diverges on the oscillator at every step size, growing exponentially; the step controls the growth rate rather than eliminating it. That is the central fact about explicit integration: it cannot stay stable on a stiff spring, only on the limit where the spring is replaced by a constraint.

Semi-implicit Euler, which integrates velocity first and then position using the new velocity, does better. Its amplification magnitude is one for $h\omega \le 2$, so it neither grows nor shrinks within that bound and diverges beyond it. The bound $\omega h < 2$, equivalently $h < 2/\omega$, is the canonical stability limit, and it says the step size must be smaller than about a third of the oscillation period of the stiffest spring. A spring that oscillates at 100 rad/s demands a step below 20 milliseconds; a stiff contact treated as a spring oscillates far faster and demands a step far smaller.

Where the stiffness comes from

Stiffness is the ratio of the fastest to the slowest time scale in the system, and in a character or cloth scene the fast scale is a constraint enforced by a spring. A distance constraint between two particles with a spring of stiffness k has natural frequency $\omega = \sqrt{k/m}$ for particle mass m, and the simulator’s step is tied to the largest such frequency in the scene. A bone-length constraint on a skeleton, a stretch constraint on cloth, and a collision penalty all raise the maximum $\omega$ and lower the allowed step.

The symptom of exceeding the limit is specific and worth recognizing. Near the limit the energy oscillates within the bound and the system jitters; past it the energy grows per step and the particles fly apart, cloth inflates, and constraints that held the previous frame release explosively. The explosion is not a bug in the force model. It is the integrator amplifying the natural oscillation of the stiffest spring, and it can be triggered by shrinking a time step’s mass or raising a stiffness during tuning without changing anything else.

Mass interacts directly. The frequency is $\sqrt{k/m}$, so a light particle on a stiff spring is the worst case: halving the mass multiplies the frequency by $\sqrt2$ and moves the stability bound. This is why reducing mass to make cloth feel lighter can destabilize a simulation that was stable, and why solvers often clamp a minimum mass rather than let it approach zero, since a zero-mass particle has infinite frequency and no stable explicit step exists.

/* Stability of a damped harmonic oscillator under semi-implicit Euler.
   The bound is h*omega < 2 for zero damping; damping enlarges it slightly. */
#include <math.h>

static void damped_oscillator_step(double *x, double *v,
                                   double omega, double zeta, double h) {
    double a = -omega * omega * *x - 2.0 * zeta * omega * *v;
    *v += h * a;           /* velocity first (semi-implicit) */
    *x += h * *v;          /* position uses the new velocity */
}

/* Caller keeps h below the stable region. For zeta=0, h*omega must be < 2;
   at omega=100 rad/s that is h < 20 ms, one step of a 50 Hz simulator. */

Damping shifts the bound. A damped oscillator has amplification that depends on the damping ratio $\zeta$, and positive damping enlarges the stable region slightly, which is why adding damping to a marginally unstable solver can stabilize it. The cure is real but local: it changes the bound, not the mechanism, and a sufficiently stiff spring still overcomes it.

Verlet and the integrators that cheat the bound

Verlet integration stores positions and implicitly carries velocity through the previous position: $x_{n+1} = 2x_n - x_{n-1} + h^2 a_n$. It has the same stability structure as semi-implicit Euler, with the same $\omega h < 2$ bound, and it has one property that made it the default for cloth and particles: it is symplectic in the undamped, constant-step case, meaning it preserves a modified energy exactly rather than drifting in one direction. Energy does not systematically grow or shrink; it oscillates around a value, which keeps a simulation alive far longer than a method that pumps energy in.

The modified energy is not the true energy. Verlet’s preserved quantity differs from the physical energy by a term of order $h^2$, so a very long simulation can still drift, and the drift is not visible as a divergence but as a slow, unphysical change in the total energy that accumulates over minutes of simulated time. The fix is a small step or a symplectic integrator of higher order, and the drift is a different failure from instability: it is bounded per step but unbounded in total.

Velocity Verlet, which carries an explicit velocity and computes the force at the midpoint, is a rearrangement of the same method with a cleaner velocity, and it inherits the bound. The choice among these three is about convenience and damping behaviour, not about stability, because they share the second-order accuracy class and the same linear stability boundary. An integrator that claims better stability at the same order is usually trading accuracy for it, and the trade shows up as increased damping of the physical oscillations.

Stiff constraints need a different class of method

When the stiff spring is replaced by a hard constraint, the natural frequency is effectively infinite and no explicit step is stable. The historical route around this was to make the spring soft enough to fit the step, which limits how stiff a constraint can be and produces visible stretching of cloth and bending of rigid bodies. The modern route is to solve the constraints implicitly, or to enforce them by projection.

Implicit integration evaluates the force at the next state, turning the step into a solve of a nonlinear system. Backward Euler is unconditionally stable for the linear oscillator at any step size, which is the property that pays for the solve: it is not more accurate, it is stable. An implicit cloth solver can take a step twenty times larger than the explicit bound and remain bounded, at the cost of damping high-frequency modes and of solving a linear system per step whose size scales with the number of constraints. The damping is a real trade: implicit methods dissipate the fast modes that an explicit small-step method would preserve, and for cloth that reads as a loss of bounce and detail.

Projection methods, associated with position-based dynamics, take a different route. They integrate position forward with a stable explicit step, then move particles to satisfy the constraints directly, iterating until the constraints hold. The projection is not integration; it is a correction applied after the fact, and its stability comes from the fact that each correction shrinks the constraint error rather than amplifying an oscillation. The step size is then bounded by the explicit part of the integration, and the constraints contribute no stiffness because they are never expressed as forces at all.

Choosing a step is a constrained optimization

The step is set by the stiffest active feature, and the stiffest feature changes with the scene. A cloth resting on the ground has a mild maximum frequency; the same cloth pulled tight against a pin has a high one; a collision that drives two particles to the same point momentarily acts as an infinite stiffness. A fixed step chosen for the worst case wastes accuracy everywhere else, and a variable step that shrinks on demand pays the cost of bookkeeping and of changing the integrator’s error behaviour.

The practical resolution is a fixed step at a stable multiple of the render rate, with the stiffnesses clamped so the bound is never violated. Sub-stepping is the tool: run several simulation steps per rendered frame, so a 60 Hz renderer can run a 240 Hz simulation and hit a stability bound four times tighter. The cost is four evaluations of the forces per frame, which is the honest price of a stiff constraint, and it is why simulation-heavy scenes are CPU-bound by force evaluation rather than by anything else.

The bound interacts with the rest of the pipeline. A skinned character driven by simulation needs its simulated pose evaluated at the render rate, and the interpolation that smooths between simulation steps is a second approximation layered on the integrator’s error. The integrator’s step and the interpolation’s segment are two time scales that must be kept apart: the simulation step is a stability choice, and the render interpolation is a smoothness choice, and conflating them either destabilizes the solver or reintroduces the jerky motion the smoothing was meant to remove.