In the Newton–Euler formulation for an \( n \)-link serial manipulator, the forward recursion of Lesson 1 provides, for each link \( i \), the angular velocity \( \boldsymbol{\omega}_i \), angular acceleration \( \boldsymbol{\alpha}_i \), and linear acceleration of the center of mass \( \mathbf{a}_i^c \), all expressed in a body-fixed frame \( \{i\} \). The backward recursion uses these kinematic quantities and inertial parameters to compute:

• The net inertial force \( \mathbf{F}_i \) and moment \( \mathbf{N}_i \) of each link.
• The internal wrench \( (\mathbf{f}_i, \mathbf{n}_i) \) transmitted across each joint.
• The generalized joint effort \( \tau_i \) for the \( i \)-th joint.

The overall structure is:

\[ \text{Forward sweep:} \quad \{q_i, \dot{q}_i, \ddot{q}_i\}_{i=1}^n \;\Rightarrow\; \{\boldsymbol{\omega}_i, \boldsymbol{\alpha}_i, \mathbf{a}_i^c\}_{i=1}^n \]

\[ \text{Backward sweep:} \quad \{\boldsymbol{\omega}_i, \boldsymbol{\alpha}_i, \mathbf{a}_i^c\}_{i=1}^n \;\Rightarrow\; \{\mathbf{F}_i, \mathbf{N}_i, \mathbf{f}_i, \mathbf{n}_i, \tau_i\}_{i=1}^n. \]

The backward recursion is essentially a subtree equilibrium computation: starting from the terminal body (possibly with external wrench at the end-effector), we propagate wrenches inward, enforcing Newton–Euler laws on each link.

Consider link \( i \) with mass \( m_i \), center of mass (CoM) located at vector \( \mathbf{r}_i^c \) from frame origin \( \{i\} \), and inertia tensor about its CoM in frame \( \{i\} \), \( \mathbf{I}_i \in \mathbb{R}^{3\times 3} \). From rigid-body mechanics, the inertial (net) force and moment acting on the CoM are:

\[ \mathbf{F}_i = m_i \,\mathbf{a}_i^c \]

\[ \mathbf{N}_i = \mathbf{I}_i \,\boldsymbol{\alpha}_i + \boldsymbol{\omega}_i \times (\mathbf{I}_i \,\boldsymbol{\omega}_i). \]

Here \( \mathbf{a}_i^c \), \( \boldsymbol{\omega}_i \), \( \boldsymbol{\alpha}_i \) are obtained from the forward recursion (Lesson 1). These quantities represent the inertial reaction of link \( i \) to its motion; they will be balanced by internal and external wrenches in the backward recursion.

Denote by \( \mathbf{f}_i \) the force transmitted from link \( i \) to link \( i-1 \) through joint \( i \), and \( \mathbf{n}_i \) the corresponding moment at the origin of frame \( \{i\} \). We adopt the sign convention that \( \mathbf{f}_i \) and \( \mathbf{n}_i \) act on link \( i \) from its parent.

If the end-effector is subject to a known external wrench \( (\mathbf{f}_{n+1}^{\text{ext}}, \mathbf{n}_{n+1}^{\text{ext}}) \), that wrench is treated as acting on a fictitious link \( n+1 \) attached to link \( n \). In many models we set \( \mathbf{f}_{n+1}^{\text{ext}} = \mathbf{0},\; \mathbf{n}_{n+1}^{\text{ext}} = \mathbf{0} \), corresponding to free end conditions.

Let \( \mathcal{T}_i \) be the subtree consisting of link \( i \) and all its descendants. The total inertial wrench of the subtree, expressed at frame \( \{i\} \), must be balanced by external wrenches and the internal reaction at joint \( i \). For a serial chain, each link has at most one child (link \( i+1 \)), so the recursion is one-dimensional along the chain.

Let \( \;^i\mathbf{R}_{i+1} \) be the rotation matrix that maps vectors from frame \( \{i+1\} \) to frame \( \{i\} \), and \( \mathbf{p}_i \) the vector from origin of frame \( \{i\} \) to origin of frame \( \{i+1\} \), both expressed in \( \{i\} \). The child wrench at joint \( i+1 \), expressed in frame \( \{i+1\} \), is \( (\mathbf{f}_{i+1}, \mathbf{n}_{i+1}) \); expressed in \( \{i\} \), it is:

