Recall from Lesson 1 that for an \(n\)-DOF serial manipulator in generalized coordinates \( \mathbf{q} = [q_1,\dots,q_n]^\top \), the Lagrange–Euler equations in joint space can be written as

\[ \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{M}(\mathbf{q}) \) is the symmetric positive definite inertia matrix, \( \mathbf{g}(\mathbf{q}) \) is the gravity vector, and \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \) collects Coriolis and centrifugal torques. The Christoffel symbols provide a coordinate-invariant way of constructing this term from \( \mathbf{M}(\mathbf{q}) \).

At the component level, the \(i\)-th joint equation can be written as

\[ \sum_{j=1}^n M_{ij}(\mathbf{q})\,\ddot{q}_j \;+\; \sum_{j=1}^n\sum_{k=1}^n c_{ijk}(\mathbf{q})\,\dot{q}_j \dot{q}_k \;+\; g_i(\mathbf{q}) \;=\; \tau_i, \]

where the coefficients \( c_{ijk}(\mathbf{q}) \) are the Christoffel symbols of the second kind associated with the kinetic energy metric defined by \( \mathbf{M}(\mathbf{q}) \). In this lesson we make the passage \( \mathbf{M}(\mathbf{q}) \Longrightarrow c_{ijk}(\mathbf{q}) \Longrightarrow \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) explicit and algorithmic.

flowchart TD
  T["Specify kinetic energy T(q, q_dot)"] --> M["Extract inertia matrix M(q)"]
  M --> G["Compute Christoffel symbols c_ijk(q)"]
  G --> C["Build Coriolis matrix C(q, q_dot)"]
  C --> E["Joint-space dynamics: M(q) q_ddot + C(q, q_dot) q_dot + g(q) = tau"]
        

For a manipulator, the kinetic energy can always be written in quadratic form

\[ T(\mathbf{q},\dot{\mathbf{q}}) = \tfrac{1}{2}\,\dot{\mathbf{q}}^\top \mathbf{M}(\mathbf{q}) \dot{\mathbf{q}} = \tfrac{1}{2}\,\sum_{j=1}^n\sum_{k=1}^n M_{jk}(\mathbf{q})\,\dot{q}_j \dot{q}_k . \]

Lagrange's equations for purely mechanical systems (neglecting dissipation) read

\[ \frac{\mathrm{d}}{\mathrm{d}t}\!\left( \frac{\partial T}{\partial \dot{q}_i} \right) \;-\; \frac{\partial T}{\partial q_i} \;+\; \frac{\partial U}{\partial q_i} \;=\; \tau_i, \qquad i = 1,\dots,n, \]

where \( U(\mathbf{q}) \) is the potential energy. The Christoffel symbols emerge when we expand the kinetic terms \( \frac{\mathrm{d}}{\mathrm{d}t}\big( \partial T / \partial \dot{q}_i \big) - \partial T / \partial q_i \).

2.1. Computing \( \partial T / \partial \dot{q}_i \) and its time derivative

First compute the partial derivative with respect to the generalized velocity:

\[ \frac{\partial T}{\partial \dot{q}_i} = \frac{1}{2}\sum_{j,k} M_{jk}(\mathbf{q})\, \frac{\partial (\dot{q}_j \dot{q}_k)}{\partial \dot{q}_i} = \sum_{j=1}^n M_{ij}(\mathbf{q})\,\dot{q}_j , \]

because \( \mathbf{M}(\mathbf{q}) \) is symmetric. Taking the total time derivative gives

\[ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial T}{\partial \dot{q}_i} \right) = \sum_{j=1}^n M_{ij}(\mathbf{q})\,\ddot{q}_j \;+\; \sum_{j=1}^n \sum_{k=1}^n \frac{\partial M_{ij}(\mathbf{q})}{\partial q_k}\,\dot{q}_k \dot{q}_j . \]

2.2. Computing \( \partial T / \partial q_i \)

The kinetic energy depends on \( \mathbf{q} \) only through the entries of \( \mathbf{M}(\mathbf{q}) \), so

\[ \frac{\partial T}{\partial q_i} = \frac{1}{2}\sum_{j=1}^n\sum_{k=1}^n \frac{\partial M_{jk}(\mathbf{q})}{\partial q_i}\, \dot{q}_j \dot{q}_k . \]

