Dynamics for Sphere Linkages

Computing a metric for which the natural physical motion of the linkages are geodesics

Third in a series on spherical configuration spaces of planar linkages, joint work with Aaron Abrams, Dave Bachmann, and Edmund Harriss. The setup is in the first post.

In the previous post we built explicit smooth coordinates on $\mathcal{M}^n_L \cong S^{n-2}$ — a homeomorphism $\Psi^n_L : S^{n-2} \to \mathcal{M}^n_L$ for every $L$ in the regime $n - 2 < L < n$. Combined with the topological result of post 1, we can now name a configuration by a point on a sphere, drag it in code, and watch the linkage flex.

What we haven’t yet asked is the question that turns this from a topological curiosity into physics: how does the linkage actually move?

Pin the endpoints to a frictionless table, give the linkage a push, and let it go. With no friction and no potential, Lagrangian mechanics tells us the trajectory is the path minimizing the action $\int T\,dt$, where $T$ is the linkage’s kinetic energy.

Here is the geometric reframing — the new content of this post. At each configuration, kinetic energy is a positive-definite symmetric bilinear form on the tangent space, quadratic in the velocity. That’s the same thing as a Riemannian metric $g$:

$$T \;=\; \tfrac{1}{2}\, g(\dot\gamma, \dot\gamma).$$

The action of a path now reads

$$\int T\, dt \;=\; \tfrac{1}{2}\int g(\dot\gamma, \dot\gamma)\, dt,$$

which up to the factor of $\tfrac{1}{2}$ is the Riemannian energy functional. Its critical curves — the paths picked out as physical trajectories — are exactly the constant-speed geodesics of $g$. Lagrange’s variational principle and the variational principle for geodesics coincide.

Theorem. Force-free trajectories of a mechanical system are constant-speed geodesics of the kinetic-energy metric on its configuration space.

So to understand how the linkage moves, we equip $\mathcal{M}^n_L$ with its kinetic metric and find the geodesics. The metric depends on the physical model — where the mass lives — which we set up in §1. Section 2 pulls it back along $\Psi^n_L$ to a metric on the sphere; §3 computes it explicitly for $\mathcal{M}^4_L \cong S^2$. The geodesics themselves we approach two ways: §4 hands the §3 metric to a standard Christoffel-based integrator, while §5 derives them directly from the embedding into $\mathbb{R}^6$, never touching the metric formulas. The two serve as independent cross-checks on each other.

The kinetic metric on $(\mathbb{R}^2)^{n-1}$

A chain has $n + 1$ vertices $\mathbf{p}_0, \mathbf{p}_1, \ldots, \mathbf{p}_n$ connected by $n$ unit-length rigid rods, with $\mathbf{p}_0 = 0$ and $\mathbf{p}_n = L$ pinned. The $n - 1$ interior hinges are the moving parts.

To set up dynamics we specify where the mass lives. The most general choice has

• a point mass $M_j$ at each interior hinge $\mathbf{p}_j$, and
• a linear mass density $\lambda_k(s)$ along each rod, parameterized by arclength $s \in [0, 1]$.

(Mass on the pinned endpoints doesn’t move, so it doesn’t enter the dynamics.)

The rod-distributed mass is determined by the motion of its hinges — that’s what rigidity buys us. A material point at arclength $s$ on rod $k$ sits at

$$\mathbf{q}_k(s) \;=\; (1 - s)\,\mathbf{p}_{k-1} \;+\; s\,\mathbf{p}_k,$$

with velocity $\dot{\mathbf{q}}_k(s) = (1-s)\dot{\mathbf{p}}_{k-1} + s\,\dot{\mathbf{p}}_k$, linear in $s$. The total kinetic energy is

$$T \;=\; \tfrac{1}{2}\sum_{j=1}^{n-1} M_j\,|\dot{\mathbf{p}}_j|^2 \;+\; \tfrac{1}{2}\sum_{k=1}^{n}\int_0^1 \lambda_k(s)\,|\dot{\mathbf{q}}_k(s)|^2\,ds.$$

Expanding $|\dot{\mathbf{q}}_k(s)|^2$ inside the integral gives a quadratic form in $\dot{\mathbf{p}}_{k-1}$ and $\dot{\mathbf{p}}_k$ — with a cross term $\langle\dot{\mathbf{p}}_{k-1}, \dot{\mathbf{p}}_k\rangle$ from the $s(1-s)$ part of the expansion. So in general the kinetic-energy metric on $(\mathbb{R}^2)^{n-1}$ has nearest-neighbor hinge correlations, and the metric matrix is block-tridiagonal rather than diagonal.

For this post we make the simplest choice: massless rods ($\lambda_k \equiv 0$) and unit hinge masses ($M_j = 1$). The cross terms vanish, and

$$T \;=\; \tfrac{1}{2}\sum_{j=1}^{n-1} |\dot{\mathbf{p}}_j|^2.$$

This is exactly half the squared length of $\dot\gamma \in T(\mathbb{R}^2)^{n-1}$ in the standard Euclidean metric on $(\mathbb{R}^2)^{n-1}$. With this physical choice, the kinetic-energy metric on the unconstrained system is just Euclidean.