\[ \mathbf{f}_{i+1}^{(i)} = \;^i\mathbf{R}_{i+1}\,\mathbf{f}_{i+1},\quad \mathbf{n}_{i+1}^{(i)} = \;^i\mathbf{R}_{i+1}\,\mathbf{n}_{i+1} + \mathbf{p}_i \times \;^i\mathbf{R}_{i+1}\,\mathbf{f}_{i+1}. \]

The equilibrium of link \( i \) (Newton–Euler) in frame \( \{i\} \) is:

\[ \mathbf{f}_i + \mathbf{f}_{i+1}^{(i)} + \mathbf{f}_i^{\text{ext}} = \mathbf{F}_i, \]

\[ \mathbf{n}_i + \mathbf{n}_{i+1}^{(i)} + \mathbf{n}_i^{\text{ext}} = \mathbf{N}_i + \mathbf{r}_i^c \times \mathbf{F}_i, \]

where \( \mathbf{f}_i^{\text{ext}}, \mathbf{n}_i^{\text{ext}} \) denote any external wrenches acting directly on link \( i \) (e.g., gravity can be modeled as inertial force via \( \mathbf{a}_i^c \), so often the only external load is at the end-effector).

Solving for \( \mathbf{f}_i \) and \( \mathbf{n}_i \) gives the backward recursion:

\[ \mathbf{f}_i = \mathbf{F}_i + \mathbf{f}_i^{\text{ext}} + \;^i\mathbf{R}_{i+1}\,\mathbf{f}_{i+1}, \]

\[ \mathbf{n}_i = \mathbf{N}_i + \mathbf{r}_i^c \times \mathbf{F}_i + \mathbf{n}_i^{\text{ext}} + \;^i\mathbf{R}_{i+1}\,\mathbf{n}_{i+1} + \mathbf{p}_i \times \;^i\mathbf{R}_{i+1}\,\mathbf{f}_{i+1}. \]

At the terminal end \( i = n \), we use the boundary condition:

\[ \mathbf{f}_{n+1} = \mathbf{f}_{n+1}^{\text{ext}}, \quad \mathbf{n}_{n+1} = \mathbf{n}_{n+1}^{\text{ext}}, \]

which is often chosen as zero for a free end-effector or as the known reaction wrench when the end-effector is in contact with the environment (but still treated as a known, not an unknown for this recursive pass).

The internal wrench \( (\mathbf{f}_i, \mathbf{n}_i) \) at joint \( i \) must be projected onto the joint’s motion subspace to obtain the generalized effort \( \tau_i \). Let \( \mathbf{z}_i \) be the unit joint axis expressed in frame \( \{i\} \).

Revolute joint. The generalized coordinate is an angle \( q_i \), and the power equality \( \tau_i \dot{q}_i = \mathbf{n}_i^{\top}\boldsymbol{\omega}_i^J \) for joint-relative angular velocity \( \boldsymbol{\omega}_i^J = \dot{q}_i\mathbf{z}_i \) yields:

\[ \tau_i = \mathbf{z}_i^{\top}\mathbf{n}_i. \]

Prismatic joint. The generalized coordinate is a displacement \( q_i \), and the power equality \( \tau_i \dot{q}_i = \mathbf{f}_i^{\top}\mathbf{v}_i^J \) for joint-relative linear velocity \( \mathbf{v}_i^J = \dot{q}_i\mathbf{z}_i \) yields:

\[ \tau_i = \mathbf{z}_i^{\top}\mathbf{f}_i. \]

In a compact notation using the joint motion subspace \( \mathbf{s}_i \in \mathbb{R}^6 \) and the spatial wrench \( \mathbf{w}_i \in \mathbb{R}^6 \), the joint effort is \( \tau_i = \mathbf{s}_i^{\top}\mathbf{w}_i \), but full spatial vector theory will only be systematically introduced in later chapters on inertia parameterization.

