[Robotics] Manipulator Dynamics: Newton-Euler and Lagrangian Formulations

Kinematics tells us where the links of a manipulator are and how they move geometrically. Dynamics, by contrast, tells us why they move—the relationship between forces and torques applied at the joints and the accelerations they produce. In every practical control problem, from gravity compensation in a teach-mode arm to optimal trajectory tracking on a torque-controlled humanoid, the equations of motion sit at the heart of the controller. This post derives those equations for an $n$-link serial manipulator, presenting both the energy-based Lagrangian formulation and the recursive Newton-Euler algorithm, and shows that the two yield the same compact form

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

The development closely follows Spong et al. (Chapter 7) and Lynch & Park (Chapter 8). The first emphasizes the Lagrangian route, which is most natural for the model-based controllers we study in later posts; the second emphasizes a modern body-frame Newton-Euler recursion built on twists and wrenches.

Forward and Inverse Dynamics

Before deriving anything we should be clear about the two distinct problems that “manipulator dynamics” can refer to.

$\color{blue}{\mathbf{Definition.}}$ Forward and inverse dynamics
Let $\mathbf{q} \in \mathbb{R}^n$ denote the vector of joint coordinates of an $n$-link serial manipulator, $\dot{\mathbf{q}}$ and $\ddot{\mathbf{q}}$ the joint velocities and accelerations, and $\boldsymbol{\tau} \in \mathbb{R}^n$ the vector of generalized forces (joint torques for revolute joints, forces for prismatic joints). - The **forward dynamics** problem is: given $(\mathbf{q}, \dot{\mathbf{q}}, \boldsymbol{\tau})$, compute the resulting acceleration $\ddot{\mathbf{q}}$. This is the equation we integrate in a physics simulator. - The **inverse dynamics** problem is: given a desired motion $(\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}})$, compute the joint torques $\boldsymbol{\tau}$ that realize it. This is the workhorse of model-based control (e.g. computed-torque control).

Both problems are answered by the same equation of motion. Whether we read it left-to-right or right-to-left determines which problem we solve.

Lagrangian Mechanics Refresher

The Lagrangian approach is an energy-based derivation that, once set up, mechanically produces the equations of motion without the bookkeeping of constraint forces that plagues a free-body analysis.

Let $\mathcal{K}(\mathbf{q}, \dot{\mathbf{q}})$ denote the total kinetic energy of the system and $\mathcal{P}(\mathbf{q})$ the total potential energy. The Lagrangian is

\[\mathcal{L}(\mathbf{q}, \dot{\mathbf{q}}) = \mathcal{K}(\mathbf{q}, \dot{\mathbf{q}}) - \mathcal{P}(\mathbf{q}).\]

For a system of $n$ generalized coordinates with no nonholonomic constraints, the Euler-Lagrange equations are

\[\frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \dot{q}_i} - \frac{\partial \mathcal{L}}{\partial q_i} = \tau_i, \qquad i = 1, \ldots, n,\]

where $\tau_i$ is the generalized force conjugate to $q_i$. For a serial manipulator the generalized coordinates are exactly the joint variables, and the generalized forces are exactly the joint actuator torques (plus any external wrenches mapped to joint space through the Jacobian).

$\color{red}{\mathbf{Theorem.}}$ Euler-Lagrange equations
Among all paths $\mathbf{q}(t)$ that connect two fixed configurations in a fixed time interval, the path actually realized by a holonomic mechanical system makes the action integral $\int_{t_0}^{t_1} \mathcal{L}\, dt$ stationary. The necessary conditions for stationarity are precisely the Euler-Lagrange equations above. When non-conservative generalized forces $\tau_i$ act on the system, they appear as a source term on the right-hand side.

Kinetic Energy of a Rigid Body

To compute the manipulator’s kinetic energy we first recall the expression for a single rigid body.

