Realistic Sphere Linkages

Redoing the calculations for configurations with mass distributed along the rods and at the vertices

Fifth 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.

The last two posts computed the dynamics and curvature of the configuration space $\mathcal{M}^4_L \cong S^2$, but with the simplest possible physical model: unit point masses at the interior hinges, massless rods. That’s the cleanest setting to set up the geometric machinery, but it’s a long way from a physical linkage you’d actually build. A real linkage will have both mass distributed along its rods (say they are made of plastic, or wood) and also mass at the hinges (metal bearings, and whatever else).

This post upgrades the physics. We pick two real parameters,

• $m_e$ — the total mass of each rod/edge,
• $m_v$ — the mass at each interior hinge/vertex,

with each rod assumed to have its mass distributed uniformly along its length. The previous posts are the special case $m_e = 0, m_v = 1$.

The plan: derive the kinetic-energy metric on the unconstrained ambient space $(\mathbb{R}^2)^{n-1}$, then use it numerically as the inner product for the post 4 pullback machinery to get first and second fundamental forms on the sphere, and finally redo the curvature plots and geodesic dynamics with the new metric. The ambient metric is the only piece that genuinely changes; everything downstream is the same code, with a different inner product plugged in.

The unconstrained metric

We work in the unconstrained ambient $(\mathbb{R}^2)^{n-1}$ with coordinates $(\mathbf{p}_1, \ldots, \mathbf{p}_{n-1})$, the moving interior hinges. Pinned endpoints $\mathbf{p}_0 = 0$ and $\mathbf{p}_n = L$ have $\dot{\mathbf{p}}_0 = \dot{\mathbf{p}}_n = 0$ and don’t enter the dynamics regardless of mass.

Total kinetic energy is the sum of the hinge contributions and the rod contributions:

$$T \;=\; \tfrac{1}{2}\,m_v\sum_{j=1}^{n-1}|\dot{\mathbf{p}}_j|^2 \;+\; \sum_{k=1}^n T_k,$$

with $T_k$ the kinetic energy of rod $k$. The hinge term is elementary; the rod term is where the work happens.

Stretchy rods

Here’s the subtlety with the rod term. In the constrained configuration space $\mathcal{M}^n_L$ every rod has length $1$, so saying “total mass $m_e$ per rod” is the same as saying “linear mass density $m_e$ per rod.” The two specifications agree on the constraint surface.

In the unconstrained ambient $(\mathbb{R}^2)^{n-1}$ they don’t. The rod from $\mathbf{p}_{k-1}$ to $\mathbf{p}_k$ has length $\ell_k = |\mathbf{p}_k - \mathbf{p}_{k-1}|$, which can be anything. Two natural conventions:

• Fixed total mass: rod $k$ has mass $m_e$ no matter what its length is. Stretching the rod thins it out; squashing concentrates the mass.
• Fixed linear density: rod $k$ has density $\lambda$, with total mass $\lambda\ell_k$ proportional to length. Stretching makes the rod heavier.

These give different ambient metrics on $(\mathbb{R}^2)^{n-1}$ — but, crucially, they agree on the constraint surface $\mathcal{M}^n_L$, where every $\ell_k = 1$ and the two specifications coincide. Since the ambient metric only enters our story through its restriction to $\mathcal{M}^n_L$, the pullback to the sphere is the same either way. We’re free to pick whichever ambient metric is more convenient to compute with.

We go with fixed total mass. Stretchy physical objects (rubber bands, springs) really do have this property — pulling on the band doesn’t conjure new material into being, it just spreads what’s there over a longer interval. Fixed linear density, by contrast, would have stretching create mass, which is harder to make sense of physically. And as we’ll see in a moment, the fixed-total-mass metric turns out to be configuration-independent — a single constant symmetric form on the whole ambient — which makes numerical pullback trivial.

So rod $k$ has uniform linear density

$$\lambda_k \;=\; \frac{m_e}{\ell_k},$$

depending on configuration through $\ell_k$, but with the configuration-dependence arranged so that the integrated mass is always $m_e$.

The rod integral

Parameterize each rod by a material coordinate $s \in [0, 1]$ — a label for points along the rod that doesn’t change as the rod stretches. The point at material parameter $s$ on rod $k$ sits in the plane at

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

with velocity (linear in $s$, regardless of the rod’s length or orientation)

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

The kinetic energy of the rod is the integral of $\tfrac{1}{2}|\dot{\mathbf{q}}|^2$ against the mass element. In arclength $\xi \in [0, \ell_k]$ the mass element is $\lambda_k\,d\xi$; switching to the material coordinate $s = \xi/\ell_k$ gives $d\xi = \ell_k\,ds$, and