Our actual configuration space $\mathcal{M}^n_L$ sits inside $(\mathbb{R}^2)^{n-1}$ as a high-codimension submanifold: an $(n-2)$-dimensional sphere in a $2(n-1)$-dimensional ambient. The kinetic metric on $\mathcal{M}^n_L$ is the restriction of the Euclidean metric to tangent vectors of $\mathcal{M}^n_L$ — equivalently, the pullback of the Euclidean metric along the inclusion $\mathcal{M}^n_L \hookrightarrow (\mathbb{R}^2)^{n-1}$.

Pulling back to $S^{n-2}$

Composing the spherical parameterization $\Psi^n_L : S^{n-2} \to \mathcal{M}^n_L$ from post 2 (position-space version) with the inclusion gives a smooth embedding

$$S^{n-2} \;\xrightarrow{\Psi^n_L}\; \mathcal{M}^n_L \;\hookrightarrow\; (\mathbb{R}^2)^{n-1},$$

sending each sphere coordinate to its tuple of interior hinge positions $(\mathbf{p}_1, \ldots, \mathbf{p}_{n-1})$. Pulling the Euclidean metric back along this composite gives the kinetic metric on the sphere.

For a Euclidean metric, the pullback is just the sum of inner products of partials with respect to the sphere coordinates:

$$\boxed{\; h_{ab} \;=\; \sum_{k=1}^{n-1} \big\langle \partial_a \mathbf{p}_k,\; \partial_b \mathbf{p}_k\big\rangle, \qquad a, b \;\text{sphere coordinates}. \;}$$

The task is now concrete: take partials of the hinge positions with respect to sphere coordinates, then sum the inner products. We do this for $n = 4$ in the next section.

The metric on $\mathcal{M}^4_L \cong S^2$

We apply the pullback at $n = 4$, working directly with the position-space hinges from post 2. The plan: write down the hinge derivatives recursively, identify the named scalar partials that appear, assemble $h_{tt}, h_{\phi t}, h_{\phi\phi}$.

Hinge positions

Citing the boxed position-space formula from post 2 directly:

$$\mathbf{p}_3 \;=\; L - e^{i\theta}, \qquad \mathbf{p}_2 \;=\; \mathbf{p}_3 - \hat{\mathbf{p}}_3\,e^{i\theta_\text{in}}, \qquad \mathbf{p}_1 \;=\; \frac{\mathbf{p}_2}{2} + \frac{i\,\nu_\text{in}}{2\,|\mathbf{p}_2|}\,\mathbf{p}_2,$$

with $\theta = \alpha\sin\phi$, $\theta_\text{in} = \alpha_\text{in}\sin t$, and the helper definitions of $\alpha, \alpha_\text{in}, \nu_\text{in}$ from post 2. Three observations from there will do work below:

• $\mathbf{u}_4 := L - \mathbf{p}_3 = e^{i\theta}$ — the last rod direction (unit vector).
• $\mathbf{u}_3 := \mathbf{p}_3 - \mathbf{p}_2 = \hat{\mathbf{p}}_3\,e^{i\theta_\text{in}} = e^{i(R + \theta_\text{in})}$ — the third rod direction (unit vector).
• $|\mathbf{p}_1| = 1$ — the first hinge lies on the unit circle around the origin.

Derivatives

Three chain-rule identities, one for each hinge formula above. Below, $a$ stands for either of the sphere coordinates $\phi$ or $t$; each $\partial_a$ acts on the formula’s helpers via their own dependence on $a$.

$\partial_a \mathbf{p}_3$. Differentiating $\mathbf{p}_3 = L - e^{i\theta}$:

$$\partial_a \mathbf{p}_3 \;=\; -i\,\mathbf{u}_4\,\partial_a\theta.$$

The factor $-i$ rotates $90°$, so $\partial_a\mathbf{p}_3$ is perpendicular to the last rod with magnitude $|\partial_a\theta|$. Geometrically: a sphere coordinate that affects $\theta$ swings the last rod about its pinned end at $L$.

$\partial_a \mathbf{p}_2$. Differentiating $\mathbf{p}_2 = \mathbf{p}_3 - \mathbf{u}_3$ uses the unit-vector identity $\partial_a\mathbf{u}_3 = i\mathbf{u}_3\,\partial_a\theta_3$ where $\theta_3 = R + \theta_\text{in}$ is the third rod’s angle in the ambient frame:

$$\partial_a \mathbf{p}_2 \;=\; \partial_a \mathbf{p}_3 \;-\; i\,\mathbf{u}_3\,\partial_a\theta_3, \qquad \partial_a\theta_3 = \partial_a R + \partial_a\theta_\text{in}.$$

Two contributions: the next hinge moves ($\partial_a\mathbf{p}_3$), plus the third rod swings about its now-moving end at $\mathbf{p}_3$.