2.3. Identification of the Christoffel symbols

The difference of the two kinetic terms is

\[ \begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial T}{\partial \dot{q}_i} \right) - \frac{\partial T}{\partial q_i} &= \sum_{j=1}^n M_{ij}(\mathbf{q})\,\ddot{q}_j \\ &\quad+ \sum_{j=1}^n\sum_{k=1}^n \left[ \frac{\partial M_{ij}}{\partial q_k} - \tfrac{1}{2}\frac{\partial M_{jk}}{\partial q_i} \right] \dot{q}_j \dot{q}_k . \end{aligned} \]

By symmetrizing in the indices \( j \) and \( k \), we define the Christoffel symbols of the second kind:

\[ c_{ijk}(\mathbf{q}) \;:=\; \tfrac{1}{2}\left( \frac{\partial M_{ij}}{\partial q_k} + \frac{\partial M_{ik}}{\partial q_j} - \frac{\partial M_{jk}}{\partial q_i} \right), \qquad c_{ijk} = c_{ikj}. \]

Then the kinetic part of Lagrange's equations becomes

\[ \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\partial T}{\partial \dot{q}_i} \right) - \frac{\partial T}{\partial q_i} = \sum_{j=1}^n M_{ij}(\mathbf{q})\,\ddot{q}_j \;+\; \sum_{j=1}^n\sum_{k=1}^n c_{ijk}(\mathbf{q})\,\dot{q}_j \dot{q}_k . \]

The Christoffel symbols are three-index objects. In computations and control, it is more convenient to work with an \(n\times n\) matrix \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) such that \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \) reproduces the quadratic velocity term. The common robotics convention is

\[ C_{ij}(\mathbf{q},\dot{\mathbf{q}}) \;:=\; \sum_{k=1}^n c_{ijk}(\mathbf{q})\,\dot{q}_k, \qquad \big[\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\big]_{ij} = C_{ij}. \]

Then the \(i\)-th component of \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \) is

\[ \begin{aligned} \big[\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}\big]_i &= \sum_{j=1}^n C_{ij}\,\dot{q}_j = \sum_{j=1}^n\sum_{k=1}^n c_{ijk}(\mathbf{q})\,\dot{q}_k \dot{q}_j , \end{aligned} \]

which matches the quadratic term obtained from Lagrange's equations. Notice that the mapping from \( c_{ijk} \) to \( \mathbf{C} \) is not unique: different conventions for distributing the Christoffel symbols among the entries of \( \mathbf{C} \) lead to different matrices, but all such matrices produce the same vector \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \).

3.1. Skew-symmetry property and energy consistency

A key structural property used in stability and passivity analysis is that the matrix \( \dot{\mathbf{M}}(\mathbf{q}) - 2\,\mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) is skew-symmetric:

\[ \big(\dot{\mathbf{M}}(\mathbf{q}) - 2\,\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\big)^\top = -\big(\dot{\mathbf{M}}(\mathbf{q}) - 2\,\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\big) . \]

Equivalently, for any velocity vector \( \mathbf{v} \in \mathbb{R}^n \),

\[ \mathbf{v}^\top\big(\dot{\mathbf{M}} - 2\mathbf{C}\big)\mathbf{v} = 0, \qquad\text{in particular}\qquad \dot{\mathbf{q}}^\top \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} = \tfrac{1}{2}\,\dot{\mathbf{q}}^\top \dot{\mathbf{M}}(\mathbf{q}) \dot{\mathbf{q}}. \]

This identity implies that the Coriolis/centrifugal term does not inject or dissipate mechanical energy: it only redistributes kinetic energy between degrees of freedom. This is central in proofs of passivity for robot manipulators.

flowchart TD
  A["Inertia matrix M(q)"] --> B["Compute time derivative M_dot(q, q_dot)"]
  A --> C["Compute Christoffel symbols c_ijk(q)"]
  C --> D["Form Coriolis matrix C(q, q_dot)"]
  B --> E["Check M_dot - 2 C is skew-symmetric"]
  D --> E
  E --> F["Use structure for stability/passivity analysis"]
        

Consider a planar 2R manipulator with joint angles \( \mathbf{q} = [q_1, q_2]^\top \). After standard rigid-body computations (using link masses, lengths and inertias), its inertia matrix can be written in the compact form

