Consider an \(n\)-DOF serial manipulator with joint coordinates \( \mathbf{q}(t) \in \mathbb{R}^n \) and a task-space description (e.g. end-effector position or pose) \( \mathbf{x}(t) \in \mathbb{R}^m \). The forward kinematics from Chapter 5 and the differential kinematics from Chapter 7 can be summarized as

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

The Jacobian is the matrix of first partial derivatives of the forward kinematics:

\[ \mathbf{J}(\mathbf{q}) = \frac{\partial f(\mathbf{q})}{\partial \mathbf{q}} = \begin{bmatrix} \dfrac{\partial f_1}{\partial q_1} & \cdots & \dfrac{\partial f_1}{\partial q_n} \\ \vdots & \ddots & \vdots \\ \dfrac{\partial f_m}{\partial q_1} & \cdots & \dfrac{\partial f_m}{\partial q_n} \end{bmatrix}. \]

Differentiating the forward kinematics with respect to time (chain rule) yields the fundamental velocity mapping

\[ \dot{\mathbf{x}}(t) = \mathbf{J}(\mathbf{q}(t)) \, \dot{\mathbf{q}}(t). \]

Component-wise, for each task coordinate \( x_i \),

\[ \dot{x}_i(t) = \sum_{j=1}^{n} \frac{\partial f_i(\mathbf{q}(t))}{\partial q_j} \, \dot{q}_j(t), \quad i = 1,\dots,m. \]

For full 6D end-effector motion in \(SE(3)\), we often group linear and angular velocity into a twist

\[ \mathbf{V}_s(t) = \begin{bmatrix} \mathbf{v}_s(t) \\[4pt] \boldsymbol{\omega}_s(t) \end{bmatrix} \in \mathbb{R}^6, \quad \mathbf{V}_s(t) = \mathbf{J}_s(\mathbf{q}(t)) \, \dot{\mathbf{q}}(t), \]

where \( \mathbf{J}_s(\mathbf{q}) \) is the spatial Jacobian (Chapter 7). An analogous relation holds for the body Jacobian \( \mathbf{J}_b(\mathbf{q}) \).

flowchart TD
  Q["Joint trajectory q(t)"] --> FK["Forward kinematics x = f(q)"]
  FK --> J["Jacobian J(q) = df/dq"]
  J --> V["Velocity mapping x_dot = J(q) * q_dot"]
  V --> A["Acceleration mapping x_ddot = J(q) * q_ddot + Jdot(q,q_dot) * q_dot"]
  A --> USE["Analysis of trajectories, feasibility, dynamic consistency"]
        

Given a sufficiently smooth trajectory \( \mathbf{q}(t) \), the task acceleration \( \ddot{\mathbf{x}}(t) \) follows from differentiating \( \dot{\mathbf{x}} = \mathbf{J}(\mathbf{q}) \dot{\mathbf{q}} \):

\[ \ddot{\mathbf{x}}(t) = \frac{\mathrm{d}}{\mathrm{d}t}\big(\mathbf{J}(\mathbf{q}(t)) \dot{\mathbf{q}}(t)\big) = \dot{\mathbf{J}}(\mathbf{q}(t),\dot{\mathbf{q}}(t)) \, \dot{\mathbf{q}}(t) + \mathbf{J}(\mathbf{q}(t)) \, \ddot{\mathbf{q}}(t). \]

This is the fundamental acceleration mapping. The term \( \dot{\mathbf{J}}(\mathbf{q},\dot{\mathbf{q}}) \dot{\mathbf{q}} \) is purely geometric: it depends only on the kinematics and joint velocities, and not on mass or forces. It is kinematically analogous to the Coriolis and centrifugal terms in dynamics.

To understand \( \dot{\mathbf{J}} \), we can differentiate each Jacobian entry using the chain rule. For each component \( J_{ij}(\mathbf{q}) \),

\[ \dot{J}_{ij}(\mathbf{q},\dot{\mathbf{q}}) = \sum_{k=1}^{n} \frac{\partial J_{ij}(\mathbf{q})}{\partial q_k} \, \dot{q}_k, \quad i=1,\dots,m,\; j=1,\dots,n. \]

In tensor notation, this can be written as \( \dot{\mathbf{J}} = \sum_{k} \frac{\partial \mathbf{J}}{\partial q_k} \dot{q}_k \). Although \( \dot{\mathbf{J}} \) is cumbersome to store as a true third-order tensor, we rarely need the tensor itself; typically we need the product \( \dot{\mathbf{J}}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \).

