Part 4: Robot Dynamics: The Physics of Force

Welcome to Part 4 of "Robotics Zero to Hero." Kinematics (Parts 1–3) described motion without asking what causes it. But robots have mass, accelerate hard, and fight gravity. To control them accurately we need Dynamics — the relationship between forces/torques and motion.

The Math: Euler–Lagrange Equations

We derive equations of motion from energy. Define kinetic energy $T$ , potential energy $V$ , and the Lagrangian

L = T - V

For each generalized coordinate $q_i$ , the Euler–Lagrange equation gives the required generalized force (torque) $\tau_i$ :

\frac{d}{dt}\!\left(\frac{\partial L}{\partial \dot q_i}\right) - \frac{\partial L}{\partial q_i} = \tau_i

Carrying this out for every joint yields the canonical manipulator equation:

M(\mathbf q)\,\ddot{\mathbf q} + C(\mathbf q,\dot{\mathbf q})\,\dot{\mathbf q} + \mathbf g(\mathbf q) = \boldsymbol\tau

$M(\mathbf q)$ — the mass/inertia matrix ( $n\times n$ , symmetric positive-definite),
$C(\mathbf q,\dot{\mathbf q})$ — Coriolis and centrifugal terms (quadratic in velocity),
$\mathbf g(\mathbf q)$ — gravity,
$\boldsymbol\tau$ — joint torques.

Structure you can exploit

This equation is not a generic black box — it has properties that good controllers lean on:

$M(\mathbf q)$ is symmetric positive-definite, so it is always invertible: forward dynamics $\ddot{\mathbf q} = M^{-1}(\boldsymbol\tau - C\dot{\mathbf q} - \mathbf g)$ is well-posed.
Skew-symmetry: with the right choice of $C$ , the matrix $\dot M - 2C$ is skew-symmetric. This is the linchpin of passivity-based and Lyapunov control proofs.
Linearity in parameters: the dynamics are linear in the mass/inertia parameters, $\;M\ddot{\mathbf q} + C\dot{\mathbf q} + \mathbf g = Y(\mathbf q,\dot{\mathbf q},\ddot{\mathbf q})\,\boldsymbol\Theta$ , which is exactly what makes adaptive control and system identification (Part 5) possible.

Python Implementation: Symbolic Derivation

Deriving $M, C, \mathbf g$ by hand beyond two links is error-prone. sympy does it symbolically.

Loading Interactive Python Environment...

Competing Formulations: Lagrange vs. Newton–Euler vs. Recursive

The Lagrangian route is elegant for deriving equations, but it is the wrong tool for computing them on a 30-DOF robot at 1 kHz. The competing formulations differ enormously in cost:

Formulation	What it's good at	Complexity	Shortcomings
Euler–Lagrange (symbolic)	Insight, small systems, control proofs	$O(n^3)$ if expanded naively	Explodes symbolically; impractical to evaluate for large $n$
Newton–Euler (force/torque balance)	Physical intuition per link	$O(n^2)$ – $O(n^3)$	Bookkeeping-heavy; constraint forces
Recursive Newton–Euler (RNEA)	Inverse dynamics ( $\boldsymbol\tau$ from $\ddot{\mathbf q}$ )	$O(n)$	None for inverse dynamics — this is the standard
Articulated-Body Algorithm (ABA)	Forward dynamics ( $\ddot{\mathbf q}$ from $\boldsymbol\tau$ )	$O(n)$	More complex to implement
Composite-Rigid-Body (CRBA)	Building $M(\mathbf q)$ explicitly	$O(n^2)$	Only the mass matrix

The headline: a naive symbolic mass matrix costs $O(n^3)$ , but Featherstone's recursive algorithms compute the same physics in $O(n)$ by propagating velocities outward and forces inward along the chain. For a 6-DOF arm the difference is academic; for a humanoid or a finely-discretized tentacle it is the difference between real-time and impossible.

High-Dimensional vs. Low-Dimensional Dynamics

This post is where the dimensional gap becomes a computational one:

Low-dim (6-DOF arm): $M(\mathbf q)$ is $6\times 6$ . You can invert it, store it, even derive it symbolically. Coriolis terms are a nuisance but bounded.
High-dim (humanoid, hand): $M$ is $30\times 30$ +. You must use $O(n)$ recursive algorithms; the $O(n^3)$ matrix inverse becomes the bottleneck. Contact and closed loops add constraint equations on top.
Continuum (octopus tentacle): there are no discrete joints at all. The dynamics become partial differential equations — Cosserat-rod PDEs in arc-length $s$ and time $t$ — not ODEs. You discretize the rod into elements (trading dimension for fidelity, exactly as in Part 2) and the "mass matrix" becomes a banded operator over the discretization. Real-time simulation here is an active research frontier.

In short: as dimension grows, the form of the dynamics changes (ODE → PDE), not just the size of the matrices.

Focus on the Octopus: Gravitational Sag and Whiplash

For industrial robots, links are rigid, so $M$ , $C$ , and $\mathbf g$ are predictable. A robot made of flexible metallic segments is a different animal.

Gravitational sag. Our continuum tentacle bends under its own weight. We cannot lump mass at discrete joints; the potential energy is an integral over the continuous mass distribution:

V = \int_0^L \rho(s)\, g\, y\big(\mathbf r(s)\big)\, ds

where $\rho(s)$ is mass per unit length and $\mathbf r(s)$ the backbone curve. Even at rest, solving for the sagged equilibrium shape is a boundary-value problem.

Whiplash. If the tentacle accelerates too quickly, the velocity-quadratic Coriolis/centrifugal term $C(\mathbf q,\dot{\mathbf q})\dot{\mathbf q}$ can drive violent, uncontrolled oscillations — the same physics that makes a cracked whip supersonic. To prevent the octopus from destroying its surroundings (or itself), our controllers evaluate the (discretized) equations of motion in real time, predicting and damping these oscillations before they grow. This two-time-scale structure — slow heavy base, fast light tip — is exactly the setting for singular perturbation, the subject of Part 6.

In Part 5 we step away from clean equations and ask what happens when these models meet real, imperfect hardware.

Further reading: Featherstone, "Rigid Body Dynamics Algorithms" (2008) for RNEA/ABA; Till, Aloi & Rucker (2019), "Real-time dynamics of soft and continuum robots based on Cosserat rod models"; continuum-robot rod survey (arXiv:2407.05886).