$$\lambda_k\,d\xi \;=\; \tfrac{m_e}{\ell_k} \cdot \ell_k\,ds \;=\; m_e\,ds.$$

The two factors of $\ell_k$ cancel — the configuration-dependence of the density and the configuration-dependence of the arclength element annihilate each other. This is the geometric content of the fixed-total-mass spec, and it’s why this convention gives a clean metric.

So the rod-$k$ kinetic energy reduces to

$$T_k \;=\; \tfrac{1}{2}\int_0^1 |\dot{\mathbf{q}}_k(s)|^2 \,m_e\,ds,$$

an integral over a fixed interval against a constant measure.

Expanding the integrand,

$$|\dot{\mathbf{q}}_k(s)|^2 \;=\; (1-s)^2\,|\dot{\mathbf{p}}_{k-1}|^2 \;+\; 2\,s(1-s)\,\langle\dot{\mathbf{p}}_{k-1}, \dot{\mathbf{p}}_k\rangle \;+\; s^2\,|\dot{\mathbf{p}}_k|^2,$$

we see that everything can be written in terms of two simple integrals!

$$\int_0^1 (1-s)^2\,ds = \int_0^1 s^2\,ds = \tfrac{1}{3}$$ $$\int_0^1 s(1-s)\,ds = \tfrac{1}{6}$$

Plugging these in,

$$\boxed{\; T_k \;=\; \tfrac{m_e}{6}\,\Big[\,|\dot{\mathbf{p}}_{k-1}|^2 \;+\; \langle\dot{\mathbf{p}}_{k-1}, \dot{\mathbf{p}}_k\rangle \;+\; |\dot{\mathbf{p}}_k|^2\,\Big]. \;}$$

Sanity checks. Two limiting cases pin this down.

Pure translation: $\dot{\mathbf{p}}_{k-1} = \dot{\mathbf{p}}_k = \mathbf{v}$. Then $T_k = \tfrac{m_e}{6}\cdot 3|\mathbf{v}|^2 = \tfrac{m_e}{2}|\mathbf{v}|^2$ — the kinetic energy of a particle of mass $m_e$ moving at speed $|\mathbf{v}|$. ✓

Pure rotation about the rod’s center: $\dot{\mathbf{p}}_{k-1} = -\dot{\mathbf{p}}_k = \mathbf{v}$. Then $T_k = \tfrac{m_e}{6}|\mathbf{v}|^2$. Cross-check with elementary mechanics: the endpoints are at distance $\ell_k/2$ from the center moving at speed $|\mathbf{v}|$, so the angular velocity is $\omega = 2|\mathbf{v}|/\ell_k$. A uniform rod has moment of inertia $I = m_e\ell_k^2/12$ about its center, giving $\tfrac{1}{2}I\omega^2 = \tfrac{m_e}{6}|\mathbf{v}|^2$. ✓

Assembling the kinetic energy

Sum the rod contributions and add the hinge term, with $\dot{\mathbf{p}}_0 = \dot{\mathbf{p}}_n = 0$ at the pinned endpoints:

$$T \;=\; \tfrac{m_e}{6}\sum_{k=1}^n\Big[|\dot{\mathbf{p}}_{k-1}|^2 + \langle\dot{\mathbf{p}}_{k-1}, \dot{\mathbf{p}}_k\rangle + |\dot{\mathbf{p}}_k|^2\Big] \;+\; \tfrac{m_v}{2}\sum_{j=1}^{n-1}|\dot{\mathbf{p}}_j|^2.$$

This is correct, but the rod sum is awkward with the same variable appearing multiple places. To write it nicely as a quadratic form on the hinge positions $(\mathbf{p}_1, \ldots, \mathbf{p}_{n-1})$ we need to collect terms. It’s easier to see the pattern by working a small case first.

The 3-rod chain

Two interior hinges $\mathbf{p}_1, \mathbf{p}_2$, three rods. The two boundary rods (1 and 3) each touch only one interior hinge — the other end is pinned and contributes nothing:

$$T_1 = \tfrac{m_e}{6}|\dot{\mathbf{p}}_1|^2, \qquad T_3 = \tfrac{m_e}{6}|\dot{\mathbf{p}}_2|^2.$$

The middle rod (2) sits between both interior hinges and gives the full three-term expression:

$$T_2 = \tfrac{m_e}{6}\Big[|\dot{\mathbf{p}}_1|^2 + \langle\dot{\mathbf{p}}_1, \dot{\mathbf{p}}_2\rangle + |\dot{\mathbf{p}}_2|^2\Big].$$