The Newton–Euler algorithm achieves \( \mathcal{O}(n) \) complexity for an \( n \)-link chain. The backward recursion uses only local operations at each link: evaluating \( \mathbf{F}_i, \mathbf{N}_i \) and propagating wrenches from child to parent. A high-level procedure is:

1. Run forward recursion (Lesson 1) to obtain \( \boldsymbol{\omega}_i, \boldsymbol{\alpha}_i, \mathbf{a}_i^c \).
2. For each link, compute \( \mathbf{F}_i, \mathbf{N}_i \) from inertial parameters.
3. Initialize terminal wrenches \( \mathbf{f}_{n+1}, \mathbf{n}_{n+1} \) from external load.
4. For \( i = n, n-1, \dots, 1 \), compute \( \mathbf{f}_i, \mathbf{n}_i, \tau_i \).

flowchart TD
  S["Inputs: q, qdot, qddot, transforms, inertial params, end-effector wrench"] --> FWD["Forward sweep: compute omega_i, alpha_i, a_i_c"]
  FWD --> INERT["For each link i: compute F_i = m_i a_i_c and N_i = I_i alpha_i + omega_i x (I_i omega_i)"]
  INERT --> INIT["Set f_(n+1), n_(n+1) from external load (often zero)"]
  INIT --> BWD["For i = n down to 1: propagate f_i, n_i using subtree equilibrium"]
  BWD --> TAU["Project wrench on joint axis to get tau_i"]
  TAU --> OUT["Outputs: joint efforts tau_1,...,tau_n and internal wrenches"]
        

Each recursion step uses only matrix–vector products and cross products. Assuming constant-time operations per link, the total cost scales linearly with the number of links, in contrast to naive Lagrangian formulations that may scale as \( \mathcal{O}(n^3) \) if implemented without exploiting sparsity.

In Python, popular robotics libraries such as roboticstoolbox-python (Corke) provide Newton–Euler algorithms on top of a SerialLink or rigid-body-tree abstraction. Here we implement the core backward recursion from scratch using numpy, assuming the forward recursion has already supplied \( \mathbf{F}_i \) and \( \mathbf{N}_i \) for each link.

import numpy as np

def newton_euler_backward(F_list, N_list, R_ip1_i_list, p_i_list,
                          z_list, joint_type_list,
                          f_tip=None, n_tip=None):
    """
    Backward recursion for a serial chain.

    Parameters
    ----------
    F_list : list of (3,) arrays
        Inertial forces F_i at link CoM in frame {i}.
    N_list : list of (3,) arrays
        Inertial moments N_i about frame origin {i}.
    R_ip1_i_list : list of (3,3) arrays
        Rotation matrices ^i R_{i+1} (from frame {i+1} to {i}).
        Length n, but R_ip1_i_list[-1] is used only if you model a tip link.
    p_i_list : list of (3,) arrays
        Vectors p_i from origin of frame {i} to origin of frame {i+1}, in frame {i}.
    z_list : list of (3,) arrays
        Joint axes z_i expressed in frame {i}.
    joint_type_list : list of str
        "R" for revolute, "P" for prismatic.
    f_tip, n_tip : (3,) arrays or None
        External wrench at the terminal body (expressed in last frame).

    Returns
    -------
    f_list, n_list, tau_list : lists of arrays/scalars
        Internal forces, moments, and generalized efforts at each joint.
    """
    n = len(F_list)
    assert len(N_list) == n
    assert len(R_ip1_i_list) == n
    assert len(p_i_list) == n
    assert len(z_list) == n
    assert len(joint_type_list) == n

    if f_tip is None:
        f_next = np.zeros(3)
    else:
        f_next = np.asarray(f_tip).reshape(3)
    if n_tip is None:
        n_next = np.zeros(3)
    else:
        n_next = np.asarray(n_tip).reshape(3)

    f_list = [np.zeros(3) for _ in range(n)]
    n_list = [np.zeros(3) for _ in range(n)]
    tau_list = [0.0 for _ in range(n)]

    # Backward sweep: i = n-1, ..., 0
    for i in range(n - 1, -1, -1):
        R_ip1_i = R_ip1_i_list[i]
        p_i = p_i_list[i]
        F_i = F_list[i]
        N_i = N_list[i]
        z_i = z_list[i]
        jt = joint_type_list[i]

        # Propagate child wrench into frame {i}
        f_child_i = R_ip1_i @ f_next
        n_child_i = R_ip1_i @ n_next + np.cross(p_i, f_child_i)

        # Internal wrench at joint i
        f_i = F_i + f_child_i
        # r_i_c is already accounted for in N_i via N_i + r_i_c x F_i,
        # or you can add that explicitly here depending on how you define N_i.
        n_i = N_i + n_child_i

        f_list[i] = f_i
        n_list[i] = n_i

        if jt == "R":
            tau_i = float(z_i.dot(n_i))
        elif jt == "P":
            tau_i = float(z_i.dot(f_i))
        else:
            raise ValueError("joint_type_list must contain only 'R' or 'P'")

        tau_list[i] = tau_i

        # Wrench at this joint becomes "child" for previous joint
        f_next = f_i
        n_next = n_i

    return f_list, n_list, tau_list

