In previous lessons of this chapter you derived the Lagrange–Euler equations of motion for general manipulators, \( \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \boldsymbol{\tau} \), and studied the structure and energy consistency of \( \mathbf{M} \) and \( \mathbf{C} \). In this lab we ask:

• How do we compute these terms in practice for a given robot?
• When is it better to use symbolic derivation (closed forms)?
• When is it better to use numeric derivation (finite differences or automatic differentiation)?

We restrict to serial manipulators with generalized coordinates \( \mathbf{q} \in \mathbb{R}^n \), and assume you already know how to express kinetic energy \( T(\mathbf{q},\dot{\mathbf{q}}) \) and potential energy \( V(\mathbf{q}) \). The Lagrangian is

\[ L(\mathbf{q},\dot{\mathbf{q}}) = T(\mathbf{q},\dot{\mathbf{q}}) - V(\mathbf{q}). \]

The Euler–Lagrange equations for coordinate \( q_i \) are

\[ \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i} = \tau_i,\quad i=1,\dots,n. \]

The two derivation paradigms are illustrated in the following flow:

flowchart TD
  A["Robot specification (links, masses, inertias)"]
    --> B["Define T(q, qdot) and V(q)"]
  B --> C1["Symbolic pipeline"]
  B --> C2["Numeric pipeline"]

  C1 --> S1["Analytic derivatives of L(q, qdot)"]
  S1 --> S2["Closed-form M(q), C(q, qdot), g(q)"]
  S2 --> S3["Code generation / simplification"]

  C2 --> N1["Numeric evaluation of T(q, qdot), V(q)"]
  N1 --> N2["Finite difference or auto diff"]
  N2 --> N3["Numeric M(q), C(q, qdot), g(q) at sampled states"]
        

Consider a planar 2-link revolute manipulator (2R) in the horizontal plane with joint angles \( \mathbf{q} = [q_1, q_2]^\top \). Let link lengths be \( l_1, l_2 \), center-of-mass distances \( l_{c1}, l_{c2} \), masses \( m_1, m_2 \), and planar moments of inertia about centers of mass \( I_1, I_2 \). For simplicity, assume gravity acts in the plane so that potential energy is well defined; the derivation is identical for horizontal operation with \( V(\mathbf{q}) = 0 \) (then only kinetic energy contributes).

2.1 Kinetic and potential energy

Using planar kinematics (derived in earlier chapters), one can express the squared linear velocity of the link centers as a quadratic form in \( \dot{\mathbf{q}} \). The total kinetic energy is

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

with inertia matrix

\[ \mathbf{M}(\mathbf{q}) = \begin{bmatrix} M_{11} & M_{12} \\ M_{21} & M_{22} \end{bmatrix}, \quad M_{12} = M_{21}. \]

A standard symbolic derivation (using homogeneous transforms or PoE) gives

\[ \begin{aligned} M_{11}(\mathbf{q}) &= I_1 + I_2 + m_1 l_{c1}^2 + m_2 \bigl(l_1^2 + l_{c2}^2 + 2 l_1 l_{c2}\cos q_2\bigr), \\ M_{12}(\mathbf{q}) &= I_2 + m_2 \bigl(l_{c2}^2 + l_1 l_{c2}\cos q_2\bigr), \\ M_{22}(\mathbf{q}) &= I_2 + m_2 l_{c2}^2. \end{aligned} \]

The potential energy (with vertical coordinate defined from previous kinematic derivation) is

\[ V(\mathbf{q}) = (m_1 l_{c1} + m_2 l_1) g \cos q_1 + m_2 l_{c2} g \cos(q_1 + q_2). \]

2.2 From energies to equations of motion

By definition, \( M_{ij}(\mathbf{q}) = \frac{\partial^2 T}{\partial \dot{q}_i \partial \dot{q}_j} \). It follows immediately that \( \mathbf{M}(\mathbf{q}) \) is symmetric:

\[ M_{ij} = \frac{\partial}{\partial \dot{q}_j} \left(\frac{\partial T}{\partial \dot{q}_i}\right) = \frac{\partial}{\partial \dot{q}_i} \left(\frac{\partial T}{\partial \dot{q}_j}\right) = M_{ji}, \]

since mixed partial derivatives of a smooth scalar field commute.

Expanding the Euler–Lagrange equation for the 2R arm yields

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