Add the rod terms together: the diagonal $|\dot{\mathbf{p}}_1|^2$ appears in both $T_1$ and $T_2$ with coefficient $\tfrac{m_e}{6}$ each, totaling $\tfrac{m_e}{3}$. Same for $|\dot{\mathbf{p}}_2|^2$ between $T_2$ and $T_3$. The cross term $\langle\dot{\mathbf{p}}_1, \dot{\mathbf{p}}_2\rangle$ comes from $T_2$ alone with coefficient $\tfrac{m_e}{6}$. Adding the hinge masses:

$$T \;=\; \tfrac{1}{2}\Big(\tfrac{2m_e}{3} + m_v\Big)\Big[|\dot{\mathbf{p}}_1|^2 + |\dot{\mathbf{p}}_2|^2\Big] \;+\; \tfrac{1}{2}\cdot\tfrac{m_e}{3}\langle\dot{\mathbf{p}}_1, \dot{\mathbf{p}}_2\rangle.$$

The cross term is contributed entirely by the middle rod — the one that’s between the two interior hinges. The boundary rods don’t couple anything, because each of them has one end pinned. And the diagonal coefficient is the same for both hinges: $\tfrac{m_e}{3}$ from the two adjacent rods plus $m_v$ for the hinge itself.

Now to read off the metric tensor. The unconstrained ambient is $(\mathbb{R}^2)^2 = \mathbb{R}^4$, with linear coordinates $\mathbf{p}_1 = (p_{1x}, p_{1y})$ and $\mathbf{p}_2 = (p_{2x}, p_{2y})$. Expanding the inner products in terms of these scalar components:

$$|\dot{\mathbf{p}}_1|^2 = \dot p_{1x}^2 + \dot p_{1y}^2, \qquad |\dot{\mathbf{p}}_2|^2 = \dot p_{2x}^2 + \dot p_{2y}^2, \qquad \langle\dot{\mathbf{p}}_1, \dot{\mathbf{p}}_2\rangle = \dot p_{1x}\dot p_{2x} + \dot p_{1y}\dot p_{2y}.$$

Substituting and writing $T = \tfrac{1}{2}\dot{\mathbf{q}}^\top g\,\dot{\mathbf{q}}$ where $\dot{\mathbf{q}} = (\dot p_{1x}, \dot p_{1y}, \dot p_{2x}, \dot p_{2y})$, the metric tensor is the $4 \times 4$ matrix

$$g_3 \;=\; \begin{pmatrix} \tfrac{2m_e}{3} + m_v & 0 & \tfrac{m_e}{6} & 0 \\[4pt] 0 & \tfrac{2m_e}{3} + m_v & 0 & \tfrac{m_e}{6} \\[4pt] \tfrac{m_e}{6} & 0 & \tfrac{2m_e}{3} + m_v & 0 \\[4pt] 0 & \tfrac{m_e}{6} & 0 & \tfrac{2m_e}{3} + m_v \end{pmatrix}.$$

(The off-diagonals are $\tfrac{m_e}{6}$ rather than $\tfrac{m_e}{3}$ because the cross term $\langle\dot{\mathbf{p}}_1, \dot{\mathbf{p}}_2\rangle$ in $T$ gets distributed symmetrically between $g_{13} = g_{31}$ and $g_{24} = g_{42}$.)

Look at where the entries actually sit. The diagonal is constant, the only off-diagonals are $g_{13} = g_{31}$ and $g_{24} = g_{42}$, and the two zero blocks $g_{12}$ and $g_{34}$ tell us the metric never mixes the $x$ and $y$ components of any single hinge. In fact every $2 \times 2$ block of $g_3$ is a scalar multiple of $I_2$ — the diagonal blocks are $aI_2$ and the off-diagonal blocks are $bI_2$. That structure is exactly the tensor product:

$$g_3 \;=\; \begin{pmatrix} a I_2 & b I_2 \\ b I_2 & a I_2 \end{pmatrix} \;=\; M_3 \otimes I_2,$$ $$M_3 \;=\; \begin{pmatrix} a & b \\ b & a \end{pmatrix}, \qquad a \;=\; \tfrac{2m_e}{3} + m_v, \quad b \;=\; \tfrac{m_e}{6}.$$

The $2 \times 2$ factor $M_3$ records the inter-hinge coupling, and $\otimes I_2$ records that each hinge has identical and decoupled $x, y$ structure. We’ll see this factorization appear at every $n$, so for the rest of the post we’ll work with the smaller coupling matrices $M$ and tensor up to the full metric only when we need to.

The 4-rod chain

