Let \( \mathbf{q}(t) \in \mathbb{R}^n \) denote the joint configuration of an \( n \)-DOF manipulator, assumed sufficiently smooth (at least three times differentiable) on a time interval \([0,T]\). We define

\[ \dot{\mathbf{q}}(t) = \frac{d\mathbf{q}(t)}{dt}, \quad \ddot{\mathbf{q}}(t) = \frac{d^2\mathbf{q}(t)}{dt^2}, \quad \mathbf{q}^{(3)}(t) = \frac{d^3\mathbf{q}(t)}{dt^3}, \]

where \( \mathbf{q}^{(3)}(t) \) is called the joint jerk vector. In task space, for an end-effector Cartesian trajectory \( \mathbf{x}(t) \in \mathbb{R}^m \), we similarly define velocity, acceleration, and jerk. High jerk manifests as rapid changes in acceleration and is associated with mechanical stress, excitation of flexible modes, and poor human comfort in physical human–robot interaction.

Jerk provides a natural building block for smoothness measures that assign a scalar cost to a trajectory. A common family of measures penalizes the squared Euclidean norm of jerk, integrated over time. For joint trajectories:

\[ J_{\text{jerk}}(\mathbf{q}) \;=\; \int_{0}^{T} \left\| \mathbf{q}^{(3)}(t) \right\|_2^2 \, dt, \]

where small values of \( J_{\text{jerk}} \) correspond to smooth trajectories. In this lesson we treat \( J_{\text{jerk}} \) and related functionals purely as quantitative descriptors of trajectories, not yet as optimization objectives for trajectory planning.

flowchart TD
  A["Joint trajectory q(t)"] --> B["Compute qdot(t), qddot(t)"]
  B --> C["Compute jerk q'''(t)"]
  C --> D["Form scalar smoothness functional"]
  D --> E["Evaluate smoothness index"]
  E --> F["Compare candidate trajectories"]
        

Consider a differentiable forward kinematics mapping \( \mathbf{x} = f(\mathbf{q}) \) from joint space to task space. Its first differential defines the manipulator Jacobian \( \mathbf{J}(\mathbf{q}) \in \mathbb{R}^{m \times n} \):

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

Differentiating once more yields the well-known acceleration relationship

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

where \( \dot{\mathbf{J}} \) is the time derivative of the Jacobian along the trajectory. Differentiating a third time provides the task-space jerk \( \mathbf{x}^{(3)}(t) \):

\[ \mathbf{x}^{(3)}(t) = \frac{d}{dt}\ddot{\mathbf{x}}(t) = \mathbf{J}\,\mathbf{q}^{(3)} + 2\,\dot{\mathbf{J}}\,\ddot{\mathbf{q}} + \ddot{\mathbf{J}}\,\dot{\mathbf{q}}, \]

where, for brevity, we suppressed the explicit time dependence and wrote \( \mathbf{J} \equiv \mathbf{J}(\mathbf{q}(t)) \), \( \dot{\mathbf{J}} \equiv \dot{\mathbf{J}}(\mathbf{q}(t),\dot{\mathbf{q}}(t)) \), and \( \ddot{\mathbf{J}} \equiv \ddot{\mathbf{J}}(\mathbf{q}(t),\dot{\mathbf{q}}(t),\ddot{\mathbf{q}}(t)) \).

The derivation follows from repeated application of the product rule and the chain rule. Let us illustrate the first step. Starting from \( \ddot{\mathbf{x}} = \mathbf{J}\,\ddot{\mathbf{q}} + \dot{\mathbf{J}}\,\dot{\mathbf{q}} \), we have

\[ \mathbf{x}^{(3)} = \frac{d}{dt}\left(\mathbf{J}\,\ddot{\mathbf{q}}\right) + \frac{d}{dt}\left(\dot{\mathbf{J}}\,\dot{\mathbf{q}}\right) = \dot{\mathbf{J}}\,\ddot{\mathbf{q}} + \mathbf{J}\,\mathbf{q}^{(3)} + \ddot{\mathbf{J}}\,\dot{\mathbf{q}} + \dot{\mathbf{J}}\,\ddot{\mathbf{q}}, \]