Let body $i$ have mass $m_i$, and let $\mathbf{v}_{c_i} \in \mathbb{R}^3$ be the linear velocity of its center of mass expressed in an inertial frame. Let $\boldsymbol{\omega}_i \in \mathbb{R}^3$ be the angular velocity of the body and $\mathbf{I}_i \in \mathbb{R}^{3 \times 3}$ its inertia tensor about the center of mass expressed in the same frame in which $\boldsymbol{\omega}_i$ is written. Then the kinetic energy decomposes cleanly as

\[\mathcal{K}_i = \tfrac{1}{2} m_i\, \mathbf{v}_{c_i}^T \mathbf{v}_{c_i} + \tfrac{1}{2}\, \boldsymbol{\omega}_i^T \mathbf{I}_i\, \boldsymbol{\omega}_i.\]

The first term is the translational kinetic energy of a point mass located at the centroid; the second is the rotational kinetic energy of the body spinning about the centroid. The inertia tensor is symmetric and positive definite. If $\mathbf{R}i \in SO(3)$ is the rotation of the body frame relative to the inertial frame and $\mathbf{I}_i^{b} = \mathrm{diag}(I{xx}, I_{yy}, I_{zz})$ is the (constant) inertia in the body frame about a principal-axis basis at the centroid, then the inertia expressed in the inertial frame is

\[\mathbf{I}_i(\mathbf{q}) = \mathbf{R}_i(\mathbf{q})\, \mathbf{I}_i^{b}\, \mathbf{R}_i(\mathbf{q})^T,\]

which is configuration-dependent because $\mathbf{R}_i$ depends on $\mathbf{q}$.

Kinetic Energy of a Link via the Jacobian

For link $i$ of the manipulator, both $\mathbf{v}_{c_i}$ and $\boldsymbol{\omega}_i$ are linear in $\dot{\mathbf{q}}$. They are obtained from the body Jacobian of the centroid of link $i$:

\[\mathbf{v}_{c_i} = \mathbf{J}_{v_i}(\mathbf{q})\, \dot{\mathbf{q}}, \qquad \boldsymbol{\omega}_i = \mathbf{J}_{\omega_i}(\mathbf{q})\, \dot{\mathbf{q}},\]

where $\mathbf{J}{v_i}, \mathbf{J}{\omega_i} \in \mathbb{R}^{3 \times n}$ are the translational and rotational parts of the geometric Jacobian of the centroid of link $i$, with zero columns beyond column $i$ (later joints cannot move earlier links). Substituting,

\[\mathcal{K}_i = \tfrac{1}{2} \dot{\mathbf{q}}^T \Big[ m_i\, \mathbf{J}_{v_i}^T \mathbf{J}_{v_i} + \mathbf{J}_{\omega_i}^T \mathbf{R}_i \mathbf{I}_i^{b} \mathbf{R}_i^T \mathbf{J}_{\omega_i} \Big]\, \dot{\mathbf{q}}.\]

The Mass Matrix M(q)

Summing the contributions of all $n$ links,

\[\mathcal{K}(\mathbf{q}, \dot{\mathbf{q}}) = \tfrac{1}{2} \dot{\mathbf{q}}^T \mathbf{M}(\mathbf{q})\, \dot{\mathbf{q}}, \qquad \mathbf{M}(\mathbf{q}) := \sum_{i=1}^n \Big[ m_i\, \mathbf{J}_{v_i}^T \mathbf{J}_{v_i} + \mathbf{J}_{\omega_i}^T \mathbf{R}_i \mathbf{I}_i^{b} \mathbf{R}_i^T \mathbf{J}_{\omega_i} \Big].\]

This $n \times n$ matrix is the mass matrix (also called the inertia matrix or generalized inertia) of the manipulator. Its properties are crucial in both analysis and control.