# Example usage: 2-link planar manipulator (simplified)
if __name__ == "__main__":
    n = 2
    F_list = [np.array([1.0, 0.0, 0.0]),
              np.array([0.5, 0.0, 0.0])]
    N_list = [np.array([0.0, 0.0, 0.2]),
              np.array([0.0, 0.0, 0.1])]
    R_ip1_i_list = [np.eye(3), np.eye(3)]
    p_i_list = [np.array([0.5, 0.0, 0.0]),
                np.array([0.3, 0.0, 0.0])]
    z_list = [np.array([0.0, 0.0, 1.0]),
              np.array([0.0, 0.0, 1.0])]
    joint_type_list = ["R", "R"]

    f_list, n_list, tau_list = newton_euler_backward(
        F_list, N_list, R_ip1_i_list, p_i_list,
        z_list, joint_type_list,
        f_tip=np.zeros(3), n_tip=np.zeros(3)
    )

    print("Joint torques:", tau_list)
      

In practice, the forward recursion provides \( \mathbf{F}_i, \mathbf{N}_i \) consistently with the modeling of gravity and any CoM offsets. A robotics toolbox can be used to validate this implementation by comparing with its built-in inverse dynamics routine.

In C++, widely used libraries such as Orocos KDL and RBDL implement Newton–Euler inverse dynamics efficiently and are integrated into many robotics frameworks (e.g., ROS). Below is a minimal illustration using Eigen for vector/matrix operations; it mirrors the Python logic at a lower level.

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

struct LinkDynamics {
    Eigen::Vector3d F;   // inertial force
    Eigen::Vector3d N;   // inertial moment
    Eigen::Matrix3d R_ip1_i; // ^i R_(i+1)
    Eigen::Vector3d p_i; // from frame i to i+1 in frame i
    Eigen::Vector3d z;   // joint axis
    char joint_type;     // 'R' or 'P'
};

void newtonEulerBackward(const std::vector<LinkDynamics>& links,
                         const Eigen::Vector3d& f_tip,
                         const Eigen::Vector3d& n_tip,
                         std::vector<Eigen::Vector3d>& f_list,
                         std::vector<Eigen::Vector3d>& n_list,
                         std::vector<double>& tau_list)
{
    int n = static_cast<int>(links.size());
    f_list.assign(n, Eigen::Vector3d::Zero());
    n_list.assign(n, Eigen::Vector3d::Zero());
    tau_list.assign(n, 0.0);

    Eigen::Vector3d f_next = f_tip;
    Eigen::Vector3d n_next = n_tip;

    for (int i = n - 1; i >= 0; --i) {
        const auto& L = links[i];

        Eigen::Vector3d f_child_i = L.R_ip1_i * f_next;
        Eigen::Vector3d n_child_i = L.R_ip1_i * n_next
                                  + L.p_i.cross(f_child_i);

        Eigen::Vector3d f_i = L.F + f_child_i;
        Eigen::Vector3d n_i = L.N + n_child_i;

        f_list[i] = f_i;
        n_list[i] = n_i;

        double tau_i;
        if (L.joint_type == 'R') {
            tau_i = L.z.dot(n_i);
        } else { // 'P'
            tau_i = L.z.dot(f_i);
        }
        tau_list[i] = tau_i;

        f_next = f_i;
        n_next = n_i;
    }
}