which, by grouping terms in \( \dot{\mathbf{J}}\,\ddot{\mathbf{q}} \), yields the expression above. This explicit representation shows that task-space jerk depends on joint jerk, but also on lower-order derivatives through the time-varying Jacobian.

We now formalize jerk-based smoothness measures. Let \( \mathbf{q}(t) \) be an admissible trajectory with prescribed boundary conditions (e.g., specified position and velocity at \( t=0 \) and \( t=T \)). A general quadratic jerk functional in joint space is

\[ J_{\text{jerk}}(\mathbf{q}) = \int_{0}^{T} \mathbf{q}^{(3)}(t)^{\top} \mathbf{W}\,\mathbf{q}^{(3)}(t)\,dt, \]

where \( \mathbf{W} \in \mathbb{R}^{n \times n} \) is symmetric positive definite and can be used to weigh different joints differently (e.g., according to gear ratios or human comfort).

Convexity. To see that \( J_{\text{jerk}} \) is convex in the trajectory, consider two admissible trajectories \( \mathbf{q}_1(t) \) and \( \mathbf{q}_2(t) \) and a scalar \( \lambda \in [0,1] \). Define \( \mathbf{q}_\lambda(t) = \lambda \mathbf{q}_1(t) + (1-\lambda)\mathbf{q}_2(t) \). Then

\[ \mathbf{q}_\lambda^{(3)}(t) = \lambda \mathbf{q}_1^{(3)}(t) + (1-\lambda)\mathbf{q}_2^{(3)}(t), \]

and by positive definiteness of \( \mathbf{W} \) and the fact that \( \mathbf{y} \mapsto \mathbf{y}^{\top}\mathbf{W}\mathbf{y} \) is convex, we obtain

\[ J_{\text{jerk}}(\mathbf{q}_\lambda) \le \lambda J_{\text{jerk}}(\mathbf{q}_1) + (1-\lambda) J_{\text{jerk}}(\mathbf{q}_2), \]

so jerk-based smoothness measures are convex in the trajectory space (under linear boundary constraints). This property is important when such functionals are later used in optimization.

Task-space smoothness. Analogously, one may define a task-space functional

\[ J_{\text{jerk}}^{\text{task}}(\mathbf{q}) = \int_{0}^{T} \mathbf{x}^{(3)}(t)^{\top} \mathbf{W}_x\,\mathbf{x}^{(3)}(t)\,dt, \]

where \( \mathbf{x}^{(3)}(t) \) is expressed in terms of joint quantities via the previous section. Note that joint-space and task-space measures are generally not equivalent because the mapping between them is configuration dependent.

Alternative smoothness measures. Besides integrated squared jerk, other commonly used functionals include

• Integrated squared acceleration:

  \[ J_{\text{acc}}(\mathbf{q}) = \int_{0}^{T} \ddot{\mathbf{q}}(t)^{\top} \mathbf{W}_a\,\ddot{\mathbf{q}}(t)\,dt, \]

• Total variation of acceleration:

  \[ J_{\text{TVa}}(\mathbf{q}) = \int_{0}^{T} \left\| \frac{d}{dt}\ddot{\mathbf{q}}(t) \right\|_1 \, dt = \int_{0}^{T} \left\| \mathbf{q}^{(3)}(t) \right\|_1 \, dt, \]

• Infinity-norm jerk measure:

  \[ J_{\infty}(\mathbf{q}) = \max_{t \in [0,T]} \left\| \mathbf{q}^{(3)}(t) \right\|_{\infty}, \]

Quadratic measures (\( J_{\text{jerk}} \), \( J_{\text{acc}} \)) are differentiable and convenient for analytic treatment, while non-quadratic measures (e.g. \( J_{\infty} \)) capture peak values that may correlate better with actuator limits.

