Core concept

One call to mujoco.mj_step(model, data) does two jobs. First it runs forward dynamics: from the current positions, velocities and inputs it computes every force and the resulting accelerations. Then it integrates: it uses the accelerations to move the state forward by one timestep. mj_forward is the first job alone, which is why it is safe to call any number of times without moving time forward.

Here is mj_step as it is written in MuJoCo 3.14.0 (src/engine/engine_forward.c), with the work grouped into the stages the rest of the course refers to:

Notice where contacts come from: collision detection runs in the position stage, from the positions at the start of the step. That is why, in Lesson 0.1, a contact appeared one step after the ball had already crossed the floor.

[!implementation] mj_step1 and mj_step2 MuJoCo also offers the step split in two: mj_step1 runs the position and velocity stages, mj_step2 runs actuation, acceleration, constraint and integration. Between the two calls, mjData holds fresh sensor readings and kinematics for the current state, so a controller can read them, write ctrl, and then finish the step. The documentation notes one catch: mj_step2 always uses a single-step integrator, so with RK4 selected the split falls back to Euler.

Visual intuition

The lab in the side panel runs the same model from the same initial state under all five integrators at once, each in its own compiled model, and plots how far each one’s total energy has drifted from its starting value. The double pendulum has no damping and no contacts, so its true energy is constant: any drift is numerical error.

```lab integrators {“dock”: true, “model”: “double_pendulum”, “key”: 0, “timestep”: 0.002}

Run it at 2 ms, then reset and try 0.5 ms and 10 ms. Then switch to the single pendulum and compare.

## Mathematics

### The integrators

MuJoCo's documentation writes the single-step methods as one implicit-in-velocity update, differing in which force derivatives they treat implicitly. With $M$ the mass matrix, $a(v_t)$ the acceleration from forward dynamics and $D$ a matrix of (negative) force derivatives with respect to velocity,

$$v_{t+h} = v_t + h\,\widehat{M}^{-1} M\, a(v_t), \qquad \widehat M = M + hD, \qquad q_{t+h} = q_t \oplus h\, v_{t+h},$$

where $\oplus$ is the position update of Lesson 1.2 (quaternion-aware).

| Integrator | $D$ contains | Notes from the 3.14.0 documentation |
|---|---|---|
| `Euler` (default) | joint damping only | semi-implicit Euler with implicit joint damping; "use for compatibility with older models" |
| `implicitfast` | derivatives of all forces except Coriolis and centripetal, symmetrized | "the recommended integrator for most models"; similar cost to `Euler`, more stable |
| `implicit` | derivatives of all forces, including Coriolis and centripetal (not constraint forces) | for coupled rotational systems with large velocity-dependent forces; needs an LU factorization |
| `discrete` (new in 3.13) | positive damping and stiffness terms, solved jointly with the constraints; also implicit in position | "choose when stiff springs, position servos or strong damping interact with contacts"; `qacc` then holds $(v^+ - v)/h$; documented as new and subject to change |
| `RK4` | (four evaluations per step) | "qualitatively better than the single-step methods" for energy-conserving systems, even at a 4 times smaller timestep for equal cost |

### Why explicit integration goes unstable

Take one joint with inertia $I$ held by a spring-damper: $I\ddot q = -k q - c \dot q$, the shape of a position servo with gains $k = k_p$ and $c = k_v$.

> [!derivation] Two stability limits for explicit integration
> Treat the damper explicitly and ignore the spring: $v_{t+h} = v_t - \tfrac{h c}{I} v_t = (1 - \tfrac{hc}{I})\, v_t$. The velocity shrinks in magnitude only if $|1 - hc/I| < 1$, that is
> $$h < \frac{2I}{c}.$$
> Treat the spring explicitly (semi-implicit Euler, no damping): the update matrix for $(q, v)$ has eigenvalues on the unit circle only if $h\omega < 2$ with $\omega = \sqrt{k/I}$, that is
> $$h < \frac{2}{\omega} = 2\sqrt{I/k}.$$
> Beyond either limit errors grow geometrically every step.

Try it before reading the measurements. The lab holds the pendulum with a position servo (target angle 0) and lets you choose the gains, the timestep and the integrator. Predict the outcome of each setting from the two limits above, then press Play.

```lab simlab
{"title": "Servo stability: gains, timestep, integrator", "model": "pendulum_servo", "height": 220,
 "camera": {"azimuth": -90, "elevation": 5, "distance": 1.9, "target": [0, 0, 0.75]},
 "setup": "data.qpos[0] = 1.0",
 "controls": [
   {"label": "stiffness kp", "unit": "N m/rad", "apply": "model.actuator_gainprm[0] = value; model.actuator_biasprm[1] = -value", "min": 10, "max": 40000, "step": 10, "value": 100, "digits": 0},
   {"label": "damping kv", "unit": "N m s/rad", "apply": "model.actuator_biasprm[2] = -value", "min": 0, "max": 100, "step": 0.5, "value": 5, "digits": 1},
   {"label": "timestep h", "unit": "s", "set": "model.opt.timestep", "min": 0.001, "max": 0.02, "step": 0.001, "value": 0.002, "digits": 3, "reset": true},
   {"type": "select", "label": "integrator", "set": "model.opt.integrator", "options": [["Euler", 0], ["implicitfast", 3], ["discrete", 4]], "value": 0, "reset": true}
 ],
 "readouts": [
   {"label": "explicit spring limit 2/w (ms)", "expr": "2000 * Math.sqrt(0.2501 / model.actuator_gainprm[0])", "digits": 2},
   {"label": "explicit damping limit 2I/kv (ms)", "expr": "model.actuator_biasprm[2] < 0 ? 2000 * 0.2501 / -model.actuator_biasprm[2] : Infinity", "digits": 2},
   {"label": "angle (rad)", "expr": "data.qpos[0]", "digits": 4}],
 "plots": [{"ylabel": "angle (rad)", "window": 3, "ymin": -1.2, "ymax": 1.2, "traces": [{"label": "q (rad)", "expr": "data.qpos[0]"}]}],
 "note": "The pendulum starts at 1 rad and the servo pulls it to 0. If MuJoCo detects divergence it prints a warning in the view and resets the state, which you will see as the pendulum snapping back to hang straight down. Changing the timestep or integrator restarts the run."}

