In robot modeling, the configuration of an \(n\)-DOF serial manipulator is represented by a joint vector \( \mathbf{q} \in \mathbb{R}^n \). A standard kinematic map sends joint coordinates to a task-space quantity (for example, end-effector position):

\[ f : \mathbb{R}^n \to \mathbb{R}^m, \quad \mathbf{q} \mapsto \mathbf{x} = f(\mathbf{q}). \]

In addition, we often define scalar fields on configuration or task space, such as potential functions, distance costs, or energy:

\[ \varphi : \mathbb{R}^n \to \mathbb{R}, \quad \mathbf{q} \mapsto \varphi(\mathbf{q}). \]

For modeling, analysis, and later control and optimization, we need:

• The gradient \( \nabla \varphi(\mathbf{q}) \): local sensitivity of a scalar field to changes in joints.
• The Jacobian \( J_f(\mathbf{q}) \): linear map relating small joint variations \( \Delta \mathbf{q} \) to task variations \( \Delta \mathbf{x} \approx J_f(\mathbf{q}) \Delta \mathbf{q} \).
• The Hessian \( H_\varphi(\mathbf{q}) \): second-order information describing curvature.

flowchart TD
  Q["q in R^n (joint space)"] --> F["Task map: x = f(q)"]
  F --> SF["Scalar field phi(x) (cost, energy)"]
  SF --> GQ["Gradient grad_q phi(q)"]
  F --> J["Jacobian J(q) = d f / d q"]
  GQ --> H["Hessian H(q) = d^2 phi / d q^2"]
  

In this lesson we work purely at the mathematical level (no control laws yet), but always keep in mind that these objects will eventually be evaluated on concrete robot models.

Let \( f : \mathbb{R}^n \to \mathbb{R} \) be a scalar field, e.g., a cost on joint space. Fix a point \( \mathbf{x}_0 \in \mathbb{R}^n \). The partial derivative with respect to coordinate \(x_i\) is

\[ \frac{\partial f}{\partial x_i}(\mathbf{x}_0) = \lim_{h \to 0} \frac{f(\mathbf{x}_0 + h \mathbf{e}_i) - f(\mathbf{x}_0)}{h}, \]

where \( \mathbf{e}_i \) is the \(i\)-th standard basis vector of \( \mathbb{R}^n \). Collecting all partial derivatives gives the gradient:

\[ \nabla f(\mathbf{x}_0) = \begin{bmatrix} \dfrac{\partial f}{\partial x_1}(\mathbf{x}_0) \\ \vdots \\ \dfrac{\partial f}{\partial x_n}(\mathbf{x}_0) \end{bmatrix} \in \mathbb{R}^n. \]

A key formal fact is that if \( f \) is differentiable at \( \mathbf{x}_0 \), then for small \( \Delta \mathbf{x} \) we have the first-order approximation

\[ f(\mathbf{x}_0 + \Delta \mathbf{x}) = f(\mathbf{x}_0) + \nabla f(\mathbf{x}_0)^{\top} \Delta \mathbf{x} + o(\|\Delta \mathbf{x}\|), \]

where \( o(\|\Delta \mathbf{x}\|) \) denotes higher-order terms small compared to \( \|\Delta \mathbf{x}\| \). In robot modeling this means that, locally, a scalar field on joint space behaves like an affine function with slope given by the gradient.

Let \( f : \mathbb{R}^n \to \mathbb{R} \) and \( \mathbf{u} \in \mathbb{R}^n \) with \( \|\mathbf{u}\|_2 = 1 \). The directional derivative of \( f \) at \( \mathbf{x}_0 \) in direction \( \mathbf{u} \) is

\[ D_{\mathbf{u}} f(\mathbf{x}_0) = \lim_{h \to 0} \frac{ f(\mathbf{x}_0 + h \mathbf{u}) - f(\mathbf{x}_0) }{h}, \]

provided the limit exists. If \( f \) is differentiable at \( \mathbf{x}_0 \), we can use the linear approximation from the previous section with \( \Delta \mathbf{x} = h \mathbf{u} \):

\[ f(\mathbf{x}_0 + h \mathbf{u}) = f(\mathbf{x}_0) + \nabla f(\mathbf{x}_0)^{\top} (h \mathbf{u}) + o(|h|). \]