An alternative expression uses kinematic Hessians. Writing \( \mathbf{x}(t) = f(\mathbf{q}(t)) \), a second time derivative yields

\[ \ddot{x}_i(t) = \sum_{j=1}^{n} \frac{\partial f_i}{\partial q_j} \ddot{q}_j(t) + \sum_{j=1}^{n}\sum_{k=1}^{n} \frac{\partial^2 f_i}{\partial q_j \partial q_k} \dot{q}_j(t) \dot{q}_k(t), \quad i=1,\dots,m. \]

Collecting terms into matrix form,

\[ \ddot{\mathbf{x}}(t) = \mathbf{J}(\mathbf{q}(t)) \ddot{\mathbf{q}}(t) + \mathbf{h}(\mathbf{q}(t),\dot{\mathbf{q}}(t)), \]

where \( \mathbf{h}(\mathbf{q},\dot{\mathbf{q}}) \) stacks all the \( \sum_{j,k} \frac{\partial^2 f_i}{\partial q_j \partial q_k} \dot{q}_j \dot{q}_k \) terms. By comparing with the previous expression, we see

\[ \mathbf{h}(\mathbf{q},\dot{\mathbf{q}}) = \dot{\mathbf{J}}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}. \]

This equivalence is extremely useful: one can compute \( \ddot{\mathbf{x}} \) either via explicit Hessians or by first forming \( \mathbf{J} \) and then differentiating it.

Consider a planar 2R arm with link lengths \( \ell_1, \ell_2 \) and joint angles \( q_1, q_2 \). The end-effector point \( P \) in the base frame has coordinates

\[ x(q_1,q_2) = \ell_1 \cos(q_1) + \ell_2 \cos(q_1 + q_2), \quad y(q_1,q_2) = \ell_1 \sin(q_1) + \ell_2 \sin(q_1 + q_2). \]

The Jacobian (mapping joint velocity to end-effector linear velocity in the plane) is

\[ \mathbf{J}(\mathbf{q}) = \begin{bmatrix} \dfrac{\partial x}{\partial q_1} & \dfrac{\partial x}{\partial q_2} \\ \dfrac{\partial y}{\partial q_1} & \dfrac{\partial y}{\partial q_2} \end{bmatrix} = \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}. \]

Let \( \dot{q}_1, \dot{q}_2 \) denote joint velocities. Then \( \dot{\mathbf{x}} = [\dot{x},\dot{y}]^\top = \mathbf{J}(\mathbf{q}) \dot{\mathbf{q}} \). Differentiating \( \mathbf{J}(\mathbf{q}) \) with respect to time gives

\[ \dot{\mathbf{J}}(\mathbf{q},\dot{\mathbf{q}}) = \begin{bmatrix} -\dot{q}_1 \ell_1 \cos(q_1) -(\dot{q}_1+\dot{q}_2)\ell_2 \cos(q_1+q_2) & -(\dot{q}_1+\dot{q}_2)\ell_2 \cos(q_1+q_2) \\[4pt] -\dot{q}_1 \ell_1 \sin(q_1) -(\dot{q}_1+\dot{q}_2)\ell_2 \sin(q_1+q_2) & -(\dot{q}_1+\dot{q}_2)\ell_2 \sin(q_1+q_2) \end{bmatrix}. \]

For given joint accelerations \( \ddot{q}_1, \ddot{q}_2 \), the end-effector acceleration is

\[ \ddot{\mathbf{x}} = \mathbf{J}(\mathbf{q}) \ddot{\mathbf{q}} + \dot{\mathbf{J}}(\mathbf{q},\dot{\mathbf{q}}) \dot{\mathbf{q}}. \]

Notice how the nonlinear coupling in \( q_1, q_2 \) and \( \dot{q}_1, \dot{q}_2 \) appears explicitly in the second term. Even though this example is planar and low-DOF, the same structure extends to full 6D end-effector motion in spatial or body frames.

flowchart TD
  IN["Input: q1, q2, q1_dot, q2_dot, q1_ddot, q2_ddot"] --> JSTEP["Compute J(q1,q2)"]
  JSTEP --> JDSTEP["Compute Jdot(q, q_dot) via chain rule"]
  JDSTEP --> VSTEP["Compute x_dot = J * q_dot"]
  VSTEP --> ASTEP["Compute x_ddot = J * q_ddot + Jdot * q_dot"]
  ASTEP --> OUT["Output: planar end-effector velocity and acceleration"]
        