$\partial_a \mathbf{p}_1$. Since $|\mathbf{p}_1| = 1$, the innermost hinge moves on the unit circle around the origin: $\partial_a\mathbf{p}_1 = i\mathbf{p}_1 \cdot \Omega_a$ for some real angular velocity $\Omega_a$. From post 2, $\arg \mathbf{p}_1 = \arg \mathbf{p}_2 + \sigma\beta_\text{in}$, so $\Omega_a = \partial_a\arg\mathbf{p}_2 + \sigma\,\partial_a\beta_\text{in}$. The smooth replacement $\sigma\,\partial_a\beta_\text{in} = -\partial_a|\mathbf{p}_2|/\nu_\text{in}$ from post 2’s $\nu$-trick eliminates $\sigma$:

$$\partial_a \mathbf{p}_1 \;=\; i\,\mathbf{p}_1\,\Omega_a, \qquad \Omega_a \;=\; \operatorname{Im}\!\left(\frac{\partial_a \mathbf{p}_2}{\mathbf{p}_2}\right) \;-\; \frac{\partial_a|\mathbf{p}_2|}{\nu_\text{in}}.$$

Here $\Omega_a$ is a linear functional of $\partial_a\mathbf{p}_2$ — one cleanly named scalar carries everything we need to know about how the elbow swings.

The recursion above bottoms out in five named scalar partials, one for each helper that depends on $(\phi, t)$:

$$\boxed{\; \begin{aligned} \partial_\phi\theta \;&=\; \alpha\cos\phi, \\ \partial_\phi R \;&=\; \frac{1 - L\cos\theta}{d^2}\,\partial_\phi\theta, \\ \partial_\phi d \;&=\; \frac{L\sin\theta}{d}\,\partial_\phi\theta, \\ \partial_\phi\theta_\text{in} \;&=\; \partial_d\alpha_\text{in}\sin t \cdot \partial_\phi d, \\ \partial_t\theta_\text{in} \;&=\; \alpha_\text{in}\cos t, \end{aligned} \qquad \partial_d\alpha_\text{in} \;=\; -\,\frac{d^2 + 3}{2\, d^2 \sin\alpha_\text{in}}. \;}$$

(With $\partial_t\theta = \partial_t R = \partial_t d = 0$ since these depend on $\phi$ alone, and $\sin\alpha_\text{in} = \sqrt{(d^2-1)(9-d^2)}/(2d)$ on the valid range $d \in (1, 3)$.)

The chain $\partial_\phi\theta_\text{in}$ inherits the $\partial_d\alpha_\text{in}$ factor — the only derivative in the problem that doesn’t have a clean geometric interpretation. Everything else propagates through it.

Hinge velocities

Substituting the recursion into the auxiliary partials:

$$\begin{aligned} \partial_t\mathbf{p}_3 \;&=\; 0, &\partial_\phi\mathbf{p}_3 \;&=\; -i\,\mathbf{u}_4\,\partial_\phi\theta, \\[4pt] \partial_t\mathbf{p}_2 \;&=\; -i\,\mathbf{u}_3\,\partial_t\theta_\text{in}, &\partial_\phi\mathbf{p}_2 \;&=\; -i\,\mathbf{u}_4\,\partial_\phi\theta \;-\; i\,\mathbf{u}_3\,(\partial_\phi R + \partial_\phi\theta_\text{in}), \\[4pt] \partial_t\mathbf{p}_1 \;&=\; i\,\mathbf{p}_1\,\Omega_t, &\partial_\phi\mathbf{p}_1 \;&=\; i\,\mathbf{p}_1\,\Omega_\phi, \end{aligned}$$

with $\Omega_a$ (and $\partial_a|\mathbf{p}_2|$ inside it) evaluated on each side from the corresponding $\partial_a\mathbf{p}_2$:

$$\partial_a|\mathbf{p}_2| \;=\; \frac{\langle\mathbf{p}_2, \partial_a\mathbf{p}_2\rangle}{|\mathbf{p}_2|}.$$

Each hinge velocity is now a linear combination of the rod-direction unit vectors $\mathbf{u}_3, \mathbf{u}_4$ (with $i$-factors making them perpendicular) — except for $\partial_a\mathbf{p}_1$, which by unit-circle confinement is just $i\mathbf{p}_1$ times the scalar $\Omega_a$.

The metric components

Plug into $h_{ab} = \sum_k \langle\partial_a\mathbf{p}_k, \partial_b\mathbf{p}_k\rangle$. Two simplifications fall out for free:

• $|\partial_\phi\mathbf{p}_3|^2 = (\partial_\phi\theta)^2 = \alpha^2\cos^2\phi$ and $|\partial_t\mathbf{p}_2|^2 = (\partial_t\theta_\text{in})^2 = \alpha_\text{in}^2\cos^2 t$, since the rod directions are unit vectors.
• $|\partial\mathbf{p}_1|^2 = \Omega^2$ on either side (the angular velocity squared) since $|\mathbf{p}_1| = 1$, and $\langle\partial_\phi\mathbf{p}_1, \partial_t\mathbf{p}_1\rangle = \Omega_\phi\Omega_t$.