Subtracting \( f(\mathbf{x}_0) \), dividing by \( h \neq 0 \), and taking the limit yields

\[ D_{\mathbf{u}} f(\mathbf{x}_0) = \nabla f(\mathbf{x}_0)^{\top} \mathbf{u}. \]

By the Cauchy–Schwarz inequality,

\[ |D_{\mathbf{u}} f(\mathbf{x}_0)| = |\nabla f(\mathbf{x}_0)^{\top} \mathbf{u}| \leq \|\nabla f(\mathbf{x}_0)\|_2 \, \|\mathbf{u}\|_2 = \|\nabla f(\mathbf{x}_0)\|_2. \]

The maximum value is attained when \( \mathbf{u} \) is parallel to \( \nabla f(\mathbf{x}_0) \). Thus:

• The gradient points in the direction of steepest ascent of \( f \).
• \( -\nabla f(\mathbf{x}_0) \) points in the direction of steepest descent (used later in optimization).

In robot configuration space, this means that if \( \varphi(\mathbf{q}) \) measures, e.g., a collision penalty, then \( -\nabla \varphi(\mathbf{q}) \) gives the joint-space direction that most rapidly reduces that penalty locally.

Let \( f : \mathbb{R}^n \to \mathbb{R}^m \) be a vector-valued map with components \( f_1, \dots, f_m \). The Jacobian matrix of \( f \) at \( \mathbf{x}_0 \) is

\[ J_f(\mathbf{x}_0) = \begin{bmatrix} \dfrac{\partial f_1}{\partial x_1}(\mathbf{x}_0) & \cdots & \dfrac{\partial f_1}{\partial x_n}(\mathbf{x}_0) \\ \vdots & \ddots & \vdots \\ \dfrac{\partial f_m}{\partial x_1}(\mathbf{x}_0) & \cdots & \dfrac{\partial f_m}{\partial x_n}(\mathbf{x}_0) \end{bmatrix} \in \mathbb{R}^{m \times n}. \]

The first-order approximation of \( f \) around \( \mathbf{x}_0 \) is

\[ f(\mathbf{x}_0 + \Delta \mathbf{x}) = f(\mathbf{x}_0) + J_f(\mathbf{x}_0) \, \Delta \mathbf{x} + o(\|\Delta \mathbf{x}\|). \]

If \( f \) is scalar-valued, \( m = 1 \), then the Jacobian is a \( 1 \times n \) row vector and relates directly to the gradient:

\[ J_f(\mathbf{x}) = \nabla f(\mathbf{x})^{\top}. \]

For compositions, the matrix chain rule is fundamental. Let \( g : \mathbb{R}^k \to \mathbb{R}^n \), \( f : \mathbb{R}^n \to \mathbb{R}^m \), and \( h = f \circ g : \mathbb{R}^k \to \mathbb{R}^m \). Then

\[ J_h(\mathbf{u}_0) = J_f(g(\mathbf{u}_0)) \, J_g(\mathbf{u}_0), \]

where dimensions are \( J_f \in \mathbb{R}^{m \times n} \), \( J_g \in \mathbb{R}^{n \times k} \), and \( J_h \in \mathbb{R}^{m \times k} \).

In manipulator kinematics, taking \( f(\mathbf{q}) = \mathbf{x} \) and differentiating with respect to time gives the classical relation \( \dot{\mathbf{x}} = J_f(\mathbf{q}) \dot{\mathbf{q}} \), which will be studied in detail in the later chapter on differential kinematics.

flowchart TD
  X["Joint config q"] --> F["x = f(q) (FK map)"]
  F --> L["Linearization: x + dx = f(q) + J(q) dq"]
  L --> V["Velocity relation: xdot = J(q) qdot"]
        

For a twice differentiable scalar field \( f : \mathbb{R}^n \to \mathbb{R} \), the Hessian matrix at \( \mathbf{x}_0 \) is defined by

\[ H_f(\mathbf{x}_0) = \begin{bmatrix} \dfrac{\partial^2 f}{\partial x_1^2}(\mathbf{x}_0) & \cdots & \dfrac{\partial^2 f}{\partial x_1 \partial x_n}(\mathbf{x}_0) \\ \vdots & \ddots & \vdots \\ \dfrac{\partial^2 f}{\partial x_n \partial x_1}(\mathbf{x}_0) & \cdots & \dfrac{\partial^2 f}{\partial x_n^2}(\mathbf{x}_0) \end{bmatrix}. \]