Three interior hinges $\mathbf{p}_1, \mathbf{p}_2, \mathbf{p}_3$, four rods. By analogy with the 3-rod case: the boundary rods (1 and 4) touch one interior hinge each, the interior rods (2 and 3) connect adjacent hinge pairs:

$$\begin{aligned} T_1 &= \tfrac{m_e}{6}|\dot{\mathbf{p}}_1|^2, \\ T_2 &= \tfrac{m_e}{6}\Big[|\dot{\mathbf{p}}_1|^2 + \langle\dot{\mathbf{p}}_1, \dot{\mathbf{p}}_2\rangle + |\dot{\mathbf{p}}_2|^2\Big], \\ T_3 &= \tfrac{m_e}{6}\Big[|\dot{\mathbf{p}}_2|^2 + \langle\dot{\mathbf{p}}_2, \dot{\mathbf{p}}_3\rangle + |\dot{\mathbf{p}}_3|^2\Big], \\ T_4 &= \tfrac{m_e}{6}|\dot{\mathbf{p}}_3|^2. \end{aligned}$$

Each interior hinge’s $|\dot{\mathbf{p}}_j|^2$ appears in exactly two of these — the rods on either side — so the diagonal coefficient is again $\tfrac{m_e}{3}$ for every hinge, including the boundary ones $j = 1$ and $j = 3$. (The boundary hinge $\mathbf{p}_1$ only has one interior rod next to it — but the boundary rod $T_1$ still picks up its motion, contributing the same $\tfrac{m_e}{6}$ as any other rod-end. The pinned endpoint is what gets killed, not the live one.) The cross terms come from the two interior rods. Collecting:

$$T \;=\; \tfrac{1}{2}\Big(\tfrac{2m_e}{3} + m_v\Big)\Big[|\dot{\mathbf{p}}_1|^2 + |\dot{\mathbf{p}}_2|^2 + |\dot{\mathbf{p}}_3|^2\Big] \;+\; \tfrac{1}{2}\cdot\tfrac{m_e}{3}\Big[\langle\dot{\mathbf{p}}_1, \dot{\mathbf{p}}_2\rangle + \langle\dot{\mathbf{p}}_2, \dot{\mathbf{p}}_3\rangle\Big].$$

The ambient is now $\mathbb{R}^6$ with coordinates $(p_{1x}, p_{1y}, p_{2x}, p_{2y}, p_{3x}, p_{3y})$. Expanding the inner products as before, the metric tensor is the $6 \times 6$ matrix

$$g_4 \;=\; \begin{pmatrix} a & 0 & b & 0 & 0 & 0 \\ 0 & a & 0 & b & 0 & 0 \\ b & 0 & a & 0 & b & 0 \\ 0 & b & 0 & a & 0 & b \\ 0 & 0 & b & 0 & a & 0 \\ 0 & 0 & 0 & b & 0 & a \end{pmatrix},$$

with $a = \tfrac{2m_e}{3} + m_v$ and $b = \tfrac{m_e}{6}$ as before. Same structure as $g_3$: zeros wherever $x$ and $y$ from the same or different hinges would mix, and a checkerboard of $a$‘s and $b$‘s on the diagonal-like positions where they don’t. The new feature is the missing $b$‘s connecting hinges 1 and 3 — at $g_{15}$ and $g_{26}$ — because the non-adjacent hinges $\mathbf{p}_1$ and $\mathbf{p}_3$ aren’t directly coupled. The matrix factors as

$$g_4 \;=\; M_4 \otimes I_2, \qquad M_4 \;=\; \begin{pmatrix} a & b & 0 \\ b & a & b \\ 0 & b & a \end{pmatrix},$$

with the same tensor-product structure as before. $M_4$ is tridiagonal — nonzero only on the main diagonal and the two adjacent off-diagonals — reflecting that no single rod touches two non-adjacent hinges.

The general formula

The pattern extends verbatim to any $n$. Every interior hinge contributes a diagonal term with coefficient $\tfrac{2m_e}{3} + m_v$, and every adjacent pair contributes a cross term with coefficient $\tfrac{m_e}{3}$ from the rod between them. The kinetic energy is

$$\boxed{\; T \;=\; \tfrac{1}{2}\bigg[a\sum_{j=1}^{n-1}|\dot{\mathbf{p}}_j|^2 \;+\; 2b\sum_{j=1}^{n-2}\langle\dot{\mathbf{p}}_j, \dot{\mathbf{p}}_{j+1}\rangle\bigg], \;}$$ $$a \;=\; \tfrac{2m_e}{3} + m_v, \qquad b \;=\; \tfrac{m_e}{6},$$

with corresponding $(n - 1) \times (n - 1)$ coupling matrix