where \(\mathbf{g}(\mathbf{q}) = \nabla_{\mathbf{q}} V(\mathbf{q})\) and \(\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\) is defined from the Christoffel symbols (previous lesson):

\[ c_{ijk}(\mathbf{q}) = \tfrac{1}{2} \left( \frac{\partial M_{ik}}{\partial q_j} + \frac{\partial M_{ij}}{\partial q_k} - \frac{\partial M_{jk}}{\partial q_i} \right), \quad \bigl[\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}\bigr]_i = \sum_{j,k} c_{ijk}(\mathbf{q}) \dot{q}_j \dot{q}_k. \]

In symbolic derivation, all these partial derivatives are computed analytically, typically by a CAS (computer algebra system), yielding exact expressions for \( \mathbf{M}, \mathbf{C}, \mathbf{g} \).

Numeric derivation replaces analytic partial derivatives with approximations. Assume you are able to compute \( T(\mathbf{q},\dot{\mathbf{q}}) \) and \( V(\mathbf{q}) \) numerically for any \( (\mathbf{q},\dot{\mathbf{q}}) \). This is often the case when you have a forward kinematics routine and rigid-body parameters but not explicit formulas.

3.1 Inertia matrix from kinetic energy

For a fixed configuration \( \mathbf{q} \), the entries of \( \mathbf{M}(\mathbf{q}) \) satisfy

\[ M_{ij}(\mathbf{q}) = \frac{\partial^2 T(\mathbf{q},\dot{\mathbf{q}})} {\partial \dot{q}_i \partial \dot{q}_j}. \]

A second-order accurate central-difference approximation for the first derivative is, for a scalar function \( f(x) \),

\[ f'(x) \approx \frac{f(x + h) - f(x - h)}{2h} + \mathcal{O}(h^2). \]

Applying this twice yields an approximation of \( M_{ij} \). For example, for a fixed \( \mathbf{q} \) and nominal \( \dot{\mathbf{q}} = \mathbf{0} \),

\[ \begin{aligned} M_{ij}(\mathbf{q}) &\approx \frac{1}{4 h_i h_j} \Bigl( T(\mathbf{q}, \mathbf{e}_i h_i + \mathbf{e}_j h_j) - T(\mathbf{q}, \mathbf{e}_i h_i - \mathbf{e}_j h_j) \\ &\quad - T(\mathbf{q}, -\mathbf{e}_i h_i + \mathbf{e}_j h_j) + T(\mathbf{q}, -\mathbf{e}_i h_i - \mathbf{e}_j h_j) \Bigr), \end{aligned} \]

where \( \mathbf{e}_i \) is the \( i \)-th unit vector and \( h_i, h_j \) are small step sizes.

3.2 Gravity and Coriolis terms

The gravity term is the gradient of potential energy:

\[ g_i(\mathbf{q}) = \frac{\partial V(\mathbf{q})}{\partial q_i} \approx \frac{V(\mathbf{q} + h_i \mathbf{e}_i) - V(\mathbf{q} - h_i \mathbf{e}_i)}{2 h_i}. \]

The Coriolis/centrifugal term can be approximated by numerically evaluating \( \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} \) and taking differences of the full left-hand side, or by directly approximating Christoffel symbols via finite differences of \( \mathbf{M} \) with respect to \( \mathbf{q} \). For example,

\[ \frac{\partial M_{ij}}{\partial q_k}(\mathbf{q}) \approx \frac{ M_{ij}(\mathbf{q} + h_k \mathbf{e}_k) - M_{ij}(\mathbf{q} - h_k \mathbf{e}_k) }{2 h_k}. \]

Numeric derivation is approximate but can be implemented for arbitrary robot models when symbolic expressions would be unwieldy. It is also closely related to automatic differentiation methods, which compute exact derivatives of the implemented program without explicit symbolic manipulation.

In a practical lab you typically:

1. Derive a symbolic model for a simple manipulator (2R) using a CAS.
2. Generate numeric functions \( \mathbf{M}_\text{sym}, \mathbf{C}_\text{sym}, \mathbf{g}_\text{sym} \).
3. Implement numeric finite-difference approximations \( \mathbf{M}_\text{num}, \mathbf{C}_\text{num}, \mathbf{g}_\text{num} \).
4. Sample many random states \( (\mathbf{q},\dot{\mathbf{q}}) \) in the joint limits.
5. Compare the two models via norms of their differences and timing.

