For an \(n\)-DOF serial manipulator, the joint-space equations of motion can be written in the standard form

\[ \mathbf{M}(\mathbf{q}) \ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \boldsymbol{\tau}, \quad \mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}} \in \mathbb{R}^n. \]

The Lagrange–Euler formulation uses energies \( T(\mathbf{q}, \dot{\mathbf{q}}) \) and \( V(\mathbf{q}) \) to derive the same equations via

\[ L(\mathbf{q}, \dot{\mathbf{q}}) = T(\mathbf{q}, \dot{\mathbf{q}}) - V(\mathbf{q}), \quad \frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} = Q_i, \quad i = 1,\dots,n, \]

where \( Q_i \) are generalized forces (joint torques/forces). This derivation is symbolic and global: each equation couples all joints through the kinetic energy and potential energy expressions.

The Newton–Euler formulation instead operates at the rigid-body level: it propagates angular and linear velocities and accelerations forward along the chain, computes link wrenches, and then propagates forces and moments backward to obtain joint \( \tau_i \). It is local and recursive.

Both methods produce identical dynamics, but with very different computational structures and costs. This lesson quantifies that difference and shows how to implement the efficient recursive mapping numerically.

Recall that for a serial manipulator with configuration \( \mathbf{q} \), the kinetic energy can be written as

\[ T(\mathbf{q}, \dot{\mathbf{q}}) = \tfrac{1}{2} \dot{\mathbf{q}}^\top \mathbf{M}(\mathbf{q}) \dot{\mathbf{q}}, \]

where \( \mathbf{M}(\mathbf{q}) \) is the positive definite inertia matrix obtained by summing link contributions. The potential energy is \( V(\mathbf{q}) = \sum_{i=1}^n m_i g^\top \mathbf{p}_{c_i}(\mathbf{q}) \), with \( \mathbf{p}_{c_i} \) the center-of-mass position of link \(i\).

Using Lagrange, the \(i\)-th equation reads

\[ \frac{d}{dt} \left( \frac{\partial T}{\partial \dot{q}_i} \right) - \frac{\partial T}{\partial q_i} + \frac{\partial V}{\partial q_i} = \tau_i. \]

To evaluate the left-hand side numerically at a given \( (\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}}) \), the naive Lagrange-based approach requires:

• Constructing (symbolically or automatically) expressions for \( \frac{\partial T}{\partial \dot{q}_i} \), \( \frac{\partial T}{\partial q_i} \), \( \frac{\partial V}{\partial q_i} \), for all \(i\).
• Differentiating \( \frac{\partial T}{\partial \dot{q}_i} \) with respect to time, which involves the chain rule and the full dependence of \( T \) on all \( \dot{q}_j \).

Even after simplification, the inertia matrix entries satisfy

\[ M_{ij}(\mathbf{q}) = \frac{\partial^2 T(\mathbf{q}, \dot{\mathbf{q}})} {\partial \dot{q}_i \, \partial \dot{q}_j}, \quad 1 \leq i, j \leq n, \]

so that computing all \( M_{ij} \) involves roughly \( \mathcal{O}(n^2) \) second derivatives. Each derivative itself is a sum over links, and naive evaluation leads to a computational cost that grows approximately like \( \mathcal{O}(n^3) \) floating-point operations per evaluation of \( \boldsymbol{\tau} \).

In practice, specialized symbolic manipulation and careful factoring can lower the constant factors, but without exploiting the chain structure directly, the asymptotic behavior for large \(n\) is still cubic in \(n\). This becomes problematic for high-DOF systems or high-rate control loops.

The Newton–Euler algorithm computes \( \boldsymbol{\tau} \) using:

1. A forward recursion (from base to end-effector) for angular velocities \( \boldsymbol{\omega}_i \), angular accelerations \( \dot{\boldsymbol{\omega}}_i \), linear accelerations \( \mathbf{a}_i \), and center-of-mass accelerations \( \mathbf{a}_{c_i} \).
2. A backward recursion (from end-effector to base) for forces \( \mathbf{f}_i \), moments \( \mathbf{n}_i \), and joint efforts \( \tau_i \).