\[ \mathbf{M}(\mathbf{q}) = \begin{bmatrix} a_1 + 2a_2\cos(q_2) & a_3 + a_2\cos(q_2) \\ a_3 + a_2\cos(q_2) & a_3 \end{bmatrix}, \]

where the constants \( a_1, a_2, a_3 \) depend on the physical parameters (masses, link lengths, inertias). Denote \( c_2 := \cos(q_2) \) and \( s_2 := \sin(q_2) \).

4.1. Derivatives of \( \mathbf{M}(\mathbf{q}) \)

Only \( q_2 \) appears in the inertia matrix, so derivatives with respect to \( q_1 \) vanish:

\[ \begin{aligned} &\frac{\partial M_{11}}{\partial q_1} = 0, &&\frac{\partial M_{11}}{\partial q_2} = -2a_2 s_2, \\ &\frac{\partial M_{12}}{\partial q_1} = 0, &&\frac{\partial M_{12}}{\partial q_2} = -a_2 s_2, \\ &\frac{\partial M_{22}}{\partial q_1} = 0, &&\frac{\partial M_{22}}{\partial q_2} = 0. \end{aligned} \]

4.2. Christoffel symbols and Coriolis matrix

Using the definition \( c_{ijk} = \tfrac{1}{2}\big(\partial M_{ij}/\partial q_k + \partial M_{ik}/\partial q_j - \partial M_{jk}/\partial q_i\big) \), and exploiting symmetry \( M_{12} = M_{21} \), the nonzero Christoffel symbols are

\[ \begin{aligned} c_{112} &= \tfrac{1}{2}\left( \frac{\partial M_{11}}{\partial q_2} + \frac{\partial M_{12}}{\partial q_1} - \frac{\partial M_{12}}{\partial q_1} \right) = \tfrac{1}{2}(-2a_2 s_2) = -a_2 s_2, \\ c_{121} &= c_{112} = -a_2 s_2, \\ c_{122} &= \tfrac{1}{2}\left( \frac{\partial M_{12}}{\partial q_2} + \frac{\partial M_{12}}{\partial q_2} - \frac{\partial M_{22}}{\partial q_1} \right) = -a_2 s_2, \\ c_{212} &= \tfrac{1}{2}\left( \frac{\partial M_{21}}{\partial q_2} + \frac{\partial M_{22}}{\partial q_1} - \frac{\partial M_{12}}{\partial q_2} \right) = 0, \\ c_{221} &= \tfrac{1}{2}\left( \frac{\partial M_{22}}{\partial q_1} + \frac{\partial M_{21}}{\partial q_2} - \frac{\partial M_{21}}{\partial q_2} \right) = 0. \end{aligned} \]

Let the joint velocities be \( \dot{\mathbf{q}} = [\dot{q}_1,\dot{q}_2]^\top \). With the convention \( C_{ij} = \sum_k c_{ijk}\dot{q}_k \), one standard choice of \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) that satisfies the skew-symmetry property is

\[ \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) = \begin{bmatrix} -a_2 s_2 \dot{q}_2 & -a_2 s_2 (\dot{q}_1 + \dot{q}_2) \\ a_2 s_2 \dot{q}_1 & 0 \end{bmatrix}. \]

Direct computation verifies that \( \dot{\mathbf{M}}(\mathbf{q}) - 2\mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \) is skew-symmetric and that \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \) reproduces the quadratic velocity terms in the 2R manipulator dynamics.

Python with sympy (symbolic) and numpy (numeric) is convenient for constructing Christoffel symbols and the Coriolis matrix directly from \( \mathbf{M}(\mathbf{q}) \).

import sympy as sp

# General n-DOF symbolic construction
n = 2
q1, q2 = sp.symbols('q1 q2')
dq1, dq2 = sp.symbols('dq1 dq2')
a1, a2, a3 = sp.symbols('a1 a2 a3')

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

# Inertia matrix for 2R example
M = sp.Matrix([
    [a1 + 2*a2*sp.cos(q2), a3 + a2*sp.cos(q2)],
    [a3 + a2*sp.cos(q2),   a3]
])

# Christoffel symbols c[i][j][k]
c = [[[0 for k in range(n)] for j in range(n)] for i in range(n)]