The following flowchart summarizes this experimental pipeline:

flowchart TD
  S["Start: choose robot (e.g., 2R)"]
    --> A1["Symbolic model: derive T(q, qdot), V(q)"]
  A1 --> A2["CAS: compute M_sym, C_sym, g_sym"]
  A2 --> A3["Generate fast numeric code"]

  S --> B1["Numeric model: \nimplement T(q, qdot), V(q) evaluator"]
  B1 --> B2["Finite diff / auto diff \nto get M_num, C_num, g_num"]

  A3 --> C["Sample states in joint limits"]
  B2 --> C
  C --> D["Compare errors: norms of differences, energy checks"]
  D --> E["Measure compute time and memory"]
  E --> F["Decide strategy for real-time control / simulation"]
        

In the rest of the lesson we construct implementations for this pipeline in several programming languages and discuss numerical issues such as step size and conditioning.

Python is a natural environment for this lab because it offers both symbolic (SymPy) and numeric (NumPy) computation. Below we sketch a workflow for the 2R arm:

import sympy as sp
import numpy as np

# 1. Symbolic definitions
q1, q2 = sp.symbols('q1 q2', real=True)
dq1, dq2 = sp.symbols('dq1 dq2', real=True)
g, l1, l2, lc1, lc2, m1, m2, I1, I2 = sp.symbols(
    'g l1 l2 lc1 lc2 m1 m2 I1 I2', positive=True, real=True
)

q = sp.Matrix([q1, q2])
dq = sp.Matrix([dq1, dq2])

# 2. Kinetic energy T and potential energy V (from Section 2)
M11 = I1 + I2 + m1*lc1**2 + m2*(l1**2 + lc2**2 + 2*l1*lc2*sp.cos(q2))
M12 = I2 + m2*(lc2**2 + l1*lc2*sp.cos(q2))
M22 = I2 + m2*lc2**2

M = sp.Matrix([[M11, M12],
               [M12, M22]])

T = sp.Rational(1, 2) * dq.T * M * dq
T = sp.simplify(T[0])  # scalar

V = (m1*lc1 + m2*l1)*g*sp.cos(q1) + m2*lc2*g*sp.cos(q1 + q2)

L = T - V

# 3. Euler-Lagrange equations to get closed-form dynamics
tau = sp.symbols('tau1 tau2', real=True)

dLd_dq = [sp.diff(L, q_i) for q_i in q]
dLd_ddq = [sp.diff(L, dq_i) for dq_i in dq]

# Time derivatives d/dt(dL/ddq_i): treat q, dq as independent and replace dq, ddq
ddq1, ddq2 = sp.symbols('ddq1 ddq2', real=True)
ddq = sp.Matrix([ddq1, ddq2])

# Here we use the known structure T = 0.5 * dq.T * M(q) * dq
# so that M(q) is already identified; we directly build M, C, g:
g_vec = sp.Matrix([sp.diff(V, q1), sp.diff(V, q2)])

# Christoffel-based C(q, dq):
C = sp.zeros(2, 2)
for i in range(2):
    for j in range(2):
        Cij = 0
        for k in range(2):
            c_ijk = 0.5 * (
                sp.diff(M[i, j], q[k]) +
                sp.diff(M[i, k], q[j]) -
                sp.diff(M[j, k], q[i])
            )
            Cij += c_ijk * dq[k]
        C[i, j] = sp.simplify(Cij)

# 4. Lambdify numeric versions M(q), C(q, dq), g(q)
M_func = sp.lambdify((q1, q2, l1, l2, lc1, lc2, m1, m2, I1, I2),
                     M, 'numpy')
C_func = sp.lambdify((q1, q2, dq1, dq2, l1, l2, lc1, lc2, m1, m2, I1, I2),
                     C, 'numpy')
g_func = sp.lambdify((q1, q2, g, l1, l2, lc1, lc2, m1, m2),
                     g_vec, 'numpy')

# 5. Numeric finite-difference approximation for M(q) as a check
def T_numeric(q_vec, dq_vec, params):
    # params = (g, l1, l2, lc1, lc2, m1, m2, I1, I2)
    g_, l1_, l2_, lc1_, lc2_, m1_, m2_, I1_, I2_ = params
    q1_, q2_ = q_vec
    dq1_, dq2_ = dq_vec
    M_num = M_func(q1_, q2_, l1_, l2_, lc1_, lc2_, m1_, m2_, I1_, I2_)
    return 0.5 * dq_vec.T @ M_num @ dq_vec