In compact 3D vector notation (for a revolute joint \(i\) with axis \( \mathbf{z}_i \), rotation \( \mathbf{R}_i \) from frame \(i-1\) to \(i\), and offset \( \mathbf{p}_i \)), the forward recursion is (assuming a fixed base and \( \mathbf{a}_0 = -\mathbf{g} \)):

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

\[ \dot{\boldsymbol{\omega}}_i = \mathbf{R}_i \dot{\boldsymbol{\omega}}_{i-1} + \mathbf{z}_i \ddot{q}_i + \boldsymbol{\omega}_i \times \mathbf{z}_i \dot{q}_i, \]

\[ \mathbf{a}_i = \mathbf{R}_i \big( \mathbf{a}_{i-1} + \dot{\boldsymbol{\omega}}_{i-1} \times \mathbf{p}_i + \boldsymbol{\omega}_{i-1} \times (\boldsymbol{\omega}_{i-1} \times \mathbf{p}_i) \big), \]

\[ \mathbf{a}_{c_i} = \mathbf{a}_i + \dot{\boldsymbol{\omega}}_i \times \mathbf{r}_{c_i} + \boldsymbol{\omega}_i \times (\boldsymbol{\omega}_i \times \mathbf{r}_{c_i}), \]

where \( \mathbf{r}_{c_i} \) is the COM offset within link \(i\). The corresponding link wrench contributions are

\[ \mathbf{f}_i^{(c)} = m_i \mathbf{a}_{c_i}, \quad \mathbf{n}_i^{(c)} = \mathbf{I}_i \dot{\boldsymbol{\omega}}_i + \boldsymbol{\omega}_i \times (\mathbf{I}_i \boldsymbol{\omega}_i), \]

with \( \mathbf{I}_i \) the inertia tensor about the link origin.

The backward recursion (assuming no external end-effector wrench) is

\[ \mathbf{f}_i = \mathbf{R}_{i+1} \mathbf{f}_{i+1} + \mathbf{f}_i^{(c)}, \]

\[ \mathbf{n}_i = \mathbf{n}_i^{(c)} + \mathbf{R}_{i+1} \mathbf{n}_{i+1} + \mathbf{r}_{c_i} \times \mathbf{f}_i^{(c)} + \mathbf{p}_{i+1} \times (\mathbf{R}_{i+1} \mathbf{f}_{i+1}), \]

\[ \tau_i = \begin{cases} \mathbf{n}_i^\top \mathbf{z}_i & \text{(revolute)} \\ \mathbf{f}_i^\top \mathbf{z}_i & \text{(prismatic)}. \end{cases} \]

For each link, the number of vector operations (matrix–vector products, cross products, and additions) is bounded by a constant independent of \(n\). Since both the forward and backward phases visit each link exactly once,

\[ \text{cost}_{\text{Newton–Euler}}(n) = c_{\text{fwd}} n + c_{\text{bwd}} n = \mathcal{O}(n). \]

flowchart TD
  A["Joint states (q, dq, ddq)"] --> B["Forward pass: velocities, accelerations"]
  B --> C["Per-link wrenches: f_i, n_i"]
  C --> D["Backward pass: accumulate f_i, n_i"]
  D --> E["Joint efforts tau_i"]
  A2["Lagrange side: T(q,dq), V(q)"] --> F["Differentiate, assemble M(q), C(q,dq), g(q)"]
  F --> E2["Joint efforts tau = M ddq + C dq + g"]
  B ---|"O(n)"| D
  F ---|"O(n^3) naive"| E2
        

To compare complexity more concretely, consider computing the inverse dynamics mapping \( \boldsymbol{\tau} = \mathbf{M}(\mathbf{q}) \ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) \) at a single configuration.