For the course pendulum, $I = mL^2 + I_c = 0.2501$ kg·m². With a servo damping $k_v$ = 60 N·m·s/rad, the explicit limit is $2I/k_v$ = 8.3 ms; with $k_p$ = 20 000 N·m/rad, $\omega$ = 283 rad/s and the limit is 7.1 ms. Now look at what MuJoCo does at 2, 10 and 20 ms:

INPUT: `double_pendulum.xml`; `pendulum.xml` with a position servo of stiffness kp and damping kv
PROCESS: energy drift over 20 s per integrator; servo stability for 3 s at three timesteps
OUTPUT: two tables

```python file=examples/l1_3_integrators.py “"”Lesson 1.3: what the integrator choice changes, measured.

INPUT double_pendulum.xml (conservative, chaotic) and pendulum.xml with its torque motor swapped for a position servo of varying stiffness and damping PROCESS (1) energy error after 20 s for each integrator at h = 2 ms; (2) does a servo-held pendulum stay stable at h = 2, 10 and 20 ms? OUTPUT two tables; “UNSTABLE” means MuJoCo raised its bad-acceleration warning (and, as it always does then, reset the state)

Run: python examples/l1_3_integrators.py “””

import mujoco

from mjcourse import model_path

INTEGRATORS = {“Euler”: 0, “RK4”: 1, “implicit”: 2, “implicitfast”: 3, “discrete”: 4}

def energy_error(integrator: int, timestep: float, seconds: float = 20.0) -> float: model = mujoco.MjModel.from_xml_path(str(model_path(“double_pendulum”))) model.opt.integrator, model.opt.timestep = integrator, timestep model.opt.enableflags |= mujoco.mjtEnableBit.mjENBL_ENERGY data = mujoco.MjData(model) mujoco.mj_resetDataKeyframe(model, data, 0) mujoco.mj_forward(model, data) e0 = data.energy.sum() for _ in range(round(seconds / timestep)): mujoco.mj_step(model, data) return 100.0 * (data.energy.sum() - e0) / abs(e0)

def servo_is_stable(kp: float, kv: float, integrator: int, timestep: float, seconds: float = 3.0) -> bool: xml = model_path(“pendulum”).read_text().replace( ‘’, f’’) model = mujoco.MjModel.from_xml_string(xml) model.opt.integrator, model.opt.timestep = integrator, timestep data = mujoco.MjData(model) data.qpos[0] = 1.0 # start 1 rad away from the target 0 for _ in range(round(seconds / timestep)): mujoco.mj_step(model, data) return data.warning[mujoco.mjtWarning.mjWARN_BADQACC].number == 0

if name == “main”: print(“double pendulum, energy error after 20 s at h = 2 ms”) for name, code in INTEGRATORS.items(): print(f” {name:<13s} {energy_error(code, 0.002):+8.2f} %”)

print("\nposition servo on the pendulum (I = 0.25 kg m^2), stable for 3 s?")
cases = [(100, 60), (20000, 5)]
print(f"  {'kp, kv':<14s}{'h (ms)':>7s}" + "".join(f"{n:>14s}" for n in ("Euler", "implicitfast", "discrete")))
for kp, kv in cases:
    for h in (0.002, 0.01, 0.02):
        cells = ["stable" if servo_is_stable(kp, kv, INTEGRATORS[n], h) else "UNSTABLE"
                 for n in ("Euler", "implicitfast", "discrete")]
        print(f"  {f'{kp}, {kv}':<14s}{1000 * h:>7.0f}" + "".join(f"{c:>14s}" for c in cells)) ``` Output with MuJoCo 3.14.0 (the warnings MuJoCo prints to standard error are omitted):