$\color{green}{\mathbf{Property.}}$ Symmetry and positive definiteness of $\mathbf{M}(\mathbf{q})$
For every configuration $\mathbf{q}$, the mass matrix $\mathbf{M}(\mathbf{q})$ is 1. **symmetric**: $\mathbf{M}(\mathbf{q}) = \mathbf{M}(\mathbf{q})^T$, and 2. **positive definite**: $\dot{\mathbf{q}}^T \mathbf{M}(\mathbf{q})\, \dot{\mathbf{q}} > 0$ for every $\dot{\mathbf{q}} \neq \mathbf{0}$. The first property follows from the symmetry of each summand. The second follows because the kinetic energy $\frac{1}{2} \dot{\mathbf{q}}^T \mathbf{M}(\mathbf{q})\, \dot{\mathbf{q}}$ is strictly positive whenever the manipulator is in motion. Consequently $\mathbf{M}(\mathbf{q})$ is always invertible, and the forward dynamics problem $\ddot{\mathbf{q}} = \mathbf{M}^{-1}(\boldsymbol{\tau} - \mathbf{C}\dot{\mathbf{q}} - \mathbf{g})$ is well-posed.

The diagonal entries $M_{ii}(\mathbf{q})$ are the effective inertias seen by joint $i$ when all other joints are locked; the off-diagonal entries $M_{ij}(\mathbf{q}) = M_{ji}(\mathbf{q})$ couple accelerations between joints $i$ and $j$.

Potential Energy and the Gravity Vector

Gravity is the only configuration-dependent potential force we typically include in a manipulator model. (Joint springs, when present, contribute additional terms.) Let $\mathbf{g}0 \in \mathbb{R}^3$ be the gravitational acceleration vector expressed in the inertial frame (e.g. $\mathbf{g}_0 = [0, 0, -9.81]^T\, \text{m/s}^2$). If $\mathbf{p}{c_i}(\mathbf{q})$ is the position of the centroid of link $i$, then the gravitational potential energy of link $i$ is

\[\mathcal{P}_i(\mathbf{q}) = -m_i\, \mathbf{g}_0^T \mathbf{p}_{c_i}(\mathbf{q}),\]

and the total potential energy is $\mathcal{P}(\mathbf{q}) = \sum_{i=1}^n \mathcal{P}_i(\mathbf{q})$. The corresponding generalized force vector (i.e. the joint torques that must be supplied to hold the manipulator stationary against gravity) is

\[\mathbf{g}(\mathbf{q}) := \frac{\partial \mathcal{P}(\mathbf{q})}{\partial \mathbf{q}} \in \mathbb{R}^n.\]

Note carefully the notation: $\mathbf{g}_0$ is the $3 \times 1$ acceleration of gravity (a physical constant), while $\mathbf{g}(\mathbf{q})$ is the $n \times 1$ joint-space gravity torque vector (a configuration-dependent quantity). Both Spong and Lynch & Park use the same symbol; context disambiguates.

Equations of Motion

We now apply the Euler-Lagrange equations to the Lagrangian $\mathcal{L} = \frac{1}{2} \dot{\mathbf{q}}^T \mathbf{M}(\mathbf{q})\, \dot{\mathbf{q}} - \mathcal{P}(\mathbf{q})$. Because $\mathcal{P}$ is independent of $\dot{\mathbf{q}}$,

\[\frac{\partial \mathcal{L}}{\partial \dot{q}_i} = \sum_{j=1}^n M_{ij}(\mathbf{q})\, \dot{q}_j,\]

and therefore

\[\frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \dot{q}_i} = \sum_{j=1}^n M_{ij}(\mathbf{q})\, \ddot{q}_j + \sum_{j=1}^n \sum_{k=1}^n \frac{\partial M_{ij}(\mathbf{q})}{\partial q_k}\, \dot{q}_k\, \dot{q}_j.\]

The remaining derivative of $\mathcal{L}$ is

\[\frac{\partial \mathcal{L}}{\partial q_i} = \tfrac{1}{2} \sum_{j=1}^n \sum_{k=1}^n \frac{\partial M_{jk}(\mathbf{q})}{\partial q_i}\, \dot{q}_j\, \dot{q}_k - \frac{\partial \mathcal{P}(\mathbf{q})}{\partial q_i}.\]