• Lagrange-based evaluation. If you store closed-form expressions for the entries of \( \mathbf{M}(\mathbf{q}) \), \( \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}}) \), and \( \mathbf{g}(\mathbf{q}) \), evaluating them typically costs on the order of \( \mathcal{O}(n^3) \) multiplications and additions due to the dense nature of \( \mathbf{M} \) and the quadratic terms in \( \mathbf{C}\dot{\mathbf{q}} \).
• Newton–Euler evaluation. The recursive algorithm directly outputs \( \boldsymbol{\tau} \) in \( \mathcal{O}(n) \) operations with a relatively small constant factor. No intermediate dense matrices are formed.

In fact, one can construct the inertia matrix numerically using the inverse dynamics oracle:

\[ \mathbf{M}(\mathbf{q}) \mathbf{e}_j = \boldsymbol{\tau}(\mathbf{q}, \dot{\mathbf{q}} = \mathbf{0}, \ddot{\mathbf{q}} = \mathbf{e}_j), \quad j = 1,\dots,n, \]

where \( \mathbf{e}_j \) is the \(j\)-th standard basis vector. Each column of \( \mathbf{M}(\mathbf{q}) \) can be obtained with a single call to inverse dynamics, so computing the full inertia matrix costs \( \mathcal{O}(n^2) \) using an underlying \( \mathcal{O}(n) \) Newton–Euler call per column.

Similarly, by evaluating inverse dynamics at \( \ddot{\mathbf{q}} = \mathbf{0} \) and nonzero \( \dot{\mathbf{q}} \) or by turning gravity on/off, one can obtain the gravity and Coriolis terms numerically in \( \mathcal{O}(n^2) \), still significantly more efficient than a naive symbolic Lagrange approach for large \(n\).

Despite its unfavorable scaling for numerical evaluation, the Lagrange–Euler formulation remains extremely valuable:

• Structural insight. The energy-based derivation directly reveals the symmetry and positive definiteness of \( \mathbf{M}(\mathbf{q}) \), and passivity properties of the dynamics.
• Symbolic simplification. For low-DOF systems (e.g., 2R or 3R arms), closed-form Lagrange-derived expressions are compact and useful for analysis, education, and controller theory.
• Validation. Newton–Euler results are often cross-checked against Lagrange-derived equations on simple examples to validate implementations.

In practice, the derivation is often done once using Lagrange, while the runtime evaluation is implemented using the Newton–Euler recursion (or its rigid-body-operator variants) to meet real-time constraints.

In Python, inverse dynamics based on Newton–Euler can be implemented with numpy for vector operations. High-quality robotics libraries such as roboticstoolbox-python (from Peter Corke) and pin (Python bindings of Pinocchio) internally use recursive algorithms for efficiency.

The following sketch emphasizes the \( \mathcal{O}(n) \) loop structure. We assume that, for a fixed configuration \( \mathbf{q} \), we already have the per-link rotations \( \mathbf{R}_i \) and offsets \( \mathbf{p}_i \) from forward kinematics (Chapter 5).

import numpy as np

Z_AXIS = np.array([0.0, 0.0, 1.0])

class Link:
    def __init__(self, mass, inertia, r_com, joint_type="R"):
        """
        mass    : scalar
        inertia : 3x3 inertia matrix about link frame origin
        r_com   : 3D vector from link frame origin to center of mass
        joint_type: "R" for revolute, "P" for prismatic
        """
        self.m = mass
        self.I = inertia
        self.r_com = r_com
        self.joint_type = joint_type