If all mixed partials are continuous in a neighborhood of \( \mathbf{x}_0 \), Clairaut's (Schwarz's) theorem ensures symmetry:

\[ \frac{\partial^2 f}{\partial x_i \partial x_j}(\mathbf{x}_0) = \frac{\partial^2 f}{\partial x_j \partial x_i}(\mathbf{x}_0), \quad \Rightarrow \quad H_f(\mathbf{x}_0) = H_f(\mathbf{x}_0)^{\top}. \]

The second-order Taylor expansion of \( f \) around \( \mathbf{x}_0 \) is

\[ f(\mathbf{x}_0 + \Delta \mathbf{x}) \approx f(\mathbf{x}_0) + \nabla f(\mathbf{x}_0)^{\top} \Delta \mathbf{x} + \tfrac{1}{2} \Delta \mathbf{x}^{\top} H_f(\mathbf{x}_0) \Delta \mathbf{x}. \]

At a critical point \( \mathbf{x}_\star \) with \( \nabla f(\mathbf{x}_\star) = \mathbf{0} \), the eigenvalues of \( H_f(\mathbf{x}_\star) \) classify the point:

• All eigenvalues positive → strict local minimum.
• All eigenvalues negative → strict local maximum.
• Mixed signs → saddle point.

In robotics, these concepts appear in potential-field methods, energy shaping, and second-order optimization for model fitting or trajectory refinement.

Many energies and costs used in robot modeling are quadratic. Consider \( f : \mathbb{R}^n \to \mathbb{R} \) of the form

\[ f(\mathbf{x}) = \tfrac{1}{2} \mathbf{x}^{\top} Q \mathbf{x} + \mathbf{b}^{\top} \mathbf{x} + c, \]

where \( Q \in \mathbb{R}^{n \times n} \), \( \mathbf{b} \in \mathbb{R}^n \), \( c \in \mathbb{R} \). Taking derivatives:

\[ \nabla f(\mathbf{x}) = \tfrac{1}{2}(Q + Q^{\top}) \mathbf{x} + \mathbf{b}, \quad H_f(\mathbf{x}) = \tfrac{1}{2}(Q + Q^{\top}). \]

If \( Q \) is symmetric, then \( \nabla f(\mathbf{x}) = Q \mathbf{x} + \mathbf{b} \) and \( H_f(\mathbf{x}) = Q \) is constant. Positive definiteness of \( Q \) guarantees a unique minimizer \( \mathbf{x}_\star = -Q^{-1} \mathbf{b} \). This structure underlies inertia matrices, least-squares identification, and many other quadratic models later in the course.

Consider a planar 2-link manipulator with joint angles \( \mathbf{q} = [q_1, q_2]^{\top} \) and link lengths \( \ell_1, \ell_2 \). From introductory robotics, the end-effector position in the plane is

\[ \begin{aligned} x(\mathbf{q}) &= \ell_1 \cos q_1 + \ell_2 \cos(q_1 + q_2), \\ y(\mathbf{q}) &= \ell_1 \sin q_1 + \ell_2 \sin(q_1 + q_2). \end{aligned} \]

We write \( f(\mathbf{q}) = [x(\mathbf{q}),\, y(\mathbf{q})]^{\top} \). The Jacobian \( J_f(\mathbf{q}) \in \mathbb{R}^{2 \times 2} \) has entries \( (J_f)_{ij} = \partial f_i / \partial q_j \):

\[ J_f(\mathbf{q}) = \begin{bmatrix} -\ell_1 \sin q_1 - \ell_2 \sin(q_1 + q_2) & -\ell_2 \sin(q_1 + q_2) \\ \ell_1 \cos q_1 + \ell_2 \cos(q_1 + q_2) & \ell_2 \cos(q_1 + q_2) \end{bmatrix}. \]

Now define a task-space tracking cost for a desired point \( \mathbf{x}_d = [x_d, y_d]^{\top} \):

\[ \varphi(\mathbf{q}) = \tfrac{1}{2} \| f(\mathbf{q}) - \mathbf{x}_d \|_2^2 = \tfrac{1}{2} (f(\mathbf{q}) - \mathbf{x}_d)^{\top} (f(\mathbf{q}) - \mathbf{x}_d). \]