Assembling the Euler-Lagrange equations component by component and grouping yields, for each $i$,

\[\sum_{j=1}^n M_{ij}(\mathbf{q})\, \ddot{q}_j + \sum_{j=1}^n \sum_{k=1}^n \Gamma_{ijk}(\mathbf{q})\, \dot{q}_j\, \dot{q}_k + g_i(\mathbf{q}) = \tau_i,\]

where the coefficients $\Gamma_{ijk}$, called the Christoffel symbols of the first kind associated with the metric $\mathbf{M}(\mathbf{q})$, are

\[\Gamma_{ijk}(\mathbf{q}) := \tfrac{1}{2} \left( \frac{\partial M_{ij}}{\partial q_k} + \frac{\partial M_{ik}}{\partial q_j} - \frac{\partial M_{jk}}{\partial q_i} \right).\]

In matrix-vector form,

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

with the Coriolis/centrifugal matrix $\mathbf{C}(\mathbf{q}, \dot{\mathbf{q}}) \in \mathbb{R}^{n \times n}$ defined entrywise by

\[C_{ij}(\mathbf{q}, \dot{\mathbf{q}}) := \sum_{k=1}^n \Gamma_{ijk}(\mathbf{q})\, \dot{q}_k = \sum_{k=1}^n \tfrac{1}{2} \left( \frac{\partial M_{ij}}{\partial q_k} + \frac{\partial M_{ik}}{\partial q_j} - \frac{\partial M_{jk}}{\partial q_i} \right) \dot{q}_k.\]

The terms involving products of distinct joint velocities $\dot{q}_j \dot{q}_k$ with $j \neq k$ are called Coriolis terms; those involving $\dot{q}_j^2$ are centrifugal terms. They vanish only when the manipulator is at rest.

Skew-Symmetry of $\dot{\mathbf{M}} - 2\mathbf{C}$

The Christoffel definition of $\mathbf{C}$ is not the only valid choice—any matrix $\mathbf{C}’$ satisfying $\mathbf{C}’(\mathbf{q}, \dot{\mathbf{q}})\, \dot{\mathbf{q}} = \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}})\, \dot{\mathbf{q}}$ would yield the same equations of motion. The Christoffel choice is privileged because it produces a deep structural property.

$\color{green}{\mathbf{Property.}}$ Skew-symmetry of $\dot{\mathbf{M}} - 2\mathbf{C}$
With $\mathbf{C}$ defined via Christoffel symbols, the matrix $$ \mathbf{N}(\mathbf{q}, \dot{\mathbf{q}}) := \dot{\mathbf{M}}(\mathbf{q}, \dot{\mathbf{q}}) - 2\, \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}}) $$ is skew-symmetric: $\mathbf{N} = -\mathbf{N}^T$, i.e. $\mathbf{x}^T \mathbf{N}(\mathbf{q}, \dot{\mathbf{q}}) \mathbf{x} = 0$ for every $\mathbf{x} \in \mathbb{R}^n$.

The proof is a direct calculation: $\dot{M}{ij} = \sum_k \frac{\partial M{ij}}{\partial q_k} \dot{q}k$, and substituting the Christoffel form of $C{ij}$ shows

\[\dot{M}_{ij} - 2 C_{ij} = \sum_k \left( \frac{\partial M_{ij}}{\partial q_k} - \frac{\partial M_{ij}}{\partial q_k} - \frac{\partial M_{ik}}{\partial q_j} + \frac{\partial M_{jk}}{\partial q_i} \right) \dot{q}_k = \sum_k \left( \frac{\partial M_{jk}}{\partial q_i} - \frac{\partial M_{ik}}{\partial q_j} \right) \dot{q}_k,\]

whose transpose is its negation.

This skew-symmetry is the key passivity property of manipulator dynamics. Differentiating the kinetic energy along trajectories,

