Section 6.4: Numerical integration and stability

"A small error today is a curiosity. The same error compounded ten thousand steps forward is a catastrophe you can no longer trace back to its cause."

Section 6.4

Every Gymnasium MuJoCo task (HalfCheetah, Hopper, Ant) and every Isaac Lab locomotion environment for ANYmal or Unitree H1 hides one decision that the reward curve never shows you: which integrator advances the state, and at what timestep. Get it wrong and the policy that scores highest in simulation is precisely the one most overfit to numerical artifact. We define that update rule, trace how its truncation error accumulates over a rollout, and run a pendulum stress test that exposes the energy drift a learning agent would otherwise silently exploit.

The practical question is sharp: for a given robot and control rate, does the integrator preserve the physical invariant the task depends on (energy in swing phase, momentum through contact), and can you prove it from a saved trace before you spend GPU-days training on it?

Theory

Every physics simulator reduces a continuous differential equation to a finite sequence of state updates. How you perform that reduction, the integrator, decides one of two outcomes. A reinforcement learning experiment either produces a policy that transfers to hardware, or one that exploits artifacts that exist only in software. The gap between those two outcomes comes down to truncation error: the small difference between the exact continuous solution and what the discrete update actually computes at each step, which the worked example below measures directly as energy drift. A controller tuned with explicit Euler at \(\Delta t = 0.05\) s experiences a subtly different dynamical system than one trained with symplectic Euler at \(\Delta t = 0.01\) s, even when every other parameter matches. Figure 6.4A makes this visible: on the same spring-mass trajectory at one fixed step size, explicit Euler spirals outward and goes unstable while the semi-implicit method stays bounded. Treat the integrator and timestep as explicit model parameters, or the resulting policy stays neither reproducible nor comparable across labs. That is why numerical integration belongs at the center of the embodied AI stack, not in a footnote about solver settings.

The practical design rule is to make the interface inspectable before optimization begins: inputs, outputs, units, latency, bounds, and failure labels should all be visible in the saved artifact.

Worked Example: Euler vs RK4 vs Symplectic Euler on a Pendulum

To make that inspectable interface concrete, the fastest way to expose what an integrator does to a saved trace is to run it on a system whose true behavior you already know exactly.

Before reading on, ask yourself: if you run the same pendulum simulation for 200 seconds with explicit Euler versus symplectic Euler, how large do you expect the energy error to be? A fraction of a percent? A few percent? Ten times worse? The answer matters because that gap is exactly what a reinforcement learning policy will learn to exploit.

The sharpest stress test for an integrator is an undamped pendulum, because its true total energy is conserved exactly. Any energy the integrator gains or loses is pure numerical artifact. The pendulum equation of motion comes from the Lagrangian \(L = T - V\) with \(T = \tfrac{1}{2}m\ell^2\dot{\theta}^2\) and \(V = -mg\ell\cos\theta\). Applying the Euler-Lagrange equation \(\tfrac{d}{dt}\tfrac{\partial L}{\partial\dot{\theta}} - \tfrac{\partial L}{\partial\theta} = 0\) gives \(\ddot{\theta} = -\tfrac{g}{\ell}\sin\theta\).

Three integrators handle the same dynamics very differently:

• Explicit (forward) Euler uses the old velocity for both updates. It systematically injects energy and the orbit spirals outward: it drifts. On a real robot this matters because a walking policy trained in a drifting simulator perceives artificially high kinetic energy at foot-strike, learns to exploit that phantom energy for propulsion, and then falls the moment real hardware absorbs the energy through leg damping instead of amplifying it. To make the scale concrete: an 18% energy surplus in simulation is equivalent to training a legged robot on a planet with 18% stronger gravity than the one it will be deployed on; the policy's learned timing, foot-placement, and push-off forces are all calibrated for conditions that never existed. The mechanical cause is a first-order truncation: the Taylor expansion of \(x(t+\Delta t)\) leaves a residual of \(\frac{1}{2}\Delta t^2 \ddot{x}\), where truncation error is the gap between the exact solution and the discrete update caused by dropping the higher-order Taylor terms. For an oscillator whose acceleration points inward, that residual always adds a small outward kick, accumulating across every step until the orbit visibly expands.
• Symplectic (semi-implicit) Euler updates velocity first, then uses the new velocity to update position. A symplectic method is one that preserves the phase-space area a trajectory sweeps out, the geometric structure Hamiltonian mechanics conserves. It does not conserve energy exactly but the error stays bounded for all time, so it is energy-stable. This is why MuJoCo and most game physics engines use it.
• RK4 takes four stage evaluations per step. Its local truncation error is \(O(\Delta t^5)\), so for smooth motion it is far more accurate, though over very long horizons it still slowly dissipates energy.

import numpy as np

g, L = 9.81, 1.0
def accel(theta):
    return -(g / L) * np.sin(theta)

def energy(theta, omega):
    return 0.5 * L**2 * omega**2 + g * L * (1 - np.cos(theta))