$$\boxed{\; M \;=\; \begin{pmatrix} a & b & & & \\ b & a & b & & \\ & b & a & \ddots & \\ & & \ddots & \ddots & b \\ & & & b & a \end{pmatrix}. \;}$$

The full $2(n-1) \times 2(n-1)$ metric tensor on $(\mathbb{R}^2)^{n-1}$ is then $g = M \otimes I_2$.

Configuration-independence. No $\ell_k$ appears anywhere in $M$. The ambient kinetic metric is constant on $(\mathbb{R}^2)^{n-1}$ — the same symmetric form at every configuration.

In code

// (n-1) x (n-1) symmetric tridiagonal coupling matrix M.
// The full ambient metric on (R^2)^{n-1} is M ⊗ I_2.
// Kinetic energy: T = (1/2) Σ_{i,j} M_{ij} <p_dot_i, p_dot_j>.
function ambientCouplingMatrix(n, m_e, m_v) {
  const a = 2*m_e/3 + m_v;
  const b = m_e/6;
  const M = Array.from({length: n-1}, () => new Array(n-1).fill(0));
  for (let i = 0; i < n-1; i++) {
    M[i][i] = a;
    if (i+1 < n-1) M[i][i+1] = M[i+1][i] = b;
  }
  return M;
}

// M-weighted inner product on (R^2)^{n-1}, viewing U, V as arrays
// of (n-1) vectors in R^2: <U, V>_M = Σ_{i,j} M_{ij} <U_i, V_j>.
function weightedDot(M, U, V) {
  let s = 0;
  for (let i = 0; i < U.length; i++) {
    for (let j = 0; j < V.length; j++) {
      const m = M[i][j];
      if (m !== 0) s += m * (U[i][0]*V[j][0] + U[i][1]*V[j][1]);
    }
  }
  return s;
}

This is the only piece of new infrastructure the upgrade requires. Everything in posts 3 and 4 — pullback to the sphere, fundamental forms, Christoffel symbols, geodesic ODEs, Brioschi, the Gauss equation — is already phrased in terms of an ambient inner product, and we just swap the Euclidean one for weightedDot(M, ·, ·). The rest of this post is the consequences.

Fundamental forms

The first and second fundamental forms on $\mathcal{M}^4_L \cong S^2$ are computed exactly as in post 4, with the only change being the ambient inner product. Where post 4 used the standard dot product on $\mathbb{R}^6$, we now use $\langle \cdot,\cdot\rangle_g$ — the $M$-weighted inner product determined by the ambient metric tensor $g = M \otimes I_2$.

The embedding itself is unchanged — it’s the same $\Psi^4_L$ from post 2, mapping $S^2$ into $(\mathbb{R}^2)^3$. We just measure inner products of vectors in the ambient differently. For convenience we’ll work with the embedding as a list of three hinges $(\mathbf{p}_1, \mathbf{p}_2, \mathbf{p}_3)$ rather than flattening to $\mathbb{R}^6$, since weightedDot is naturally indexed that way.

The functions below all take the rod and hinge masses $(m_e, m_v)$ as direct arguments and build the coupling matrix internally on each call. This is overkill for performance but makes the functions self-contained, which is exactly what we want for hooking sliders up to interactive demos: changing a slider just changes the next argument value.

// Embedding S^2 -> (R^2)^3 as a list of three hinges.
function embedHinges(phi, t, L) {
  const { p1, p2, p3 } = psi4Position(phi, t, L);
  return [p1, p2, p3];
}

The first fundamental form

Pulling back the ambient metric along the embedding gives the metric on $S^2$:

$$h_{ab}(\phi, t) \;=\; \langle \partial_a\mathbf{p}, \partial_b\mathbf{p}\rangle_g, \qquad a, b \in \{\phi, t\}.$$

In code, this is exactly post 4’s firstFundamentalForm with the dot product replaced by weightedDot:

function firstFundamentalForm(phi, t, L, m_e, m_v) {
  const M = ambientCouplingMatrix(4, m_e, m_v);
  const eps = 1e-4;
  const pPh = embedHinges(phi + eps, t, L);
  const pPl = embedHinges(phi - eps, t, L);
  const pTh = embedHinges(phi, t + eps, L);
  const pTl = embedHinges(phi, t - eps, L);

  // Central differences: ∂_φ p, ∂_t p as lists of 3 hinges
  const dPp = pPh.map((p, k) => [
    (p[0] - pPl[k][0]) / (2*eps),
    (p[1] - pPl[k][1]) / (2*eps),
  ]);
  const dTp = pTh.map((p, k) => [
    (p[0] - pTl[k][0]) / (2*eps),
    (p[1] - pTl[k][1]) / (2*eps),
  ]);

  return {
    E: weightedDot(M, dPp, dPp),
    F: weightedDot(M, dPp, dTp),
    G: weightedDot(M, dTp, dTp),
  };
}