We now implement the 2R velocity and acceleration mappings in Python. First we code them directly, then we indicate how a robotics library can be used to obtain Jacobians for more complex robots.

import numpy as np

def fk_2R(q, l1, l2):
    """
    Planar 2R forward kinematics (position only).
    q: array-like of shape (2,) with [q1, q2]
    """
    q1, q2 = q
    x = l1 * np.cos(q1) + l2 * np.cos(q1 + q2)
    y = l1 * np.sin(q1) + l2 * np.sin(q1 + q2)
    return np.array([x, y])

def jacobian_2R(q, l1, l2):
    q1, q2 = q
    s1 = np.sin(q1)
    c1 = np.cos(q1)
    s12 = np.sin(q1 + q2)
    c12 = np.cos(q1 + q2)
    J = np.array([
        [-l1 * s1 - l2 * s12, -l2 * s12],
        [ l1 * c1 + l2 * c12,  l2 * c12]
    ])
    return J

def jdot_2R(q, dq, l1, l2):
    q1, q2 = q
    dq1, dq2 = dq
    c1 = np.cos(q1)
    s1 = np.sin(q1)
    c12 = np.cos(q1 + q2)
    s12 = np.sin(q1 + q2)
    # Derived symbolically in Section 3
    Jdot = np.array([
        [-dq1 * l1 * c1 - (dq1 + dq2) * l2 * c12,
         -(dq1 + dq2) * l2 * c12],
        [-dq1 * l1 * s1 - (dq1 + dq2) * l2 * s12,
         -(dq1 + dq2) * l2 * s12]
    ])
    return Jdot

def task_velocity_acceleration_2R(q, dq, ddq, l1, l2):
    J = jacobian_2R(q, l1, l2)
    Jdot = jdot_2R(q, dq, l1, l2)
    dq = np.asarray(dq).reshape(2,)
    ddq = np.asarray(ddq).reshape(2,)
    xdot = J @ dq
    xddot = J @ ddq + Jdot @ dq
    return xdot, xddot

if __name__ == "__main__":
    l1, l2 = 1.0, 0.7
    q = np.deg2rad([30.0, 20.0])
    dq = np.array([0.5, -0.3])
    ddq = np.array([0.2, 0.1])

    x = fk_2R(q, l1, l2)
    xdot, xddot = task_velocity_acceleration_2R(q, dq, ddq, l1, l2)

    print("End-effector position:", x)
    print("End-effector velocity:", xdot)
    print("End-effector acceleration:", xddot)
      

For more complex robots, one can use a Python robotics library such as roboticstoolbox-python. For example, with a DH-based model:

from roboticstoolbox import DHRobot, RevoluteDH
import numpy as np

# Define a simple 2R planar arm in the xy-plane
l1, l2 = 1.0, 0.7
robot = DHRobot([
    RevoluteDH(a=l1, alpha=0.0, d=0.0),
    RevoluteDH(a=l2, alpha=0.0, d=0.0)
], name="planar2R")

q = np.deg2rad([30.0, 20.0])
dq = np.array([0.5, -0.3])
ddq = np.array([0.2, 0.1])

# Spatial Jacobian at the end-effector
J6 = robot.jacob0(q)  # 6x2: linear and angular velocity
xdot6 = J6 @ dq
# For pure planar motion, use the first two rows of J6 for xy-linear velocity.
      

The same library can provide geometric Jacobians in different frames, enabling a direct link between symbolic theory and numerical implementation.

In C++, Eigen provides efficient matrix primitives, while the Kinematics and Dynamics Library (KDL) in the ROS ecosystem offers high-level Jacobian solvers. Below is a minimal Eigen-based implementation for the 2R example.

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

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

Matrix2d jacobian2R(double q1, double q2, double l1, double l2) {
    double s1  = std::sin(q1);
    double c1  = std::cos(q1);
    double s12 = std::sin(q1 + q2);
    double c12 = std::cos(q1 + q2);
    Matrix2d 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;
}

Matrix2d jdot2R(double q1, double q2,
                double dq1, double dq2,
                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);
    Matrix2d Jdot;
    Jdot(0,0) = -dq1 * l1 * c1 - (dq1 + dq2) * l2 * c12;
    Jdot(0,1) = -(dq1 + dq2) * l2 * c12;
    Jdot(1,0) = -dq1 * l1 * s1 - (dq1 + dq2) * l2 * s12;
    Jdot(1,1) = -(dq1 + dq2) * l2 * s12;
    return Jdot;
}