flowchart TD
  Q["Trajectory q(t)"] --> J1["Compute jerk q'''(t)"]
  J1 --> S1["Quadratic measure: integral of ||q'''(t)||^2"]
  J1 --> S2["Nonquadratic measure: max or L1 norm"]
  S1 --> R1["Smoothness index (convex)"]
  S2 --> R2["Smoothness index (may be nonsmooth)"]
        

Real robots are commanded at discrete sampling times \( t_k = k\Delta t \) for \( k=0,\dots,N \), where \( \Delta t > 0 \) is the sampling period. Let \( \mathbf{q}_k \approx \mathbf{q}(t_k) \) denote the sampled joint positions. A standard third-order finite difference approximation of jerk is

\[ \mathbf{j}_k \;\approx\; \frac{\mathbf{q}_{k+2} - 3\mathbf{q}_{k+1} + 3\mathbf{q}_{k} - \mathbf{q}_{k-1}} {\Delta t^3}, \quad k = 1,\dots,N-2. \]

Stack the samples into a vector \( \mathbf{q}_{\text{stack}} \in \mathbb{R}^{(N+1)n} \) and define a block Toeplitz matrix \( \mathbf{D}_3 \) implementing the above difference operator. Then the discrete jerk vector \( \mathbf{j}_{\text{stack}} = \mathbf{D}_3 \mathbf{q}_{\text{stack}} \) and a discrete approximation of the continuous functional can be written as

\[ J_{\text{jerk}}^{\text{disc}}(\mathbf{q}_{\text{stack}}) = \Delta t \sum_{k=1}^{N-2} \mathbf{j}_k^{\top}\mathbf{W}\,\mathbf{j}_k = \mathbf{q}_{\text{stack}}^{\top} \underbrace{\left( \Delta t \, \mathbf{D}_3^{\top} (\mathbf{I} \otimes \mathbf{W}) \mathbf{D}_3 \right)}_{\mathbf{H} \succeq 0} \mathbf{q}_{\text{stack}}. \]

Here \( \otimes \) denotes the Kronecker product. We see that \( J_{\text{jerk}}^{\text{disc}} \) is a quadratic form with symmetric positive semidefinite matrix \( \mathbf{H} \), since it can be written as \( J_{\text{jerk}}^{\text{disc}} = \left\| \mathbf{A}\mathbf{q}_{\text{stack}} \right\|_2^2 \) for \( \mathbf{A} = (\Delta t^{1/2})(\mathbf{I} \otimes \mathbf{W}^{1/2})\mathbf{D}_3 \).

This matrix representation is useful when jerk-based smoothness is combined with other quadratic quantities (e.g., in least-squares estimation or as regularization terms) but here we simply regard \( J_{\text{jerk}}^{\text{disc}} \) as a computable index.

The following Python example computes discrete joint jerk and a scalar jerk-based smoothness index for a given sampled trajectory. We use numpy for array manipulation. In a larger robotics project, this can be integrated with toolboxes such as roboticstoolbox (to generate trajectories) or modern_robotics.

import numpy as np

def finite_difference_jerk(q_samples, dt):
    """
    Compute discrete joint jerk from sampled joint positions.
    Parameters
    ----------
    q_samples : ndarray, shape (N+1, n_joints)
        Sampled joint positions at times t_k = k*dt, k = 0,...,N.
    dt : float
        Sampling period.

    Returns
    -------
    jerk : ndarray, shape (N-2, n_joints)
        Approximate jerk at times t_1,...,t_{N-2}.
    """
    q_samples = np.asarray(q_samples)
    if q_samples.ndim == 1:
        q_samples = q_samples[:, None]

    N_plus_1, n_joints = q_samples.shape
    if N_plus_1 < 4:
        raise ValueError("Need at least 4 samples to compute third derivative")

    N = N_plus_1 - 1
    jerk = np.zeros((N - 2, n_joints))
    dt3 = dt ** 3

    for k in range(1, N - 1):
        jerk[k - 1, :] = (
            q_samples[k + 2, :]
            - 3.0 * q_samples[k + 1, :]
            + 3.0 * q_samples[k, :]
            - q_samples[k - 1, :]
        ) / dt3

    return jerk