double pendulum, energy error after 20 s at h = 2 ms
  Euler           -10.05 %
  RK4              -0.00 %
  implicit        +21.64 %
  implicitfast    -10.05 %
  discrete        -10.05 %

position servo on the pendulum (I = 0.25 kg m^2), stable for 3 s?
  kp, kv         h (ms)         Euler  implicitfast      discrete
  100, 60             2        stable        stable        stable
  100, 60            10      UNSTABLE        stable        stable
  100, 60            20      UNSTABLE        stable        stable
  20000, 5            2        stable        stable        stable
  20000, 5           10      UNSTABLE      UNSTABLE        stable
  20000, 5           20      UNSTABLE      UNSTABLE        stable

The second table is the derivation, confirmed. With heavy servo damping, Euler fails above the 8.3 ms limit because it treats actuator damping explicitly (its $D$ contains only joint damping), while implicitfast and discrete treat it implicitly and stay stable. With a stiff servo, both Euler and implicitfast fail above the 7.1 ms limit because the spring term is explicit in both; only discrete, implicit in position, stays stable.

The first table needs more care. All the single-step methods except implicit agree here because the double pendulum has no damping, so their $D$ matrices are empty and they reduce to the same update. They lose 10% of the energy in 20 s, and implicit gains 22%. RK4 conserves it to the printed precision.

[!research] Reading the energy table honestly The double pendulum is chaotic, so the per-integrator numbers jump around as the timestep changes: in this course’s runs, implicit gained between 5% and 27% over 20 s at timesteps from 0.5 to 2 ms, and at 5 ms on this point-mass model its energy grew without bound (velocities reached 45 rad/s) with no warning raised. Treat the pattern (RK4 conserves, the single-step methods drift by percent-level amounts on an undamped chaotic system) as the finding, not the individual percentages. The documentation recommends implicit for coupled rotational systems with large velocity-dependent forces; this measurement is a reminder that its behaviour on an undamped system is a different question, worth checking on your own model before relying on it.

Choosing timestep and integrator

[!recommendation] A procedure, not a rule

1. Start from implicitfast at 2 ms, as the documentation recommends for most models.
2. If the model is energy-conserving and long-horizon accuracy matters (orbits, undamped pendula), compare against RK4 at a quarter of the timestep, which costs the same.
3. If stiff servos, springs or strong damping are the reason you cannot raise the timestep, try discrete, and read its documented limitations first (it is new in 3.13 and marked as subject to change).
4. Raise the timestep until a quantity you care about (a contact force, a tracking error, a success rate) changes by more than your tolerance. Use the largest timestep below that point. The trade-off is direct: halving the timestep doubles the cost of every experiment.

Debugging

“The simulation suddenly snaps back to the starting pose.” MuJoCo detected NaN or huge accelerations, printed WARNING: Nan, Inf or huge value in QACC at DOF ... The simulation is unstable, and reset the state. Check data.warning[mujoco.mjtWarning.mjWARN_BADQACC].number after every rollout in an experiment; a reset in the middle of an episode silently corrupts any metric computed from it.

“My controller reads sensors that are one step old.” With mj_step, the controller you run before the call sees sensordata computed in the previous step. Use mj_step1, read the fresh values, set ctrl, then mj_step2, or install a control callback.

“RK4 made my contact scene worse.” RK4 evaluates forward dynamics four times per step and treats damping and constraint stabilization explicitly; the documentation notes that with large velocity-dependent forces, a single-step method that integrates them implicitly “can be significantly more stable than RK4”.