def inverse_dynamics_newton_euler(q, qd, qdd, g, R_list, p_list, links):
    """
    q, qd, qdd : arrays of shape (n,)
    g          : 3D gravity vector (e.g. [0, 0, -9.81])
    R_list[i]  : 3x3 rotation from frame (i-1) to i
    p_list[i]  : 3D vector from frame (i-1) origin to frame i origin
    links[i]   : Link object for joint i
    returns    : tau of shape (n,)
    """
    n = len(links)

    # Forward recursion: velocities and accelerations
    w = [np.zeros(3) for _ in range(n + 1)]
    wd = [np.zeros(3) for _ in range(n + 1)]
    a = [np.zeros(3) for _ in range(n + 1)]

    # Base acceleration includes gravity
    a[0] = -np.asarray(g, dtype=float)

    for i in range(1, n + 1):
        R_i = R_list[i - 1]
        p_i = p_list[i - 1]
        link = links[i - 1]

        if link.joint_type == "R":
            # revolute joint
            w[i] = R_i @ w[i - 1] + Z_AXIS * qd[i - 1]
            wd[i] = (R_i @ wd[i - 1]
                     + Z_AXIS * qdd[i - 1]
                     + np.cross(w[i], Z_AXIS * qd[i - 1]))
            a[i] = (R_i @ (a[i - 1]
                           + np.cross(wd[i - 1], p_i)
                           + np.cross(w[i - 1],
                                      np.cross(w[i - 1], p_i))))
        else:
            # prismatic joint (simplified)
            w[i] = R_i @ w[i - 1]
            wd[i] = R_i @ wd[i - 1]
            a[i] = (R_i @ (a[i - 1]
                           + np.cross(wd[i - 1], p_i)
                           + np.cross(w[i - 1],
                                      np.cross(w[i - 1], p_i)))
                    + Z_AXIS * qdd[i - 1]
                    + 2.0 * np.cross(w[i], Z_AXIS * qd[i - 1]))

    # Center-of-mass accelerations, forces and moments
    F = [np.zeros(3) for _ in range(n)]
    N = [np.zeros(3) for _ in range(n)]
    for i in range(n):
        link = links[i]
        r_c = link.r_com
        a_c = (a[i + 1]
               + np.cross(wd[i + 1], r_c)
               + np.cross(w[i + 1], np.cross(w[i + 1], r_c)))
        F[i] = link.m * a_c
        N[i] = (link.I @ wd[i + 1]
                + np.cross(w[i + 1], link.I @ w[i + 1]))

    # Backward recursion: forces, moments, and joint efforts
    f_next = np.zeros(3)
    n_next = np.zeros(3)
    tau = np.zeros(n)

    for i in reversed(range(n)):
        R_ip1 = R_list[i] if i < n - 1 else np.eye(3)
        p_ip1 = p_list[i] if i < n - 1 else np.zeros(3)
        r_c = links[i].r_com

        f_i = R_ip1 @ f_next + F[i]
        n_i = (N[i]
               + R_ip1 @ n_next
               + np.cross(r_c, F[i])
               + np.cross(p_ip1, R_ip1 @ f_next))

        if links[i].joint_type == "R":
            tau[i] = n_i.dot(Z_AXIS)
        else:
            tau[i] = f_i.dot(Z_AXIS)

        f_next, n_next = f_i, n_i

    return tau
      

Every loop over i is linear in \(n\). The constant factor is dominated by 3D vector operations (matrix–vector products and cross products), which are highly optimized in numerical libraries.

In C++, highly optimized dynamics libraries such as RBDL (Rigid Body Dynamics Library), KDL (Kinematics and Dynamics Library), and Pinocchio implement Newton–Euler and related algorithms (composite rigid-body, articulated-body) using Eigen for linear algebra.

Below is a stripped-down skeleton of a Newton–Euler inverse dynamics routine using Eigen. The emphasis is again on the linear pass structure, not on all geometric details.

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

struct Link {
    double m;
    Eigen::Matrix3d I;
    Eigen::Vector3d r_com;
    char joint_type; // 'R' or 'P'
};