Let \( \mathbf{e}(\mathbf{q}) = f(\mathbf{q}) - \mathbf{x}_d \). Using the chain rule,

\[ \nabla_{\mathbf{q}} \varphi(\mathbf{q}) = J_f(\mathbf{q})^{\top} \, \mathbf{e}(\mathbf{q}). \]

Sketch of derivation: View \( \varphi = g \circ f \) with \( g(\mathbf{x}) = \tfrac{1}{2} \|\mathbf{x} - \mathbf{x}_d\|_2^2 \). Then \( \nabla_{\mathbf{x}} g(\mathbf{x}) = \mathbf{x} - \mathbf{x}_d \). The chain rule for gradients of compositions gives

\[ \nabla_{\mathbf{q}} \varphi(\mathbf{q}) = J_f(\mathbf{q})^{\top} \, \nabla_{\mathbf{x}} g(f(\mathbf{q})) = J_f(\mathbf{q})^{\top} (f(\mathbf{q}) - \mathbf{x}_d). \]

The Hessian of \( \varphi \) is

\[ H_{\varphi}(\mathbf{q}) = J_f(\mathbf{q})^{\top} J_f(\mathbf{q}) + \sum_{k=1}^{2} e_k(\mathbf{q}) \, H_{f_k}(\mathbf{q}), \]

where \( H_{f_k} \) is the Hessian of the \(k\)-th component \( f_k \). In practice, a common Gauss–Newton approximation discards the second term and uses \( H_{\varphi}(\mathbf{q}) \approx J_f(\mathbf{q})^{\top} J_f(\mathbf{q}) \), which is always symmetric positive semidefinite.

Python offers symbolic and numeric tools that are convenient for experimenting with gradients, Jacobians, and Hessians. Here we use sympy for exact expressions and numpy for numerical evaluation. Robotic modeling libraries like robopy or the Python port of the Robotics Toolbox can later supply f(\( \mathbf{q} \) ) while you still differentiate symbolically.

import sympy as sp
import numpy as np

# Symbolic variables
q1, q2, l1, l2, xd, yd = sp.symbols('q1 q2 l1 l2 xd yd', real=True)

# Forward kinematics for planar 2R
x = l1*sp.cos(q1) + l2*sp.cos(q1 + q2)
y = l1*sp.sin(q1) + l2*sp.sin(q1 + q2)

q = sp.Matrix([q1, q2])
f = sp.Matrix([x, y])

# Jacobian J(q)
J = f.jacobian(q)

# Task-space error and cost
e = f - sp.Matrix([xd, yd])
phi = sp.Rational(1, 2) * (e.dot(e))

# Gradient and Hessian in joint space
grad_phi = sp.Matrix([sp.diff(phi, q1), sp.diff(phi, q2)])
H_phi = sp.hessian(phi, q)

print("J(q) =")
sp.pprint(J)
print("grad_phi(q) =")
sp.pprint(grad_phi)
print("H_phi(q) =")
sp.pprint(H_phi)

# Numeric evaluation
subs_vals = {
    l1: 1.0, l2: 0.8,
    xd: 1.2, yd: 0.3,
    q1: 0.5, q2: -0.3
}
J_num = np.array(J.evalf(subs=subs_vals), dtype=float)
grad_num = np.array(grad_phi.evalf(subs=subs_vals), dtype=float).reshape(-1)
H_num = np.array(H_phi.evalf(subs=subs_vals), dtype=float)

print("J_num =", J_num)
print("grad_num =", grad_num)
print("H_num =", H_num)
      

Later, once you define more complex kinematic chains via a robotics library, you can either differentiate symbolically (for small models) or use automatic differentiation frameworks that integrate with NumPy.

In C++, the Eigen library is standard for linear algebra in robotics. Below is an analytic implementation of the Jacobian and the gradient of the 2R cost \( \varphi(\mathbf{q}) = \tfrac{1}{2} \| f(\mathbf{q}) - \mathbf{x}_d \|_2^2 \). More advanced projects often combine Eigen with model libraries such as KDL, RBDL, or Pinocchio.

#include <Eigen/Dense>
#include <cmath>

using Vec2 = Eigen::Vector2d;
using Mat2 = Eigen::Matrix2d;

