Chapter 12: Newton–Euler Recursive Dynamics

Lesson 1: Forward Recursion (velocities/accelerations)

In this lesson we derive and implement the forward recursion of the Newton–Euler algorithm for serial manipulators. Given joint positions, velocities, and accelerations, we propagate angular and linear velocities and accelerations link-by-link, including center-of-mass accelerations, which will be used in the backward recursion (forces and joint torques) in the next lesson.

1. Conceptual Overview of Forward Recursion

The Newton–Euler method views a manipulator as a serial chain of rigid bodies and applies Newton's and Euler's equations to each link. The algorithm has two passes:

  • Forward recursion: propagate kinematic quantities (angular velocity, angular acceleration, linear acceleration, and center-of-mass acceleration) from the base to the end-effector.
  • Backward recursion: propagate forces and moments from the end-effector back to the base to obtain joint forces/torques (next lesson).

Let the generalized coordinates be \( \mathbf{q} = [q_1,\dots,q_n]^\top \) with velocities \( \dot{\mathbf{q}} \) and accelerations \( \ddot{\mathbf{q}} \). The forward recursion computes, for each link \( i=1,\dots,n \),

\[ \boldsymbol{\omega}_i,\quad \boldsymbol{\alpha}_i,\quad \mathbf{a}_i,\quad \mathbf{a}_{c_i}, \]

where \( \boldsymbol{\omega}_i \) and \( \boldsymbol{\alpha}_i \) are the angular velocity and angular acceleration of link \( i \), \( \mathbf{a}_i \) is the linear acceleration of the frame origin, and \( \mathbf{a}_{c_i} \) is the linear acceleration of the center of mass of link \( i \), all expressed in frame \( i \). 