Eigen::VectorXd inverseDynamicsNewtonEuler(
    const Eigen::VectorXd& q,
    const Eigen::VectorXd& qd,
    const Eigen::VectorXd& qdd,
    const Eigen::Vector3d& g,
    const std::vector<Eigen::Matrix3d>& R,
    const std::vector<Eigen::Vector3d>& p,
    const std::vector<Link>& links)
{
    const int n = static_cast<int>(links.size());
    const Eigen::Vector3d Z_AXIS(0.0, 0.0, 1.0);

    std::vector<Eigen::Vector3d> w(n + 1, Eigen::Vector3d::Zero());
    std::vector<Eigen::Vector3d> wd(n + 1, Eigen::Vector3d::Zero());
    std::vector<Eigen::Vector3d> a(n + 1, Eigen::Vector3d::Zero());
    a[0] = -g;

    // Forward recursion
    for (int i = 1; i <= n; ++i) {
        const Eigen::Matrix3d& R_i = R[i - 1];
        const Eigen::Vector3d& p_i = p[i - 1];
        const Link& link = links[i - 1];

        if (link.joint_type == 'R') {
            w[i] = R_i * w[i - 1] + Z_AXIS * qd[i - 1];
            wd[i] = R_i * wd[i - 1]
                    + Z_AXIS * qdd[i - 1]
                    + w[i].cross(Z_AXIS * qd[i - 1]);
            a[i] = R_i * (a[i - 1]
                          + wd[i - 1].cross(p_i)
                          + w[i - 1].cross(w[i - 1].cross(p_i)));
        } else {
            w[i] = R_i * w[i - 1];
            wd[i] = R_i * wd[i - 1];
            a[i] = R_i * (a[i - 1]
                          + wd[i - 1].cross(p_i)
                          + w[i - 1].cross(w[i - 1].cross(p_i)))
                   + Z_AXIS * qdd[i - 1]
                   + 2.0 * w[i].cross(Z_AXIS * qd[i - 1]);
        }
    }

    // Per-link forces and moments
    std::vector<Eigen::Vector3d> F(n), N(n);
    for (int i = 0; i < n; ++i) {
        const Link& link = links[i];
        Eigen::Vector3d r_c = link.r_com;
        Eigen::Vector3d a_c = a[i + 1]
                              + wd[i + 1].cross(r_c)
                              + w[i + 1].cross(w[i + 1].cross(r_c));
        F[i] = link.m * a_c;
        N[i] = link.I * wd[i + 1]
               + w[i + 1].cross(link.I * w[i + 1]);
    }

    // Backward recursion
    Eigen::Vector3d f_next = Eigen::Vector3d::Zero();
    Eigen::Vector3d n_next = Eigen::Vector3d::Zero();
    Eigen::VectorXd tau(n);
    tau.setZero();

    for (int i = n - 1; i >= 0; --i) {
        Eigen::Matrix3d R_ip1 = (i < n - 1) ? R[i] : Eigen::Matrix3d::Identity();
        Eigen::Vector3d p_ip1 = (i < n - 1) ? p[i] : Eigen::Vector3d::Zero();
        Eigen::Vector3d r_c = links[i].r_com;

        Eigen::Vector3d f_i = R_ip1 * f_next + F[i];
        Eigen::Vector3d n_i = N[i]
                              + R_ip1 * n_next
                              + r_c.cross(F[i])
                              + p_ip1.cross(R_ip1 * f_next);

        if (links[i].joint_type == 'R')
            tau[i] = n_i.dot(Z_AXIS);
        else
            tau[i] = f_i.dot(Z_AXIS);

        f_next = f_i;
        n_next = n_i;
    }
    return tau;
}
      

The outer loops are both linear in \(n\). Modern C++ robotics libraries further exploit cache locality, template metaprogramming, and spatial vector algebra to reduce constant factors.

Java is less common in low-level robot dynamics, but the same recursive algorithm structure applies. Often, dynamics are implemented in project-specific classes. Below is a minimalistic sketch in plain Java with arrays and small vector helpers.

class Link {
    double m;
    double[][] I;      // 3x3 inertia matrix
    double[] rCom;     // length-3 vector
    char jointType;    // 'R' or 'P'

    Link(double mass, double[][] inertia, double[] rCom, char jointType) {
        this.m = mass;
        this.I = inertia;
        this.rCom = rCom;
        this.jointType = jointType;
    }
}

public class NewtonEuler {