// Forward kinematics
Vec2 fk2R(double q1, double q2, double l1, double l2) {
    double c1 = std::cos(q1);
    double s1 = std::sin(q1);
    double c12 = std::cos(q1 + q2);
    double s12 = std::sin(q1 + q2);

    double x = l1 * c1 + l2 * c12;
    double y = l1 * s1 + l2 * s12;
    return Vec2(x, y);
}

// Analytic Jacobian
Mat2 jac2R(double q1, double q2, double l1, double l2) {
    double c1 = std::cos(q1);
    double s1 = std::sin(q1);
    double c12 = std::cos(q1 + q2);
    double s12 = std::sin(q1 + q2);

    Mat2 J;
    J(0,0) = -l1 * s1 - l2 * s12;
    J(0,1) = -l2 * s12;
    J(1,0) =  l1 * c1 + l2 * c12;
    J(1,1) =  l2 * c12;
    return J;
}

// Gradient of phi(q) = 0.5 * ||f(q) - x_d||^2
Vec2 gradPhi2R(const Vec2& q, const Vec2& x_d,
               double l1, double l2) {
    Vec2 x = fk2R(q(0), q(1), l1, l2);
    Vec2 e = x - x_d;
    Mat2 J = jac2R(q(0), q(1), l1, l2);
    // grad_phi = J(q)^T * e
    return J.transpose() * e;
}
      

For second-order methods, one can implement finite-difference approximations of the Hessian using repeated evaluations of gradPhi2R. Libraries such as CppAD and CasADi provide automatic differentiation back-ends that compute Jacobians and Hessians for complex robot models.

In Java, the EJML library offers efficient matrix operations. Below is an example using EJML's SimpleMatrix interface to compute the Jacobian and gradient for the same 2R example.

import org.ejml.simple.SimpleMatrix;

public class TwoRGrad {

    public static SimpleMatrix fk2R(double q1, double q2,
                                    double l1, double l2) {
        double c1 = Math.cos(q1);
        double s1 = Math.sin(q1);
        double c12 = Math.cos(q1 + q2);
        double s12 = Math.sin(q1 + q2);

        double x = l1 * c1 + l2 * c12;
        double y = l1 * s1 + l2 * s12;

        return new SimpleMatrix(2, 1, true, new double[]{x, y});
    }

    public static SimpleMatrix jac2R(double q1, double q2,
                                     double l1, double l2) {
        double c1 = Math.cos(q1);
        double s1 = Math.sin(q1);
        double c12 = Math.cos(q1 + q2);
        double s12 = Math.sin(q1 + q2);

        double[] data = new double[]{
            -l1 * s1 - l2 * s12,  -l2 * s12,
             l1 * c1 + l2 * c12,   l2 * c12
        };
        return new SimpleMatrix(2, 2, true, data);
    }

    // grad_phi = J(q)^T * (f(q) - x_d)
    public static SimpleMatrix gradPhi2R(SimpleMatrix q,
                                         SimpleMatrix x_d,
                                         double l1, double l2) {
        double q1 = q.get(0);
        double q2 = q.get(1);

        SimpleMatrix x = fk2R(q1, q2, l1, l2);
        SimpleMatrix e = x.minus(x_d);
        SimpleMatrix J = jac2R(q1, q2, l1, l2);
        return J.transpose().mult(e);
    }

    public static void main(String[] args) {
        double l1 = 1.0, l2 = 0.8;
        SimpleMatrix q = new SimpleMatrix(2, 1, true, new double[]{0.5, -0.3});
        SimpleMatrix x_d = new SimpleMatrix(2, 1, true, new double[]{1.2, 0.3});

        SimpleMatrix g = gradPhi2R(q, x_d, l1, l2);
        System.out.println("grad_phi(q) = ");
        g.print();
    }
}
      

For larger robots, Java-based simulation environments can use the same pattern: symbolic or analytic expressions for kinematics and then EJML for evaluating Jacobians and gradients numerically.

MATLAB is ubiquitous in control and robotics. Symbolic derivatives can be obtained with the Symbolic Math Toolbox, while the Robotics System Toolbox provides kinematic maps for standard manipulators. Below is a small script computing the Jacobian, gradient, and Gauss–Newton Hessian approximation for the planar 2R example.

syms q1 q2 l1 l2 xd yd real

% Forward kinematics
x = l1*cos(q1) + l2*cos(q1 + q2);
y = l1*sin(q1) + l2*sin(q1 + q2);