int main() {
    std::vector<LinkDynamics> links(2);
    // Initialize with some toy data (in a real application, fill from model)
    links[0].F = Eigen::Vector3d(1.0, 0.0, 0.0);
    links[0].N = Eigen::Vector3d(0.0, 0.0, 0.2);
    links[0].R_ip1_i = Eigen::Matrix3d::Identity();
    links[0].p_i = Eigen::Vector3d(0.5, 0.0, 0.0);
    links[0].z = Eigen::Vector3d(0.0, 0.0, 1.0);
    links[0].joint_type = 'R';

    links[1].F = Eigen::Vector3d(0.5, 0.0, 0.0);
    links[1].N = Eigen::Vector3d(0.0, 0.0, 0.1);
    links[1].R_ip1_i = Eigen::Matrix3d::Identity();
    links[1].p_i = Eigen::Vector3d(0.3, 0.0, 0.0);
    links[1].z = Eigen::Vector3d(0.0, 0.0, 1.0);
    links[1].joint_type = 'R';

    Eigen::Vector3d f_tip = Eigen::Vector3d::Zero();
    Eigen::Vector3d n_tip = Eigen::Vector3d::Zero();

    std::vector<Eigen::Vector3d> f_list, n_list;
    std::vector<double> tau_list;

    newtonEulerBackward(links, f_tip, n_tip, f_list, n_list, tau_list);

    std::cout << "Joint torques: ";
    for (double tau : tau_list) {
        std::cout << tau << " ";
    }
    std::cout << std::endl;
    return 0;
}
      

Libraries like RBDL encapsulate this pattern in functions such as InverseDynamics(model, q, qd, qdd, tau), which internally perform both forward and backward recursions on a tree-structured model.

8.1 Java Skeleton (with EJML for Linear Algebra)

Java does not have a canonical robotics library, but matrix packages such as EJML simplify linear algebra. Below is a conceptual skeleton (only the core loop is shown):

import org.ejml.simple.SimpleMatrix;

class LinkDynamics {
    SimpleMatrix F;        // 3x1
    SimpleMatrix N;        // 3x1
    SimpleMatrix R_ip1_i;  // 3x3
    SimpleMatrix p_i;      // 3x1
    SimpleMatrix z;        // 3x1
    char jointType;        // 'R' or 'P'
}

public class NewtonEulerBackward {

    static class Result {
        SimpleMatrix[] fList;
        SimpleMatrix[] nList;
        double[] tauList;
    }

    public static Result compute(LinkDynamics[] links,
                                 SimpleMatrix fTip,
                                 SimpleMatrix nTip) {
        int n = links.length;
        SimpleMatrix[] fList = new SimpleMatrix[n];
        SimpleMatrix[] nList = new SimpleMatrix[n];
        double[] tauList = new double[n];

        SimpleMatrix fNext = fTip.copy();
        SimpleMatrix nNext = nTip.copy();

        for (int i = n - 1; i >= 0; --i) {
            LinkDynamics L = links[i];

            SimpleMatrix fChild_i = L.R_ip1_i.mult(fNext);
            SimpleMatrix nChild_i = L.R_ip1_i.mult(nNext)
                    .plus(L.p_i.cross(fChild_i)); // cross: helper returning 3x1

            SimpleMatrix f_i = L.F.plus(fChild_i);
            SimpleMatrix n_i = L.N.plus(nChild_i);

            fList[i] = f_i;
            nList[i] = n_i;

            double tau_i;
            if (L.jointType == 'R') {
                tau_i = L.z.dot(n_i);
            } else {
                tau_i = L.z.dot(f_i);
            }
            tauList[i] = tau_i;

            fNext = f_i;
            nNext = n_i;
        }

        Result res = new Result();
        res.fList = fList;
        res.nList = nList;
        res.tauList = tauList;
        return res;
    }
}
      

