Section 6.4: Numerical integration and stability

A Careful Control Loop

This section develops the technical contract for Numerical integration and stability into a usable mental model. First we define the object of study, then we connect it to the agent loop, then we test it with a compact implementation.

The key question in Numerical integration and stability is practical: what must the agent know, what can it observe, what action is available, and what evidence shows that the action worked under the stated conditions?

Theory

For Numerical integration and stability, 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

The cleanest 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.
• Symplectic (semi-implicit) Euler updates velocity first, then uses the new velocity to update position. 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: explicit Euler gains a large positive percentage (the pendulum swings ever higher, which is physically impossible), symplectic Euler keeps the energy oscillating within a small bounded band around $E_0$, and RK4 holds energy to within a tiny fraction of a percent. 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.

Practical Recipe

1. Write the observation, action, and success metric before choosing a model.
2. Build a baseline that is simple enough to debug by inspection.
3. Add the library implementation only after the baseline behavior is understood.
4. Record failures as structured cases: perception error, state error, planning error, control error, or evaluation error.
5. Run at least one perturbation test before trusting the result.

Production Pattern

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.

Tie numerical integration choices to stability, energy behavior, timestep, and controller update rate. This makes the section useful to students, builders, and researchers at the same time: the idea has an intuitive role, a formal interface, a runnable check, and a failure mode that can be reproduced.

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.

For Numerical integration and stability, 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.