def jerk_smoothness_index(q_samples, dt, W=None):
    """
    Compute integrated squared jerk index via discrete approximation.
    J = dt * sum_k j_k^T W j_k.
    """
    q_samples = np.asarray(q_samples)
    jerk = finite_difference_jerk(q_samples, dt)

    n_joints = jerk.shape[1]
    if W is None:
        W = np.eye(n_joints)
    W = np.asarray(W)
    if W.shape != (n_joints, n_joints):
        raise ValueError("Weight matrix W must have shape (n_joints, n_joints)")

    # compute sum over time
    J = 0.0
    for j_k in jerk:
        J += float(j_k.T @ W @ j_k)
    J *= dt
    return J


if __name__ == "__main__":
    # Example: quintic interpolation between two joint configurations
    # q(t) = a0 + a1 t + ... + a5 t^5 chosen to satisfy boundary conditions.
    # Here we simply construct a smooth scalar trajectory using numpy.
    T = 2.0
    dt = 0.01
    t = np.arange(0.0, T + dt, dt)
    N_plus_1 = t.size

    # Single joint example: smooth "S"-shaped motion from 0 to 1 rad
    s = t / T
    q = 10 * s**3 - 15 * s**4 + 6 * s**5  # C^2 quintic polynomial
    q_samples = q[:, None]  # shape (N+1, 1)

    j = finite_difference_jerk(q_samples, dt)
    J_index = jerk_smoothness_index(q_samples, dt)

    print("Discrete jerk shape:", j.shape)
    print("Jerk-based smoothness index J =", J_index)
      

The polynomial used above, \( q(t) = 10 (t/T)^3 - 15 (t/T)^4 + 6 (t/T)^5 \), is a standard smooth trajectory with vanishing velocity and acceleration at the boundaries and relatively small jerk throughout the interval.

In C++, the Eigen library is widely used for linear algebra in robotics. The following function computes discrete jerk and a scalar smoothness index from a matrix of samples q whose rows correspond to time steps and columns to joints.

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

Eigen::MatrixXd finiteDifferenceJerk(const Eigen::MatrixXd& q_samples, double dt) {
    // q_samples: (N+1) x n_joints
    const std::ptrdiff_t N_plus_1 = q_samples.rows();
    const std::ptrdiff_t n_joints = q_samples.cols();

    if (N_plus_1 < 4) {
        throw std::runtime_error("Need at least 4 samples to compute third derivative");
    }

    const std::ptrdiff_t N = N_plus_1 - 1;
    Eigen::MatrixXd jerk(N - 2, n_joints);
    double dt3 = dt * dt * dt;

    for (std::ptrdiff_t k = 1; k < N - 1; ++k) {
        jerk.row(k - 1) =
            (q_samples.row(k + 2)
            - 3.0 * q_samples.row(k + 1)
            + 3.0 * q_samples.row(k)
            - q_samples.row(k - 1)) / dt3;
    }
    return jerk;
}

double jerkSmoothnessIndex(const Eigen::MatrixXd& q_samples, double dt,
                           const Eigen::MatrixXd* W_ptr = nullptr) {
    Eigen::MatrixXd jerk = finiteDifferenceJerk(q_samples, dt);
    const std::ptrdiff_t n_joints = jerk.cols();

    Eigen::MatrixXd W;
    if (W_ptr == nullptr) {
        W = Eigen::MatrixXd::Identity(n_joints, n_joints);
    } else {
        W = *W_ptr;
        if (W.rows() != n_joints || W.cols() != n_joints) {
            throw std::runtime_error("Weight matrix W has wrong dimensions");
        }
    }

    double J = 0.0;
    for (std::ptrdiff_t k = 0; k < jerk.rows(); ++k) {
        Eigen::VectorXd j_k = jerk.row(k).transpose();
        J += j_k.transpose() * W * j_k;
    }
    J *= dt;
    return J;
}

