Every MuJoCo step solves one equation for the joint accelerations:
\[M(q)\,\ddot q + c(q, \dot q) = \tau + \tau_\text{passive} + J_c^\top f_c .\]$M(q)$ is the mass matrix (joint-space inertia): symmetric, positive definite, and dependent on the configuration. $c(q, \dot q)$ is the bias force: gravity plus the Coriolis and centrifugal forces that appear when joints move. On the right are the forces that cause motion: actuators and applied forces $\tau$, passive springs and dampers, and constraint forces such as contacts and joint limits, which act through their Jacobians.
MuJoCo names each term in mjData: qM (the mass matrix in a sparse factorized form; mj_fullM turns it into a dense matrix), qfrc_bias, qfrc_actuator, qfrc_applied, qfrc_passive, qfrc_constraint, and the result qacc. This lesson derives the left-hand side for the two-link arm by hand, matches it to MuJoCo term by term, and then checks the whole equation on the 7-DOF arm.
[!established] How MuJoCo computes the terms
qfrc_biascomes from the recursive Newton-Euler algorithm (mj_rne) run with zero acceleration, so it is “the torque needed to keep the joints from accelerating”. The mass matrix comes from the composite rigid-body algorithm (mj_crb), then armature is added to its diagonal. Joint damping and springs go intoqfrc_passive, not into the bias. All of this happens inmj_forwardand in everymj_step.
The lab swings the two-link arm with no motor torque and prints the two mass-matrix entries that change: MuJoCo’s value and the Lagrangian formula below, side by side.
```lab simlab
{“dock”: true, “title”: “The mass matrix of a swinging arm”, “model”: “arm2”, “height”: 260,
“camera”: {“azimuth”: -90, “elevation”: 0, “distance”: 2.3, “target”: [0.15, 0, 1.0]}, “trail”: “ee”,
“setup”: “data.qpos[0] = 0; data.qpos[1] = 0”,
“controls”: [
{“type”: “button”, “label”: “Release horizontal”, “run”: “data.qpos[0] = 0; data.qpos[1] = 0; data.qvel.fill(0)”},
{“type”: “button”, “label”: “Spin the elbow”, “run”: “data.qpos[0] = 0.6; data.qpos[1] = -1.2; data.qvel[0] = 0; data.qvel[1] = 8”},
{“type”: “select”, “label”: “armature, both joints”, “apply”: “model.dof_armature.fill(value)”, “options”: [[“0 (as modelled)”, 0], [“0.1 kg m²”, 0.1]], “value”: 0}],
“readouts”: [
{“label”: “M11, mj_fullM”, “expr”: “lib.massMatrix()[0][0]”, “digits”: 5},
{“label”: “M11, Lagrangian”, “expr”: “0.0208333 + 0.0106667 + 0.0625 + 0.8 * (0.25 + 0.04 + 0.2 * Math.cos(data.qpos[1])) + model.dof_armature[0]”, “digits”: 5},
{“label”: “M12, mj_fullM”, “expr”: “lib.massMatrix()[0][1]”, “digits”: 5},
{“label”: “M12, Lagrangian”, “expr”: “0.0106667 + 0.8 * (0.04 + 0.1 * Math.cos(data.qpos[1]))”, “digits”: 5},
{“label”: “Coriolis + centrifugal, shoulder (N m)”, “expr”: “data.qfrc_bias[0] - 9.81 * (0.65 * Math.cos(data.qpos[0]) + 0.16 * Math.cos(data.qpos[0] + data.qpos[1]))”, “digits”: 3},
{“label”: “kinetic energy ½ q̇ᵀMq̇ (J)”, “expr”: “(() => { const M = lib.massMatrix(), v = data.qvel; return 0.5 * (M[0][0] * v[0] * v[0] + 2 * M[0][1] * v[0] * v[1] + M[1][1] * v[1] * v[1]); })()”, “digits”: 4},
{“label”: “kinetic energy, MuJoCo (J)”, “expr”: “lib.energy().kinetic”, “digits”: 4}],
“plots”: [{“ylabel”: “inertia (kg m²)”, “window”: 8, “traces”: [
{“label”: “M11”, “expr”: “lib.massMatrix()[0][0]”}, {“label”: “M12”, “expr”: “lib.massMatrix()[0][1]”, “dash”: true}]}],
“note”: “The readouts call mj_forward before reading, so MuJoCo’s numbers and the formulas refer to the same state. Gravity is subtracted from qfrc_bias with the formula of Lesson 5.1; what remains depends only on the joint velocities.”}
## Mathematics
### The two-link Lagrangian
Measure both joint angles as in Lesson 5.1, with link masses $m_i$, centre-of-mass distances $c_i$ from each joint, and inertias $I_i$ about each centre of mass. The centre of link 2 moves with the elbow plus its own rotation, so its speed squared is
$$\|v_{c2}\|^2 = L_1^2\dot q_1^2 + c_2^2(\dot q_1 + \dot q_2)^2 + 2L_1c_2\,\dot q_1(\dot q_1+\dot q_2)\cos q_2 .$$
The kinetic energy $T = \tfrac12 m_1 c_1^2\dot q_1^2 + \tfrac12 I_1\dot q_1^2 + \tfrac12 m_2\|v_{c2}\|^2 + \tfrac12 I_2(\dot q_1+\dot q_2)^2$ is a quadratic form, $T = \tfrac12\dot q^\top M(q)\dot q$, with
$$M(q) = \begin{bmatrix} I_1 + I_2 + m_1c_1^2 + m_2(L_1^2 + c_2^2 + 2L_1c_2\cos q_2) & I_2 + m_2(c_2^2 + L_1c_2\cos q_2) \\ I_2 + m_2(c_2^2 + L_1 c_2\cos q_2) & I_2 + m_2c_2^2 \end{bmatrix}.$$
> [!derivation] The bias forces from Lagrange's equations
> Lagrange's equations $\frac{d}{dt}\frac{\partial T}{\partial\dot q} - \frac{\partial T}{\partial q} + \frac{\partial V}{\partial q} = \tau$ give $M\ddot q + \dot M\dot q - \tfrac12\frac{\partial}{\partial q}(\dot q^\top M\dot q) + g(q) = \tau$. Only $\cos q_2$ appears in $M$; with $h = m_2L_1c_2\sin q_2$, $\partial M_{11}/\partial q_2 = -2h$ and $\partial M_{12}/\partial q_2 = -h$. Collecting terms,
> $$c(q,\dot q) = \begin{bmatrix} -h\,(2\dot q_1\dot q_2 + \dot q_2^2) \\ h\,\dot q_1^2 \end{bmatrix} + g(q),$$
> with $g(q)$ the gravity torque of Lesson 5.1. The $\dot q_1\dot q_2$ term is the Coriolis force, the squared terms are centrifugal.
### Properties worth knowing
$M$ is symmetric and positive definite, because the kinetic energy is positive for any motion. Its entries depend on $q$ but never on $\dot q$, and the bias is quadratic in $\dot q$ plus a configuration-dependent gravity term. MuJoCo never forms $M^{-1}$: it factorizes $M = L^\top D L$ in a sparse form that follows the kinematic tree, which is why the dense matrix is something you ask for (`mj_fullM`) rather than a field you read.
### Where armature and damping enter
Armature is reflected rotor inertia (Lesson 5.2). It moves with the joint, not with any link, so it adds to the diagonal of $M$ only: $M \leftarrow M + \mathrm{diag}(\text{armature})$. Joint damping is a force, $-b\,\dot q$, in `qfrc_passive`. Both change the dynamics, in different ways: armature makes a joint slower to accelerate, damping makes it slower to move.
### Conditioning
The ratio of $M$'s largest to smallest eigenvalue measures how differently the joints respond to the same torque. A badly conditioned $M$ makes explicit integration and torque control harder, because the light joints set the stable timestep (Lesson 5.2). Output (4) shows armature turning `arm7`'s condition number from 7372 into 35.
## Implementation
```io
INPUT: arm2.xml, arm7.xml
PROCESS: compare M(q) and the bias with mj_fullM and qfrc_bias; add armature and read qfrc_passive; check the full equation of motion on arm7; compare kinetic energy and M's spectrum
OUTPUT: printed comparisons
```python file=examples/l7_1_equations_of_motion.py “"”Lesson 7.1: the equations of motion MuJoCo solves, derived and read out of the engine.
INPUT arm2.xml (planar two-link arm, explicit inertias), arm7.xml PROCESS (1) the two-link mass matrix and bias forces from the Lagrangian, against mj_fullM and qfrc_bias at random states; (2) where armature and joint damping enter: M’s diagonal and qfrc_passive; (3) the full balance M qacc + qfrc_bias = qfrc_actuator + qfrc_passive + qfrc_constraint + qfrc_applied on arm7, with random states and controls; (4) kinetic energy as 1/2 qdot^T M qdot against MuJoCo’s energy, and the spread of M’s eigenvalues with and without armature OUTPUT printed comparisons
Run: python examples/l7_1_equations_of_motion.py “””
import os
import mujoco import numpy as np
from mjcourse import model_path
N = 200 if os.environ.get(“MJC_FAST”) == “1” else 1000 G = 9.81 M1, M2, L1, C1, C2 = 1.0, 0.8, 0.5, 0.25, 0.2 # masses (kg), link 1 length, COM distances (m) I1, I2 = 0.0208333, 0.0106667 # COM inertias as written in arm2.xml (kg m^2) I1_ROD, I2_ROD = 1.0 * 0.52 / 12, 0.8 * 0.42 / 12 # the exact thin-rod values they round
def mass_matrix_2link(q: np.ndarray, i1: float = I1, i2: float = I2) -> np.ndarray: c2 = np.cos(q[1]) m11 = i1 + i2 + M1 * C12 + M2 * (L12 + C22 + 2 * L1 * C2 * c2) m12 = i2 + M2 * (C22 + L1 * C2 * c2) m22 = i2 + M2 * C2**2 return np.array([[m11, m12], [m12, m22]])
def bias_2link(q: np.ndarray, qd: np.ndarray) -> np.ndarray: “"”Coriolis, centrifugal and gravity terms: the torque needed when qacc = 0.””” h = M2 * L1 * C2 * np.sin(q[1]) coriolis = np.array([-h * (2 * qd[0] * qd[1] + qd[1] ** 2), h * qd[0] ** 2]) g2 = G * M2 * C2 * np.cos(q[0] + q[1]) gravity = np.array([G * (M1 * C1 + M2 * L1) * np.cos(q[0]) + g2, g2]) return coriolis + gravity
def full_mass(model, data) -> np.ndarray: m = np.zeros((model.nv, model.nv)) mujoco.mj_fullM(model, data, m) return m
def two_link() -> None: model = mujoco.MjModel.from_xml_path(str(model_path(“arm2”))) data = mujoco.MjData(model) rng = np.random.default_rng(0) e_m = e_rod = e_b = 0.0 for _ in range(N): data.qpos[:] = rng.uniform(model.jnt_range[:, 0], model.jnt_range[:, 1]) data.qvel[:] = rng.normal(scale=3.0, size=2) mujoco.mj_forward(model, data) m = full_mass(model, data) e_m = max(e_m, np.abs(m - mass_matrix_2link(data.qpos)).max()) e_rod = max(e_rod, np.abs(m - mass_matrix_2link(data.qpos, I1_ROD, I2_ROD)).max()) e_b = max(e_b, np.abs(data.qfrc_bias - bias_2link(data.qpos, data.qvel)).max()) print(f” {N} random states (q in range, qdot ~ N(0, 3 rad/s)):”) print(f” max |mj_fullM - M(q)| = {e_m:.1e} kg m^2 with the inertias as written in the XML, “ f”{e_rod:.1e} with the exact m L^2 / 12 they round”) print(f” max |qfrc_bias - c(q, qdot) - g(q)| = {e_b:.1e} N m”) data.qpos[:] = [0.0, 0.0] mujoco.mj_forward(model, data) print(f” M at q = (0, 0): {np.round(full_mass(model, data), 5).tolist()} kg m^2”) data.qpos[:] = [0.0, np.pi / 2] mujoco.mj_forward(model, data) print(f” M at q = (0, pi/2): {np.round(full_mass(model, data), 5).tolist()} kg m^2 (the off-diagonal shrinks as the elbow bends)”)
def armature_and_damping() -> None: model = mujoco.MjModel.from_xml_path(str(model_path(“arm2”))) data = mujoco.MjData(model) data.qpos[:] = [0.3, -0.8] data.qvel[:] = [1.5, -2.0] mujoco.mj_forward(model, data) bare = full_mass(model, data) model.dof_armature[:] = 0.1 mujoco.mj_forward(model, data) print(f” armature 0.1 kg m^2 on both joints: M changes by {np.round(full_mass(model, data) - bare, 6).tolist()}”) print(f” qfrc_passive = {np.round(data.qfrc_passive, 6)} N m; -damping * qvel = {np.round(-model.dof_damping * data.qvel, 6)} N m “ f”(damping is not part of qfrc_bias)”)
def balance() -> None: model = mujoco.MjModel.from_xml_path(str(model_path(“arm7”))) data = mujoco.MjData(model) rng = np.random.default_rng(1) worst = 0.0 for _ in range(N): data.qpos[:] = rng.uniform(model.jnt_range[:, 0], model.jnt_range[:, 1]) data.qvel[:] = rng.normal(scale=1.0, size=model.nv) data.ctrl[:] = rng.uniform(model.actuator_ctrlrange[:, 0], model.actuator_ctrlrange[:, 1]) data.qfrc_applied[:] = rng.normal(scale=2.0, size=model.nv) mujoco.mj_forward(model, data) lhs = full_mass(model, data) @ data.qacc + data.qfrc_bias rhs = data.qfrc_actuator + data.qfrc_passive + data.qfrc_constraint + data.qfrc_applied worst = max(worst, np.abs(lhs - rhs).max()) print(f” {N} random states, controls and applied forces: max |M qacc + bias - (actuator + passive + “ f”constraint + applied)| = {worst:.1e} N m”)
def energy_and_spectrum() -> None: model = mujoco.MjModel.from_xml_path(str(model_path(“arm7”))) model.opt.enableflags |= mujoco.mjtEnableBit.mjENBL_ENERGY data = mujoco.MjData(model) data.qpos[:] = model.key_qpos[model.key(“home”).id] data.qvel[:] = np.random.default_rng(2).normal(size=model.nv) mujoco.mj_forward(model, data) m = full_mass(model, data) print(f” 1/2 qdot^T M qdot = {0.5 * data.qvel @ m @ data.qvel:.10f} J, MuJoCo energy[1] = {data.energy[1]:.10f} J”) for arm in (0.05, 0.0): model.dof_armature[:] = arm mujoco.mj_forward(model, data) eig = np.linalg.eigvalsh(full_mass(model, data)) print(f” armature {arm:<4}: M symmetric {np.allclose(full_mass(model, data), full_mass(model, data).T)}, “ f”eigenvalues from {eig[0]:.2e} to {eig[-1]:.3f} kg m^2, condition number {eig[-1] / eig[0]:.0f}”)
if name == “main”: print(“(1) two-link arm: Lagrangian against MuJoCo”) two_link() print(“(2) armature and damping”) armature_and_damping() print(“(3) the balance MuJoCo solves, arm7”) balance() print(“(4) kinetic energy and the conditioning of M, arm7 at home”) energy_and_spectrum()
Output:
```text
(1) two-link arm: Lagrangian against MuJoCo
1000 random states (q in range, qdot ~ N(0, 3 rad/s)):
max |mj_fullM - M(q)| = 3.9e-16 kg m^2 with the inertias as written in the XML, 3.3e-08 with the exact m L^2 / 12 they round
max |qfrc_bias - c(q, qdot) - g(q)| = 5.3e-15 N m
M at q = (0, 0): [[0.486, 0.12267], [0.12267, 0.04267]] kg m^2
M at q = (0, pi/2): [[0.326, 0.04267], [0.04267, 0.04267]] kg m^2 (the off-diagonal shrinks as the elbow bends)
(2) armature and damping
armature 0.1 kg m^2 on both joints: M changes by [[0.1, 0.0], [0.0, 0.1]]
qfrc_passive = [-0.03 0.04] N m; -damping * qvel = [-0.03 0.04] N m (damping is not part of qfrc_bias)
(3) the balance MuJoCo solves, arm7
1000 random states, controls and applied forces: max |M qacc + bias - (actuator + passive + constraint + applied)| = 9.2e-12 N m
(4) kinetic energy and the conditioning of M, arm7 at home
1/2 qdot^T M qdot = 2.2154066574 J, MuJoCo energy[1] = 2.2154066574 J
armature 0.05: M symmetric True, eigenvalues from 5.02e-02 to 1.781 kg m^2, condition number 35
armature 0.0 : M symmetric True, eigenvalues from 2.35e-04 to 1.731 kg m^2, condition number 7372
Reading the output:
arm2.xml rounds $1/48$ to 0.0208333. Machine-precision agreement is a test of the model as written, not of the physics you intended.qfrc_passive as $-b\dot q$ and not in the bias.arm7’s smallest eigenvalue is $2.35 \times 10^{-4}$ kg m², the bare flange spinning about its axis, and $M$ is 200 times worse conditioned.mjData has no field qM in your code, or mj_fullM complains about its arguments. mj_fullM takes (model, data, dst) since MuJoCo 3.10 and reads the factorized matrix from data; data.qM itself was removed in 3.11, so older examples that pass it fail.
Your hand-derived dynamics disagree with MuJoCo by $10^{-8}$. Check the inertias in the XML before suspecting the derivation: rounded constants are the usual cause (output 1).
Gravity compensation with qfrc_bias makes a moving arm drift. qfrc_bias includes the Coriolis and centrifugal terms of the current velocity, which is correct for full compensation; if you meant gravity only, compute it at zero velocity: copy the state into a scratch MjData, set its qvel to 0, run mj_forward and read its qfrc_bias.
Damping “does not show up” in qfrc_bias. It is in qfrc_passive (output 2).
Add a 0.5 kg point mass at the tip of arm2 (a small sphere welded to link 2). Update $M(q)$ and $g(q)$ by hand, then check them against mj_fullM and qfrc_bias at five poses.
Verify the passivity property of the equations of motion numerically: along any trajectory with zero actuation, zero damping and no gravity, $\frac{d}{dt}\left(\tfrac12\dot q^\top M\dot q\right) = 0$. Then show that this implies $\dot q^\top\left(\dot M - 2C\right)\dot q = 0$ for the matrix $C$ with $C\dot q$ equal to the Coriolis and centrifugal terms, and check the skew-symmetry of $\dot M - 2C$ for the two-link arm with the Christoffel-symbol choice of $C$.
[!research] Learned and structured dynamics Lagrangian and Hamiltonian neural networks build exactly this structure into a learned model: a positive-definite $M(q)$, a potential, and the bias forces that follow from them. The structure guarantees energy behaviour that an unstructured network does not, at the cost of assuming rigid-body mechanics. MuJoCo’s analytic terms are the reference such models are usually tested against; Level 17 compares learned world models with this kind of ground truth.
{"id": "7.1-check", "title": "Knowledge check", "questions": [
{"kind": "numeric", "q": "For the two-link arm, what is M11 (kg m²) with the arm straight (q2 = 0)? Use I1 = 0.0208333, I2 = 0.0106667, m1 = 1, c1 = 0.25, m2 = 0.8, L1 = 0.5, c2 = 0.2.",
"answer": 0.486, "tol": 0.0005, "unit": "kg m^2",
"explain": "<p>$0.0208333 + 0.0106667 + 0.0625 + 0.8 \\times (0.25 + 0.04 + 0.2)$ = 0.486 kg m², as output (1) prints.</p>"},
{"kind": "mcq", "q": "Which of these does armature change?",
"options": ["The off-diagonal entries of M", "The diagonal entries of M only", "qfrc_bias", "qfrc_passive"],
"answer": 1,
"explain": "<p>Armature is inertia that turns with the joint itself, so it adds to that joint's diagonal entry and nothing else (output 2).</p>"},
{"kind": "mcq", "q": "The elbow of the two-link arm spins while the shoulder is held still. Which term makes the shoulder motor feel a torque that depends on the elbow speed squared?",
"options": ["Gravity", "The Coriolis term 2 q̇1 q̇2", "The centrifugal term in q̇2²", "Armature"],
"answer": 2,
"explain": "<p>With $\\dot q_1 = 0$, the shoulder bias is $-h\\dot q_2^2$ plus gravity: centrifugal. The Coriolis term needs both joints moving.</p>"},
{"kind": "mcq", "q": "After mj_forward, which identity holds for a model without contacts?",
"options": ["M qacc = qfrc_actuator", "M qacc + qfrc_bias = qfrc_actuator + qfrc_passive + qfrc_constraint + qfrc_applied", "qfrc_bias = qfrc_passive", "qacc = M qfrc_bias"],
"answer": 1,
"explain": "<p>Output (3) checks it to $10^{-11}$ N m on <code>arm7</code>. Joint limits enter through <code>qfrc_constraint</code>.</p>"}
]}
Lesson 7.2 runs the equation in both directions: forward dynamics (torques to accelerations, mj_forward) and inverse dynamics (accelerations to torques, mj_inverse), and finds out what it takes to make them agree exactly through a timestep.