For the $k = 2$ contributions involving $\partial_\phi\mathbf{p}_2$, the cross between the two rod-direction terms produces $\langle\mathbf{u}_3, \mathbf{u}_4\rangle = \cos(\theta_3 - \theta) = \cos(R + \theta_\text{in} - \theta)$ — the (signed) cosine of the angle between the third and fourth rods. Carrying this through:

$$\boxed{\; \begin{aligned} h_{tt} \;&=\; \Omega_t^2 \;+\; (\partial_t\theta_\text{in})^2, \\[4pt] h_{\phi t} \;&=\; \Omega_t\Omega_\phi \;+\; \partial_t\theta_\text{in}\,\big[\partial_\phi\theta\,\langle\mathbf{u}_3, \mathbf{u}_4\rangle + \partial_\phi R + \partial_\phi\theta_\text{in}\big], \\[4pt] h_{\phi\phi} \;&=\; \Omega_\phi^2 \;+\; 2\,(\partial_\phi\theta)^2 \;+\; (\partial_\phi R + \partial_\phi\theta_\text{in})^2 \;+\; 2\,\partial_\phi\theta\,(\partial_\phi R + \partial_\phi\theta_\text{in})\,\langle\mathbf{u}_3, \mathbf{u}_4\rangle. \end{aligned} \;}$$

(In $h_{\phi\phi}$, the two factors of $(\partial_\phi\theta)^2$ come from $|\partial_\phi\mathbf{p}_3|^2$ and the $\mathbf{u}_4$-piece of $|\partial_\phi\mathbf{p}_2|^2$ — varying $\phi$ moves $\mathbf{p}_3$ and also drags $\mathbf{p}_2$ via the dependency $\mathbf{p}_2 = \mathbf{p}_3 - \mathbf{u}_3$.)

A clean structural fact falls out: $h_{tt}$ depends only on inner-level helpers ($\alpha_\text{in}, \theta_\text{in}, \mathbf{p}_2, \nu_\text{in}$) — no $\partial_\phi R$, no $\partial_\phi d$, no last-rod direction. The reason is geometric: varying $t$ at fixed $\phi$ leaves the outer rotation $\hat{\mathbf{p}}_3$ frozen and the outer endpoint $\mathbf{p}_3$ stationary, so the inner config simply executes its own $\Psi^3_d$ motion at the inherited sub-base length $d = d(\phi)$. Rotation invariance of the Euclidean inner product gives $h_{tt}^{(4)}(\phi, t) = h_{tt}^{(3)}\big(t;\, d(\phi)\big)$ — the circle metric, evaluated at the (slowly $\phi$-varying) sub-base length.

The components above feed on a stack of helper formulas. Most appeared in the derivation, but laying them out together — angular-velocity functionals, velocity vectors, named scalar partials, and the helpers from posts 1–2 — gives a self-contained reference for what metric4 consumes.

The angular-velocity functionals $\Omega_a$ and the velocity vectors $\partial_a\mathbf{p}_2$ they consume:

$$\begin{aligned} \Omega_a \;&=\; \operatorname{Im}\!\big(\partial_a\mathbf{p}_2/\mathbf{p}_2\big) \;-\; \partial_a|\mathbf{p}_2|\,/\,\nu_\text{in}, \qquad a \in \{\phi, t\}, \\[4pt] \partial_t\mathbf{p}_2 \;&=\; -i\,\mathbf{u}_3\,\partial_t\theta_\text{in}, \\[2pt] \partial_\phi\mathbf{p}_2 \;&=\; -i\,\mathbf{u}_4\,\partial_\phi\theta \;-\; i\,\mathbf{u}_3\,(\partial_\phi R + \partial_\phi\theta_\text{in}), \\[2pt] \partial_a|\mathbf{p}_2| \;&=\; \langle\mathbf{p}_2, \partial_a\mathbf{p}_2\rangle\,/\,|\mathbf{p}_2|. \end{aligned}$$

The named scalar partials:

$$\begin{aligned} \partial_\phi\theta \;&=\; \alpha\cos\phi, & \partial_t\theta_\text{in} \;&=\; \alpha_\text{in}\cos t, \\[2pt] \partial_\phi R \;&=\; \frac{1 - L\cos\theta}{d^2}\,\partial_\phi\theta, & \partial_\phi\theta_\text{in} \;&=\; \partial_d\alpha_\text{in}\,\sin t \cdot \partial_\phi d, \\[2pt] \partial_\phi d \;&=\; \frac{L\sin\theta}{d}\,\partial_\phi\theta, & \partial_d\alpha_\text{in} \;&=\; -\frac{d^2 + 3}{2\,d^2\,\sin\alpha_\text{in}}. \end{aligned}$$

And the helpers from posts 1–2:

$$\begin{aligned} \alpha \;&=\; \arccos\!\tfrac{L^2 - 8}{2L}, & \alpha_\text{in} \;&=\; \arccos\!\tfrac{d^2 - 3}{2d}, \\[2pt] \theta \;&=\; \alpha\sin\phi, & \theta_\text{in} \;&=\; \alpha_\text{in}\sin t, \\[2pt] d \;&=\; \sqrt{L^2 - 2L\cos\theta + 1}, & \sin\alpha_\text{in} \;&=\; \sqrt{(d^2 - 1)(9 - d^2)}\,/\,(2d), \\[2pt] R \;&=\; \arg\!\big(L - e^{i\theta}\big), & \nu_\text{in} \;&=\; \cos t\,\sqrt{2\,d\,(\cos\theta_\text{in} - \cos\alpha_\text{in})\,/\,\cos^2 t}, \\[6pt] \mathbf{p}_3 \;&=\; L - e^{i\theta}, & \mathbf{u}_4 \;&=\; e^{i\theta}, \\[2pt] \mathbf{p}_2 \;&=\; \mathbf{p}_3 - \mathbf{u}_3, & \mathbf{u}_3 \;&=\; e^{i(R + \theta_\text{in})}, \\[6pt] \langle\mathbf{u}_3,\mathbf{u}_4\rangle \;&=\; \cos(R + \theta_\text{in} - \theta). \end{aligned}$$

In code

function metric4(phi, t, L) {
  // Outer-level helpers and their phi-partials
  const alpha   = Math.acos((L*L - 8) / (2*L));
  const cphi    = Math.cos(phi), sphi = Math.sin(phi);
  const theta   = alpha * sphi;
  const ctheta  = Math.cos(theta), stheta = Math.sin(theta);
  const d       = Math.sqrt(L*L - 2*L*ctheta + 1);
  const R       = Math.atan2(-stheta, L - ctheta);
  const dthetaP = alpha * cphi;
  const ddP     = (L * stheta / d) * dthetaP;
  const dRP     = ((1 - L * ctheta) / (d * d)) * dthetaP;

  // Inner-level helpers
  const alpha_in = Math.acos((d*d - 3) / (2*d));
  const sin_ai   = Math.sqrt((d*d - 1) * (9 - d*d)) / (2*d);
  const ct = Math.cos(t), st = Math.sin(t);
  const theta_in = alpha_in * st;
  const ctin = Math.cos(theta_in), stin = Math.sin(theta_in);

  const dthetainT = alpha_in * ct;
  const dalphainP = -((d*d + 3) / (2 * d * d * sin_ai)) * ddP;
  const dthetainP = dalphainP * st;

  // Hinges (as (re, im) pairs)
  const p3 = [L - ctheta, -stheta];
  const u4 = [ctheta, stheta];                                   // last rod
  const u3_arg = R + theta_in;
  const u3 = [Math.cos(u3_arg), Math.sin(u3_arg)];               // third rod
  const p2 = [p3[0] - u3[0], p3[1] - u3[1]];
  const p2_abs = Math.hypot(p2[0], p2[1]);
  // nu_in: smooth signed amplitude. The ratio (4 - r^2)/cos^2(t) extends
  // to d * alpha_in * sin(alpha_in) at cos t = 0; use a Taylor branch there.
  const nu_in = Math.abs(ct) < 1e-6
              ? ct * Math.sqrt(d * alpha_in * sin_ai)
              : ct * Math.sqrt((4 - p2_abs * p2_abs) / (ct * ct));

  // Hinge velocities (as (re, im) pairs); -i * u * s = (u_im * s, -u_re * s)
  const dt_p2 = [u3[1] * dthetainT, -u3[0] * dthetainT];
  const c1 = dthetaP, c2 = dRP + dthetainP;
  const dphi_p2 = [u4[1]*c1 + u3[1]*c2, -(u4[0]*c1 + u3[0]*c2)];

  // Omega functional applied to a velocity v of p2
  const Omega = (v) => {
    // Im(v / p2) = (p2_re * v_im - p2_im * v_re) / |p2|^2
    const im = (p2[0]*v[1] - p2[1]*v[0]) / (p2_abs * p2_abs);
    // <p2, v> / (|p2| * nu_in)
    const re = (p2[0]*v[0] + p2[1]*v[1]) / (p2_abs * nu_in);
    return im - re;
  };
  const Omega_t   = Omega(dt_p2);
  const Omega_phi = Omega(dphi_p2);

  // Inner products
  const u34_dot = u3[0]*u4[0] + u3[1]*u4[1];      // <u3, u4>

  const h_tt = Omega_t * Omega_t + dthetainT * dthetainT;
  const h_pp = Omega_phi * Omega_phi
             + 2 * dthetaP * dthetaP
             + (dRP + dthetainP) * (dRP + dthetainP)
             + 2 * dthetaP * (dRP + dthetainP) * u34_dot;
  const h_pt = Omega_t * Omega_phi
             + dthetainT * (dthetaP * u34_dot + dRP + dthetainP);

  return [[h_pp, h_pt], [h_pt, h_tt]];
}

The metric matrix is what we hand to a geodesic solver in §4.

Geodesics

Physical trajectories of the linkage are constant-speed geodesics of $g$. Writing $\dot{}$ for derivative with respect to the geodesic parameter (distinct from our coordinate $t$), and following the 2D geodesic cheatsheet, the geodesic ODEs in our chart $(\phi, t)$ are

