Consider an \( n \)-link serial manipulator with generalized coordinates \( \mathbf{q} \in \mathbb{R}^n \) and joint velocities \( \dot{\mathbf{q}} \in \mathbb{R}^n \). Each link \( i \) carries a frame \(\{i\}\) whose pose with respect to a base frame \(\{0\}\) is \( T_i(\mathbf{q}) \in \mathrm{SE}(3) \). Differential kinematics asks:

• How do we compute the twist (angular and linear velocity) of each link from \( \dot{\mathbf{q}} \) and possibly the base velocity?
• How do contributions of different joints add up to yield end-effector velocity?
• How can this computation be organized recursively along the chain so that it is numerically efficient and structurally clear?

Let \( V_i \in \mathbb{R}^6 \) denote the twist of link \( i \), decomposed as \( V_i = \begin{bmatrix} \boldsymbol{\omega}_i \\ \mathbf{v}_i \end{bmatrix} \), where \( \boldsymbol{\omega}_i \in \mathbb{R}^3 \) is angular velocity and \( \mathbf{v}_i \in \mathbb{R}^3 \) is linear velocity of the frame origin. A central result of this lesson is the velocity propagation recursion

\[ V_i^{(i)} = \mathrm{Ad}_{T_{i-1,i}(\mathbf{q})^{-1}} \, V_{i-1}^{(i-1)} + S_i \, \dot{q}_i , \quad i = 1,\dots,n , \]

where \( S_i \in \mathbb{R}^6 \) is the joint screw of joint \( i \) expressed in frame \(\{i\}\), \( T_{i-1,i}(\mathbf{q}) \) is the relative transform from frame \(\{i\}\) to \(\{i-1\}\), and \( \mathrm{Ad}_{T} \) is the adjoint representation of \( T \in \mathrm{SE}(3) \). If the base is fixed, \( V_0^{(0)} = \mathbf{0} \).

This recursion is purely kinematic: it involves only positions, orientations, and velocities, and will later be reused in Newton–Euler dynamics algorithms. In the next lessons, stacking the columns \( \mathrm{Ad}_{T_{0}^{i-1}} S_i \) will yield the geometric Jacobian, but here we restrict ourselves to the propagation view.

Let \( T(t) \in \mathrm{SE}(3) \) describe the pose of a rigid body as a function of time. We write

\[ T(t) = \begin{bmatrix} R(t) & p(t) \\ 0 & 1 \end{bmatrix}, \quad R(t) \in \mathrm{SO}(3), \; p(t) \in \mathbb{R}^3 . \]

The spatial twist of the body (twist expressed in the spatial frame) is defined by

\[ \hat{V}_s(t) \triangleq \dot{T}(t) \, T(t)^{-1} = \begin{bmatrix} \hat{\boldsymbol{\omega}}(t) & \mathbf{v}(t) \\ 0 & 0 \end{bmatrix}, \]

where \( \boldsymbol{\omega}(t) \in \mathbb{R}^3 \) is angular velocity, \( \mathbf{v}(t) \in \mathbb{R}^3 \) a linear velocity, and \( \hat{\boldsymbol{\omega}} \) is the skew-symmetric matrix