    private static final double[] Z_AXIS = new double[]{0.0, 0.0, 1.0};

    // Vector and matrix helpers (3D) omitted for brevity:
    // add(a,b), sub(a,b), cross(a,b), matVec(R,v), etc.
    // Assume they allocate small arrays and perform 3x3 ops.

    public static double[] inverseDynamics(
            double[] q,
            double[] qd,
            double[] qdd,
            double[] g,
            double[][][] R,
            double[][] p,
            Link[] links) {

        int n = links.length;

        double[][] w = new double[n + 1][3];
        double[][] wd = new double[n + 1][3];
        double[][] a = new double[n + 1][3];

        // a[0] = -g
        a[0][0] = -g[0];
        a[0][1] = -g[1];
        a[0][2] = -g[2];

        // forward recursion
        for (int i = 1; i <= n; ++i) {
            double[][] R_i = R[i - 1];
            double[] p_i = p[i - 1];
            Link link = links[i - 1];

            if (link.jointType == 'R') {
                double[] wPrev = w[i - 1];
                w[i] = add(matVec(R_i, wPrev), scale(Z_AXIS, qd[i - 1]));
                wd[i] = add(
                        add(matVec(R_i, wd[i - 1]),
                            scale(Z_AXIS, qdd[i - 1])),
                        cross(w[i], scale(Z_AXIS, qd[i - 1])));
                a[i] = matVec(R_i,
                        add(a[i - 1],
                            add(cross(wd[i - 1], p_i),
                                cross(w[i - 1], cross(w[i - 1], p_i)))));
            } else {
                // prismatic (simplified)
                w[i] = matVec(R_i, w[i - 1]);
                wd[i] = matVec(R_i, wd[i - 1]);
                a[i] = add(
                        matVec(R_i,
                               add(a[i - 1],
                                   add(cross(wd[i - 1], p_i),
                                       cross(w[i - 1], cross(w[i - 1], p_i))))),
                        add(scale(Z_AXIS, qdd[i - 1]),
                            scale(cross(w[i], Z_AXIS), 2.0 * qd[i - 1])));
            }
        }

        // forces and moments
        double[][] F = new double[n][3];
        double[][] N = new double[n][3];
        for (int i = 0; i < n; ++i) {
            Link link = links[i];
            double[] rC = link.rCom;
            double[] aC = add(a[i + 1],
                              add(cross(wd[i + 1], rC),
                                  cross(w[i + 1], cross(w[i + 1], rC))));
            F[i] = scale(aC, link.m);
            N[i] = add(matVec(link.I, wd[i + 1]),
                       cross(w[i + 1], matVec(link.I, w[i + 1])));
        }

        double[] fNext = new double[3];
        double[] nNext = new double[3];
        double[] tau = new double[n];

        for (int i = n - 1; i >= 0; --i) {
            double[][] R_ip1 = (i < n - 1) ? R[i] : identity3();
            double[] p_ip1 = (i < n - 1) ? p[i] : new double[3];
            double[] rC = links[i].rCom;

            double[] f_i = add(matVec(R_ip1, fNext), F[i]);
            double[] n_i = add(
                    add(N[i], matVec(R_ip1, nNext)),
                    add(cross(rC, F[i]),
                        cross(p_ip1, matVec(R_ip1, fNext))));

            if (links[i].jointType == 'R') {
                tau[i] = dot(n_i, Z_AXIS);
            } else {
                tau[i] = dot(f_i, Z_AXIS);
            }

            fNext = f_i;
            nNext = n_i;
        }
        return tau;
    }
}
      

Although Java lacks the specialized robotics libraries of C++, the Newton–Euler structure is language-agnostic: forward and backward linear passes with constant per-link work.

MATLAB supports efficient Newton–Euler dynamics in several ways:

• Robotics System Toolbox. The function inverseDynamics applied to a rigidBodyTree object uses recursive algorithms internally.
• Simscape Multibody. Provides block-diagram modeling of multi-body systems; its internal solvers exploit sparse and recursive structures.
• Custom scripts. You can implement your own Newton–Euler recursion when interfacing with special hardware.