\[\frac{d \mathcal{K}}{dt} = \dot{\mathbf{q}}^T \mathbf{M} \ddot{\mathbf{q}} + \tfrac{1}{2} \dot{\mathbf{q}}^T \dot{\mathbf{M}} \dot{\mathbf{q}} = \dot{\mathbf{q}}^T (\boldsymbol{\tau} - \mathbf{g}) - \dot{\mathbf{q}}^T \mathbf{C} \dot{\mathbf{q}} + \tfrac{1}{2} \dot{\mathbf{q}}^T \dot{\mathbf{M}} \dot{\mathbf{q}} = \dot{\mathbf{q}}^T (\boldsymbol{\tau} - \mathbf{g}) + \tfrac{1}{2} \dot{\mathbf{q}}^T \mathbf{N} \dot{\mathbf{q}} = \dot{\mathbf{q}}^T (\boldsymbol{\tau} - \mathbf{g}),\]

which is just the work-energy theorem. Many adaptive and robust manipulator controllers exploit this property in their Lyapunov stability arguments.

Worked Example: Planar 2R Manipulator

To make the abstractions concrete, consider a two-link planar manipulator lying in a vertical plane under gravity. Let the link lengths be $\ell_1, \ell_2$, the link masses $m_1, m_2$ (assumed to be point masses concentrated at the distal end of each link, for simplicity), and the joint angles $q_1, q_2$ measured from the horizontal positive $x$-axis (for joint 1) and from link 1 (for joint 2). Gravity acts along $-\hat{y}$.

The centroid positions are

\[\mathbf{p}_{c_1} = \begin{bmatrix} \ell_1 \cos q_1 \\\\ \ell_1 \sin q_1 \end{bmatrix}, \qquad \mathbf{p}_{c_2} = \begin{bmatrix} \ell_1 \cos q_1 + \ell_2 \cos(q_1 + q_2) \\\\ \ell_1 \sin q_1 + \ell_2 \sin(q_1 + q_2) \end{bmatrix}.\]

Differentiating,

\[\mathbf{v}_{c_1} = \ell_1 \dot{q}_1 \begin{bmatrix} -\sin q_1 \\\\ \cos q_1 \end{bmatrix}, \quad \|\mathbf{v}_{c_1}\|^2 = \ell_1^2 \dot{q}_1^2,\]

and

\[\|\mathbf{v}_{c_2}\|^2 = \ell_1^2 \dot{q}_1^2 + \ell_2^2 (\dot{q}_1 + \dot{q}_2)^2 + 2 \ell_1 \ell_2 \dot{q}_1 (\dot{q}_1 + \dot{q}_2) \cos q_2.\]

Therefore the kinetic energy is

\[\mathcal{K} = \tfrac{1}{2} (m_1 + m_2) \ell_1^2 \dot{q}_1^2 + \tfrac{1}{2} m_2 \ell_2^2 (\dot{q}_1 + \dot{q}_2)^2 + m_2 \ell_1 \ell_2 \dot{q}_1 (\dot{q}_1 + \dot{q}_2) \cos q_2.\]

Identifying terms with $\frac{1}{2} \dot{\mathbf{q}}^T \mathbf{M}(\mathbf{q})\, \dot{\mathbf{q}}$ yields

\[\mathbf{M}(\mathbf{q}) = \begin{bmatrix} (m_1 + m_2) \ell_1^2 + m_2 \ell_2^2 + 2 m_2 \ell_1 \ell_2 \cos q_2 & m_2 \ell_2^2 + m_2 \ell_1 \ell_2 \cos q_2 \\\\ m_2 \ell_2^2 + m_2 \ell_1 \ell_2 \cos q_2 & m_2 \ell_2^2 \end{bmatrix}.\]

Clearly $\mathbf{M}$ is symmetric, and one checks $\det \mathbf{M} = m_2 \ell_2^2 \big[(m_1 + m_2) \ell_1^2 + m_2 \ell_1^2 \ell_2^2 \sin^2 q_2 / \ell_2^2 \cdot \ldots\big] > 0$, confirming positive definiteness.