\[ \hat{\boldsymbol{\omega}} = \begin{bmatrix} 0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0 \end{bmatrix}, \quad \boldsymbol{\omega} = \begin{bmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{bmatrix}. \]

The body twist (twist expressed in the body-attached frame) is similarly defined by

\[ \hat{V}_b(t) \triangleq T(t)^{-1} \dot{T}(t). \]

Spatial and body twists are related by the adjoint operator:

\[ V_s = \mathrm{Ad}_{T} \, V_b, \quad \mathrm{Ad}_{T} = \begin{bmatrix} R & 0 \\ \hat{p} R & R \end{bmatrix}, \quad T = \begin{bmatrix} R & p \\ 0 & 1 \end{bmatrix}. \]

A key identity, used later, is

\[ T \, \hat{S} \, T^{-1} = \widehat{\mathrm{Ad}_{T} S} \]

for any twist \( S \in \mathbb{R}^6 \). This follows from the definition of \( \mathrm{Ad}_{T} \) and the way homogeneous transforms act on points; it encodes the change of coordinates of twists between frames.

Finally, the instantaneous velocity of a point rigidly attached to the body at spatial position \( x(t) = R(t) \, r + p(t) \), where \( r \) is the fixed body-frame coordinate, can be shown to satisfy

\[ \dot{x}(t) = \boldsymbol{\omega}(t) \times x(t) + \mathbf{v}(t) , \]

so the twist \( V_s = \begin{bmatrix} \boldsymbol{\omega} \\ \mathbf{v} \end{bmatrix} \) compactly encodes all velocities of points on the body.

For a serial manipulator, each joint contributes an instantaneous motion described by a joint screw \( S_i \in \mathbb{R}^6 \). Recall from the Product-of-Exponentials (PoE) formulation that for a single joint, the relative motion is

\[ T(\,q_i\,) = \exp\!\big(\hat{S}_i \, q_i \big), \]

where \( q_i \) is the scalar joint variable (angle or displacement). The screw \( S_i \) has the structure

• Revolute joint: axis direction \( \mathbf{\omega}_i \) (unit vector) and a point \( q^{(p)}_i \) on the axis give

  \[ S_i = \begin{bmatrix} \mathbf{\omega}_i \\ -\mathbf{\omega}_i \times q^{(p)}_i \end{bmatrix}, \quad \|\mathbf{\omega}_i\| = 1 . \]

• Prismatic joint: axis direction \( \mathbf{v}_i \) (unit vector), pure translation:

  \[ S_i = \begin{bmatrix} \mathbf{0} \\ \mathbf{v}_i \end{bmatrix}, \quad \|\mathbf{v}_i\| = 1. \]

For small \( \Delta q_i \), the relative transform satisfies

\[ T(q_i + \Delta q_i) \approx \big( I + \hat{S}_i \, \Delta q_i \big) \, T(q_i), \]

hence the instantaneous joint twist contribution is \( S_i \dot{q}_i \). The central question of velocity propagation is: how do these local contributions combine along a chain when expressed in different link frames?

Let \( T_i(\mathbf{q}) \) denote the pose of frame \(\{i\}\) with respect to \(\{0\}\), and denote the relative transform from \(\{i\}\) to \(\{i-1\}\) by

\[ T_{i-1,i}(q_i) = \begin{bmatrix} R_{i-1,i}(q_i) & p_{i-1,i}(q_i) \\ 0 & 1 \end{bmatrix}. \]

Then \( T_i = T_{i-1} \, T_{i-1,i}(q_i) \). Differentiating with respect to time,

\[ \dot{T}_i = \dot{T}_{i-1} \, T_{i-1,i} + T_{i-1} \, \dot{T}_{i-1,i}. \]

Multiply on the right by \( T_i^{-1} = T_{i-1,i}^{-1} T_{i-1}^{-1} \):

\[ \dot{T}_i T_i^{-1} = \dot{T}_{i-1} T_{i-1}^{-1} + T_{i-1} \dot{T}_{i-1,i} T_{i-1,i}^{-1} T_{i-1}^{-1}. \]

By definition of spatial twists, \( \hat{V}_i^{(s)} = \dot{T}_i T_i^{-1} \) and \( \hat{V}_{i-1}^{(s)} = \dot{T}_{i-1} T_{i-1}^{-1} \), so

\[ \hat{V}_i^{(s)} = \hat{V}_{i-1}^{(s)} + T_{i-1} \dot{T}_{i-1,i} T_{i-1,i}^{-1} T_{i-1}^{-1}. \]

For a single joint, \( T_{i-1,i}(q_i) = \exp(\hat{S}_i q_i) \), hence

\[ \dot{T}_{i-1,i} = \hat{S}_i \dot{q}_i \exp(\hat{S}_i q_i) = \hat{S}_i \dot{q}_i T_{i-1,i}. \]

Substituting,

\[ T_{i-1} \dot{T}_{i-1,i} T_{i-1,i}^{-1} T_{i-1}^{-1} = T_{i-1} \hat{S}_i \dot{q}_i T_{i-1}^{-1} = \widehat{\mathrm{Ad}_{T_{i-1}} S_i} \, \dot{q}_i , \]

where we used \( T \hat{S} T^{-1} = \widehat{\mathrm{Ad}_{T} S} \). Thus,

\[ \hat{V}_i^{(s)} = \hat{V}_{i-1}^{(s)} + \widehat{\mathrm{Ad}_{T_{i-1}} S_i} \, \dot{q}_i . \]

Removing the hat operator yields the spatial velocity recursion

\[ V_i^{(s)} = V_{i-1}^{(s)} + \mathrm{Ad}_{T_{i-1}} S_i \, \dot{q}_i . \]

Using instead twists expressed in the local link frames, one obtains

\[ V_i^{(i)} = \mathrm{Ad}_{T_{i-1,i}(q_i)^{-1}} \, V_{i-1}^{(i-1)} + S_i \, \dot{q}_i. \]

If the base is fixed, \( V_0^{(0)} = \mathbf{0} \), and we can unroll the recursion to see that \( V_n \) is a linear function of \( \dot{\mathbf{q}} \). In the next lesson this linear map will be identified explicitly as the manipulator Jacobian.

flowchart TD
  A["Given q, qdot, base velocity V0"] --> B["i = 1"]
  B --> C["Compute relative transform T_{i-1,i}(q_i)"]
  C --> D["Transform previous twist: V_prev = Ad_{T_{i-1,i}^{-1}}(V_{i-1})"]
  D --> E["Add joint contribution: V_i = V_prev + S_i * qdot_i"]
  E --> F{"i < n ?"}
  F -->|yes| G["i = i + 1"] --> C
  F -->|no| H["Collect all V_i and end-effector twist"]
        

To illustrate velocity propagation concretely, consider a planar 2R arm in the \( x\text{-}y \) plane with link lengths \( L_1, L_2 \) and joint angles \( q_1, q_2 \). The end-effector position is

\[ p(q_1, q_2) = \begin{bmatrix} L_1 \cos q_1 + L_2 \cos(q_1+q_2) \\ L_1 \sin q_1 + L_2 \sin(q_1+q_2) \\ 0 \end{bmatrix}. \]

Differentiating with respect to time gives the end-effector linear velocity components

\[ \begin{aligned} \dot{p}_x &= -L_1 \sin q_1 \, \dot{q}_1 - L_2 \sin(q_1+q_2) \, (\dot{q}_1 + \dot{q}_2), \\ \dot{p}_y &= L_1 \cos q_1 \, \dot{q}_1 + L_2 \cos(q_1+q_2) \, (\dot{q}_1 + \dot{q}_2). \end{aligned} \]

The angular velocity of the end-effector about the out-of-plane \( z \)-axis is

\[ \omega_z = \dot{q}_1 + \dot{q}_2. \]

We can derive the same result using the twist recursion. If both joints rotate about the same \( z \)-axis and the base is fixed, the joint screws in spatial coordinates are

\[ S_1 = \begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \quad S_2 = \begin{bmatrix} 0 \\ 0 \\ 1 \\ -L_1 \sin q_1 \\ L_1 \cos q_1 \\ 0 \end{bmatrix}, \]

evaluated at the current configuration. With \( V_0^{(s)} = 0 \),

\[ V_1^{(s)} = S_1 \dot{q}_1, \quad V_2^{(s)} = V_1^{(s)} + \mathrm{Ad}_{T_1(q_1)} S_2 \dot{q}_2, \]

and the translational part of \( V_2^{(s)} \) coincides with \( \dot{p} \) above. This illustrates how joint contributions accumulate through the adjoint transforms rather than by re-deriving trigonometric expressions from scratch at each step.

We now implement velocity propagation using NumPy, assuming that rigid-body motions are represented by \( 4 \times 4 \) homogeneous matrices and screw axes by \( 6 \times 1 \) vectors. We work with twists expressed in local link frames and use the recursion \( V_i^{(i)} = \mathrm{Ad}_{T_{i-1,i}^{-1}} V_{i-1}^{(i-1)} + S_i \dot{q}_i \).

import numpy as np

def skew(w):
    """
    w: shape (3,), returns 3x3 skew matrix.
    """
    wx, wy, wz = w
    return np.array([
        [0.0, -wz,  wy],
        [wz,  0.0, -wx],
        [-wy, wx,  0.0]
    ])

def adjoint(T):
    """
    T: 4x4 homogeneous transform.
    returns 6x6 adjoint matrix Ad_T.
    """
    R = T[0:3, 0:3]
    p = T[0:3, 3]
    p_hat = skew(p)
    Ad = np.zeros((6, 6))
    Ad[0:3, 0:3] = R
    Ad[3:6, 0:3] = p_hat @ R
    Ad[3:6, 3:6] = R
    return Ad

def forward_velocity_chain(T_list, S_list, qdot, V0=None):
    """
    Velocity propagation for a serial chain.

    Parameters
    ----------
    T_list : list of 4x4 numpy arrays
        T_list[i] is T_{i-1,i}(q_i), the transform from frame i to i-1.
    S_list : list of 6x1 numpy arrays
        Joint screw axes S_i expressed in frame i.
    qdot : numpy array of shape (n,)
        Joint velocities qdot_i.
    V0 : 6x1 numpy array or None
        Twist of the base frame expressed in frame 0. If None, assumes zero.

    Returns
    -------
    V_links : list of 6x1 numpy arrays
        V_links[i] is the twist of link i+1 expressed in frame i+1.
    """
    n = len(S_list)
    assert len(T_list) == n
    assert qdot.shape[0] == n

    if V0 is None:
        V_prev = np.zeros((6, 1))
    else:
        V_prev = np.asarray(V0).reshape(6, 1)

    V_links = []
    for i in range(n):
        T_i = T_list[i]
        Ad_inv = adjoint(np.linalg.inv(T_i))
        V_prev_in_i = Ad_inv @ V_prev
        S_i = S_list[i].reshape(6, 1)
        V_i = V_prev_in_i + S_i * qdot[i]
        V_links.append(V_i)
        V_prev = V_i
    return V_links

# Example usage for planar 2R (both joints about z axis)
L1, L2 = 1.0, 0.8
q = np.array([0.3, 0.5])
qdot = np.array([0.4, 0.2])

# Define relative transforms T_{0,1}(q1), T_{1,2}(q2) numerically
def rot_z(theta):
    c = np.cos(theta)
    s = np.sin(theta)
    T = np.eye(4)
    T[0:3, 0:3] = np.array([
        [c, -s, 0.0],
        [s,  c, 0.0],
        [0.0, 0.0, 1.0]
    ])
    return T

T01 = rot_z(q[0])
T01[0, 3] = 0.0
T01[1, 3] = 0.0
T01[2, 3] = 0.0

T12 = rot_z(q[1])
T12[0, 3] = L1
T12[1, 3] = 0.0
T12[2, 3] = 0.0

# Joint screws expressed in local joint frames
S1 = np.array([0.0, 0.0, 1.0, 0.0, 0.0, 0.0])
S2 = np.array([0.0, 0.0, 1.0, 0.0, 0.0, 0.0])

V_links = forward_velocity_chain([T01, T12], [S1, S2], qdot)
V1, V2 = V_links

print("Link 1 twist:", V1.flatten())
print("End-effector twist:", V2.flatten())
      

In practice, the relative transforms \( T_{i-1,i}(q_i) \) and screw axes \( S_i \) would come from a prior kinematics description (PoE or DH-derived). This function exposes a clean separation between geometric modeling and differential propagation.

In C++, the Eigen library is a standard choice for linear algebra in robotics. Below is a minimal implementation of the same recursion using Eigen::Matrix types.

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

using Matrix4d = Eigen::Matrix4d;
using Matrix6d = Eigen::Matrix<double, 6, 6>;
using Vector6d = Eigen::Matrix<double, 6, 1>;

Eigen::Matrix3d skew(const Eigen::Vector3d& w) {
    Eigen::Matrix3d W;
    W << 0.0,    -w.z(),  w.y(),
          w.z(),  0.0,    -w.x(),
         -w.y(),  w.x(),   0.0;
    return W;
}

Matrix6d adjoint(const Matrix4d& T) {
    Matrix6d Ad = Matrix6d::Zero();
    Eigen::Matrix3d R = T.block<3,3>(0,0);
    Eigen::Vector3d p = T.block<3,1>(0,3);
    Eigen::Matrix3d p_hat = skew(p);

    Ad.block<3,3>(0,0) = R;
    Ad.block<3,3>(3,0) = p_hat * R;
    Ad.block<3,3>(3,3) = R;
    return Ad;
}

std::vector<Vector6d> forwardVelocityChain(
    const std::vector<Matrix4d>& T_list,
    const std::vector<Vector6d>& S_list,
    const Eigen::VectorXd& qdot,
    const Vector6d& V0 = Vector6d::Zero())
{
    const std::size_t n = S_list.size();
    std::vector<Vector6d> V_links;
    V_links.reserve(n);

    Vector6d V_prev = V0;

    for (std::size_t i = 0; i < n; ++i) {
        const Matrix4d& T_i = T_list[i];
        Matrix4d T_inv = T_i.inverse();
        Matrix6d Ad_inv = adjoint(T_inv);
        Vector6d V_prev_in_i = Ad_inv * V_prev;
        Vector6d V_i = V_prev_in_i + S_list[i] * qdot(static_cast<int>(i));
        V_links.push_back(V_i);
        V_prev = V_i;
    }
    return V_links;
}
      

In more elaborate C++ robotics toolkits, this recursion is embedded into articulated body algorithms and templated on scalar types (e.g., for automatic differentiation). The structure, however, is the same: adjoint transforms and local screw contributions.

Java has a variety of numerical libraries (e.g., EJML) that can be used for robotics. Here we illustrate the algorithm using primitive arrays to emphasize the recursion logic; in practice, wrapping these operations into matrix types is recommended.

public class VelocityPropagation {

    // Multiply 4x4 matrices: C = A * B
    public static double[][] mul4(double[][] A, double[][] B) {
        double[][] C = new double[4][4];
        for (int i = 0; i < 4; ++i) {
            for (int j = 0; j < 4; ++j) {
                double sum = 0.0;
                for (int k = 0; k < 4; ++k) {
                    sum += A[i][k] * B[k][j];
                }
                C[i][j] = sum;
            }
        }
        return C;
    }

    // Compute 3x3 skew matrix of w
    public static double[][] skew(double[] w) {
        double wx = w[0], wy = w[1], wz = w[2];
        return new double[][] {
            { 0.0, -wz,  wy },
            { wz,  0.0, -wx },
            { -wy, wx,  0.0 }
        };
    }

    // Compute 6x6 adjoint of T (4x4)
    public static double[][] adjoint(double[][] T) {
        double[][] R = new double[3][3];
        double[] p = new double[3];
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                R[i][j] = T[i][j];
            }
            p[i] = T[i][3];
        }
        double[][] pHat = skew(p);

        double[][] Ad = new double[6][6];
        // Top-left (R)
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                Ad[i][j] = R[i][j];
            }
        }
        // Bottom-left (pHat * R)
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                double sum = 0.0;
                for (int k = 0; k < 3; ++k) {
                    sum += pHat[i][k] * R[k][j];
                }
                Ad[3 + i][j] = sum;
            }
        }
        // Bottom-right (R)
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                Ad[3 + i][3 + j] = R[i][j];
            }
        }
        return Ad;
    }

    // Multiply 6x6 matrix by 6x1 vector
    public static double[] mul6(double[][] A, double[] x) {
        double[] y = new double[6];
        for (int i = 0; i < 6; ++i) {
            double sum = 0.0;
            for (int j = 0; j < 6; ++j) {
                sum += A[i][j] * x[j];
            }
            y[i] = sum;
        }
        return y;
    }

    // Add two 6x1 vectors
    public static double[] add6(double[] a, double[] b) {
        double[] c = new double[6];
        for (int i = 0; i < 6; ++i) {
            c[i] = a[i] + b[i];
        }
        return c;
    }

    // Scale 6x1 vector
    public static double[] scale6(double[] a, double s) {
        double[] c = new double[6];
        for (int i = 0; i < 6; ++i) {
            c[i] = s * a[i];
        }
        return c;
    }

    // Forward velocity propagation
    public static double[][] forwardVelocityChain(
            double[][][] Tlist, // [n][4][4]
            double[][] Slist,   // [n][6]
            double[] qdot,      // [n]
            double[] V0         // [6], may be null
    ) {
        int n = Slist.length;
        double[] Vprev = new double[6];
        if (V0 != null) {
            System.arraycopy(V0, 0, Vprev, 0, 6);
        }

        double[][] Vlinks = new double[n][6];

        for (int i = 0; i < n; ++i) {
            double[][] T_i = Tlist[i];
            // For brevity, assume T_i is already the inverse transform T_{i-1,i}^{-1}
            double[][] Ad_inv = adjoint(T_i);
            double[] VprevInI = mul6(Ad_inv, Vprev);
            double[] Si = Slist[i];
            double[] jointTwist = scale6(Si, qdot[i]);
            double[] Vi = add6(VprevInI, jointTwist);
            Vlinks[i] = Vi;
            Vprev = Vi;
        }
        return Vlinks;
    }
}
      