int main() {
    // Example: single joint, smooth scalar trajectory, as in the Python example.
    double T = 2.0;
    double dt = 0.01;
    int N_plus_1 = static_cast<int>(T / dt) + 1;
    Eigen::MatrixXd q_samples(N_plus_1, 1);

    for (int k = 0; k < N_plus_1; ++k) {
        double t = k * dt;
        double s = t / T;
        double q = 10.0 * s * s * s
                 - 15.0 * s * s * s * s
                 + 6.0 * s * s * s * s * s;
        q_samples(k, 0) = q;
    }

    Eigen::MatrixXd jerk = finiteDifferenceJerk(q_samples, dt);
    double J = jerkSmoothnessIndex(q_samples, dt);

    std::cout << "Jerk samples: " << jerk.rows() << " x "
              << jerk.cols() << std::endl;
    std::cout << "Smoothness index J = " << J << std::endl;
    return 0;
}
      

This implementation is easily embedded into larger kinematics or dynamics code, for example after generating a sequence of joint configurations via a time-parameterized path.

Robotics in Java often represents trajectories as raw arrays. The following pure Java code assumes a 2D array qSamples of shape [Nplus1][nJoints].

public final class JerkUtils {

    public static double[][] finiteDifferenceJerk(double[][] qSamples, double dt) {
        int Nplus1 = qSamples.length;
        if (Nplus1 < 4) {
            throw new IllegalArgumentException("Need at least 4 samples");
        }
        int nJoints = qSamples[0].length;
        int N = Nplus1 - 1;
        double[][] jerk = new double[N - 2][nJoints];

        double dt3 = dt * dt * dt;

        for (int k = 1; k < N - 1; ++k) {
            for (int j = 0; j < nJoints; ++j) {
                double qkp2 = qSamples[k + 2][j];
                double qkp1 = qSamples[k + 1][j];
                double qk   = qSamples[k][j];
                double qkm1 = qSamples[k - 1][j];
                jerk[k - 1][j] = (qkp2 - 3.0 * qkp1 + 3.0 * qk - qkm1) / dt3;
            }
        }
        return jerk;
    }

    public static double jerkSmoothnessIndex(double[][] qSamples, double dt, double[][] W) {
        double[][] jerk = finiteDifferenceJerk(qSamples, dt);
        int Tsteps = jerk.length;
        int nJoints = jerk[0].length;

        double[][] Wlocal;
        if (W == null) {
            Wlocal = new double[nJoints][nJoints];
            for (int i = 0; i < nJoints; ++i) {
                Wlocal[i][i] = 1.0;
            }
        } else {
            Wlocal = W;
        }

        double J = 0.0;
        for (int k = 0; k < Tsteps; ++k) {
            for (int i = 0; i < nJoints; ++i) {
                for (int j = 0; j < nJoints; ++j) {
                    J += jerk[k][i] * Wlocal[i][j] * jerk[k][j];
                }
            }
        }
        J *= dt;
        return J;
    }

    public static void main(String[] args) {
        double T = 2.0;
        double dt = 0.01;
        int Nplus1 = (int) Math.round(T / dt) + 1;
        double[][] qSamples = new double[Nplus1][1];

        for (int k = 0; k < Nplus1; ++k) {
            double t = k * dt;
            double s = t / T;
            double q = 10.0 * Math.pow(s, 3)
                     - 15.0 * Math.pow(s, 4)
                     + 6.0 * Math.pow(s, 5);
            qSamples[k][0] = q;
        }

        double J = jerkSmoothnessIndex(qSamples, dt, null);
        System.out.println("Jerk-based smoothness index J = " + J);
    }
}
      

In larger software systems (e.g. Java ROS clients or simulation environments), this utility can be invoked on joint trajectories produced by planners or interpolators to rank their smoothness.

MATLAB is a standard environment for robotics. For a sampled joint trajectory Q of size (N+1) × nJoints, the code below computes discrete jerk and an integrated squared jerk index.