The Christoffel symbols depend only on $\partial M_{ij} / \partial q_2$, since none of the $M_{ij}$ depend on $q_1$. A short calculation gives

\[\mathbf{C}(\mathbf{q}, \dot{\mathbf{q}}) = m_2 \ell_1 \ell_2 \sin q_2 \begin{bmatrix} -\dot{q}_2 & -(\dot{q}_1 + \dot{q}_2) \\\\ \dot{q}_1 & 0 \end{bmatrix}.\]

The potential energy is

\[\mathcal{P}(\mathbf{q}) = m_1 g \ell_1 \sin q_1 + m_2 g \big[\ell_1 \sin q_1 + \ell_2 \sin(q_1 + q_2)\big],\]

with $g = 9.81$ m/s${}^2$, so

\[\mathbf{g}(\mathbf{q}) = \begin{bmatrix} (m_1 + m_2) g \ell_1 \cos q_1 + m_2 g \ell_2 \cos(q_1 + q_2) \\\\ m_2 g \ell_2 \cos(q_1 + q_2) \end{bmatrix}.\]

Putting these together gives a complete, explicit dynamics model

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

suitable for simulation or computed-torque control. Notice already that the model has terms of three qualitatively different types: products like $\dot{q}_1 \dot{q}_2$ (Coriolis), squares like $\dot{q}_2^2$ (centrifugal), and configuration-dependent gravity terms involving $\cos q_1, \cos(q_1 + q_2)$. Even a two-link arm exhibits the full menagerie of nonlinear couplings that make manipulator control interesting.

The Newton-Euler Recursive Algorithm

The Lagrangian derivation above scales poorly: forming $\mathbf{M}(\mathbf{q})$ symbolically and differentiating to obtain Christoffel symbols becomes intractable for, say, a 7-DOF arm. The Newton-Euler approach instead applies the linear momentum balance $\mathbf{F} = m \dot{\mathbf{v}}_c$ and the angular momentum balance $\boldsymbol{\tau} = \mathbf{I} \dot{\boldsymbol{\omega}} + \boldsymbol{\omega} \times \mathbf{I} \boldsymbol{\omega}$ to each link individually, then propagates the resulting equations along the kinematic chain. The classical recursive algorithm consists of two sweeps.

Outward Recursion: Velocities and Accelerations

Starting from the base (link 0) with known $\boldsymbol{\omega}0 = \mathbf{0}$, $\dot{\boldsymbol{\omega}}_0 = \mathbf{0}$, $\mathbf{v}_0 = \mathbf{0}$, $\dot{\mathbf{v}}_0 = -\mathbf{g}_0$ (so that gravity is bundled into the base acceleration—a useful trick due to Luh, Walker, and Paul), we propagate the angular and linear motion of each link outward toward the end-effector. For a revolute joint $i$ with axis $\hat{\mathbf{z}}{i-1}$, expressed in the frame of link $i$,

\[\boldsymbol{\omega}_i = \mathbf{R}_{i-1}^{i} \boldsymbol{\omega}_{i-1} + \dot{q}_i\, \hat{\mathbf{z}}_i,\] \[\dot{\boldsymbol{\omega}}_i = \mathbf{R}_{i-1}^{i} \dot{\boldsymbol{\omega}}_{i-1} + \mathbf{R}_{i-1}^{i} \boldsymbol{\omega}_{i-1} \times \dot{q}_i\, \hat{\mathbf{z}}_i + \ddot{q}_i\, \hat{\mathbf{z}}_i,\] \[\dot{\mathbf{v}}_i = \mathbf{R}_{i-1}^{i} \dot{\mathbf{v}}_{i-1} + \dot{\boldsymbol{\omega}}_i \times \mathbf{r}_{i-1, i} + \boldsymbol{\omega}_i \times (\boldsymbol{\omega}_i \times \mathbf{r}_{i-1, i}),\] \[\dot{\mathbf{v}}_{c_i} = \dot{\mathbf{v}}_i + \dot{\boldsymbol{\omega}}_i \times \mathbf{r}_{i, c_i} + \boldsymbol{\omega}_i \times (\boldsymbol{\omega}_i \times \mathbf{r}_{i, c_i}),\]