For serious applications, one should use a tested linear algebra library (such as EJML) and provide robust routines for matrix inversion and motion transforms. The core recursion, however, is unchanged.

MATLAB is widely used in control and robotics. Below is a vectorized implementation of velocity propagation. Homogeneous matrices are stored as a 3D array Tlist(:,:,i), and screw axes as columns of Slist.

function Vlinks = forwardVelocityChain(Tlist, Slist, qdot, V0)
% FORWARDVELOCITYCHAIN Propagate twists along a serial chain.
%   Tlist : 4x4xn array of relative transforms T_{i-1,i}(q_i)
%   Slist : 6xn matrix, column i is screw axis S_i in frame i
%   qdot  : nx1 vector of joint velocities
%   V0    : 6x1 base twist (optional, default zero)
%
%   Vlinks: 6xn matrix, column i is twist of link i in frame i

if nargin < 4
    Vprev = zeros(6,1);
else
    Vprev = V0(:);
end

n = size(Slist, 2);
Vlinks = zeros(6, n);

for i = 1:n
    Ti = Tlist(:,:,i);
    Ad_inv = adjoint(inv(Ti));
    Vprev_in_i = Ad_inv * Vprev;
    Si = Slist(:,i);
    Vi = Vprev_in_i + Si * qdot(i);
    Vlinks(:,i) = Vi;
    Vprev = Vi;