The second fundamental form

The second fundamental form $II_{ab}$ is the normal component of $\partial_a\partial_b\mathbf{p}$ in the ambient. “Normal” is defined relative to the ambient inner product — orthogonal complement of the tangent space — so the normal projection step also uses weightedDot.

The tangential component of an ambient vector $\mathbf{v}$ has the form $c^\phi\partial_\phi\mathbf{p} + c^t\partial_t\mathbf{p}$, with the coefficients $(c^\phi, c^t)$ satisfying

$$\begin{pmatrix} E & F \\ F & G \end{pmatrix}\begin{pmatrix} c^\phi \\ c^t \end{pmatrix} \;=\; \begin{pmatrix} \langle\mathbf{v}, \partial_\phi\mathbf{p}\rangle_g \\ \langle\mathbf{v}, \partial_t\mathbf{p}\rangle_g \end{pmatrix}$$

— same equation as post 4, but with all inner products $g$-weighted. Subtracting the tangential piece from $\mathbf{v}$ gives its normal component:

function secondFundamentalForm(phi, t, L, m_e, m_v) {
  const M = ambientCouplingMatrix(4, m_e, m_v);
  const eps = 1e-4, e2 = eps*eps;
  const p   = embedHinges(phi,       t,       L);
  const pPh = embedHinges(phi + eps, t,       L);
  const pPl = embedHinges(phi - eps, t,       L);
  const pTh = embedHinges(phi,       t + eps, L);
  const pTl = embedHinges(phi,       t - eps, L);
  const pPP = embedHinges(phi + eps, t + eps, L);
  const pPN = embedHinges(phi + eps, t - eps, L);
  const pNP = embedHinges(phi - eps, t + eps, L);
  const pNN = embedHinges(phi - eps, t - eps, L);

  // Hinge-list central differences (each entry a 2-vector in R^2)
  const fd = (a, b, h) => a.map((p, k) => [
    (p[0] - b[k][0]) / h,
    (p[1] - b[k][1]) / h,
  ]);
  const fd2 = (a, c, b, h) => a.map((p, k) => [
    (p[0] - 2*c[k][0] + b[k][0]) / h,
    (p[1] - 2*c[k][1] + b[k][1]) / h,
  ]);

  const dPp   = fd(pPh, pPl, 2*eps);
  const dTp   = fd(pTh, pTl, 2*eps);
  const dPPp  = fd2(pPh, p, pPl, e2);
  const dTTp  = fd2(pTh, p, pTl, e2);
  const dPTp  = p.map((_, k) => [
    (pPP[k][0] - pPN[k][0] - pNP[k][0] + pNN[k][0]) / (4*e2),
    (pPP[k][1] - pPN[k][1] - pNP[k][1] + pNN[k][1]) / (4*e2),
  ]);

  // 2x2 system from <q_a, q_b>_g
  const E = weightedDot(M, dPp, dPp);
  const F = weightedDot(M, dPp, dTp);
  const G = weightedDot(M, dTp, dTp);
  const D = E*G - F*F;

  // Subtract the g-tangent part from v to get its g-normal component
  const projectNormal = v => {
    const bP = weightedDot(M, v, dPp);
    const bT = weightedDot(M, v, dTp);
    const cP = (G*bP - F*bT) / D;
    const cT = (E*bT - F*bP) / D;
    return v.map((vk, k) => [
      vk[0] - cP*dPp[k][0] - cT*dTp[k][0],
      vk[1] - cP*dPp[k][1] - cT*dTp[k][1],
    ]);
  };

  return {
    II_pp: projectNormal(dPPp),
    II_pt: projectNormal(dPTp),
    II_tt: projectNormal(dTTp),
  };
}

Note the asymmetry between the form of the metric components and the form of the normal vectors. The metric components $E, F, G$ are scalars — direct outputs of weightedDot. The normal vectors $II_{ab}$ are still represented as lists of three hinges in $(\mathbb{R}^2)^3$ — ambient vectors — and their lengths and inner products will need to be computed via weightedDot whenever they appear downstream.

Curvature

With the fundamental forms in hand, we can now compute the curvature scalars from post 4 — the Gaussian curvature $K$ (intrinsic, governing how geodesics diverge) and the total bending $|II|^2$ (extrinsic, the magnitude of the constraint force). Both formulas are unchanged from post 4; the only difference is that all inner products on ambient vectors now go through weightedDot.

Gaussian curvature

