A simulation “explodes” when the numerical solution grows without bound while the physics it represents would not. In MuJoCo that almost always means one of three things: something in the model is stiffer than the integrator can follow at the chosen timestep; a force is computed explicitly on a body too light to take it; or the state starts somewhere the contact model must violently correct. A fourth failure looks different but has the same roots: objects passing through each other, because the contact was never seen or was too soft to stop them.

This lesson measures each one, shows how MuJoCo reports it (and how it hides it), and orders the fixes by what they cost. Earlier lessons supplied pieces: the explicit damping limit $h < 2I/c$ (Lesson 1.3), armature for light joints (Lesson 5.2), and the time-constant rule for contacts (Lesson 9.2).

Visual intuition

The lab fires a 1 cm ball at a 2 mm wall and reports whether the ball stopped, pushed through a contact that was detected, or tunnelled past without any contact being seen.

```lab simlab {“dock”: true, “title”: “Stopped, pushed through, or tunnelled”, “height”: 220, “xml”: “<body name='ball' pos='0 0 0'></body>”, “camera”: {“azimuth”: -90, “elevation”: 35, “distance”: 0.7, “target”: [0.2, 0, 0]}, “setup”: “ctx.touched = false”, “controls”: [ {“label”: “speed”, “unit”: “m/s”, “param”: “v”, “min”: 0.5, “max”: 100, “step”: 0.5, “value”: 5, “digits”: 1}, {“type”: “button”, “label”: “Fire”, “run”: “mj.mj_resetData(model, data); data.qvel[0] = ctx.v; ctx.touched = false; sim.play()”}, {“label”: “contact time constant”, “unit”: “s”, “apply”: “model.geom_solref[0] = value; model.geom_solref[2] = value”, “min”: 0.001, “max”: 0.05, “step”: 0.001, “value”: 0.02, “digits”: 3}, {“type”: “select”, “label”: “timestep”, “set”: “model.opt.timestep”, “options”: [[“2 ms”, 0.002], [“0.5 ms”, 0.0005]], “value”: 0.002}], “controller”: “ctx.touched = ctx.touched || data.ncon > 0”, “readouts”: [ {“label”: “outcome”, “expr”: “data.qpos[0] < 0.25 ? (ctx.touched ? ‘stopped’ : (data.qvel[0] > 0 ? ‘on its way’ : ‘-‘)) : (ctx.touched ? ‘push-through’ : ‘tunnelled’)”}, {“label”: “contact seen”, “expr”: “ctx.touched ? ‘yes’ : ‘no’”}, {“label”: “ball x (m)”, “expr”: “data.qpos[0]”, “digits”: 3}, {“label”: “time (s)”, “expr”: “data.time”, “digits”: 3}], “stopAt”: 0.5, “note”: “The wall is at x = 0.25 m and 2 mm thick; the ball starts at x = 0 with radius 1 cm. refsafe raises time constants below twice the timestep to 2h, so 1 ms only takes effect at the 0.5 ms step.”}


## Mathematics

### Stiffness and the stability region

For a mass $m$ on a spring of stiffness $k$, $\omega = \sqrt{k/m}$. One step of each integrator multiplies the oscillation by a factor whose magnitude decides stability:

- **Semi-implicit Euler** (`Euler`), stable for $\omega h < 2$ and energy-preserving inside it. Near the limit the oscillation's position amplitude grows to $1/\sqrt{1 - (\omega h/2)^2}$ times its true value (3.2 at $\omega h$ = 1.9 in output 1) without diverging. `implicitfast` and `implicit` treat only velocity-dependent forces implicitly, so a spring has exactly the same limit under them.
- **RK4**, stable up to $\omega h = 2\sqrt2 \approx 2.83$, and numerically damped inside: at $\omega h$ = 1.41 each step multiplies the amplitude by 0.958, so the oscillation is gone within a second.
- **discrete**, implicit in position as well, never diverges, but its factor $1/\sqrt{1 + (\omega h)^2}$ damps every oscillation; even at $\omega h$ = 0.2 it removes 2% per step. Stability bought this way changes the physics.

Contacts, joint limits and servos are springs too: a contact's time constant sets its effective $\omega$, which is why `refsafe` keeps it at least $2h$.

### Push-through and tunnelling

A body approaching a contact at speed $v$ is decelerated by the constraint's reference dynamics, which act on a time scale of the time constant: it needs a penetration of the order of $v \cdot \text{timeconst}$ to stop (an estimate). If that exceeds the obstacle's thickness, it **pushes through**: the contact is detected, but its normal flips once the ball's centre passes the wall's mid-plane, and the contact then helps it out of the far side. If instead $v\,h$ exceeds the thickness plus the body's size, no step ends with the two overlapping and the contact is never seen: **tunnelling**. MuJoCo detects collisions at discrete instants; there is no sweep between them. Output (2) shows both, and that a stiffer contact at a smaller step stops even 100 m/s.

### Starting inside something

A body that starts overlapping by $\delta$ is pushed out by the constraint's reference acceleration $-k\,r$, with $k \propto 1/\text{timeconst}^2$. With the default 20 ms time constant a 1 kg box that starts 20 mm deep rises gently (0.38 m/s, no overshoot); with 4 ms it is ejected at 2.6 m/s and flies 33 cm. Settle initial states in simulation, or generate them without overlap, before you start a stiff scene.

### Divergence is reset, not reported

Output (4) shows MuJoCo's behaviour on a diverging step: one warning, the state reset to `qpos0`, the clock restarted, and the simulation continuing as if nothing had happened. Here the reset state is a spring at rest, which never diverges again, so a run of 1000 steps finishes "successfully" with one warning. Disable `autoreset` while developing to make divergence visible (the state becomes NaN), and check the warning counters after every evaluation run.

### Fixes, in order of cost

The cheapest fixes change the model, not the simulation: realistic masses and inertias, armature on geared joints (Lesson 5.2), damping in the model rather than in the controller (Lessons 5.2 and 8.1), initial states without overlap, thicker collision geoms. Next, the integrator: `implicitfast` costs about the same as `Euler` and removes the damping limit for modelled damping and actuator `kv`; `implicit` costs about 30% more here; `discrete` was even cheaper per step on this scene (we have not isolated why; it used fewer solver iterations) but damps oscillations. Then the contact time constant, free but floored at $2h$. Last, the timestep itself: halving it doubles the cost of every simulated second, and RK4 costs nearly three times `Euler` per step while remaining explicit.

## Implementation

```io
INPUT: small scenes in the script; pick_place.xml
PROCESS: map the spring stability limit for five integrators; fire a ball at a wall for four settings and four speeds; start a box inside the floor; force a divergence and watch the reset; time pick_place under each integrator
OUTPUT: printed tables

```python file=examples/l9_4_stability.py “"”Lesson 9.4: why simulations explode, how to tell, and what the fixes cost.

INPUT small scenes written below; pick_place.xml for the cost of fixes PROCESS (1) a 1 kg mass on a spring: largest stable omega * h for five integrators, and how much each damps the oscillation; (2) a small ball fired at a thin wall: tunnelling (no contact ever detected) versus soft push-through (detected but too soft to stop it), and the fixes; (3) a box that starts inside the floor: ejection speed for two contact stiffnesses; (4) what MuJoCo does when a simulation diverges: the warning, the automatic reset; (5) the cost of the fixes: microseconds per step of pick_place.xml for each integrator and timestep OUTPUT printed tables

Core concept

Visual intuition

Debugging

Exercise

Challenge

Research connection

Next