for i in range(n):
    for j in range(n):
        for k in range(n):
            c[i][j][k] = sp.Rational(1, 2) * (
                sp.diff(M[i, j], q[k]) +
                sp.diff(M[i, k], q[j]) -
                sp.diff(M[j, k], q[i])
            )

# Coriolis matrix C(q, dq) with convention C_ij = sum_k c_ijk dq_k
C = sp.Matrix.zeros(n, n)
for i in range(n):
    for j in range(n):
        C[i, j] = sum(c[i][j][k] * dq[k] for k in range(n))

print("M(q) =")
sp.pprint(M)
print("\nC(q, dq) =")
sp.pprint(C)

# Optional: generate fast numerical functions (e.g., for real-time control)
M_fun = sp.lambdify((q1, q2, a1, a2, a3), M, "numpy")
C_fun = sp.lambdify((q1, q2, dq1, dq2, a2), C, "numpy")
      

In a numerical control loop, one would call M_fun and C_fun with the current state to compute the inertia and Coriolis matrices. For larger DOF, the nested loops generalize directly.

In C++, dynamic libraries such as RBDL, Pinocchio, or KDL internally use Christoffel symbols or equivalent formulations. For educational purposes, we implement the 2R Coriolis matrix explicitly. This style of implementation is typical in embedded control software.

#include <array>
#include <cmath>

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double, 2>, 2>;

// Compute Coriolis matrix C(q, q_dot) for 2R elbow manipulator
// with inertia matrix
// M11 = a1 + 2 a2 cos(q2)
// M12 = a3 + a2 cos(q2)
// M22 = a3
Mat2 coriolis2R(const Vec2& q, const Vec2& dq,
                double a1, double a2, double a3)
{
    (void)a1; // unused in Coriolis
    (void)a3; // unused in Coriolis

    double q2  = q[1];
    double dq1 = dq[0];
    double dq2 = dq[1];

    double s2  = std::sin(q2);
    double h   = -a2 * s2;  // scalar h(q) often used in robotics

    Mat2 C{};
    // One consistent convention:
    // C11 = h * dq2
    // C12 = h * (dq1 + dq2)
    // C21 = -h * dq1
    // C22 = 0.0
    C[0][0] = h * dq2;
    C[0][1] = h * (dq1 + dq2);
    C[1][0] = -h * dq1;
    C[1][1] = 0.0;

    return C;
}

// Example usage in a torque computation:
// tau = M(q) q_ddot + C(q, q_dot) q_dot + g(q)
Vec2 coriolisTorque2R(const Vec2& q, const Vec2& dq,
                      double a1, double a2, double a3)
{
    Mat2 C = coriolis2R(q, dq, a1, a2, a3);
    Vec2 result{0.0, 0.0};
    result[0] = C[0][0] * dq[0] + C[0][1] * dq[1];
    result[1] = C[1][0] * dq[0] + C[1][1] * dq[1];
    return result;
}
      

In practice, coriolis2R can be generated symbolically (Python, MATLAB, or Mathematica) and then embedded as optimized C++ code to guarantee both correctness and real-time performance.

Java is often used in educational simulators and some robotics middleware. The following example computes the Coriolis matrix and corresponding torque term for the same 2R model.

public class TwoRCoriolis {

    // Compute C(q, qDot) for 2R elbow manipulator
    // q[0] = q1, q[1] = q2
    public static double[][] coriolis(double[] q, double[] qDot,
                                      double a1, double a2, double a3) {
        // a1, a3 unused in Coriolis but included for interface symmetry
        double q2  = q[1];
        double dq1 = qDot[0];
        double dq2 = qDot[1];

        double s2  = Math.sin(q2);
        double h   = -a2 * s2;

        double[][] C = new double[2][2];
        C[0][0] = h * dq2;
        C[0][1] = h * (dq1 + dq2);
        C[1][0] = -h * dq1;
        C[1][1] = 0.0;
        return C;
    }

    // Compute vector C(q, qDot) qDot
    public static double[] coriolisTorque(double[] q, double[] qDot,
                                          double a1, double a2, double a3) {
        double[][] C = coriolis(q, qDot, a1, a2, a3);
        double[] tauC = new double[2];
        tauC[0] = C[0][0] * qDot[0] + C[0][1] * qDot[1];
        tauC[1] = C[1][0] * qDot[0] + C[1][1] * qDot[1];
        return tauC;
    }
}
      