$$\boxed{\; \begin{aligned} \ddot\phi \;+\; \dot\phi^2\,\Gamma^\phi_{\phi\phi} \;+\; 2\dot\phi\dot t\,\Gamma^\phi_{\phi t} \;+\; \dot t^2\,\Gamma^\phi_{tt} \;&=\; 0, \\[4pt] \ddot t \;+\; \dot\phi^2\,\Gamma^t_{\phi\phi} \;+\; 2\dot\phi\dot t\,\Gamma^t_{\phi t} \;+\; \dot t^2\,\Gamma^t_{tt} \;&=\; 0. \end{aligned} \;}$$

The Christoffel symbols are determined by $E = h_{\phi\phi}, F = h_{\phi t}, G = h_{tt}$ together with the inverse-metric components $\mathcal{E} = G/D$, $\mathcal{F} = -F/D$, $\mathcal{G} = E/D$ (where $D = EG - F^2$):

$$\begin{aligned} \Gamma^\phi_{\phi\phi} \;&=\; \tfrac{1}{2}\big[\mathcal{E}\,\partial_\phi E + \mathcal{F}\,(2\partial_\phi F - \partial_t E)\big], &\Gamma^t_{\phi\phi} \;&=\; \tfrac{1}{2}\big[\mathcal{F}\,\partial_\phi E + \mathcal{G}\,(2\partial_\phi F - \partial_t E)\big], \\[2pt] \Gamma^\phi_{tt} \;&=\; \tfrac{1}{2}\big[\mathcal{E}\,(2\partial_t F - \partial_\phi G) + \mathcal{F}\,\partial_t G\big], &\Gamma^t_{tt} \;&=\; \tfrac{1}{2}\big[\mathcal{F}\,(2\partial_t F - \partial_\phi G) + \mathcal{G}\,\partial_t G\big], \\[2pt] \Gamma^\phi_{\phi t} \;&=\; \tfrac{1}{2}\big[\mathcal{E}\,\partial_t E + \mathcal{F}\,\partial_\phi G\big], &\Gamma^t_{\phi t} \;&=\; \tfrac{1}{2}\big[\mathcal{F}\,\partial_t E + \mathcal{G}\,\partial_\phi G\big]. \end{aligned}$$

Six metric partials are needed: $\partial_\phi E, \partial_t E, \partial_\phi F, \partial_t F, \partial_\phi G, \partial_t G$. These could be hand-derived by chain-ruling through the §3 helpers, but the result is six unwieldy symbolic expressions; we instead compute them numerically by central differences of the metric4 function.

function geodesicVectorField([phi, t, phidot, tdot], L) {
  const eps = 1e-5;
  const m   = metric4(phi,       t,       L);
  const mPh = metric4(phi + eps, t,       L);
  const mPl = metric4(phi - eps, t,       L);
  const mTh = metric4(phi,       t + eps, L);
  const mTl = metric4(phi,       t - eps, L);

  const E = m[0][0], F = m[0][1], G = m[1][1];
  const Ep = (mPh[0][0] - mPl[0][0]) / (2*eps);
  const Fp = (mPh[0][1] - mPl[0][1]) / (2*eps);
  const Gp = (mPh[1][1] - mPl[1][1]) / (2*eps);
  const Et = (mTh[0][0] - mTl[0][0]) / (2*eps);
  const Ft = (mTh[0][1] - mTl[0][1]) / (2*eps);
  const Gt = (mTh[1][1] - mTl[1][1]) / (2*eps);

  // Inverse metric components (script E, F, G in cheatsheet)
  const D  = E*G - F*F;
  const Ei = G / D;
  const Fi = -F / D;
  const Gi = E / D;

  // Christoffel symbols: Gam_<upper>_<lower lower>
  const Gam_p_pp = 0.5 * (Ei*Ep + Fi*(2*Fp - Et));
  const Gam_t_pp = 0.5 * (Fi*Ep + Gi*(2*Fp - Et));
  const Gam_p_tt = 0.5 * (Ei*(2*Ft - Gp) + Fi*Gt);
  const Gam_t_tt = 0.5 * (Fi*(2*Ft - Gp) + Gi*Gt);
  const Gam_p_pt = 0.5 * (Ei*Et + Fi*Gp);
  const Gam_t_pt = 0.5 * (Fi*Et + Gi*Gp);

  // Geodesic accelerations
  const phiddot = -(phidot*phidot * Gam_p_pp
                  + 2*phidot*tdot * Gam_p_pt
                  + tdot*tdot     * Gam_p_tt);
  const tddot   = -(phidot*phidot * Gam_t_pp
                  + 2*phidot*tdot * Gam_t_pt
                  + tdot*tdot     * Gam_t_tt);

  return [phidot, tdot, phiddot, tddot];
}

Hand this to any standard solver with initial state $(\phi_0, t_0, \dot\phi_0, \dot t_0)$, and the integrated trajectory is the linkage’s force-free motion. One coordinate caveat: the chart $(\phi, t)$ is singular at $\phi = \pm\pi/2$, where $\sin\alpha_\text{in} \to 0$ and several helpers diverge — a coordinate-not-geometric artifact, like polar coordinates at the origin. Trajectories that approach a pole need to be integrated in a different chart and transferred back across the overlap; but we do not do that here. Instead, if your trajectory hits exactly the north or south pole, the simulation just breaks…but this is of course rare.