flowchart TD
  Q["Joint data: q, q_dot, q_ddot"] --> T["Kinematic model: R_i, p_i, r_ci, type_i"]
  T --> F1["Forward recursion:
compute omega_i, alpha_i,
a_i, a_ci for i = 1..n"]
  F1 --> B1["Backward recursion:
forces and joint torques
(next lesson)"]
  B1 --> EM["Joint-space EoM:
M(q) q_ddot + C(q,q_dot) q_dot + g(q)"]

Because each link depends only on its parent, the computational cost of the forward recursion scales linearly with the number of links, i.e. order \( \mathcal{O}(n) \), which is crucial for real-time control of high-degree-of-freedom robots.

2. Kinematic Preliminaries and Notation

Consider a fixed-base serial manipulator with frames \( \{0\},\dots,\{n\} \) attached to each link in the usual Denavit–Hartenberg or PoE-consistent way. For each joint \( i \) we define:

\[ \begin{aligned} R_i^{i-1} &\in \mathrm{SO}(3) && \text{rotation from frame } i-1 \text{ to frame } i, \\ \mathbf{p}_i &\in \mathbb{R}^3 && \text{vector from origin of } i-1 \text{ to origin of } i \text{ expressed in frame } i-1, \\ \mathbf{z}_i &= \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} && \text{joint axis in frame } i, \\ q_i,\ \dot{q}_i,\ \ddot{q}_i &&& \text{joint position, velocity, acceleration}, \\ \boldsymbol{\omega}_i,\ \boldsymbol{\alpha}_i &&& \text{angular velocity and acceleration of link } i, \\ \mathbf{a}_i &&& \text{linear acceleration of the origin of frame } i, \\ \mathbf{r}_{c_i} &\in \mathbb{R}^3 && \text{vector from origin of frame } i \text{ to the COM of link } i \text{ in frame } i, \\ \mathbf{a}_{c_i} &&& \text{linear acceleration of the COM of link } i. \end{aligned} \]

Joint type is encoded by \( \sigma_i \in \{0,1\} \), where \( \sigma_i=0 \) denotes a revolute joint and \( \sigma_i=1 \) denotes a prismatic joint.

The base frame \( \{0\} \) is assumed inertially fixed (ground). For a fixed-base manipulator, we initialize:

\[ \boldsymbol{\omega}_0 = \mathbf{0},\qquad \boldsymbol{\alpha}_0 = \mathbf{0},\qquad \mathbf{a}_0 = -\mathbf{g}, \]

where \( \mathbf{g} \) is the gravitational acceleration vector expressed in frame 0 (e.g. \( \mathbf{g} = [0\;0\;9.81]^\top \) if the \( z \)-axis points upward). The choice \( \mathbf{a}_0 = -\mathbf{g} \) incorporates gravity into the forward recursion and simplifies the later force balance.

We use the standard cross-product matrix (hat map) for \( \mathbf{x}\in\mathbb{R}^3 \),

\[ \widehat{\mathbf{x}} = \begin{bmatrix} 0 & -x_3 & x_2 \\ x_3 & 0 & -x_1 \\ -x_2 & x_1 & 0 \end{bmatrix}, \quad \text{so that } \widehat{\mathbf{x}}\,\mathbf{y} = \mathbf{x} \times \mathbf{y}. \]

3. Derivation of Angular Velocity Recursion

Let \( \boldsymbol{\omega}_i^0 \) denote the angular velocity of link \( i \) expressed in the inertial frame \( \{0\} \). For a revolute joint, the angular velocity of link \( i \) is the sum of the parent's angular velocity and the relative joint motion about the axis \( \mathbf{z}_i \):

\[ \boldsymbol{\omega}_i^0 = \boldsymbol{\omega}_{i-1}^0 + \dot{q}_i\,{}^0R_i\,\mathbf{z}_i \quad\text{(revolute joint)}. \]

Expressing angular velocity in the local frame \( i \) via \( \boldsymbol{\omega}_i = R_i^{0}\boldsymbol{\omega}_i^0 \) and using \( R_i^{0} = R_i^{i-1}R_{i-1}^{0} \), we obtain for a revolute joint:

\[ \boldsymbol{\omega}_i = R_i^{i-1}\boldsymbol{\omega}_{i-1} + \dot{q}_i\,\mathbf{z}_i. \]

For a prismatic joint, the relative motion is a pure translation, so there is no additional angular velocity contribution beyond that of the parent:

\[ \boldsymbol{\omega}_i = R_i^{i-1}\boldsymbol{\omega}_{i-1} \quad\text{(prismatic joint)}. \]

Both cases can be combined using the joint-type indicator \( \sigma_i \):

\[ \boldsymbol{\omega}_i = R_i^{i-1}\boldsymbol{\omega}_{i-1} + (1-\sigma_i)\,\dot{q}_i\,\mathbf{z}_i. \]

This is the angular velocity recursion, and it is the first step of the forward pass.

4. Angular Acceleration Recursion

Angular acceleration is the material time derivative of angular velocity in the rotating frame: \( \boldsymbol{\alpha}_i = \dot{\boldsymbol{\omega}}_i \). For a revolute joint we differentiate \( \boldsymbol{\omega}_i = R_i^{i-1}\boldsymbol{\omega}_{i-1} + \dot{q}_i\mathbf{z}_i \):

\[ \begin{aligned} \boldsymbol{\alpha}_i &= \frac{d}{dt}\big(R_i^{i-1}\boldsymbol{\omega}_{i-1}\big) + \ddot{q}_i\mathbf{z}_i \\ &= \dot{R}_i^{i-1}\boldsymbol{\omega}_{i-1} + R_i^{i-1}\boldsymbol{\alpha}_{i-1} + \ddot{q}_i\mathbf{z}_i. \end{aligned} \]

Using the kinematic identity \( \dot{R}_i^{i-1} = \widehat{\boldsymbol{\omega}}_i R_i^{i-1} \), we have

\[ \begin{aligned} \dot{R}_i^{i-1}\boldsymbol{\omega}_{i-1} &= \widehat{\boldsymbol{\omega}}_i R_i^{i-1}\boldsymbol{\omega}_{i-1} \\ &= \boldsymbol{\omega}_i \times \big(R_i^{i-1}\boldsymbol{\omega}_{i-1}\big) = \boldsymbol{\omega}_i \times \big(\boldsymbol{\omega}_i - \dot{q}_i\mathbf{z}_i\big) = \boldsymbol{\omega}_i \times \dot{q}_i\mathbf{z}_i, \end{aligned} \]

where we used \( \boldsymbol{\omega}_i = R_i^{i-1}\boldsymbol{\omega}_{i-1} + \dot{q}_i\mathbf{z}_i \) and the fact that \( \boldsymbol{\omega}_i \times \boldsymbol{\omega}_i = \mathbf{0} \). Hence for a revolute joint:

\[ \boldsymbol{\alpha}_i = R_i^{i-1}\boldsymbol{\alpha}_{i-1} + \ddot{q}_i\mathbf{z}_i + \boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i). \]

For a prismatic joint, since \( \boldsymbol{\omega}_i = R_i^{i-1}\boldsymbol{\omega}_{i-1} \), differentiating yields simply

\[ \boldsymbol{\alpha}_i = R_i^{i-1}\boldsymbol{\alpha}_{i-1} \quad\text{(prismatic joint)}. \]

Again combining both cases:

\[ \boldsymbol{\alpha}_i = R_i^{i-1}\boldsymbol{\alpha}_{i-1} + (1-\sigma_i)\big( \ddot{q}_i\mathbf{z}_i + \boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \big). \]

This recursion captures how a joint acceleration and joint velocity (through the gyroscopic term \( \boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \)) contribute to angular acceleration of the downstream link.

5. Linear and Center-of-Mass Acceleration Recursions

For any rigid body with angular velocity \( \boldsymbol{\omega} \), angular acceleration \( \boldsymbol{\alpha} \), and a point \( B \) at fixed position \( \mathbf{p} \) relative to a point \( A \), the acceleration of \( B \) is

\[ \mathbf{a}_B = \mathbf{a}_A + \boldsymbol{\alpha} \times \mathbf{p} + \boldsymbol{\omega} \times \big(\boldsymbol{\omega} \times \mathbf{p}\big). \]

In the manipulator chain, the origin of frame \( i \) is at position \( \mathbf{p}_i \) relative to frame \( i-1 \). We first compute the acceleration of the origin of frame \( i \) (still expressed in frame \( i \)) by transforming the parent's acceleration and then adding the rotational contributions.

For both joint types we use the expression

\[ \mathbf{a}_i^{\text{base}} = R_i^{i-1}\big( \mathbf{a}_{i-1} + \boldsymbol{\alpha}_{i-1} \times \mathbf{p}_i + \boldsymbol{\omega}_{i-1} \times (\boldsymbol{\omega}_{i-1} \times \mathbf{p}_i) \big), \]

which accounts for the motion of frame \( i \) due to the parent's rigid motion. For revolute joints this is the complete expression:

\[ \mathbf{a}_i = \mathbf{a}_i^{\text{base}} \quad\text{(revolute joint)}. \]

For prismatic joints, there is additional relative translational acceleration along the joint axis. If the slider displacement is \( q_i \) along \( \mathbf{z}_i \), then the relative motion introduces:

  • a direct term \( \ddot{q}_i\mathbf{z}_i \),
  • a Coriolis term \( 2\,\boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \).

\[ \mathbf{a}_i = \mathbf{a}_i^{\text{base}} + \ddot{q}_i\mathbf{z}_i + 2\,\boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \quad\text{(prismatic joint)}. \]

Hence, using \( \sigma_i \) we can write the unified recursion:

\[ \mathbf{a}_i = R_i^{i-1}\big( \mathbf{a}_{i-1} + \boldsymbol{\alpha}_{i-1} \times \mathbf{p}_i + \boldsymbol{\omega}_{i-1} \times (\boldsymbol{\omega}_{i-1} \times \mathbf{p}_i) \big) + \sigma_i\big( \ddot{q}_i\mathbf{z}_i + 2\,\boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \big). \]

Once \( \mathbf{a}_i \) is known, the acceleration of the center of mass of link \( i \) is

\[ \mathbf{a}_{c_i} = \mathbf{a}_i + \boldsymbol{\alpha}_i \times \mathbf{r}_{c_i} + \boldsymbol{\omega}_i \times \big(\boldsymbol{\omega}_i \times \mathbf{r}_{c_i}\big). \]

The quantities \( \mathbf{a}_{c_i} \) will directly appear in the Newton (translational) equations in the backward recursion.

6. Summary of Forward Recursion Equations and Data Flow

For a serial manipulator with base initialization \( \boldsymbol{\omega}_0 = \mathbf{0} \), \( \boldsymbol{\alpha}_0 = \mathbf{0} \), \( \mathbf{a}_0 = -\mathbf{g} \), the forward recursion for each \( i=1,\dots,n \) proceeds as:

\[ \begin{aligned} \boldsymbol{\omega}_i &= R_i^{i-1}\boldsymbol{\omega}_{i-1} + (1-\sigma_i)\,\dot{q}_i\mathbf{z}_i, \\ \boldsymbol{\alpha}_i &= R_i^{i-1}\boldsymbol{\alpha}_{i-1} + (1-\sigma_i)\big( \ddot{q}_i\mathbf{z}_i + \boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \big), \\ \mathbf{a}_i &= R_i^{i-1}\big( \mathbf{a}_{i-1} + \boldsymbol{\alpha}_{i-1} \times \mathbf{p}_i + \boldsymbol{\omega}_{i-1} \times (\boldsymbol{\omega}_{i-1} \times \mathbf{p}_i) \big) \\ &\quad + \sigma_i\big( \ddot{q}_i\mathbf{z}_i + 2\,\boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \big), \\ \mathbf{a}_{c_i} &= \mathbf{a}_i + \boldsymbol{\alpha}_i \times \mathbf{r}_{c_i} + \boldsymbol{\omega}_i \times \big(\boldsymbol{\omega}_i \times \mathbf{r}_{c_i}\big). \end{aligned} \]

Note that the recursion is strictly causal along the kinematic tree: link \( i \) depends only on link \( i-1 \). This leads to the simple data-flow structure illustrated below: 

flowchart TD
  B["Base:
omega_0=0,
alpha_0=0,
a_0=-g"] --> L1["Link 1:
R_1, p_1, r_c1,
type_1, q_1, q_dot1, q_ddot1"]
  L1 --> L2["Link 2:
R_2, p_2, r_c2,
type_2, q_2, q_dot2, q_ddot2"]
  L2 --> Ln["Link n:
...
R_n, p_n, r_cn,
type_n, q_n, q_dotn, q_ddotn"]
  Ln --> OUT["Outputs:
omega_i, alpha_i,
a_i, a_ci for all i"]

7. Compact Twist-Based Formulation

Using the twist representation from earlier chapters, we can pack angular and linear velocity into a six-dimensional spatial velocity (twist):

\[ \mathbf{V}_i = \begin{bmatrix} \boldsymbol{\omega}_i \\[4pt] \mathbf{v}_i \end{bmatrix}, \]

where \( \mathbf{v}_i \) is the linear velocity of the origin of frame \( i \). For joint \( i \) we have a joint screw axis \( \mathbf{S}_i \in \mathbb{R}^6 \) (expressed in a consistent spatial frame). Using the adjoint operator \( \operatorname{Ad}_{i-1}^{i} \) associated with the homogeneous transform from frame \( i-1 \) to frame \( i \), the twist recursion takes the compact form

\[ \mathbf{V}_i = \operatorname{Ad}_{i-1}^{i}\,\mathbf{V}_{i-1} + \mathbf{S}_i\,\dot{q}_i. \]

Differentiating (and using the spatial cross-product operator \( \times_s \) for twists) gives the acceleration recursion:

\[ \dot{\mathbf{V}}_i = \operatorname{Ad}_{i-1}^{i}\,\dot{\mathbf{V}}_{i-1} + \mathbf{S}_i\,\ddot{q}_i + \mathbf{S}_i\,\dot{q}_i \times_s \mathbf{V}_i. \]

Expanding these 6-dimensional equations into angular and linear parts reproduces exactly the recursions derived in the previous sections. In modern geometric mechanics and advanced robotics libraries (e.g. Pinocchio, RBDL), this twist-based formulation is the standard internal representation for kinematic and dynamic algorithms.

8. Python Implementation of Forward Recursion

We now implement the forward recursion in Python using numpy for linear algebra. In practice, robotics libraries such as roboticstoolbox-python provide highly optimized implementations, but it is important to understand the core computations. 


import numpy as np

# Each link is represented as a dict with:
#   "R": 3x3 rotation matrix R_i^{i-1}
#   "p": 3x1 vector p_i (origin i-1 to origin i, in frame i-1)
#   "rc": 3x1 vector r_ci (origin i to COM_i, in frame i)
#   "joint_type": "R" for revolute, "P" for prismatic
#
# The function returns arrays of shape (n, 3)
#   omega[i], alpha[i], a[i], ac[i]

def forward_recursion(links, q, qd, qdd, g=np.array([0.0, 0.0, 9.81])):
    n = len(links)
    z = np.array([0.0, 0.0, 1.0])

    omega = np.zeros((n, 3))
    alpha = np.zeros((n, 3))
    a = np.zeros((n, 3))
    ac = np.zeros((n, 3))

    # Base initialization (fixed base)
    omega_prev = np.zeros(3)
    alpha_prev = np.zeros(3)
    a_prev = -g  # incorporate gravity

    for i, link in enumerate(links):
        R = link["R"]
        p = link["p"]
        rc = link["rc"]
        jt = link["joint_type"]

        qi = q[i]
        qdi = qd[i]
        qddi = qdd[i]

        if jt == "R":
            # Revolute joint
            omega_i = R @ omega_prev + z * qdi
            alpha_i = R @ alpha_prev + z * qddi + np.cross(omega_i, z * qdi)
            a_base = (a_prev
                      + np.cross(alpha_prev, p)
                      + np.cross(omega_prev, np.cross(omega_prev, p)))
            a_i = R @ a_base
        elif jt == "P":
            # Prismatic joint
            omega_i = R @ omega_prev
            alpha_i = R @ alpha_prev
            a_base = (a_prev
                      + np.cross(alpha_prev, p)
                      + np.cross(omega_prev, np.cross(omega_prev, p)))
            a_i = (R @ a_base
                   + z * qddi
                   + np.cross(2.0 * omega_i, z * qdi))
        else:
            raise ValueError("joint_type must be 'R' or 'P'")

        ac_i = (a_i
                + np.cross(alpha_i, rc)
                + np.cross(omega_i, np.cross(omega_i, rc)))

        omega[i, :] = omega_i
        alpha[i, :] = alpha_i
        a[i, :] = a_i
        ac[i, :] = ac_i

        omega_prev = omega_i
        alpha_prev = alpha_i
        a_prev = a_i

    return omega, alpha, a, ac


# Example: tiny 2R planar arm in the x-y plane (z axis out of plane)
def make_planar_2R_example():
    l1, l2 = 1.0, 1.0
    # Simple planar rotations about z with links along x
    def Rz(theta):
        c, s = np.cos(theta), np.sin(theta)
        return np.array([[c, -s, 0.0],
                         [s,  c, 0.0],
                         [0.0, 0.0, 1.0]])

    # Here we assume identity R_i^{i-1} for simplicity and
    # store link lengths in p vectors in frame i-1
    links = []
    links.append({
        "R": np.eye(3),
        "p": np.array([0.0, 0.0, 0.0]),  # base to joint 1
        "rc": np.array([l1 / 2.0, 0.0, 0.0]),
        "joint_type": "R",
    })
    links.append({
        "R": np.eye(3),
        "p": np.array([l1, 0.0, 0.0]),   # joint 1 to joint 2
        "rc": np.array([l2 / 2.0, 0.0, 0.0]),
        "joint_type": "R",
    })
    return links

if __name__ == "__main__":
    links = make_planar_2R_example()
    q = np.array([0.5, 0.3])
    qd = np.array([0.2, 0.1])
    qdd = np.array([0.0, 0.05])

    omega, alpha, a, ac = forward_recursion(links, q, qd, qdd)
    print("omega:\n", omega)
    print("alpha:\n", alpha)
    print("a:\n", a)
    print("ac:\n", ac)

This implementation can be wrapped inside an inverse-dynamics routine that performs backward recursion to compute joint torques, but here we focus on the kinematic forward pass only.

9. C++ Implementation with Eigen

In C++, the Eigen library is widely used in robotics for matrix computations. Lower-level dynamics libraries such as RBDL and Pinocchio exploit similar recursive formulas under the hood. 


#include <vector>
#include <Eigen/Dense>

struct Link {
    Eigen::Matrix3d R;      // R_i^{i-1}
    Eigen::Vector3d p;      // p_i (in frame i-1)
    Eigen::Vector3d rc;     // r_ci (in frame i)
    char joint_type;        // 'R' or 'P'
};

struct ForwardKinematicsResult {
    std::vector<Eigen::Vector3d> omega;
    std::vector<Eigen::Vector3d> alpha;
    std::vector<Eigen::Vector3d> a;
    std::vector<Eigen::Vector3d> ac;
};

ForwardKinematicsResult forwardRecursion(
    const std::vector<Link>& links,
    const Eigen::VectorXd& q,
    const Eigen::VectorXd& qd,
    const Eigen::VectorXd& qdd,
    const Eigen::Vector3d& g = Eigen::Vector3d(0.0, 0.0, 9.81))
{
    const int n = static_cast<int>(links.size());
    Eigen::Vector3d z(0.0, 0.0, 1.0);

    ForwardKinematicsResult res;
    res.omega.assign(n, Eigen::Vector3d::Zero());
    res.alpha.assign(n, Eigen::Vector3d::Zero());
    res.a.assign(n, Eigen::Vector3d::Zero());
    res.ac.assign(n, Eigen::Vector3d::Zero());

    Eigen::Vector3d omega_prev = Eigen::Vector3d::Zero();
    Eigen::Vector3d alpha_prev = Eigen::Vector3d::Zero();
    Eigen::Vector3d a_prev = -g;  // base acceleration

    for (int i = 0; i < n; ++i) {
        const Link& link = links[i];
        const double qi = q(i);
        const double qdi = qd(i);
        const double qddi = qdd(i);

        Eigen::Vector3d omega_i, alpha_i, a_i;

        Eigen::Vector3d a_base =
            a_prev
            + alpha_prev.cross(link.p)
            + omega_prev.cross(omega_prev.cross(link.p));

        if (link.joint_type == 'R') {
            omega_i = link.R * omega_prev + z * qdi;
            alpha_i = link.R * alpha_prev
                      + z * qddi
                      + omega_i.cross(z * qdi);
            a_i = link.R * a_base;
        } else if (link.joint_type == 'P') {
            omega_i = link.R * omega_prev;
            alpha_i = link.R * alpha_prev;
            a_i = link.R * a_base
                  + z * qddi
                  + (2.0 * omega_i).cross(z * qdi);
        } else {
            throw std::runtime_error("joint_type must be 'R' or 'P'");
        }

        Eigen::Vector3d ac_i =
            a_i
            + alpha_i.cross(link.rc)
            + omega_i.cross(omega_i.cross(link.rc));

        res.omega[i] = omega_i;
        res.alpha[i] = alpha_i;
        res.a[i] = a_i;
        res.ac[i] = ac_i;

        omega_prev = omega_i;
        alpha_prev = alpha_i;
        a_prev = a_i;
    }

    return res;
}

This function can be integrated into inverse-dynamics routines, simulation code, or real-time control loops. Robotics middleware (e.g. ROS-based controllers) often wrap such routines inside higher-level interfaces.

10. Java Implementation (Array-Based)

Java is less common than C++/Python for robotics dynamics, but matrix libraries such as EJML can be used. Below is a minimal array-based implementation of the forward recursion, suitable for educational use or integration into a custom Java robotics framework. 


public class Link {
    public double[][] R;   // 3x3 rotation matrix
    public double[] p;     // length-3 vector
    public double[] rc;    // length-3 vector
    public char jointType; // 'R' or 'P'
}

public class ForwardRecursionResult {
    public double[][] omega; // n x 3
    public double[][] alpha; // n x 3
    public double[][] a;     // n x 3
    public double[][] ac;    // n x 3
}

public class NewtonEulerForward {

    private static double[] vec(double x, double y, double z) {
        return new double[]{x, y, z};
    }

    private static double[] add(double[] a, double[] b) {
        return new double[]{a[0] + b[0], a[1] + b[1], a[2] + b[2]};
    }

    private static double[] sub(double[] a, double[] b) {
        return new double[]{a[0] - b[0], a[1] - b[1], a[2] - b[2]};
    }

    private static double[] scale(double[] a, double s) {
        return new double[]{s * a[0], s * a[1], s * a[2]};
    }

    private static double[] cross(double[] a, double[] b) {
        return new double[]{
            a[1] * b[2] - a[2] * b[1],
            a[2] * b[0] - a[0] * b[2],
            a[0] * b[1] - a[1] * b[0]
        };
    }

    private static double[] matVec(double[][] M, double[] v) {
        return new double[]{
            M[0][0] * v[0] + M[0][1] * v[1] + M[0][2] * v[2],
            M[1][0] * v[0] + M[1][1] * v[1] + M[1][2] * v[2],
            M[2][0] * v[0] + M[2][1] * v[1] + M[2][2] * v[2]
        };
    }

    public static ForwardRecursionResult forwardRecursion(
            Link[] links, double[] q, double[] qd, double[] qdd,
            double[] g) {

        int n = links.length;
        double[] z = new double[]{0.0, 0.0, 1.0};

        ForwardRecursionResult res = new ForwardRecursionResult();
        res.omega = new double[n][3];
        res.alpha = new double[n][3];
        res.a = new double[n][3];
        res.ac = new double[n][3];

        double[] omegaPrev = vec(0.0, 0.0, 0.0);
        double[] alphaPrev = vec(0.0, 0.0, 0.0);
        double[] aPrev = scale(g, -1.0); // base acceleration

        for (int i = 0; i < n; ++i) {
            Link L = links[i];
            double qi = q[i];
            double qdi = qd[i];
            double qddi = qdd[i];

            double[] aBase = add(
                add(aPrev, cross(alphaPrev, L.p)),
                cross(omegaPrev, cross(omegaPrev, L.p))
            );

            double[] omega_i, alpha_i, a_i;

            if (L.jointType == 'R') {
                omega_i = add(matVec(L.R, omegaPrev), scale(z, qdi));
                alpha_i = add(
                    matVec(L.R, alphaPrev),
                    add(scale(z, qddi), cross(omega_i, scale(z, qdi)))
                );
                a_i = matVec(L.R, aBase);
            } else if (L.jointType == 'P') {
                omega_i = matVec(L.R, omegaPrev);
                alpha_i = matVec(L.R, alphaPrev);
                a_i = add(
                    matVec(L.R, aBase),
                    add(scale(z, qddi), cross(scale(2.0, omega_i), scale(z, qdi)))
                );
            } else {
                throw new IllegalArgumentException("jointType must be 'R' or 'P'");
            }

            double[] ac_i = add(
                add(a_i, cross(alpha_i, L.rc)),
                cross(omega_i, cross(omega_i, L.rc))
            );

            res.omega[i] = omega_i;
            res.alpha[i] = alpha_i;
            res.a[i] = a_i;
            res.ac[i] = ac_i;

            omegaPrev = omega_i;
            alphaPrev = alpha_i;
            aPrev = a_i;
        }

        return res;
    }
}

This implementation uses only primitive arrays and can be adapted to use matrix libraries (e.g. EJML) if desired.

11. MATLAB/Simulink Implementation

MATLAB, together with the Robotics System Toolbox, is a standard environment for manipulator dynamics. Below is a basic function implementing the forward recursion using MATLAB structs. This function can be wrapped inside a Simulink “MATLAB Function” block for real-time simulation. 


function [omega, alpha, a, ac] = forward_recursion(links, q, qd, qdd, g)
% FORWARD_RECURSION  Newton-Euler forward pass (velocities/accelerations)
%   links(i) has fields:
%     R  (3x3), p (3x1), rc (3x1), joint_type ('R' or 'P')

if nargin < 5
    g = [0; 0; 9.81];
end

n = numel(links);
z = [0; 0; 1];

omega = zeros(3, n);
alpha = zeros(3, n);
a = zeros(3, n);
ac = zeros(3, n);

omega_prev = zeros(3, 1);
alpha_prev = zeros(3, 1);
a_prev = -g;  % base acceleration includes gravity

for i = 1:n
    L = links(i);
    qi = q(i);
    qdi = qd(i);
    qddi = qdd(i);

    a_base = a_prev ...
           + cross(alpha_prev, L.p) ...
           + cross(omega_prev, cross(omega_prev, L.p));

    if L.joint_type == 'R'
        omega_i = L.R * omega_prev + z * qdi;
        alpha_i = L.R * alpha_prev + z * qddi + cross(omega_i, z * qdi);
        a_i = L.R * a_base;
    elseif L.joint_type == 'P'
        omega_i = L.R * omega_prev;
        alpha_i = L.R * alpha_prev;
        a_i = L.R * a_base ...
            + z * qddi ...
            + cross(2 * omega_i, z * qdi);
    else
        error('joint_type must be ''R'' or ''P''');
    end

    rc = L.rc;
    ac_i = a_i ...
         + cross(alpha_i, rc) ...
         + cross(omega_i, cross(omega_i, rc));

    omega(:, i) = omega_i;
    alpha(:, i) = alpha_i;
    a(:, i) = a_i;
    ac(:, i) = ac_i;

    omega_prev = omega_i;
    alpha_prev = alpha_i;
    a_prev = a_i;
end
end

Simulink integration: create a MATLAB Function block, paste the above function body inside, and pass in vectors \( \mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}} \) and a structure array of link parameters (which can be constructed from a URDF or from Robotics System Toolbox rigid body tree).

12. Wolfram Mathematica Implementation

Wolfram Mathematica is well suited for symbolic as well as numeric explorations of Newton–Euler dynamics. Below a simple implementation of the forward recursion is given in a functional style. 


Clear[ForwardRecursion];

(* Each link is an association with keys:
   "R", "p", "rc", "JointType" *)
ForwardRecursion[links_List, q_List, qd_List, qdd_List, 
                 g_: {0., 0., 9.81}] := Module[
  {n = Length[links], z, omega, alpha, a, ac,
   omegaPrev, alphaPrev, aPrev, R, p, rc, jt,
   qi, qdi, qddi, aBase, omegaI, alphaI, aI, acI},
  
  z = {0., 0., 1.};
  omega = ConstantArray[{0., 0., 0.}, n];
  alpha = ConstantArray[{0., 0., 0.}, n];
  a = ConstantArray[{0., 0., 0.}, n];
  ac = ConstantArray[{0., 0., 0.}, n];

  omegaPrev = {0., 0., 0.};
  alphaPrev = {0., 0., 0.};
  aPrev = -g;

  Do[
    R = links[[i, "R"]];
    p = links[[i, "p"]];
    rc = links[[i, "rc"]];
    jt = links[[i, "JointType"]];
    qi = q[[i]];
    qdi = qd[[i]];
    qddi = qdd[[i]];

    aBase = aPrev
      + Cross[alphaPrev, p]
      + Cross[omegaPrev, Cross[omegaPrev, p]];

    If[jt === "R",
      omegaI = R.omegaPrev + qdi z;
      alphaI = R.alphaPrev + qddi z + Cross[omegaI, qdi z];
      aI = R.aBase;,
      (* jt === "P" *)
      omegaI = R.omegaPrev;
      alphaI = R.alphaPrev;
      aI = R.aBase + qddi z + Cross[2. omegaI, qdi z];
    ];

    acI = aI + Cross[alphaI, rc] + Cross[omegaI, Cross[omegaI, rc]];

    omega[[i]] = omegaI;
    alpha[[i]] = alphaI;
    a[[i]] = aI;
    ac[[i]] = acI;

    omegaPrev = omegaI;
    alphaPrev = alphaI;
    aPrev = aI;
    ,
    {i, 1, n}
  ];

  <|"Omega" -> omega, "Alpha" -> alpha, "A" -> a, "Ac" -> ac|>
];

Symbolic parameters for \( R_i^{i-1}, \mathbf{p}_i, \mathbf{r}_{c_i} \) can be used to derive closed-form expressions for accelerations that match those obtained with Lagrangian methods, providing a cross-check at the level of algebraic expressions.

13. Problems and Solutions

Problem 1 (Angular velocity recursion for a 2R planar arm): Consider a planar 2R manipulator in the \( x\text{–}y \) plane with both joints revolute about the \( z \)-axis. The base is fixed with \( \boldsymbol{\omega}_0 = \mathbf{0} \), \( \boldsymbol{\alpha}_0 = \mathbf{0} \). Show that the angular velocities of link 1 and link 2 are \( \boldsymbol{\omega}_1 = \dot{q}_1\mathbf{k} \) and \( \boldsymbol{\omega}_2 = (\dot{q}_1+\dot{q}_2)\mathbf{k} \), where \( \mathbf{k} = [0\;0\;1]^\top \).

Solution:

For link 1, the parent angular velocity is zero, and \( \mathbf{z}_1 = \mathbf{k} \). The recursion gives \( \boldsymbol{\omega}_1 = R_1^{0}\boldsymbol{\omega}_0 + \dot{q}_1\mathbf{z}_1 = \dot{q}_1\mathbf{k} \). For link 2, again \( \mathbf{z}_2 = \mathbf{k} \), and the parent angular velocity is \( \boldsymbol{\omega}_1 = \dot{q}_1\mathbf{k} \). Hence

\[ \boldsymbol{\omega}_2 = R_2^{1}\boldsymbol{\omega}_1 + \dot{q}_2\mathbf{z}_2 = \dot{q}_1\mathbf{k} + \dot{q}_2\mathbf{k} = (\dot{q}_1+\dot{q}_2)\mathbf{k}, \]

which matches the intuitive fact that link 2 rotates about the base with the sum of the two joint rates.

Problem 2 (Linear acceleration of link 1 origin): For the same 2R planar manipulator, assume link 1 has length \( \ell_1 \), its frame origin is at the first joint, and \( \mathbf{p}_1 = \mathbf{0} \). Let \( \mathbf{a}_0 = -\mathbf{g} \) with \( \mathbf{g} = [0\;0\;g]^\top \). Compute \( \mathbf{a}_1 \) explicitly and interpret it physically.

Solution:

Since \( \mathbf{p}_1 = \mathbf{0} \), the recursion simplifies to \( \mathbf{a}_1 = R_1^{0}\mathbf{a}_0 \). If the base frame is chosen so that \( R_1^{0} \) is a pure rotation in the \( x\text{–}y \) plane and gravity is along \( +z \), then \( R_1^{0}\mathbf{g} = \mathbf{g} \). Hence

\[ \mathbf{a}_1 = R_1^{0}(-\mathbf{g}) = -\mathbf{g}, \]

i.e. the origin of link 1 has the same gravitational acceleration as the base. This is consistent with the fact that the joint axis passes through the origin and link rotation does not translate that point.

Problem 3 (COM acceleration and dynamic consistency): For a single revolute link with COM offset \( \mathbf{r}_c = [\ell\;0\;0]^\top \) in its local frame and joint variables \( q,\dot{q},\ddot{q} \), express the COM acceleration \( \mathbf{a}_c \) in terms of \( \boldsymbol{\omega}, \boldsymbol{\alpha}, \mathbf{a} \) and interpret the contributions.

Solution:

Using the general formula, \( \mathbf{a}_c = \mathbf{a} + \boldsymbol{\alpha} \times \mathbf{r}_c + \boldsymbol{\omega} \times (\boldsymbol{\omega} \times \mathbf{r}_c) \), we have

\[ \mathbf{a}_c = \underbrace{\mathbf{a}}_{\text{translation of frame}} + \underbrace{\boldsymbol{\alpha} \times \mathbf{r}_c}_{\text{tangential (due to angular acceleration)}} + \underbrace{\boldsymbol{\omega} \times (\boldsymbol{\omega} \times \mathbf{r}_c)}_{\text{centripetal (due to angular velocity)}}. \]

The second term is tangent to the circular trajectory of the COM and proportional to \( \ddot{q} \), while the third term points toward the joint axis and is proportional to \( \dot{q}^2 \). This decomposition is crucial in interpreting inertial forces in the backward recursion.

Problem 4 (Prismatic Coriolis term): Show that for a prismatic joint, the relative sliding motion yields the Coriolis term \( 2\,\boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \) in the acceleration recursion.

Solution:

Let the relative position of the slider be \( \mathbf{s}_i = q_i\mathbf{z}_i \) in the local joint frame. The relative velocity and acceleration are \( \dot{\mathbf{s}}_i = \dot{q}_i\mathbf{z}_i \) and \( \ddot{\mathbf{s}}_i = \ddot{q}_i\mathbf{z}_i \). The total acceleration of the slider point includes:

  • the parent's rigid-body acceleration acting at the joint location;
  • the linear acceleration \( \ddot{\mathbf{s}}_i \);
  • an additional term \( 2\,\boldsymbol{\omega}_i \times \dot{\mathbf{s}}_i \) from differentiation in a rotating frame.

Thus the extra piece relative to the revolute case is \( \ddot{q}_i\mathbf{z}_i + 2\,\boldsymbol{\omega}_i \times (\dot{q}_i\mathbf{z}_i) \), which is precisely the prismatic contribution in the unified acceleration recursion.

Problem 5 (Algorithmic structure for forward recursion): Sketch a high-level algorithm for the forward recursion that shows which quantities must be known before each step and what is produced at each step.

Solution (algorithmic flow):

flowchart TD
  S["Start:
given robot parameters
and joint trajectories"] --> B["Initialize base:
omega_0, alpha_0, a_0"]
  B --> LOOP["For i = 1..n:
  read R_i, p_i, r_ci, type_i,
  q_i, q_dot_i, q_ddot_i"]
  LOOP --> KIN["Compute:
omega_i, alpha_i,
a_i, a_ci"]
  KIN --> STORE["Store kinematic
results for link i"]
  STORE --> NEXT["i = i + 1
until n"]

This flow emphasizes that all dynamic information in the next lesson (forces, moments, joint torques) will be built on top of the kinematic quantities produced here.

14. Summary

In this lesson we developed the Newton–Euler forward recursion for serial robots. Starting from the base and given joint trajectories \( \mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}} \), we derived recursive formulas for angular velocities and accelerations, linear accelerations of link frames, and center-of-mass accelerations. We showed how revolute and prismatic joints are unified through a joint-type indicator, and we connected the vector formulas to a twist-based geometric formulation. Finally, we implemented the recursion in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica to emphasize both theoretical understanding and practical computation. These kinematic quantities are the foundation for the backward Newton–Euler pass that yields joint forces and torques in the next lesson.

15. References

  1. Luh, J.Y.S., Walker, M.W., & Paul, R.P.C. (1980). On-line computational scheme for mechanical manipulators. Journal of Dynamic Systems, Measurement, and Control, 102(2), 69–76.
  2. Hollerbach, J.M. (1980). A recursive Lagrangian formulation of manipulator dynamics and a comparative study of dynamics formulation complexity. IEEE Transactions on Systems, Man, and Cybernetics, 10(11), 730–736.
  3. Featherstone, R. (1983). The calculation of robot dynamics using articulated-body inertias. International Journal of Robotics Research, 2(1), 13–30.
  4. Khalil, W., & Dombre, E. (1990). General formulation of robot dynamics using recursive Newton–Euler techniques. International Journal of Robotics Research, 9(5), 40–54.
  5. Anigstein, M., & Ebert-Uphoff, I. (1999). About the Newton–Euler method for the calculation of manipulator dynamics. Journal of Robotic Systems, 16(9), 499–516.
  6. Saha, S.K. (1997). A computationally efficient method for the dynamics of serial multibody systems using the Newton–Euler formulation. Multibody System Dynamics, 1(3), 223–242.
  7. Craig, J.J. (2005). Recursive Newton–Euler techniques for robot dynamics. Various journal and book chapters related to robotic manipulators.
  8. Murray, R.M., Li, Z., & Sastry, S.S. (1994). A Mathematical Introduction to Robotic Manipulation. CRC Press (chapters on twists and recursive dynamics).
  9. Müller, A. (2023). Screw and Lie group theory in multibody dynamics – recursive algorithms and equations of motion of tree-topology systems. Multibody System Dynamics, 57(4), 359–390.

Previous

(previous chapter) Lesson 5: Symbolic vs. Numeric Derivation Lab

Next

Lesson 2: Backward Recursion (forces/torques)