q = [q1; q2];
f = [x; y];

% Jacobian J(q)
J = jacobian(f, q);

% Cost phi(q) = 0.5 * ||f(q) - x_d||^2
xd_sym = sym('xd');
yd_sym = sym('yd');
e = f - [xd_sym; yd_sym];
phi = sym(0.5) * (e.' * e);

% Gradient and Hessian
grad_phi = gradient(phi, q);
H_phi = hessian(phi, q);

disp('J(q) ='); disp(J);
disp('grad_phi(q) ='); disp(grad_phi);
disp('H_phi(q) ='); disp(H_phi);

% Numeric evaluation
l1_val = 1.0; l2_val = 0.8;
xd_val = 1.2; yd_val = 0.3;
q_val = [0.5; -0.3];

J_num = double(subs(J, {l1, l2, xd_sym, yd_sym, q1, q2}, ...
                      {l1_val, l2_val, xd_val, yd_val, q_val(1), q_val(2)}));
grad_num = double(subs(grad_phi, {l1, l2, xd_sym, yd_sym, q1, q2}, ...
                                {l1_val, l2_val, xd_val, yd_val, q_val(1), q_val(2)}));

disp('J_num ='); disp(J_num);
disp('grad_num ='); disp(grad_num);
      

Simulink realization. In Simulink, one can represent the kinematic map \( f(\mathbf{q}) \) as a subsystem that takes \( q_1, q_2 \) (and parameters \( \ell_1, \ell_2 \)) as inputs and outputs \( x, y \). A separate subsystem can implement the analytic Jacobian J(q) and compute grad_phi by matrix multiplication, which is useful later when embedding these quantities in dynamic simulation and control diagrams.

Mathematica provides high-level operators for gradients, Jacobians, and Hessians that are convenient for symbolic exploration and verification of robot models.

(* Variables *)
Clear[q1, q2, l1, l2, xd, yd];
f[q1_, q2_, l1_, l2_] := {
  l1 Cos[q1] + l2 Cos[q1 + q2],
  l1 Sin[q1] + l2 Sin[q1 + q2]
};

(* Forward kinematics as a function of q *)
xvec[q1_, q2_] := f[q1, q2, l1, l2];

(* Jacobian with respect to q = {q1, q2} *)
J[q1_, q2_] := D[xvec[q1, q2], {{q1, q2}}];

(* Cost phi(q) = 0.5 ||f(q) - x_d||^2 *)
phi[q1_, q2_] := Module[{x, e},
  x = xvec[q1, q2];
  e = x - {xd, yd};
  1/2 e.e
];

(* Gradient and Hessian *)
gradPhi[q1_, q2_] := Grad[phi[q1, q2], {q1, q2}];
HessianPhi[q1_, q2_] := D[phi[q1, q2], { {q1, q2}, 2 }];

(* Example numeric evaluation *)
l1 = 1.0; l2 = 0.8; xd = 1.2; yd = 0.3;
q1val = 0.5; q2val = -0.3;

Jnum = J[q1val, q2val] /. {l1 -> 1.0, l2 -> 0.8};
gradnum = gradPhi[q1val, q2val] /. {l1 -> 1.0, l2 -> 0.8, xd -> 1.2, yd -> 0.3};

Jnum // MatrixForm
gradnum
      

Such symbolic computations are especially helpful when deriving and checking analytic Jacobians and Hessians for more complex manipulators before implementing numeric versions.

In this lesson we rigorously introduced multivariable calculus objects that will reappear throughout the course:

• Gradients as stacked partial derivatives and as the unique vectors encoding all directional derivatives via \( D_{\mathbf{u}} f = \nabla f^{\top} \mathbf{u} \).
• Jacobians of vector-valued maps as linearizations \( f(\mathbf{x} + \Delta \mathbf{x}) \approx f(\mathbf{x}) + J_f(\mathbf{x}) \Delta \mathbf{x} \), with the matrix chain rule.
• Hessians as matrices of second derivatives, governing curvature and second-order Taylor expansions.
• Concrete manipulator examples, including the planar 2R arm, where gradients, Jacobians, and Hessians of a simple task-space tracking cost were computed.
• Basic implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica, using widely deployed numerical and symbolic libraries relevant for robotics.

These tools are the analytic backbone for later chapters on differential kinematics, energy-based modeling, and, eventually, dynamics and control design.