def M_fd(q_vec, params, h=1e-6):
    n = len(q_vec)
    M_fd = np.zeros((n, n))
    dq0 = np.zeros(n)
    # use central differences on dq only, around dq0
    for i in range(n):
        for j in range(n):
            ei = np.zeros(n); ei[i] = 1.0
            ej = np.zeros(n); ej[j] = 1.0
            Tpp = T_numeric(q_vec, dq0 + h*ei + h*ej, params)
            Tpm = T_numeric(q_vec, dq0 + h*ei - h*ej, params)
            Tmp = T_numeric(q_vec, dq0 - h*ei + h*ej, params)
            Tmm = T_numeric(q_vec, dq0 - h*ei - h*ej, params)
            M_fd[i, j] = (Tpp - Tpm - Tmp + Tmm) / (4*h*h)
    return 0.5*(M_fd + M_fd.T)  # symmetrize

# 6. Random test
params_num = (9.81, 1.0, 1.0, 0.5, 0.5, 2.0, 1.0, 0.1, 0.1)
q_vec = np.array([0.3, -0.7])

M_sym = M_func(q_vec[0], q_vec[1], *params_num[1:])
M_fd_ = M_fd(q_vec, params_num)

print("M_sym =\n", M_sym)
print("M_fd  =\n", M_fd_)
print("diff norm:", np.linalg.norm(M_sym - M_fd_))
      

This code realizes both symbolic and numeric derivation of \( \mathbf{M}(\mathbf{q}) \) and provides a quantitative comparison of their difference for random joint configurations.

In C++, symbolic derivation is usually performed offline (Python/Mathematica) and the result is pasted as hard-coded expressions, or generated via code generation tools. Here we show a manually coded numeric evaluator for the 2R inertia matrix and gravity vector using Eigen. These represent the symbolic expressions numerically.

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

using Eigen::Matrix2d;
using Eigen::Vector2d;

struct TwoRParams {
  double g;
  double l1, l2;
  double lc1, lc2;
  double m1, m2;
  double I1, I2;
};

Matrix2d M_sym(const Vector2d& q, const TwoRParams& p) {
  const double q2 = q(1);
  const double c2 = std::cos(q2);

  const double M11 =
      p.I1 + p.I2 +
      p.m1 * p.lc1 * p.lc1 +
      p.m2 * (p.l1 * p.l1 + p.lc2 * p.lc2 + 2.0 * p.l1 * p.lc2 * c2);

  const double M12 =
      p.I2 + p.m2 * (p.lc2 * p.lc2 + p.l1 * p.lc2 * c2);

  const double M22 = p.I2 + p.m2 * p.lc2 * p.lc2;

  Matrix2d M;
  M(0, 0) = M11;
  M(0, 1) = M12;
  M(1, 0) = M12;
  M(1, 1) = M22;
  return M;
}

Vector2d g_sym(const Vector2d& q, const TwoRParams& p) {
  const double q1 = q(0);
  const double q2 = q(1);

  const double c1    = std::cos(q1);
  const double c12   = std::cos(q1 + q2);

  const double g1 =
      (p.m1 * p.lc1 + p.m2 * p.l1) * p.g * c1 +
      p.m2 * p.lc2 * p.g * c12;

  const double g2 =
      p.m2 * p.lc2 * p.g * c12;

  return Vector2d(g1, g2);
}

int main() {
  TwoRParams p{9.81, 1.0, 1.0, 0.5, 0.5, 2.0, 1.0, 0.1, 0.1};
  Vector2d q;
  q << 0.3, -0.7;

  Matrix2d M = M_sym(q, p);
  Vector2d g = g_sym(q, p);

  std::cout << "M(q):\n" << M << "\n";
  std::cout << "g(q): " << g.transpose() << "\n";
  return 0;
}
      

Numeric finite-difference approximations analogous to Section 3 can also be implemented in C++ by wrapping a kinetic-energy routine and perturbing \( \dot{\mathbf{q}} \). The hard-coded symbolic model is typically used in real-time control, while finite differences are used for validation.

In Java, the EJML library provides dense linear algebra. We mirror the C++ symbolic numeric evaluator in Java, and sketch a finite-difference routine for \( \mathbf{M}(\mathbf{q}) \).