function [jerk, J] = jerk_smoothness(Q, dt, W)
% Q: (N+1)-by-nJoints matrix of joint positions.
% dt: sampling period.
% W:  nJoints-by-nJoints weighting matrix (optional, defaults to eye).

[Nplus1, nJoints] = size(Q);
if Nplus1 < 4
    error('Need at least 4 samples');
end
N = Nplus1 - 1;
jerk = zeros(N-2, nJoints);

if nargin < 3
    W = eye(nJoints);
end

dt3 = dt^3;

for k = 2:(N-1)
    jerk(k-1,:) = ( Q(k+2,:) - 3*Q(k+1,:) + 3*Q(k,:) - Q(k-1,:) ) / dt3;
end

J = 0.0;
for k = 1:(N-2)
    j_k = jerk(k,:).';
    J = J + j_k.' * W * j_k;
end
J = J * dt;
end

% Example usage:
T = 2.0;
dt = 0.01;
t = 0:dt:T;
s = t / T;
q = 10*s.^3 - 15*s.^4 + 6*s.^5;   % 1-DOF smooth trajectory
Q = q(:);                          % column vector (N+1)-by-1

[jerk, J] = jerk_smoothness(Q, dt);
fprintf('Smoothness index J = %.6f\n', J);
      

Simulink realization. The same computation can be represented in Simulink as a block diagram:

• A signal line carrying joint positions \( q(t_k) \) feeds three cascaded Discrete Derivative blocks (approximating first, second, and third derivatives).
• The output of the third derivative block (jerk) is fed into a Gain block implementing the weight matrix \( \mathbf{W} \) and then to a Dot Product block to form \( \mathbf{j}_k^{\top}\mathbf{W}\,\mathbf{j}_k \).
• A Discrete-Time Integrator accumulates this quantity over time, yielding a running estimate of the jerk-based smoothness index.

This Simulink structure can be dropped into a larger model that already simulates manipulator kinematics or dynamics, allowing on-line monitoring of trajectory smoothness.

Mathematica is convenient for symbolic manipulations involving jerk as well as numerical evaluation of smoothness functionals. The example below illustrates both.

(* Symbolic definition of a quintic joint trajectory q(t) *)
Clear[q];
q[t_] := a0 + a1 t + a2 t^2 + a3 t^3 + a4 t^4 + a5 t^5;

(* Compute derivatives *)
qd[t_]  = D[q[t], t];
qdd[t_] = D[q[t], {t, 2}];
qddd[t_] = D[q[t], {t, 3}];

(* Jerk-based smoothness functional over [0, T] *)
J[T_] := Integrate[qddd[t]^2, {t, 0, T},
                  Assumptions -> {T > 0}];

(* For numerical evaluation, choose coefficients to match boundary conditions.
   Here we directly construct a specific quintic with zero vel/acc at endpoints. *)
Tval = 2.0;
s[t_] := (t/Tval);
qSmooth[t_] := 10 s[t]^3 - 15 s[t]^4 + 6 s[t]^5;
jSmooth[t_] := D[qSmooth[t], {t, 3}];

Jnum = NIntegrate[jSmooth[t]^2, {t, 0, Tval}];
Print["Smoothness index J = ", Jnum];
      

Symbolic expressions for \( J(T) \) can be insightful when comparing different families of trajectories (e.g., polynomials vs. trigonometric parametrizations) under the same boundary conditions.

In this lesson we introduced jerk as the third time derivative of joint or task-space coordinates and derived explicit formulas connecting task-space jerk to joint-space derivatives through the Jacobian and its time derivatives. We defined continuous and discrete smoothness functionals, with emphasis on integrated squared jerk, and discussed their convexity and basic properties.

We then provided implementations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica to compute discrete jerk and jerk-based smoothness indices along sampled trajectories. These constructions will serve as building blocks for subsequent lessons on time parameterization and numerical differentiation, where such measures can be used to evaluate or regularize trajectory generation algorithms.