Exercise

Using examples/l1_3_integrators.py as a starting point, measure the largest stable timestep for the servo-held pendulum with $k_p$ = 2000 N·m/rad and $k_v$ = 5 N·m·s/rad under Euler, implicitfast and discrete, by bisection to 0.5 ms. Compare the Euler result with the bound $2\sqrt{I/k_p}$ and explain the difference.

Challenge

Implement your own fourth-order Runge-Kutta step for a model without quaternions, using only mj_forward to evaluate accelerations, and check it against MuJoCo’s RK4 on the double pendulum to $10^{-10}$ over 1 s. Then extend it to free joints with mj_integratePos. The documentation notes that users “can easily implement other integrators by calling mj_forward and integrating accelerations themselves”; this is how.

Research connection

Sampling-based model-predictive control (Level 21.4) rolls out thousands of short trajectories per decision and gains directly from larger timesteps, so it is where implicitfast and discrete pay most. Reinforcement-learning environments usually fix both timestep and integrator; changing either changes the task, and results with different settings are not comparable. Report both, with the MuJoCo version.

{"id": "1.3-check", "title": "Knowledge check", "questions": [
  {"kind": "mcq", "q": "In which stage of <code>mj_step</code> is collision detection performed?",
   "options": ["Velocity stage", "Position stage", "Constraint stage", "Integration"],
   "answer": 1,
   "explain": "<p><code>mj_collision</code> runs inside <code>mj_fwdPosition</code>, from the positions at the start of the step.</p>"},
  {"kind": "numeric", "q": "A joint with inertia 0.05 kg m² is held by a position servo with <code>kv</code> = 20 N m s/rad. What is the largest timestep (in ms) at which explicit treatment of that damping is stable?",
   "answer": 5, "tol": 0.05, "unit": "ms",
   "explain": "<p>$h < 2I/c = 2 \\times 0.05 / 20$ = 0.005 s = 5 ms. Under <code>Euler</code>, actuator damping is explicit, so this bound applies; <code>implicitfast</code> removes it.</p>"},
  {"kind": "predict", "q": "The pendulum is held by a stiff servo (kp = 20000, kv = 5) at h = 10 ms. Which integrator keeps it stable?",
   "options": ["Euler", "implicitfast", "discrete", "None of them"],
   "answer": 2,
   "explain": "<p>Only <code>discrete</code> is implicit in position (the $h^2K$ term), so only it is stable above the explicit spring limit of 7.1 ms. The lab below runs the case under <code>implicitfast</code>: watch the warning appear and the state reset, then switch the integrator to discrete.</p>",
   "sim": {"model": "pendulum_servo", "height": 200, "camera": {"azimuth": -90, "elevation": 5, "distance": 1.9, "target": [0, 0, 0.75]},
           "setup": "data.qpos[0] = 1.0",
           "controls": [
             {"label": "stiffness kp", "unit": "N m/rad", "apply": "model.actuator_gainprm[0] = value; model.actuator_biasprm[1] = -value", "min": 10, "max": 40000, "step": 10, "value": 20000, "digits": 0},
             {"label": "timestep h", "unit": "s", "set": "model.opt.timestep", "min": 0.001, "max": 0.02, "step": 0.001, "value": 0.01, "digits": 3, "reset": true},
             {"type": "select", "label": "integrator", "set": "model.opt.integrator", "options": [["Euler", 0], ["implicitfast", 3], ["discrete", 4]], "value": 3, "reset": true}],
           "readouts": [{"label": "time (s)", "expr": "data.time"}, {"label": "angle (rad)", "expr": "data.qpos[0]", "digits": 4}]}},
  {"kind": "open", "q": "Your RL environment uses Euler at 2 ms and training is slow. A colleague proposes switching to 10 ms. What do you check before accepting the change, and what would make you reject it?",
   "reference": "<p>Check (1) stability: run the existing scripted or trained policy at 10 ms and count MuJoCo's bad-acceleration warnings; (2) fidelity: compare task-relevant quantities between 2 ms and 10 ms from identical states (contact forces during grasps, object trajectories, success of a fixed policy); (3) whether switching to <code>implicitfast</code> recovers stability at 10 ms. Reject if success of a fixed policy or contact behaviour changes beyond your tolerance, because then you have changed the task, not just sped it up. Report the final timestep and integrator either way.</p>"}
]}

Next

Level 2 teaches MJCF, the language every model in this course is written in, starting with the kinematic tree: bodies, geoms, joints and sites.