import org.ejml.simple.SimpleMatrix;

public class TwoRDynamics {

    public static class Params {
        public double g, l1, l2, lc1, lc2, m1, m2, I1, I2;
    }

    public static SimpleMatrix M_sym(SimpleMatrix q, Params p) {
        double q2 = q.get(1);
        double c2 = Math.cos(q2);

        double M11 = p.I1 + p.I2
                + p.m1 * p.lc1 * p.lc1
                + p.m2 * (p.l1 * p.l1 + p.lc2 * p.lc2
                + 2.0 * p.l1 * p.lc2 * c2);

        double M12 = p.I2
                + p.m2 * (p.lc2 * p.lc2 + p.l1 * p.lc2 * c2);

        double M22 = p.I2 + p.m2 * p.lc2 * p.lc2;

        SimpleMatrix M = new SimpleMatrix(2, 2);
        M.set(0, 0, M11);
        M.set(0, 1, M12);
        M.set(1, 0, M12);
        M.set(1, 1, M22);
        return M;
    }

    public static SimpleMatrix g_sym(SimpleMatrix q, Params p) {
        double q1 = q.get(0);
        double q2 = q.get(1);
        double c1 = Math.cos(q1);
        double c12 = Math.cos(q1 + q2);

        double g1 = (p.m1 * p.lc1 + p.m2 * p.l1) * p.g * c1
                + p.m2 * p.lc2 * p.g * c12;
        double g2 = p.m2 * p.lc2 * p.g * c12;

        SimpleMatrix gVec = new SimpleMatrix(2, 1);
        gVec.set(0, 0, g1);
        gVec.set(1, 0, g2);
        return gVec;
    }

    public static void main(String[] args) {
        Params p = new Params();
        p.g = 9.81;
        p.l1 = 1.0; p.l2 = 1.0;
        p.lc1 = 0.5; p.lc2 = 0.5;
        p.m1 = 2.0; p.m2 = 1.0;
        p.I1 = 0.1; p.I2 = 0.1;

        SimpleMatrix q = new SimpleMatrix(2, 1, true, new double[]{0.3, -0.7});
        SimpleMatrix M = M_sym(q, p);
        SimpleMatrix gVec = g_sym(q, p);

        System.out.println("M(q) = ");
        M.print();
        System.out.println("g(q) = ");
        gVec.print();
    }
}
      

As in Python and C++, numeric finite-difference checks can be added by repeatedly calling a kinetic-energy function and using central differences.

MATLAB offers a symbolic toolbox and tight integration with Simulink. A common workflow:

1. Derive symbolic \( \mathbf{M}(\mathbf{q}) \), \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \), \( \mathbf{g}(\mathbf{q}) \).
2. Generate a MATLAB function file using matlabFunction.
3. Use the generated function inside a Simulink MATLAB Function block to build a dynamics subsystem.

syms q1 q2 dq1 dq2 real
syms g l1 l2 lc1 lc2 m1 m2 I1 I2 real

q  = [q1; q2];
dq = [dq1; dq2];

M11 = I1 + I2 + m1*lc1^2 + m2*(l1^2 + lc2^2 + 2*l1*lc2*cos(q2));
M12 = I2 + m2*(lc2^2 + l1*lc2*cos(q2));
M22 = I2 + m2*lc2^2;

M = [M11, M12;
     M12, M22];

T = sym(0.5) * dq.' * M * dq;
V = (m1*lc1 + m2*l1)*g*cos(q1) + m2*lc2*g*cos(q1 + q2);

% Gravity vector
g_vec = jacobian(V, q).';

% Christoffel-based C(q, dq)
C = sym(zeros(2,2));
for i = 1:2
  for j = 1:2
    Cij = sym(0);
    for k = 1:2
      c_ijk = 0.5 * ( diff(M(i,j), q(k)) + diff(M(i,k), q(j)) - diff(M(j,k), q(i)) );
      Cij = Cij + c_ijk * dq(k);
    end
    C(i,j) = simplify(Cij);
  end
end

% Generate numeric functions
matlabFunction(M, 'File', 'M_2R', ...
    'Vars', {q1, q2, l1, l2, lc1, lc2, m1, m2, I1, I2});
matlabFunction(C, 'File', 'C_2R', ...
    'Vars', {q1, q2, dq1, dq2, l1, l2, lc1, lc2, m1, m2, I1, I2});