Demo

Here’s a small program running these dynamics. Click and drag-to-release the point on configuraiton space: this sets an initial position and velocity, which is then integrated to produce the geodesic of the kinetic energy metric.

Geodesics from the embedding

The pullback metric and geodesic equations are quite complicated, leaving plenty of room to make a small error. Even after going through it all carefully, I was still worried an errant minus sign could easily have snuck through, and decided a good test would be to compute an independent derivation of the dynamics, and compare the two. Happily, the extrinsic geometry of our configuration space provides exactly such a means.

From §1, with our physical choice of unit point masses at the hinges and massless rods, the kinetic metric on the unconstrained system $(\mathbb{R}^2)^3$ is the standard Euclidean metric on $\mathbb{R}^6$. So $\mathcal{M}^4_L$ sits inside a 6-dimensional Euclidean space as a 2-dimensional surface, with metric induced from the ambient. That setting opens a second route to the geodesic equations — one that bypasses the §3 metric formulas entirely and gives us an independent computation we can cross-check against.

There’s a fact from undergraduate differential geometry that’s tailor-made for this situation. Take any surface $\Sigma$ in $\mathbb{R}^3$ — a sphere, a torus — and ask which curves on it are geodesics. The answer: a curve $\gamma : I \to \Sigma$ is a geodesic iff at every point its acceleration $\ddot\gamma$ is normal to $\Sigma$. The canonical picture is a great circle on the sphere: at each point the velocity bends toward the center, and the centripetal acceleration points radially inward — perpendicular to the sphere.

The same statement works in any Euclidean dimension. Our setting just has two cosmetic differences from the textbook one. The ambient is $\mathbb{R}^6$ rather than $\mathbb{R}^3$, so points have six coordinates instead of three. And the surface is 2-dimensional in 6 (codimension 4) rather than 2-dimensional in 3 (codimension 1), so the normal space at each point is 4-dimensional rather than 1-dimensional. Neither matters for the principle. “Normal to the surface” still means “no tangential component,” and the surface still has two tangent directions, so the normality condition produces exactly two scalar equations — the right number for a second-order ODE on a 2-manifold.

So the plan is concrete: compute the acceleration of a candidate curve in $\mathbb{R}^6$, project it onto the two tangent directions of the surface, and demand both projections vanish.

Step 1: the acceleration

Pick a curve in $S^2$ in our chart coordinates, $(\phi(s), t(s))$. Pushing through the parameterization $\Psi^4_L$ — which we’ll just write $q(\phi, t) = (\mathbf{p}_1, \mathbf{p}_2, \mathbf{p}_3) \in \mathbb{R}^6$, the §3 hinge positions stacked into one vector — gives the corresponding curve in $\mathbb{R}^6$:

$$q(s) \;=\; q\big(\phi(s),\, t(s)\big).$$

The acceleration is just the second derivative with respect to $s$. Chain-ruling once gives $\dot q = q_\phi\,\dot\phi + q_t\,\dot t$, and differentiating again:

$$\ddot q \;=\; q_\phi\,\ddot\phi \;+\; q_t\,\ddot t \;+\; q_{\phi\phi}\,\dot\phi^2 \;+\; 2\,q_{\phi t}\,\dot\phi\dot t \;+\; q_{tt}\,\dot t^2.$$

Step 2: impose normality

A vector in $\mathbb{R}^6$ is normal to the surface at $q(\phi, t)$ when it’s orthogonal to both basis tangent vectors $q_\phi$ and $q_t$. So our geodesic condition $\ddot q \perp T_q\Sigma$ unpacks as two scalar equations:

$$\langle q_\phi,\, \ddot q\rangle \;=\; 0, \qquad \langle q_t,\, \ddot q\rangle \;=\; 0.$$

Plug in the expression for $\ddot q$ from Step 1 and distribute the inner products:

$$\langle q_\phi, q_\phi\rangle\,\ddot\phi \;+\; \langle q_\phi, q_t\rangle\,\ddot t \;+\; \langle q_\phi, q_{\phi\phi}\rangle\,\dot\phi^2 \;+\; 2\langle q_\phi, q_{\phi t}\rangle\,\dot\phi\dot t \;+\; \langle q_\phi, q_{tt}\rangle\,\dot t^2 \;=\; 0,$$ $$\langle q_t, q_\phi\rangle\,\ddot\phi \;+\; \langle q_t, q_t\rangle\,\ddot t \;+\; \langle q_t, q_{\phi\phi}\rangle\,\dot\phi^2 \;+\; 2\langle q_t, q_{\phi t}\rangle\,\dot\phi\dot t \;+\; \langle q_t, q_{tt}\rangle\,\dot t^2 \;=\; 0.$$

Each equation is linear in the unknown second derivatives $\ddot\phi, \ddot t$ — everything else is determined by the curve’s current position and velocity. Moving the velocity-quadratic pieces to the right and stacking the two equations:

$$\boxed{\; \begin{pmatrix} \langle q_\phi, q_\phi\rangle & \langle q_\phi, q_t\rangle \\ \langle q_t, q_\phi\rangle & \langle q_t, q_t\rangle \end{pmatrix} \begin{pmatrix} \ddot\phi \\ \ddot t \end{pmatrix} \;=\; -\begin{pmatrix} \langle q_\phi, q_{\phi\phi}\rangle\,\dot\phi^2 + 2\langle q_\phi, q_{\phi t}\rangle\,\dot\phi\dot t + \langle q_\phi, q_{tt}\rangle\,\dot t^2 \\ \langle q_t, q_{\phi\phi}\rangle\,\dot\phi^2 + 2\langle q_t, q_{\phi t}\rangle\,\dot\phi\dot t + \langle q_t, q_{tt}\rangle\,\dot t^2 \end{pmatrix}. \;}$$

A 2×2 linear system. Inverting solves for $(\ddot\phi, \ddot t)$ — the second-order ODE on $S^2$ that integrates to a geodesic.

The coefficient matrix on the left is, of course, the metric tensor $g$ from §3 — its entries $\langle q_a, q_b\rangle$ are exactly the components $E, F, G$ we worked out there. We don’t need that observation to solve the system, but it’s worth noting the connection.

In code

The whole right-hand side is finite-difference operations on psi4Position from post 2.

function geodesicVectorFieldEmbedded([phi, t, phidot, tdot], L) {
  const eps = 1e-4;

  // q : (phi, t) -> R^6, using psi4Position from post 2.
  const Q = (ph, tt) => {
    const {p1, p2, p3} = psi4Position(ph, tt, L);
    return [p1[0], p1[1], p2[0], p2[1], p3[0], p3[1]];
  };

  // R^6 vector helpers
  const sub   = (a, b) => a.map((x, i) => x - b[i]);
  const add   = (...vs) => vs.reduce((s, v) => s.map((x, i) => x + v[i]));
  const scale = (a, c) => a.map(x => c * x);
  const dot   = (a, b) => a.reduce((s, x, i) => s + x * b[i], 0);

  // 3x3 stencil of position samples
  const Q00 = Q(phi, t);
  const Qpp = Q(phi + eps, t),  Qpm = Q(phi - eps, t);
  const Qtp = Q(phi, t + eps),  Qtm = Q(phi, t - eps);
  const Qpt_pp = Q(phi + eps, t + eps), Qpt_pm = Q(phi + eps, t - eps);
  const Qpt_mp = Q(phi - eps, t + eps), Qpt_mm = Q(phi - eps, t - eps);

  // First and second partials by central differences
  const q_phi    = scale(sub(Qpp, Qpm), 1/(2*eps));
  const q_t      = scale(sub(Qtp, Qtm), 1/(2*eps));
  const q_phiphi = scale(add(Qpp, scale(Q00, -2), Qpm), 1/(eps*eps));
  const q_tt     = scale(add(Qtp, scale(Q00, -2), Qtm), 1/(eps*eps));
  const q_phit   = scale(sub(sub(Qpt_pp, Qpt_pm), sub(Qpt_mp, Qpt_mm)), 1/(4*eps*eps));

  // Velocity-quadratic part of qddot
  const qdd_vel = add(
    scale(q_phiphi, phidot * phidot),
    scale(q_phit,   2 * phidot * tdot),
    scale(q_tt,     tdot * tdot),
  );

  // 2x2 system [<q_a, q_b>] (phiddot, tddot)^T = -[<q_a, qdd_vel>]
  const g_pp = dot(q_phi, q_phi);
  const g_pt = dot(q_phi, q_t);
  const g_tt = dot(q_t, q_t);
  const det  = g_pp * g_tt - g_pt * g_pt;
  const bP   = -dot(q_phi, qdd_vel);
  const bT   = -dot(q_t,   qdd_vel);
  const phiddot = (g_tt * bP - g_pt * bT) / det;
  const tddot   = (g_pp * bT - g_pt * bP) / det;

  return [phidot, tdot, phiddot, tddot];
}

This shares no code path with §4’s geodesicVectorField past psi4Position. Run both integrators from the same initial state: agreement validates the §3 pipeline (named partials, the $\Omega$ functional, the metric components, the Christoffel formulas) all at once. Disagreement diagnoses a bug somewhere along that pipeline — by elimination, since the embedded version exercises none of it.

One practical note. The $\nu_\text{in}$ singularity at $t = \pm\pi/2$ never enters here, provided psi4Position uses the smooth (sinc-form) expression for $\nu_\text{in}$ from post 2. The whole geodesic integration then sees a smooth $q$ with no removable singularity to worry about.

Comparison demo

The two integrators in this post — geodesicVectorField from §4 and geodesicVectorFieldEmbedded from §5 — should produce the same trajectories from the same initial state. Below, both run side-by-side: the trajectories drawn over each other on $S^2$, and $\Delta$ records their coordinate divergence over integration time.