where $\mathbf{R}{i-1}^{i}$ is the rotation from frame $i{-}1$ to frame $i$, $\mathbf{r}{i-1, i}$ is the position of frame $i$’s origin relative to frame $i{-}1$’s origin, and $\mathbf{r}_{i, c_i}$ is the position of link $i$’s centroid relative to frame $i$’s origin. For a prismatic joint the formulas are obtained by replacing $\dot{q}_i \hat{\mathbf{z}}_i$ with the corresponding translational term.

Inward Recursion: Forces and Torques

After the outward sweep, we know $\boldsymbol{\omega}i$, $\dot{\boldsymbol{\omega}}_i$, $\dot{\mathbf{v}}{c_i}$ for every link. Newton’s and Euler’s equations then give the net force and torque acting on each link’s centroid:

\[\mathbf{F}_i = m_i\, \dot{\mathbf{v}}_{c_i}, \qquad \mathbf{N}_i = \mathbf{I}_i \dot{\boldsymbol{\omega}}_i + \boldsymbol{\omega}_i \times \mathbf{I}_i \boldsymbol{\omega}_i.\]

Now we propagate the constraint forces $\mathbf{f}i$ and moments $\mathbf{n}_i$ that link $i{-}1$ exerts on link $i$ *inward*, from the end-effector back to the base. Starting with the end-effector load $\mathbf{f}{n+1}, \mathbf{n}_{n+1}$ (often zero or a known environmental wrench),

\[\mathbf{f}_i = \mathbf{R}_{i+1}^{i} \mathbf{f}_{i+1} + \mathbf{F}_i,\] \[\mathbf{n}_i = \mathbf{R}_{i+1}^{i} \mathbf{n}_{i+1} + \mathbf{r}_{i, c_i} \times \mathbf{F}_i + \mathbf{r}_{i, i+1} \times \mathbf{R}_{i+1}^{i} \mathbf{f}_{i+1} + \mathbf{N}_i.\]

The required joint actuator torque is the projection of $\mathbf{n}_i$ onto the joint axis,

\[\tau_i = \mathbf{n}_i^T \hat{\mathbf{z}}_i \qquad \text{(revolute)}, \qquad \tau_i = \mathbf{f}_i^T \hat{\mathbf{z}}_i \qquad \text{(prismatic)}.\]

This algorithm computes the inverse dynamics in $O(n)$ time per evaluation, in contrast with $O(n^2)$ (or worse, naively $O(n^3)$) for forming and using the closed-form Lagrangian model. It is the standard real-time inverse-dynamics method used in industrial controllers.

Body-Frame Newton-Euler with Twists and Wrenches

Lynch & Park present the same recursion in the modern language of $SE(3)$, using body-frame twists $\mathcal{V}_i \in \mathbb{R}^6$, wrenches $\mathcal{F}_i \in \mathbb{R}^6$, the spatial inertia matrix

\[\mathcal{G}_i = \begin{bmatrix} \mathbf{I}_i^b & \mathbf{0} \\\\ \mathbf{0} & m_i \mathbf{I}_3 \end{bmatrix} \in \mathbb{R}^{6 \times 6},\]

and adjoint maps $\mathrm{Ad}T, \mathrm{ad}{\mathcal{V}}$. The outward step propagates body twists along the chain,

\[\mathcal{V}_i = \mathrm{Ad}_{T_{i, i-1}}(\mathcal{V}_{i-1}) + \mathcal{A}_i\, \dot{q}_i,\] \[\dot{\mathcal{V}}_i = \mathrm{Ad}_{T_{i, i-1}}(\dot{\mathcal{V}}_{i-1}) + \mathrm{ad}_{\mathcal{V}_i}(\mathcal{A}_i \dot{q}_i) + \mathcal{A}_i\, \ddot{q}_i,\]