8.2 MATLAB/Simulink Implementation

MATLAB’s Robotics System Toolbox offers rigidBodyTree and inverseDynamics, which encapsulate the Newton–Euler recursions. A simple from-scratch backward recursion for an \( n \)-link chain is:

function [f_list, n_list, tau] = ne_backward(F_list, N_list, R_ip1_i_list, ...
                                             p_i_list, z_list, joint_type_list, ...
                                             f_tip, n_tip)
% F_list, N_list: 3xN
% R_ip1_i_list: 3x3xN
% p_i_list, z_list: 3xN
% joint_type_list: 1xN char array: 'R' or 'P'
N = size(F_list, 2);
f_list = zeros(3, N);
n_list = zeros(3, N);
tau = zeros(1, N);

f_next = f_tip(:);
n_next = n_tip(:);

for i = N:-1:1
    R_ip1_i = R_ip1_i_list(:, :, i);
    p_i = p_i_list(:, i);
    F_i = F_list(:, i);
    N_i = N_list(:, i);
    z_i = z_list(:, i);
    jt = joint_type_list(i);

    f_child_i = R_ip1_i * f_next;
    n_child_i = R_ip1_i * n_next + cross(p_i, f_child_i);

    f_i = F_i + f_child_i;
    n_i = N_i + n_child_i;

    f_list(:, i) = f_i;
    n_list(:, i) = n_i;

    if jt == 'R'
        tau(i) = z_i.' * n_i;
    else
        tau(i) = z_i.' * f_i;
    end

    f_next = f_i;
    n_next = n_i;
end
end
      

In Simulink, one can encapsulate this function in a MATLAB Function block, driven by the joint states \( q, \dot{q}, \ddot{q} \) and model parameters (possibly precomputed in the workspace), and use the outputs as inputs to actuator or control blocks.

8.3 Wolfram Mathematica (Symbolic 2-Link Example)

Mathematica is particularly convenient for symbolic verification of the Newton–Euler algorithm. For a planar 2R arm, the backward recursion can be formulated as:

(* Symbolic variables *)
ClearAll["Global`*"];
n = 2;
m1 = Symbol["m1"]; m2 = Symbol["m2"];
I1 = Symbol["I1"]; I2 = Symbol["I2"];

(* Assume planar motion about z-axis, so we keep only scalar moments about z *)
F1 = {Symbol["Fx1"], Symbol["Fy1"], 0};
F2 = {Symbol["Fx2"], Symbol["Fy2"], 0};
N1 = {0, 0, I1*Symbol["alpha1"] + Symbol["omega1"]*I1*Symbol["omega1"]};
N2 = {0, 0, I2*Symbol["alpha2"] + Symbol["omega2"]*I2*Symbol["omega2"]};

R21 = IdentityMatrix[3];
p1 = {Symbol["l1"], 0, 0};
z1 = {0, 0, 1}; z2 = {0, 0, 1};

f3 = {0, 0, 0};
n3 = {0, 0, 0};

(* Backward recursion *)
f2 = F2 + R21.f3;
n2 = N2 + R21.n3 + Cross[p1, R21.f3];

f1 = F1 + R21.f2;
n1 = N1 + R21.n2 + Cross[p1, R21.f2];

tau2 = z2.n2;
tau1 = z1.n1;

Simplify[{tau1, tau2}]
      

Comparing the resulting symbolic expressions for \( \tau_1, \tau_2 \) with those obtained via a Lagrangian derivation gives a strong consistency check for both modeling and implementation.

In this lesson, we developed the backward recursion of the Newton–Euler algorithm. Starting from rigid-body inertial forces and moments, we derived subtree equilibrium relations that propagate wrenches from the end-effector to the base in \( \mathcal{O}(n) \) time. We showed how to project internal wrenches onto joint motion subspaces to obtain the generalized efforts \( \tau_i \), and we implemented the algorithm in Python, C++, Java, MATLAB/Simulink, and Mathematica. These results provide a computationally efficient and physically transparent method for inverse dynamics, which will be compared to Lagrange–Euler formulations and extended in later chapters to more complex multibody topologies.