int main() {
    double l1 = 1.0, l2 = 0.7;
    double q1 = 30.0 * M_PI / 180.0;
    double q2 = 20.0 * M_PI / 180.0;
    double dq1 = 0.5, dq2 = -0.3;
    double ddq1 = 0.2, ddq2 = 0.1;

    Matrix2d J = jacobian2R(q1, q2, l1, l2);
    Matrix2d Jdot = jdot2R(q1, q2, dq1, dq2, l1, l2);

    Vector2d dq(dq1, dq2);
    Vector2d ddq(ddq1, ddq2);

    Vector2d xdot = J * dq;
    Vector2d xddot = J * ddq + Jdot * dq;

    std::cout << "xdot = " << xdot.transpose() << std::endl;
    std::cout << "xddot = " << xddot.transpose() << std::endl;
    return 0;
}
      

To use KDL for more general manipulators, one typically constructs a Chain from URDF or DH parameters and calls a Jacobian solver:

#include <kdl/chain.hpp>
#include <kdl/chainjnttojacsolver.hpp>
#include <kdl/jacobian.hpp>

// ...

KDL::Chain chain; // filled from robot description
// ...
KDL::ChainJntToJacSolver jac_solver(chain);
KDL::JntArray q(chain.getNrOfJoints());
// set q(i) from current joint configuration
KDL::Jacobian J(chain.getNrOfJoints());
jac_solver.JntToJac(q, J); // 6 x n geometric Jacobian at the end-effector
      

Finite-difference approximations of \( \dot{\mathbf{J}} \) can be implemented by numerically differentiating \( \mathbf{J}(\mathbf{q}(t)) \) along a sampled trajectory, but analytic differentiation as shown is more accurate when available.

Java does not have a standard robotics library in the core language, but numerical libraries such as EJML provide dense matrix operations. The following code illustrates the 2R mapping using EJML.

import org.ejml.simple.SimpleMatrix;

public class Planar2R {

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

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

    public static SimpleMatrix jdot2R(double q1, double q2,
                                      double dq1, double dq2,
                                      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 = {
            { -dq1 * l1 * c1 - (dq1 + dq2) * l2 * c12,
              -(dq1 + dq2) * l2 * c12 },
            { -dq1 * l1 * s1 - (dq1 + dq2) * l2 * s12,
              -(dq1 + dq2) * l2 * s12 }
        };
        return new SimpleMatrix(data);
    }

    public static void main(String[] args) {
        double l1 = 1.0, l2 = 0.7;
        double q1 = Math.toRadians(30.0);
        double q2 = Math.toRadians(20.0);
        double dq1 = 0.5, dq2 = -0.3;
        double ddq1 = 0.2, ddq2 = 0.1;

        SimpleMatrix J = jacobian2R(q1, q2, l1, l2);
        SimpleMatrix Jdot = jdot2R(q1, q2, dq1, dq2, l1, l2);

        SimpleMatrix dq = new SimpleMatrix(2, 1, true, new double[]{dq1, dq2});
        SimpleMatrix ddq = new SimpleMatrix(2, 1, true, new double[]{ddq1, ddq2});

        SimpleMatrix xdot = J.mult(dq);
        SimpleMatrix xddot = J.mult(ddq).plus(Jdot.mult(dq));

        System.out.println("xdot = " + xdot.transpose());
        System.out.println("xddot = " + xddot.transpose());
    }
}
      

For more sophisticated Java robotics stacks (e.g. ROSJava-based systems), the same mapping structure is preserved; only the Jacobian computation is delegated to a library instead of coded explicitly.

MATLAB is widely used in robotics for both symbolic and numeric kinematics. A direct script implementation of the 2R mappings is:

function [x, xdot, xddot, J, Jdot] = planar2R_kin(q, dq, ddq, l1, l2)
% q, dq, ddq are 2x1 vectors

q1  = q(1);   q2  = q(2);
dq1 = dq(1);  dq2 = dq(2);

x = [l1*cos(q1) + l2*cos(q1 + q2);
     l1*sin(q1) + l2*sin(q1 + q2)];

J = [ -l1*sin(q1) - l2*sin(q1 + q2),   -l2*sin(q1 + q2);
       l1*cos(q1) + l2*cos(q1 + q2),    l2*cos(q1 + q2) ];