end
end

function Ad = adjoint(T)
R = T(1:3,1:3);
p = T(1:3,4);
p_hat = skew(p);
Ad = [R, zeros(3,3); p_hat*R, R];
end

function W = skew(w)
wx = w(1); wy = w(2); wz = w(3);
W = [   0, -wz,  wy;
      wz,   0, -wx;
     -wy,  wx,   0];
end
      

A basic Simulink implementation can treat joint velocities and transforms as signals flowing through a for-loop subsystem or a MATLAB Function block:

• Use a MATLAB Function block that encapsulates the code above.
• Feed joint velocities \( \dot{\mathbf{q}} \) and the current transforms \( T_{i-1,i}(q_i) \) (computed, e.g., from a kinematics block) into the function.
• Output per-link twists to other subsystems (e.g., for control law design or for computing contact point velocities).

Wolfram Mathematica is well suited for symbolic differential kinematics. Below we compute symbolic twists for a planar 2R arm and verify the mapping from joint velocities to end-effector velocity.

ClearAll["Global`*"];

(* Skew matrix for 3-vector *)
Skew[w_List] := {
  {0, -w[[3]],  w[[2]]},
  {w[[3]], 0, -w[[1]]},
  {-w[[2]], w[[1]], 0}
};

(* Adjoint of homogeneous transform *)
Adjoint[T_] := Module[{R, p, ph},
  R = T[[1 ;; 3, 1 ;; 3]];
  p = T[[1 ;; 3, 4]];
  ph = Skew[p];
  ArrayFlatten[{
    {R, 0 IdentityMatrix[3]},
    {ph.R, R}
  }]
];