matlabFunction(g_vec, 'File', 'g_2R', ...
    'Vars', {q1, q2, g, l1, l2, lc1, lc2, m1, m2});
      

In Simulink, create a subsystem with states \( \mathbf{q}, \dot{\mathbf{q}} \), call M_2R, C_2R, and g_2R in a MATLAB Function block, and implement the state update \( \ddot{\mathbf{q}} = \mathbf{M}^{-1}(\boldsymbol{\tau} - \mathbf{C}\dot{\mathbf{q}} - \mathbf{g}) \). Numeric finite-difference checks can be implemented in separate scripts.

Mathematica is well suited for concise symbolic derivation and structural analysis (e.g., verifying symmetry and positive definiteness of \( \mathbf{M}(\mathbf{q}) \)). The 2R example can be implemented as:

ClearAll["Global`*"];

(* Parameters and variables *)
{q1, q2, dq1, dq2, g, l1, l2, lc1, lc2, m1, m2, I1, I2} =
  Array[Symbol, 14];

q  = {q1, q2};
dq = {dq1, dq2};

M11 = I1 + I2 + m1 lc1^2 + m2 (l1^2 + lc2^2 + 2 l1 lc2 Cos[q2]);
M12 = I2 + m2 (lc2^2 + l1 lc2 Cos[q2]);
M22 = I2 + m2 lc2^2;

M = { {M11, M12},
     {M12, M22} };

T = 1/2 dq . M . dq;
V = (m1 lc1 + m2 l1) g Cos[q1] + m2 lc2 g Cos[q1 + q2];

(* Gravity vector *)
gVec = Grad[V, q];

(* Christoffel symbols and C(q, dq) *)
c[i_, j_, k_] :=
  1/2 (D[M[[i, j]], q[[k]]] +
       D[M[[i, k]], q[[j]]] -
       D[M[[j, k]], q[[i]]]);

Cmat =
  Table[
    Sum[c[i, j, k] dq[[k]], {k, 1, 2}],
    {i, 1, 2}, {j, 1, 2}
  ] // Simplify;

(* Check structural properties *)
Simplify[Transpose[M] - M]         (* symmetry check *)
FullSimplify[dq . (M - Transpose[M]) . dq]  (* should be 0 *)

(* Example numeric substitution *)
sub = {g -> 9.81, l1 -> 1, l2 -> 1, lc1 -> 0.5, lc2 -> 0.5,
       m1 -> 2, m2 -> 1, I1 -> 0.1, I2 -> 0.1,
       q1 -> 0.3, q2 -> -0.7, dq1 -> 0.0, dq2 -> 0.0};

N[M /. sub]
N[gVec /. sub]
      

Here, symbolic checks like Simplify[Transpose[M] - M] and FullSimplify[dq . (M - Transpose[M]) . dq] verify the expected symmetry property. Similar code can be extended to higher-DOF manipulators.

Finite-difference approximations balance two sources of error:

• Truncation error: for central differences it scales as \( \mathcal{O}(h^2) \).
• Roundoff error: subtracting nearly equal quantities amplifies floating-point noise, growing like \( \mathcal{O}(\varepsilon / h) \), where \( \varepsilon \) is machine precision.

For a scalar derivative, a classical heuristic is

\[ h^\star \approx \sqrt{\varepsilon} \, \max(1,|x|), \]

which minimizes a model of the sum of truncation and roundoff errors. In our context, a component-wise step choice \( h_i^\star \approx \sqrt{\varepsilon}\max(1,|q_i|) \) works well for \( \mathbf{M}(\mathbf{q}) \) and \( \mathbf{g}(\mathbf{q}) \).

For well-conditioned manipulators, symbolic and numeric derivatives will match up to numerical precision. For near-singular configurations, numeric differentiation can be noisy, and symbolic simplifications become more valuable.

In this lab you implemented both symbolic and numeric derivations of Lagrange–Euler dynamics for a simple 2R manipulator. Symbolic methods (SymPy, MATLAB, Mathematica) provide exact expressions, enable structural verification (symmetry, positive definiteness, passivity), and allow code generation for efficient numeric evaluation. Numeric methods, based on finite differences or automatic differentiation, trade some accuracy for generality and ease of implementation on complex models where explicit expressions become unwieldy. In real robotic systems, symbolic and numeric approaches are often combined: symbolic models for controller cores, and numeric derivatives for validation, identification, and rapid prototyping.