This pattern generalizes to larger manipulators by using higher-dimensional arrays or matrix libraries (e.g., EJML), with \( c_{ijk} \) precomputed offline.

8.1 MATLAB symbolic derivation and Simulink integration

MATLAB's Symbolic Math Toolbox provides an interface similar to sympy. The code below derives Christoffel symbols and exports a numeric function suitable for use in a Simulink block.

syms q1 q2 dq1 dq2 a1 a2 a3 real
q  = [q1; q2];
dq = [dq1; dq2];

% Inertia matrix
M = [a1 + 2*a2*cos(q2), a3 + a2*cos(q2);
     a3 + a2*cos(q2),   a3];

n = 2;
c = sym(zeros(n, n, n));

for i = 1:n
  for j = 1:n
    for k = 1:n
      c(i,j,k) = 1/2 * ( diff(M(i,j), q(k)) ...
                       + diff(M(i,k), q(j)) ...
                       - diff(M(j,k), q(i)) );
    end
  end
end

C = sym(zeros(n, n));
for i = 1:n
  for j = 1:n
    C(i,j) = c(i,j,1)*dq1 + c(i,j,2)*dq2;
  end
end

C = simplify(C)

% Export to a numeric MATLAB function for Simulink
matlabFunction(C, 'File', 'Coriolis2R', ...
    'Vars', [q1, q2, dq1, dq2, a1, a2, a3]);
      

In Simulink, one can place a MATLAB Function block that calls Coriolis2R and multiplies the result by \(\dot{\mathbf{q}}\) to obtain the Coriolis/centrifugal torque contribution inside a joint-space dynamics subsystem.

8.2 Wolfram Mathematica symbolic construction

Mathematica's symbolic capabilities allow very concise definitions of Christoffel symbols and Coriolis matrix:

(* General symbolic setup *)
ClearAll["Global`*"];
n   = 2;
q   = Array[q, n];
dq  = Array[dq, n];
{a1, a2, a3} = Array[a, 3];

M = { {a1 + 2 a2 Cos[q[[2]]], a3 + a2 Cos[q[[2]]]},
     {a3 + a2 Cos[q[[2]]],   a3} };

(* Christoffel symbols of the second kind *)
c[i_, j_, k_] :=
  1/2 (D[M[[i, j]], q[[k]]] +
       D[M[[i, k]], q[[j]]] -
       D[M[[j, k]], q[[i]]]);

C = Array[
      Function[{i, j},
        Sum[c[i, j, k] dq[[k]], {k, 1, n}]
      ],
      {n, n}
    ] // Simplify;

MatrixForm[M]
MatrixForm[C]
      

For higher-DOF manipulators, this pattern generalizes by increasing n and inserting the appropriate \( \mathbf{M}(\mathbf{q}) \). The resulting expressions can be exported as C code via CCodeGenerate for embedding in low-level controllers.

The following conceptual flow summarizes a generic algorithm for computing the Coriolis matrix from the inertia matrix in software, independent of programming language.

flowchart TD
  S["Start: inertia matrix M(q)"] --> D1["For each i, j, k compute c_ijk(q)"]
  D1 --> D2["Form C_ij(q, q_dot) = sum_k c_ijk(q) q_dot_k"]
  D2 --> D3["Evaluate C(q, q_dot) at current state"]
  D3 --> D4["Compute torque term: tau_C = C(q, q_dot) q_dot"]
  D4 --> E["Use in dynamics: M q_ddot + tau_C + g(q) = tau"]
        

In this lesson we connected the configuration-dependent inertia matrix \( \mathbf{M}(\mathbf{q}) \) to the Christoffel symbols of the second kind \( c_{ijk}(\mathbf{q}) \), and from there to the Coriolis/centrifugal matrix \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \). We derived the standard expression for \( c_{ijk} \) directly from the kinetic energy and Lagrange's equations, established the skew-symmetry of \( \dot{\mathbf{M}} - 2\mathbf{C} \), and interpreted this structure in terms of energy conservation.

On the computational side, we showed how to implement these constructions in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica for a 2R example and in general. These results form the theoretical and algorithmic foundation upon which subsequent lessons will build, including energy-based analysis, passivity proofs, and multi-DOF manipulator examples.