(* Planar 2R kinematics *)
L1 = Symbol["L1"]; L2 = Symbol["L2"];
q1 = Symbol["q1"]; q2 = Symbol["q2"];

RotZ[theta_] := {
  {Cos[theta], -Sin[theta], 0},
  {Sin[theta],  Cos[theta], 0},
  {0,           0,          1}
};

T01 = ArrayFlatten[{
  {RotZ[q1], { {0}, {0}, {0} } },
  { {0, 0, 0, 1} }
}];

T12 = ArrayFlatten[{
  {RotZ[q2], { {L1}, {0}, {0} } },
  {{0, 0, 0, 1}}
}];

T02 = Simplify[T01.T12];

(* Joint screws in local frames *)
S1 = {0, 0, 1, 0, 0, 0};
S2 = {0, 0, 1, 0, 0, 0};

qdot1 = Symbol["qdot1"];
qdot2 = Symbol["qdot2"];

(* Spatial twists via adjoint recursion *)
Ad01 = Adjoint[T01];
V1s = S1*qdot1;
V2s = V1s + Ad01.S2*qdot2 // Simplify;

(* Extract translational velocity component of end-effector origin *)
p2 = T02[[1 ;; 3, 4]];
px = p2[[1]]; py = p2[[2]];

vx = D[px, q1]*qdot1 + D[px, q2]*qdot2 // Simplify;
vy = D[py, q1]*qdot1 + D[py, q2]*qdot2 // Simplify;

{vx, vy}
      

The symbolic expressions for vx and vy match those derived analytically in Section 5, confirming that the twist-based recursion is consistent with classical planar kinematics.

In this lesson we:

• Reviewed spatial and body twists as concise representations of rigid-body velocities in \( \mathrm{SE}(3) \), linked by the adjoint operator.
• Expressed revolute and prismatic joints as screw axes whose instantaneous contributions to motion are \( S_i \dot{q}_i \).
• Derived the velocity propagation recursion \( V_i = \mathrm{Ad}_{T_{i-1,i}^{-1}} V_{i-1} + S_i \dot{q}_i \) for link twists in local frames, and the equivalent spatial-frame recursion.
• Demonstrated the method on a planar 2R manipulator and implemented the recursion in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.
• Showed that the end-effector twist is linear in the joint velocities, preparing the ground for the formal introduction of Jacobians in the next lesson.