The Gauss equation expresses $K$ in terms of both fundamental forms:

$$K \;=\; \frac{\langle II_{\phi\phi}, II_{tt}\rangle_g \;-\; |II_{\phi t}|^2_g}{EG - F^2}.$$

Same expression as post 4, with the inner products now $g$-weighted:

function gaussianCurvature(phi, t, L, m_e, m_v) {
  const M = ambientCouplingMatrix(4, m_e, m_v);
  const { E, F, G }            = firstFundamentalForm(phi, t, L, m_e, m_v);
  const { II_pp, II_pt, II_tt } = secondFundamentalForm(phi, t, L, m_e, m_v);

  const numerator = weightedDot(M, II_pp, II_tt) - weightedDot(M, II_pt, II_pt);
  return numerator / (E*G - F*F);
}

(As in post 4, this duplicates work between the two helper calls — both functions internally compute the first-derivative stencils. A merged version would run the stencil once and return both forms together. We’re keeping them separate here for clarity.)

Here’s a small simulation computing this curvature on the configuration space, in terms of sliders for $L$, as well as $m_e$ and $m_v$.

Total bending

Total bending is the squared norm of $II$ as a bilinear form. In an orthonormal tangent frame $\mathbf{e}_1, \mathbf{e}_2$,

$$|II|^2 \;=\; |II_{11}|^2_g \;+\; 2\,|II_{12}|^2_g \;+\; |II_{22}|^2_g.$$

Gram–Schmidt on the coordinate tangents gives the orthonormal frame, and changing the second fundamental form into this frame uses the same coefficients as post 4:

$$\begin{aligned} II_{11} \;&=\; \tfrac{1}{E}\,II_{\phi\phi}, \\ II_{12} \;&=\; -\tfrac{F}{E\sqrt{D}}\,II_{\phi\phi} \;+\; \tfrac{1}{\sqrt{D}}\,II_{\phi t}, \\ II_{22} \;&=\; \tfrac{F^2}{ED}\,II_{\phi\phi} \;-\; \tfrac{2F}{D}\,II_{\phi t} \;+\; \tfrac{E}{D}\,II_{tt}. \end{aligned}$$

These are vectors in the ambient — linear combinations of the $II_{ab}$‘s with scalar coefficients depending on $E, F, G$. Their squared norms are again $g$-weighted:

function totalBending(phi, t, L, m_e, m_v) {
  const M = ambientCouplingMatrix(4, m_e, m_v);
  const { E, F, G }             = firstFundamentalForm(phi, t, L, m_e, m_v);
  const { II_pp, II_pt, II_tt } = secondFundamentalForm(phi, t, L, m_e, m_v);
  const D = E*G - F*F, sqrtD = Math.sqrt(D);

  // Linear combination of hinge-list vectors with scalar coefficients
  const combine = (cs, vs) => vs[0].map((_, k) => [
    cs.reduce((s, c, i) => s + c * vs[i][k][0], 0),
    cs.reduce((s, c, i) => s + c * vs[i][k][1], 0),
  ]);

  const II_11 = combine([1/E],                       [II_pp]);
  const II_12 = combine([-F/(E*sqrtD),  1/sqrtD],     [II_pp, II_pt]);
  const II_22 = combine([F*F/(E*D), -2*F/D, E/D],     [II_pp, II_pt, II_tt]);

  return weightedDot(M, II_11, II_11)
       + 2*weightedDot(M, II_12, II_12)
       + weightedDot(M, II_22, II_22);
}

Here’s the corresponding demo for the extrinsic curvature, again with sliders for both distance between the endpoints and the masses.

Dynamics

The force-free dynamics of the linkage are geodesics of the kinetic-energy metric on $\mathcal{M}^4_L$. Post 3 derived these two ways: the standard route via Christoffel symbols of the pulled-back metric, and a more elementary embedded route exploiting the fact that the geodesic acceleration must be normal to the surface in the ambient. The two were used as cross-checks on each other.

For our mass-weighted setting we’ll go straight to the embedded approach, which fits naturally with the machinery we just built. The condition is unchanged from post 3, except that “normal” is now defined relative to the ambient inner product:

A curve $(\phi(s), t(s))$ on $\mathcal{M}^4_L$ is a geodesic iff its acceleration $\ddot{\mathbf{p}}$ in the ambient is $g$-normal to the surface — that is, $\langle\ddot{\mathbf{p}}, \partial_\phi\mathbf{p}\rangle_g = \langle\ddot{\mathbf{p}}, \partial_t\mathbf{p}\rangle_g = 0$.

Following post 3 §5: differentiate $\mathbf{p}(\phi(s), t(s))$ twice in $s$ to get