def integrate(method, theta0, omega0, dt, steps):
    theta, omega = theta0, omega0
    E = [energy(theta, omega)]
    for _ in range(steps):
        if method == "euler":                 # explicit Euler
            theta_new = theta + dt * omega
            omega_new = omega + dt * accel(theta)
            theta, omega = theta_new, omega_new
        elif method == "symplectic":          # semi-implicit Euler
            omega = omega + dt * accel(theta)
            theta = theta + dt * omega
        elif method == "rk4":
            def f(s):
                th, om = s
                return np.array([om, accel(th)])
            s = np.array([theta, omega])
            k1 = f(s)
            k2 = f(s + 0.5 * dt * k1)
            k3 = f(s + 0.5 * dt * k2)
            k4 = f(s + dt * k3)
            s = s + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4)
            theta, omega = s
        E.append(energy(theta, omega))
    return np.array(E)

theta0, omega0 = 0.5, 0.0       # release from 0.5 rad at rest
dt, steps = 0.05, 4000          # 200 s of simulated time
E0 = energy(theta0, omega0)

for method in ["euler", "symplectic", "rk4"]:
    E = integrate(method, theta0, omega0, dt, steps)
    drift = (E[-1] - E0) / E0 * 100
    print(f"{method:<11} energy drift over 200 s: {drift:+8.2f}%   "
          f"(min={E.min():.4f}, max={E.max():.4f})")

The printed drift makes the contrast concrete. With \(\Delta t = 0.05\) s over 200 simulated seconds, explicit Euler gains roughly +18% energy: the pendulum swings ever higher, which is physically impossible. Symplectic Euler keeps the energy oscillating within a small bounded band of about ±0.3% around \(E_0\), and RK4 holds energy to within 0.01% of the initial value. That 60-fold gap in energy error between explicit and symplectic Euler is the cost of one wrong integrator setting. A policy trained inside it learns to exploit the drift rather than solve the task. The lesson for embodied AI: choose the integrator for the property the task needs. Use symplectic methods for long-horizon stability in a learning loop, and a high-order method like RK4 when short-horizon trajectory accuracy matters more than energy bookkeeping. Always log \(\Delta t\) and the integrator, because the policy learns the discrete system you actually simulate, not the continuous one you intended.

A simulator that stays bounded is not automatically honest: numerical stability and physical fidelity are separate guarantees, and only one of them transfers to hardware.

Practical Recipe

1. Fix the integrator and timestep before any policy training run. For legged locomotion with MuJoCo (e.g., Unitree H1, ANYmal C), use integrator="implicitfast" at \(\Delta t \leq 0.005\) s; for manipulation on a Franka Panda arm use \(\Delta t \leq 0.002\) s to keep joint-torque stiffness artifacts below 1% of commanded torque. Log both values in every checkpoint name.
2. Run an energy-conservation diagnostic on a free-swinging link before contact-rich training. If explicit Euler at your chosen \(\Delta t\) drifts more than 2% over 100 s of simulated time, the policy will learn to exploit that drift in contact transitions rather than solve the task.
3. Separate free-motion from contact-rich rollouts in your benchmark set. A common failure in sim-to-real transfer for Boston Dynamics Spot gait policies is that the integrator looks stable in swing-phase but injects energy during foot-strike contact phases, which the real leg damping absorbs but the simulation does not.
4. Record the integrator type, substep count, solver tolerance, and control update rate in every saved artifact. A Franka Panda policy trained at 500 Hz control with MuJoCo substeps=4 is a different discrete system than the same policy at 1000 Hz with substeps=2, even when wall-clock physics time is identical.
5. Before scaling to GPU-parallel training in Isaac Lab or MJX, run a single-environment convergence check: sweep \(\Delta t\) by a factor of two and confirm that the task metric does not change by more than 5%. If it does, the policy is learning timestep artifacts, not the physical task.

Production Pattern

Once those diagnostics tell you which integrator and timestep keep the task honest, the next job is to fit that choice into the larger robotics pipeline that surrounds it.

Numerical integration and stability sits inside the Part II robotics contract: geometry defines where things are, kinematics defines what motion is possible, dynamics defines what motion costs, control defines how errors are corrected, and sensing defines what the agent can know on time.

Where the integrator sits in the pipeline

Tie every integration choice to four things you can read off a trace: stability, energy behavior, timestep, and controller update rate. Each choice then carries its own runnable check and a reproducible failure mode, not just an intuition.

An integrator is the contract that turns continuous dynamics into a sequence of states. The danger is that the sequence can be numerically stable even when it is physically wrong, or physically plausible for one timestep and unstable for another. For robot learning, that means a policy can learn simulator artifacts such as extra damping, energy gain, or contact jitter.

Use this recipe when turning Numerical integration and stability into code, a simulator experiment, or a robot diagnostic. The point is not to use every library. The point is to keep the hand-built baseline and the maintained-tool path comparable.

1. Specify mass, inertia, actuator limits, contact model, timestep, and solver tolerance before running a rollout.
2. Run one free-motion test and one contact test with logged energy, constraint violation, and penetration depth.
3. Compare the hand calculation with MuJoCo, Drake, Pinocchio, or MJX on the same model and timestep.
4. Store solver settings, random seed, initial state, trajectory, and failure labels in one artifact.
5. Scale to Isaac Lab or GPU-parallel simulation only after a small model passes deterministic checks.

Distrust smooth simulation until the section-specific physical assumption has been stress-tested: timestep, contact stiffness, damping, friction, actuation, and energy behavior should each have a small diagnostic.