Example: calling the toolbox implementation for a robot already defined as a rigidBodyTree:

% Assume 'robot' is a rigidBodyTree from Robotics System Toolbox

q   = rand(robot.NumBodies, 1);   % joint positions
qd  = rand(robot.NumBodies, 1);   % joint velocities
qdd = rand(robot.NumBodies, 1);   % joint accelerations

tau = inverseDynamics(robot, q, qd, qdd);  % uses Newton-Euler internally

% Extract gravity, Coriolis, and inertia numerically:
g_vec = inverseDynamics(robot, q, zeros(size(q)), zeros(size(q)));
Cqd   = inverseDynamics(robot, q, qd, zeros(size(q))) - g_vec;
M_col1 = inverseDynamics(robot, q, zeros(size(q)), [1; zeros(robot.NumBodies-1, 1)]) - g_vec;
      

To deploy in Simulink, one can wrap inverseDynamics in a MATLAB Function block and feed it joint trajectories, or rely directly on Simscape Multibody for dynamics. In all cases, the underlying computational cost is linear in the number of joints using recursive algorithms, not cubic.

Wolfram Mathematica is well suited for symbolic Lagrange derivations for low-DOF arms and for comparing them with numeric Newton–Euler implementations.

Consider a simple 2R planar manipulator. Symbolically, we can derive its equations of motion via Lagrange:

(* Joint coordinates and parameters *)
q1[t_]; q2[t_];
m1 = Symbol["m1"]; m2 = Symbol["m2"];
l1 = Symbol["l1"]; l2 = Symbol["l2"];
g  = Symbol["g"];

(* Positions of centers of mass in plane *)
x1[t_] := (l1/2) Cos[q1[t]];
y1[t_] := (l1/2) Sin[q1[t]];

x2[t_] := l1 Cos[q1[t]] + (l2/2) Cos[q1[t] + q2[t]];
y2[t_] := l1 Sin[q1[t]] + (l2/2) Sin[q1[t] + q2[t]];

(* Kinetic and potential energy *)
vx1[t_] := D[x1[t], t];
vy1[t_] := D[y1[t], t];
vx2[t_] := D[x2[t], t];
vy2[t_] := D[y2[t], t];

T = Simplify[
   1/2 m1 (vx1[t]^2 + vy1[t]^2) +
   1/2 m2 (vx2[t]^2 + vy2[t]^2)
];

V = Simplify[
   m1 g y1[t] + m2 g y2[t]
];

L = T - V;

(* Lagrange equations *)
coords = {q1[t], q2[t]};
tauSym =
  Table[
    FullSimplify[
      D[D[L, D[coords[[i]], t]], t] - D[L, coords[[i]]]
    ],
    {i, 1, 2}
  ];
      

The resulting symbolic expressions are algebraically identical to those obtainable from Newton–Euler. For numerical evaluation at a given configuration, however, we can convert the Newton–Euler recursion to a Mathematica function that performs only numeric 2D vector operations, which will scale linearly with the number of links. Timing experiments using AbsoluteTiming reveal the expected advantage of recursive evaluation as the number of DOFs increases.

In this lesson we contrasted the Lagrange–Euler and Newton–Euler formulations of robot dynamics from a computational viewpoint. While Lagrange offers a clear energy-based derivation and structural insight, its naive numerical evaluation scales poorly with the number of joints. The Newton–Euler recursion, by exploiting the tree structure of serial manipulators, evaluates the inverse dynamics mapping in linear time, and serves as the backbone of modern robotics dynamics libraries in Python, C++, MATLAB, and other languages.

You also saw how to construct the inertia matrix and other dynamic quantities numerically using inverse dynamics calls, and how symbolic tools (e.g., Mathematica) complement recursive numeric algorithms. In the next lesson, we will focus on a more detailed implementation lab for the Newton–Euler algorithm and its validation against Lagrange results.