where $\mathcal{A}_i \in \mathbb{R}^6$ is the screw axis of joint $i$ expressed in the link-$i$ body frame. The inward step propagates wrenches,

\[\mathcal{F}_i = \mathrm{Ad}_{T_{i+1, i}}^T (\mathcal{F}_{i+1}) + \mathcal{G}_i \dot{\mathcal{V}}_i - \mathrm{ad}_{\mathcal{V}_i}^T (\mathcal{G}_i \mathcal{V}_i),\]

and the joint torque is $\tau_i = \mathcal{F}_i^T \mathcal{A}_i$. The combinatorial structure is identical to the classical Luh-Walker-Paul recursion; the screw-theoretic formulation simply replaces $3 \times 3$ rotations and cross products with $6 \times 6$ adjoints, yielding a uniform treatment of revolute, prismatic, helical, and other one-DOF joints.

Equivalence and Practical Choice

The Lagrangian and Newton-Euler formulations describe the same physics and produce the same trajectories when applied to the same manipulator. Their formal equivalence can be proven by recognizing that the Newton-Euler equations are precisely the component-by-component statement of $\frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \dot{q}_i} - \frac{\partial \mathcal{L}}{\partial q_i} = \tau_i$ after the constraint forces between adjacent links are eliminated by the inward recursion.

$\color{red}{\mathbf{Theorem.}}$ Equivalence of formulations
For an $n$-link serial manipulator with conservative external forces, the closed-form equations of motion produced by the Euler-Lagrange equations, $$ \mathbf{M}(\mathbf{q})\, \ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}})\, \dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \boldsymbol{\tau}, $$ and the joint torques $\tau_i$ produced by the Newton-Euler recursive algorithm coincide for every $(\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}}) \in \mathbb{R}^{3n}$.

In practice the two formulations serve different purposes:

• The Lagrangian form gives a symbolic, structured model $(\mathbf{M}, \mathbf{C}, \mathbf{g})$ that is ideal for analysis, control design (computed torque, adaptive control, passivity-based control), and gaining insight into properties such as skew-symmetry, decoupling, and parameter linearity.
• The Newton-Euler algorithm gives an efficient numerical procedure for evaluating inverse dynamics in $O(n)$ time, and underlies most modern simulators (the related Articulated-Body Algorithm of Featherstone extends it to $O(n)$ forward dynamics).

Modern robot dynamics libraries (RBDL, Pinocchio, Drake, MuJoCo, etc.) implement variants of these recursions internally, so that high-level code can call inverse_dynamics(q, qd, qdd) and mass_matrix(q) without ever forming closed-form expressions. Nevertheless, an understanding of the structure $\mathbf{M}\ddot{\mathbf{q}} + \mathbf{C}\dot{\mathbf{q}} + \mathbf{g} = \boldsymbol{\tau}$, together with the properties of $\mathbf{M}(\mathbf{q})$ and the skew-symmetry of $\dot{\mathbf{M}} - 2\mathbf{C}$, is essential for designing the model-based controllers that turn these equations into useful motion.

Reference

[1] Mark W. Spong, Seth Hutchinson, and M. Vidyasagar, Robot Modeling and Control, John Wiley & Sons, 2nd ed., 2020. Chapter 7, “Dynamics.”

[2] Kevin M. Lynch and Frank C. Park, Modern Robotics: Mechanics, Planning, and Control, Cambridge University Press, 2017. Chapter 8, “Dynamics of Open Chains.”

[3] John J. Craig, Introduction to Robotics: Mechanics and Control, Pearson, 3rd ed., 2005. Chapter 6, “Manipulator Dynamics.”

[4] Roy Featherstone, Rigid Body Dynamics Algorithms, Springer, 2008.