$$\ddot{\mathbf{p}} \;=\; \mathbf{p}_\phi\,\ddot\phi \;+\; \mathbf{p}_t\,\ddot t \;+\; \mathbf{p}_{\phi\phi}\,\dot\phi^2 \;+\; 2\,\mathbf{p}_{\phi t}\,\dot\phi\dot t \;+\; \mathbf{p}_{tt}\,\dot t^2,$$

then enforce the two $g$-normality conditions. They give a 2×2 linear system for the unknown second derivatives:

$$\begin{pmatrix} \langle\mathbf{p}_\phi, \mathbf{p}_\phi\rangle_g & \langle\mathbf{p}_\phi, \mathbf{p}_t\rangle_g \\ \langle\mathbf{p}_t, \mathbf{p}_\phi\rangle_g & \langle\mathbf{p}_t, \mathbf{p}_t\rangle_g \end{pmatrix} \begin{pmatrix} \ddot\phi \\ \ddot t \end{pmatrix} \;=\; -\begin{pmatrix} \langle\mathbf{p}_\phi, \mathbf{p}_{vv}\rangle_g \\ \langle\mathbf{p}_t, \mathbf{p}_{vv}\rangle_g \end{pmatrix},$$

where $\mathbf{p}_{vv} := \mathbf{p}_{\phi\phi}\dot\phi^2 + 2\mathbf{p}_{\phi t}\dot\phi\dot t + \mathbf{p}_{tt}\dot t^2$ is the velocity-quadratic part of $\ddot{\mathbf{p}}$. The 2×2 coefficient matrix is exactly $\begin{pmatrix} E & F \\ F & G \end{pmatrix}$ — the first fundamental form, computed inline from the same partials.

The code re-uses the 9-point stencil from secondFundamentalForm. We keep the math entirely in the hinge-list representation:

function geodesicVectorField([phi, t, phidot, tdot], L, m_e, m_v) {
  const M   = ambientCouplingMatrix(4, m_e, m_v);
  const eps = 1e-4, e2 = eps*eps;

  // 9-point stencil
  const p   = embedHinges(phi,       t,       L);
  const pPh = embedHinges(phi + eps, t,       L);
  const pPl = embedHinges(phi - eps, t,       L);
  const pTh = embedHinges(phi,       t + eps, L);
  const pTl = embedHinges(phi,       t - eps, L);
  const pPP = embedHinges(phi + eps, t + eps, L);
  const pPN = embedHinges(phi + eps, t - eps, L);
  const pNP = embedHinges(phi - eps, t + eps, L);
  const pNN = embedHinges(phi - eps, t - eps, L);

  // Hinge-list arithmetic helpers
  const fd  = (a, b, h)    => a.map((p, k) => [(p[0]-b[k][0])/h, (p[1]-b[k][1])/h]);
  const fd2 = (a, c, b, h) => a.map((p, k) => [(p[0]-2*c[k][0]+b[k][0])/h, (p[1]-2*c[k][1]+b[k][1])/h]);
  const lc  = (cs, vs)     => vs[0].map((_, k) => [
    cs.reduce((s, c, i) => s + c*vs[i][k][0], 0),
    cs.reduce((s, c, i) => s + c*vs[i][k][1], 0),
  ]);

  // First and second partials
  const p_phi    = fd(pPh, pPl, 2*eps);
  const p_t      = fd(pTh, pTl, 2*eps);
  const p_phiphi = fd2(pPh, p, pPl, e2);
  const p_tt     = fd2(pTh, p, pTl, e2);
  const p_phit   = p.map((_, k) => [
    (pPP[k][0] - pPN[k][0] - pNP[k][0] + pNN[k][0]) / (4*e2),
    (pPP[k][1] - pPN[k][1] - pNP[k][1] + pNN[k][1]) / (4*e2),
  ]);

  // Velocity-quadratic part of acceleration
  const p_vv = lc(
    [phidot*phidot, 2*phidot*tdot, tdot*tdot],
    [p_phiphi, p_phit, p_tt],
  );

  // 2x2 system [<p_a, p_b>_g] (phiddot, tddot)^T = -[<p_a, p_vv>_g]
  const E = weightedDot(M, p_phi, p_phi);
  const F = weightedDot(M, p_phi, p_t);
  const G = weightedDot(M, p_t,   p_t);
  const D = E*G - F*F;
  const bP = -weightedDot(M, p_phi, p_vv);
  const bT = -weightedDot(M, p_t,   p_vv);
  const phiddot = (G*bP - F*bT) / D;
  const tddot   = (E*bT - F*bP) / D;

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

Here’s the updated demo!