Jdot = [ -dq1*l1*cos(q1) - (dq1 + dq2)*l2*cos(q1 + q2),  -(dq1 + dq2)*l2*cos(q1 + q2);
         -dq1*l1*sin(q1) - (dq1 + dq2)*l2*sin(q1 + q2),  -(dq1 + dq2)*l2*sin(q1 + q2) ];

xdot  = J * dq;
xddot = J * ddq + Jdot * dq;
end
      

Using the Robotics System Toolbox, one can define a rigidBodyTree and use geometricJacobian to compute full 6D Jacobians:

% Define 2R planar arm as a rigidBodyTree
robot = rigidBodyTree('DataFormat','row');

body1 = rigidBody('link1');
jnt1  = rigidBodyJoint('j1','revolute');
setFixedTransform(jnt1, trvec2tform([0 0 0]));
jnt1.JointAxis = [0 0 1];
body1.Joint = jnt1;
addBody(robot, body1, 'base');

body2 = rigidBody('link2');
jnt2  = rigidBodyJoint('j2','revolute');
setFixedTransform(jnt2, trvec2tform([l1 0 0]));
jnt2.JointAxis = [0 0 1];
body2.Joint = jnt2;
addBody(robot, body2, 'link1');

tool = rigidBody('tool');
setFixedTransform(tool.Joint, trvec2tform([l2 0 0]));
addBody(robot, tool, 'link2');

q = [q1, q2]; % row vector
[J6, ~] = geometricJacobian(robot, q, 'tool'); % 6 x 2 Jacobian
      

In Simulink, the same mapping can be placed inside a MATLAB Function block. The block inputs are \( q_1, q_2, \dot{q}_1, \dot{q}_2, \ddot{q}_1, \ddot{q}_2 \), and its outputs are \( \dot{\mathbf{x}}, \ddot{\mathbf{x}} \). The function body simply calls planar2R_kin. This allows joint-space trajectory generators to feed end-effector velocity and acceleration to downstream modeling blocks.

Mathematica is convenient for symbolic derivation of Jacobians and accelerations. The snippet below demonstrates both instantaneous Jacobians and time-differentiated mappings.

(* Symbolic Jacobian for planar 2R *)
Clear[l1, l2, q1, q2, x, y];
x[q1_, q2_] := l1 Cos[q1] + l2 Cos[q1 + q2];
y[q1_, q2_] := l1 Sin[q1] + l2 Sin[q1 + q2];

p[q1_, q2_] := {x[q1, q2], y[q1, q2]};
J[q1_, q2_] := D[p[q1, q2], {{q1, q2}}]; (* 2 x 2 Jacobian *)

(* Time-parametrized variables *)
Clear[t];
q1t[t_] := q1[t];
q2t[t_] := q2[t];

pt[t_] := p[q1t[t], q2t[t]];
vt[t_] := D[pt[t], t];
at[t_] := D[pt[t], {t, 2}];

(* Substitute specific trajectories, e.g. sinusoidal joints *)
q1[t_] := q10 + a1 Sin(w1 t);
q2[t_] := q20 + a2 Sin(w2 t);

Simplify[vt[t]]
Simplify[at[t]]
      

The functions vt[t] and at[t] give symbolic expressions for the end-effector velocity and acceleration in terms of \( q_1(t), q_2(t) \) and their derivatives. By using Grad and JacobianMatrix, similar procedures generalize to spatial robots defined via homogeneous transforms.

In this lesson we formalized how joint-space trajectories induce task-space velocities and accelerations through the Jacobian and its time derivative. Starting from the forward kinematics \( \mathbf{x} = f(\mathbf{q}) \), we derived \( \dot{\mathbf{x}} = \mathbf{J}(\mathbf{q}) \dot{\mathbf{q}} \) and \( \ddot{\mathbf{x}} = \mathbf{J}(\mathbf{q}) \ddot{\mathbf{q}} + \dot{\mathbf{J}}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \), and related the latter to kinematic Hessians. The planar 2R example made the structure explicit, revealing the geometric coupling in the convective term \( \dot{\mathbf{J}}\dot{\mathbf{q}} \).

We then implemented these mappings in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, linking the abstract expressions to concrete numerical computation. These mappings are essential building blocks for trajectory feasibility checks, motion analysis, and, in later chapters, consistent coupling between kinematics and dynamics.