Multibody Dynamics: Drift, Control, and the Geometry of Physical Law

Dieter Olson

Multibody Dynamics: Drift, Control, and the Geometry of Physical Law

theory-core

Derive modeled motion through energy, body forces, twists, and task coordinates, applying those connections to investigate golf and humanoid dynamics.

Author

Dieter Olson

Published

March 6, 2026

1 One Mechanism, Several Questions

A golfer’s pelvis can slow while the clubhead gains speed. That observation raises several questions: what forces act between segments, where does energy go, how do current torques affect acceleration, and what would a different torque history do to impact? These questions are related, but they are not interchangeable. A segment-speed graph cannot answer all of them.

Multibody dynamics gives us several ways to examine the same modeled mechanism. An energy derivation exposes the inertia and generalized efforts. Body force balances recover joint reactions. Spatial twists organize three-dimensional motion and wrenches. Task coordinates reveal what the clubhead or foot can do at its current posture. Geometry clarifies which conclusions survive a change of coordinates.

The value of using several formulations is that their predictions can be compared under the same model, state, loads, constraints and observation convention. Agreement is a strong implementation check. It does not establish that the chosen model accurately represents a golfer. That further step requires independent measurements.

1.1 The Shared Drift-Control Structure

For a declared input vector $u$, many models have the form

\[ \dot x=f_0(x,t)+G(x,t)u. \tag{1}\]

At a fixed state and time, $f_0$ is the result with those input channels set to zero; $Gu$ is the incremental state derivative. This is control affinity. It is linear in input at that fixed state, even when the state dependence and the resulting trajectories are nonlinear.

The river analogy can help if used carefully: paddling now changes where the kayaker goes and therefore which current they encounter later. Similarly, a torque changes future configuration and velocity, hence future drift and input effectiveness. A pointwise sum of derivatives is not a sum of independently generated finite trajectories.

The zero reference must be declared. Zero commanded motor voltage, zero joint torque and zero neural excitation are different interventions. Existing activation, passive tissue forces, shaft deformation and ground reactions can persist when a chosen command is zero. Nor does the split divide “physics” from “intent”: actuator forces are physical, and many mechanisms can produce the same net joint torque.

1.2 Scope of the Comparison

We first use independent local coordinates $q$, velocity $v=\dot q$, a symmetric positive-definite inertia $M(q)$, and generalized effort $\tau=B(q)u$. We then add regular stationary ideal bilateral constraints. The worked mechanism has two massless rigid rods with positive point masses at their distal ends, a fixed pivot, gravity and ideal planar revolute joints.

Damping or elastic coordinates can be included without necessarily losing input affinity. Actuator dynamics require an expanded state; affinity must then be checked with respect to the actual chosen input. Contact switching, impacts, friction and hysteresis need additional rules. The hybrid reference develops the event and contact distinctions.

For three-dimensional floating bodies, a quaternion configuration may have more components than its generalized velocity. Then $\dot q=N(q)v$ and the equations use $\dot v$, not an unqualified identification of $v$ with $\dot q$. A humanoid also needs its base motion and a declared support mode; a joint-angle vector alone may not be a complete state.

2 Configuration, Effort, and Inertia

2.1 State Is Sufficient Only for the Declared Model

Configuration specifies the arrangement of bodies. Velocity specifies how that arrangement is changing. Under a well-posed second-order rigid model, $x=(q,v)$, the specified inputs and parameters determine the local evolution. A muscle activation or viscoelastic memory variable cannot be omitted if it affects that evolution independently of $q,v$.

State summarizes the information a model retains; it does not uniquely reconstruct the historical causes of that state. Different torque histories can arrive at the same configuration and velocity.

2.2 Generalized Force Is Defined by Work

A generalized effort $Q$ satisfies

\[ \delta W=Q^{\mathsf T}\delta q,\qquad P=Q^{\mathsf T}v. \]

An angle coordinate pairs with a moment; a translation pairs with a force. For a point load $F$ at $p(q)$, $\delta p=J\delta q$ gives $Q=J^{\mathsf T}F$. This work relation, rather than “whatever makes this coordinate accelerate,” defines generalized force. Coupling means one effort can affect several accelerations.

The kinetic energy is

\[ T(q,v)=\tfrac12v^{\mathsf T}M(q)v. \]

$M$ maps an acceleration increment to an inertial generalized-force increment. Its inverse maps a force increment to acceleration when other state-dependent terms and constraints are treated consistently. A surface that adds rolling resistance changes a dissipative load, not the cart’s inertia merely by being rougher.

Positive definiteness means every nonzero admissible independent velocity has positive kinetic energy. It requires the corresponding inertial assumptions. Redundant coordinates, massless degrees of freedom or a singular chart need separate treatment.

2.3 The Manipulator Equation and Its Signs

Write

\[ M(q)\ddot q+b(q,v)=B(q)u+Q_{\rm ext},\qquad b=c(q,v)+g_V(q),\quad g_V=\frac{\partial V}{\partial q}. \tag{2}\]

Here $c=C(q,v)v$ contains the velocity-product terms of this coordinate description. The physical generalized gravitational force is $-g_V$; $g_V$ has the positive sign because it is on the left. $Q_{\rm ext}$ denotes additional loads, so gravity is not counted twice.

The partition of velocity products into “Coriolis” and “centripetal” terms depends on coordinates and on the nonunique choice of $C$. Centripetal acceleration points inward. Water following a circular path needs an inward net force; an outward centrifugal term belongs to a rotating-frame description. Neither term is a new source of energy.

3 A Double Pendulum Derived From Cartesian Motion

3.1 Geometry and Velocity

Let $q=(\theta_1,\theta_2)$, with $\theta_1$ measured counterclockwise from downward vertical and $\theta_2$ the second rod’s angle relative to the first. Positive vertical position points upward. Define

\[ r_1=l_1\begin{bmatrix}\sin\theta_1\\-\cos\theta_1\end{bmatrix}, \qquad r_2=l_2\begin{bmatrix}\sin(\theta_1+\theta_2)\\-\cos(\theta_1+\theta_2)\end{bmatrix}. \tag{3}\]

The mass positions are $p_1=r_1$, $p_2=r_1+r_2$. The absolute angular velocities are $\omega_1=v_1$, $\omega_2=v_1+v_2$. Confusing the second relative rate with its absolute rate changes both kinetic energy and joint power.

\[ J_1=\begin{bmatrix}l_1\cos\theta_1&0\\l_1\sin\theta_1&0\end{bmatrix}, \quad J_2=\begin{bmatrix} l_1\cos\theta_1+l_2\cos(\theta_1+\theta_2)&l_2\cos(\theta_1+\theta_2)\\ l_1\sin\theta_1+l_2\sin(\theta_1+\theta_2)&l_2\sin(\theta_1+\theta_2) \end{bmatrix}. \tag{4}\]

Then $\dot p_i=J_i v$. Differentiating once more gives

\[ \ddot p_1=J_1\ddot q-\omega_1^2r_1,\qquad \ddot p_2=J_2\ddot q-\omega_1^2r_1-\omega_2^2r_2. \tag{5}\]

The last terms are the Cartesian curvature or bias accelerations. They remain even if both relative joint accelerations are zero.

3.2 Inertia, Bias, and Potential

Start from $T=\tfrac12m_1\|\dot p_1\|^2+\tfrac12m_2\|\dot p_2\|^2$. Substitution gives $M=m_1J_1^{\mathsf T}J_1+m_2J_2^{\mathsf T}J_2$, hence

\[ M= \begin{bmatrix} (m_1+m_2)l_1^2+m_2l_2^2+2m_2l_1l_2\cos\theta_2& m_2l_2^2+m_2l_1l_2\cos\theta_2\\ m_2l_2^2+m_2l_1l_2\cos\theta_2&m_2l_2^2 \end{bmatrix}. \tag{6}\]

Its determinant is $m_2l_1^2l_2^2(m_1+m_2\sin^2\theta_2)>0$. Thus straight or folded alignment does not singularize this joint inertia when all masses and lengths are positive. The off-diagonal entry can vanish or change sign for some length ratios; “weakest coupling when folded” is not a general statement about all inertial responses.

The Euler–Lagrange operation $d(\partial T/\partial v)/dt-\partial T/\partial q$ yields

\[ c=\beta \begin{bmatrix}-2v_1v_2-v_2^2\\v_1^2\end{bmatrix}, \qquad \beta=m_2l_1l_2\sin\theta_2. \tag{7}\]

At $\sin\theta_2=0$, this vector vanishes at that instant, but $M$ is generally still coupled. It does not follow that the joints move independently.

With $V=m_1g(p_1)_y+m_2g(p_2)_y$,

\[ g_V=g\begin{bmatrix} (m_1+m_2)l_1\sin\theta_1+m_2l_2\sin(\theta_1+\theta_2)\\ m_2l_2\sin(\theta_1+\theta_2) \end{bmatrix}. \tag{8}\]

These expressions can also be recovered from Cartesian inertial loads, as the executable example below checks. MIT’s multibody notes give the standard manipulator and constraint framework; signs must be translated consistently when gravity is placed on the opposite side.

4 Pointwise Attribution and Future Control

4.1 The Instantaneous Split

With retained additional loads independent of $u$,

\[ a_0=M^{-1}(Q_{\rm ext}-b),\qquad \Delta a=M^{-1}Bu,\qquad \ddot q=a_0+\Delta a. \tag{9}\]

The pointwise Zero-Torque Counterfactual (ZTCF) sample uses the copied state and the declared zeroed channels. A forward ZTCF trajectory repeatedly solves the zero-input dynamics while evolving its own state. It is a different object from a series of pointwise samples evaluated along the original controlled motion.

High speed does not necessarily imply large $a_0$. A free particle in Cartesian coordinates can move arbitrarily fast with zero drift acceleration. In this pendulum, aligned rods have $c=0$ regardless of the rates. Posture, direction, gravity and cancellation all matter.

A norm ratio of drift and input increments also needs units, coordinate scales and a zero-denominator convention. It does not determine energy contribution, neural effort or whether a finite torque intervention improves a shot.

4.2 Linearize the Trajectory Actually Being Used

For a nominal solution $\bar x(t),\bar u(t)$, the first variation is

\[ \delta\dot x=A(t)\delta x+G(\bar x,t)\delta u,\qquad A=\left.\frac{\partial f_0}{\partial x}\right|_{\bar x} +\sum_j\bar u_j\left.\frac{\partial G_j}{\partial x}\right|_{\bar x}. \tag{10}\]

The second term is needed even though the dynamics are affine in input. If a reference path is not a solution, the error equation also contains $f_0(\bar x,t)+G(\bar x,t)\bar u-\dot{\bar x}$. Sampling only the zero-input Jacobian on a controlled path answers a different question.

Frozen eigenvalues of $A(t)$ do not establish stability of the time-varying system. For example, alternating for one time unit each between $A_1=\begin{bmatrix}-1&4\\0&-1\end{bmatrix}$ and $A_2=A_1^{\mathsf T}$ gives two matrices whose eigenvalues are all $-1$, but the two-step transition $e^{A_2}e^{A_1}$ has an eigenvalue greater than two. Stability concerns the transition over time, or a suitable Lyapunov/contraction bound.

4.3 Why a Torque Change Can Matter Elsewhere

A shoulder torque can affect clubhead acceleration immediately through coupled inertia. Over time it also changes posture, velocity, shaft state and when contact occurs. These changes alter subsequent control effectiveness. To investigate a proposed “release,” specify its torque or activation history, copied state, support mode and output. Then compare the resulting trajectories at a declared clock time or at each trajectory’s own impact event.

5 Constraints Redirect Acceleration and Change Reactions

5.1 Stationary Ideal Bilateral Constraints

Let $\phi(q)=0$ have $k$ independent rows locally, $J_c=\partial\phi/\partial q$, and $\kappa=\dot J_c v$. An admissible state satisfies $\phi=0$, $J_cv=0$. Its acceleration satisfies

\[ M\ddot q+b=Bu+Q_{\rm ext}+J_c^{\mathsf T}\lambda,\qquad J_c\ddot q=-\kappa. \tag{11}\]

The allowed accelerations form an affine subspace $\mathcal H=\{a:J_ca=-\kappa\}$ of dimension $n-k$, a hyperplane only when the rank is one. They do not generally lie in the velocity tangent space: a train following a curved track needs normal acceleration.

Define $W=J_cM^{-1}J_c^{\mathsf T}$, which is positive definite under these assumptions. If $a_{\rm free}=M^{-1}(Bu+Q_{\rm ext}-b)$, substitution gives

\[ \lambda=-W^{-1}(J_ca_{\rm free}+\kappa),\qquad \ddot q=a_{\rm free}-M^{-1}J_c^{\mathsf T}W^{-1}(J_ca_{\rm free}+\kappa). \tag{12}\]

This is equivalently the solution of the block system

\[ \begin{bmatrix}M&-J_c^{\mathsf T}\\J_c&0\end{bmatrix} \begin{bmatrix}\ddot q\\\lambda\end{bmatrix} =\begin{bmatrix}Bu+Q_{\rm ext}-b\\-\kappa\end{bmatrix}. \tag{13}\]

Use linear solves rather than explicitly forming inverses in numerical code. Redundant constraints need a rank-aware treatment; nonunique multipliers need not imply nonunique physical acceleration.

5.2 The Constrained Input Map

Let $P=I-M^{-1}J_c^{\mathsf T}W^{-1}J_c$. This projects acceleration increments onto $\ker J_c$ orthogonally in the $M$ metric: $P^2=P$, $P^{\mathsf T}M=MP$. Then

\[ a_0^c=P M^{-1}(Q_{\rm ext}-b)-M^{-1}J_c^{\mathsf T}W^{-1}\kappa, \qquad B_c=P M^{-1}B, \tag{14}\]

\[ \ddot q=a_0^c+B_cu,\qquad \lambda(u)-\lambda(0)=-W^{-1}J_cM^{-1}Bu. \tag{15}\]

The input changes the reaction as well as the acceleration. $B_cu$ lies in the homogeneous tangent directions of the affine acceleration set; the curvature offset belongs to the anchor. It is incorrect to project the entire acceleration onto a zero-offset tangent plane.

For a manufactured particle of mass $2$ kg at $r=(2,0)$ m with velocity $(0,3)$ m/s, constrained by $\phi=(r^{\mathsf T}r-4)/2=0$, $J_c=(2,0)$ and $\kappa=9$. With no other load, $a_0^c=(-4.5,0)$ m/s². An applied force $(4,6)$ N changes acceleration to $(-4.5,3)$ m/s². The reaction changes from $(-9,0)$ to $(-13,0)$ N: it cancels the new radial input while retaining the inward curvature acceleration.

5.3 What the Zero-Work Statement Requires

For these stationary ideal constraints,

\[ P_c=(J_c^{\mathsf T}\lambda)^{\mathsf T}v =\lambda^{\mathsf T}J_cv=0. \tag{16}\]

Zero total power does not mean zero force or zero influence on motion. For a moving constraint $\phi(q,t)=0$, $J_cv=-\phi_t$ and $P_c=-\lambda^{\mathsf T}\phi_t$, which can be nonzero. Stationary ideal Pfaffian constraints $A(q)v=0$ likewise have zero reaction power, although their velocity distribution need not integrate to a configuration submanifold.

A real foot-ground constraint is unilateral and friction limited. A computed reaction requiring tension from the floor or excessive tangential force invalidates that assumed contact mode. Removing torque may cause release or slip; then the same fixed-mode affine formula no longer describes the complete intervention. These conditions matter to humanoid support and to a golfer’s use of the ground.

For the pendulum endpoint constrained to a circle centered at $c_0$, choose $\phi=\|p_2-c_0\|^2-R^2$. Then $J_c=2(p_2-c_0)^{\mathsf T}J_2$. Differentiating the full expression, including $J_2$, gives the required curvature term. The same block solve applies wherever the constraint row is regular and the copied state is admissible.

6 Body Forces and Spatial Motion

6.1 Newton–Euler Reveals the Loads

For body $i$, take all forces acting on that body, expressed in a common frame. Let $r_{ik}$ run from its center of mass to the point where force $F_{ik}$ acts, and $N_{ik}$ denote applied pure moments. In a body frame fixed to the center of mass, the balances are

\[ m_i a_{C,i}=\sum_kF_{ik},\qquad I_{C,i}\dot\omega_i+\omega_i\times(I_{C,i}\omega_i) =\sum_k\bigl(r_{ik}\times F_{ik}+N_{ik}\bigr). \tag{17}\]

The first equation’s acceleration is the physical center-of-mass acceleration expressed in that frame, including the appropriate transport terms if obtained by differentiating body components. The moment equation includes the force lever arms. One cannot simply add joint moments about different origins.

In the standard planar reduction, the out-of-plane axis is principal and the scalar rotational equation uses its moment of inertia. Parallel angular velocities alone would not prove $\omega\times I\omega=0$ for an arbitrary three-dimensional inertia tensor.

The recursive Newton–Euler inverse-dynamics algorithm takes configuration, velocity and prescribed acceleration, propagates motion outward, and computes required loads and joint efforts inward. A tree recursion has linear cost in the number of links for this problem.

To obtain zero-torque motion, solve a forward-dynamics problem. One approach obtains the bias and mass-matrix columns from inverse-dynamics evaluations, then solves $Ma=-b+Q_{\rm ext}$. An articulated-body algorithm can instead solve tree forward dynamics directly. The forward-dynamics tutorial explains the former construction. Setting an inverse-dynamics output torque to zero without solving for acceleration is not an algorithm.

Once acceleration is known, body balance recovers reactions. Additional constraints require the corresponding coupled solve. Thus Euler–Lagrange and Newton–Euler are not restricted to disjoint physical questions; their different representations make particular calculations easier to inspect.

6.2 Spatial Twists Need an Origin Convention

Let $g_{sb}=[R,p;0,1]$ map body coordinates to a fixed space frame. With angular-first ordering, the space twist satisfies $\widehat{\mathcal V}_s=\dot g_{sb}g_{sb}^{-1}$, so

\[ \mathcal V_s=\begin{bmatrix}\omega_s\\v_s\end{bmatrix}, \qquad v_s=\dot p-\omega_s\times p,\qquad \dot p=v_s+\omega_s\times p. \tag{18}\]

The linear block $v_s$ is generally not the velocity of the body-frame origin. The derivative of the space twist is likewise not just ordinary endpoint acceleration:

\[ \ddot p=\dot v_s+\dot\omega_s\times p+\omega_s\times\dot p. \]

The twist tutorial supplies the adjoint and matrix conventions. A product of exponentials (PoE) is a kinematic representation; inertia and loads are still needed for dynamics.

For our downward-vertical home posture, use

\[ S_1=(0,0,1,0,0,0)^{\mathsf T},\qquad S_2=(0,0,1,-l_1,0,0)^{\mathsf T}, \tag{19}\]

and home translation $p(0)=(0,-l_1-l_2,0)$. Then $g_{sb}(q)=e^{\widehat S_1\theta_1}e^{\widehat S_2\theta_2}g_{sb}(0)$ and $J_s=[S_1,\operatorname{Ad}_{e^{\widehat S_1\theta_1}}S_2]$. Recover $\dot p$ using the cross-position term before comparing with $J_2v$. A screw through $(l_1,0,0)$ would describe a different home geometry.

For dynamics in a body center-of-mass frame, write $\mathcal V_b=(\omega_b,v_{C,b})$, $\mathcal I_b=\operatorname{diag}(I_C,mI_3)$, and

\[ \operatorname{ad}_{\mathcal V_b} =\begin{bmatrix}[\omega_b]_\times&0\\{}[v_{C,b}]_\times&[\omega_b]_\times\end{bmatrix}, \qquad \mathcal W_b=\mathcal I_b\dot{\mathcal V}_b -\operatorname{ad}_{\mathcal V_b}^{\mathsf T}\mathcal I_b\mathcal V_b. \tag{20}\]

This declares the sign rather than relying on an ambiguous coadjoint-star notation. The rigid-body dynamics tutorial uses this convention. Wrenches at other points must include their translated moments.

If $\mathcal V_s=\operatorname{Ad}_g\mathcal V_b$, then $\mathcal W_s=\operatorname{Ad}_g^{-\mathsf T}\mathcal W_b$, and $\mathcal W_s^{\mathsf T}\mathcal V_s=\mathcal W_b^{\mathsf T}\mathcal V_b$. This is invariance under consistent representation of the same wrench and motion. It does not say that mechanical power is unchanged for an observer moving with a different velocity.

7 What the Clubhead or Foot Can Do

7.1 Task Dynamics as an Induced Description

Let $y=\psi(q)$, $J=\partial\psi/\partial q$. Direct differentiation gives

\[ \ddot y=J M^{-1}(Bu+Q_{\rm ext}+J_c^{\mathsf T}\lambda-b)+\dot Jv. \tag{21}\]

This acceleration equation remains useful when a task inertia cannot be inverted. Where $J$ has full row rank, define mobility $H=JM^{-1}J^{\mathsf T}$ and task inertia $\Lambda=H^{-1}$. Then

\[ \Lambda\ddot y+\mu+p_V =F_{\rm ctrl}^{\rm eff}+F_{\rm ext}^{\rm eff}+F_c^{\rm eff}, \tag{22}\]

where $\mu=\Lambda(JM^{-1}c-\dot Jv)$, $p_V=\Lambda JM^{-1}g_V$, and $F^{\rm eff}=\Lambda JM^{-1}Q$ for each corresponding generalized load. These quantities depend on the full joint state; a redundant robot’s task position alone does not determine them. The task-dynamics tutorial illustrates the invertible-Jacobian case.

A specified generalized torque has a unique instantaneous task effect. Redundancy concerns the opposite allocation problem: many torques may produce the same task effect. In a fully actuated unconstrained model,

\[ \tau=J^{\mathsf T}F^{\rm eff}+\tau_{\rm null}, \qquad JM^{-1}\tau_{\rm null}=0. \]

This dynamically weighted null space differs from an arbitrary Euclidean projection. Its instantaneous lack of task acceleration does not guarantee no future task effect: it can change posture and later dynamics. Underactuation, effort bounds and support constraints further restrict which allocations are possible.

7.2 A Task Singularity Is Not an Inertia Singularity

For the pendulum endpoint, $\det J_2=l_1l_2\sin\theta_2$. At straight or folded alignment $J_2$ loses rank, $H$ is singular, and $\Lambda=H^{-1}$ is undefined. Describing $\Lambda$ as a finite singular matrix is wrong. A diverging effective inertia may describe a limit approaching that posture, but it does not supply a finite force that creates an unavailable instantaneous increment.

At the straight-down posture, the vertical row of $J_2$ is zero. Joint acceleration therefore supplies no vertical task-acceleration increment at that instant. Nevertheless, $(\dot J_2v)_y=l_1v_1^2+l_2(v_1+v_2)^2$ can be positive: the endpoint has inward curvature acceleration. Both joint degrees of freedom and the positive-definite $M$ remain.

For a supported mechanism, use the constrained map $JB_c$ for input authority. A clubhead Jacobian can be well conditioned while input authority is restricted by contact or missing actuators. Conversely, inability to move in one direction at one instant does not settle finite-time reachability.

8 Geometry Without Changing the Physical Question

8.1 The Kinetic Metric

For unconstrained pendulum angles with full revolutions allowed, configuration space is $Q=S^1\times S^1$. Joint limits restrict that domain. The inertia defines a Riemannian kinetic metric when it is positive definite: $\langle v,w\rangle_M=v^{\mathsf T}Mw$. It measures squared speed through kinetic energy, not the energy consumed by a displacement without a duration or force history.

With the Levi-Civita connection of this metric, the intrinsic equation is

\[ \nabla_vv=-\operatorname{grad}_M V +(Bu+Q_{\rm ext}+Q_c)^\sharp. \tag{23}\]

Here $\sharp$ raises a force covector using $M^{-1}$, and $\operatorname{grad}_M V=M^{-1}g_V$ in coordinates. The coordinate expansion is $\ddot q^i+\Gamma^i_{jk}v^jv^k$, with $M\Gamma(v,v)=c$. Nonzero Christoffel symbols can arise in flat space, as in polar coordinates; they are not by themselves a curvature measurement.

8.2 The Whole State Field Transforms Intrinsically

In a local chart the zero-input second-order field has components

\[ f_0(q,v)= \begin{bmatrix} v\\-\Gamma(v,v)-\operatorname{grad}_M V+Q_{\rm ext}^{\sharp} \end{bmatrix}. \tag{24}\]

This displays an intrinsic field in coordinates; the individual lower block is not an ordinary vector under arbitrary configuration changes. If $z=\chi(q)$, then

\[ \dot z=D\chi\,v,\qquad \ddot z=D\chi\,\ddot q+D^2\chi[v,v]. \]

For a free Cartesian particle with coordinate $q>0$, zero acceleration becomes $\ddot z=2v^2$ in the chart $z=q^2$. The physical motion has not acquired a force. Comparing two acceleration solutions at the same $q,v$, however, cancels this Hessian term: $\Delta\ddot z=D\chi\,\Delta\ddot q$. This is why the affine acceleration-space interpretation and its tangent increments require care.

The switch from relative to absolute pendulum angles is linear: $\alpha=(\theta_1,\theta_1+\theta_2)$. If $v=A^{-1}\dot\alpha$, the metric becomes $M_\alpha=A^{-\mathsf T}MA^{-1}$ and effort becomes $Q_\alpha=A^{-\mathsf T}Q$. The paired work and kinetic energy agree.

Zero input describes geodesic motion only when all nonconstraint forces vanish, including potential gradients, damping and additional loads. Ideal stationary holonomic constraints then give geodesics of the induced metric on the regular constraint manifold. Nonholonomic constraints require their own constrained dynamics.

9 Loads and Energy Transfer Through a Joint

9.1 A Body Balance Is More Informative Than a Speed Sequence

For a rigid body measured in an inertial frame, $T_i=\tfrac12m_i\|v_{C,i}\|^2+\tfrac12\omega_i^{\mathsf T}I_{C,i}\omega_i$. A force $F$ at point $P$ contributes $F\cdot v_P$; a pure moment $N$ contributes $N\cdot\omega_i$. Their sum is $\dot T_i$.

At an ideal common joint point $P$, force $F$ on the distal body and $-F$ on the proximal body contribute $F\cdot v_P$ and $-F\cdot v_P$. They cancel globally while transferring energy between bodies. A joint actuator moment $N$ contributes combined power $N\cdot(\omega_{\rm distal}-\omega_{\rm proximal})$, which need not vanish.

If contact points slip, their velocities differ and the summed contact-force power need not be zero. If a support moves, it can supply work. Forces and moments also have different units: “reaction forces are larger than torques” is not a meaningful magnitude comparison without a lever arm.

Reduced coordinates eliminate ideal joint reactions from virtual work, not from the physical bodies. Recover them from body balance or a consistent redundant-coordinate formulation. Per-coordinate effort powers are coordinate-dependent bookkeeping; they are not interchangeable with segment energy balances.

9.2 Work the Pendulum Transfer Explicitly

Let $F_{12}$ be the force of rod 1 on rod 2 at $p_1$, and $\mathbf g=(0,-g)$. For the point-mass model,

\[ F_{12}=m_2(\ddot p_2-\mathbf g),\qquad F_{01}=m_1(\ddot p_1-\mathbf g)+F_{12}. \]

The absolute rod rates are $\omega_1=v_1$, $\omega_2=v_1+v_2$. Since this model concentrates each rod’s mass at its endpoint, its body energy balances become

\[ \begin{aligned} \dot T_1&=\tau_1\omega_1-\tau_2\omega_1 -F_{12}\cdot\dot p_1+m_1\mathbf g\cdot\dot p_1,\\ \dot T_2&=\tau_2\omega_2 +F_{12}\cdot\dot p_1+m_2\mathbf g\cdot\dot p_2. \end{aligned} \tag{25}\]

The force-transfer terms cancel. The two joint-2 moment terms sum to $\tau_2v_2$. Adding potential-energy rates gives $\frac{d}{dt}(T_1+T_2+V)=\tau_1v_1+\tau_2v_2$. The executable checks verify both body equations, including nonzero applied torques.

A slowing first rod need not be supplying energy through $F_{12}$: gravity and actuator moments also enter its balance. For a general human segment, a decreasing angular speed does not alone establish decreasing total kinetic energy because its center-of-mass speed can change. Even a measured loss of segment energy does not identify its destination without the intervening power balances.

This is the stronger version of the whip-chain argument: demonstrate the signed force and moment powers and reconcile each segment’s energy change. Sequential speed peaks motivate that investigation; they do not prove it. A flexible shaft additionally stores strain energy and dissipates some energy. Its phase and boundary work must be included if the argument concerns transfer through the shaft.

9.3 There Is No Universal Bias-Power Channel

For the ideal rigid model, the kinetic-energy identity is

\[ \dot T=v^{\mathsf T}(Bu+Q_{\rm ext}+Q_c-g_V),\qquad \dot E=v^{\mathsf T}(Bu+Q_{\rm ext}+Q_c),\quad E=T+V. \tag{26}\]

The velocity products cancel the $\tfrac12v^{\mathsf T}\dot Mv$ term in differentiating $T$. Calling $(c+g_V)^{\mathsf T}v$ a positive “drift power” ignores that cancellation, the gravity sign and coordinate dependence. Gravity can convert potential to kinetic energy; inertial coupling changes how motion and energy are distributed. Neither is a new energy source.

10 Symmetry, Representation, and Conservation

Noether’s theorem relates a continuous variational symmetry of the action to a conserved quantity, with appropriate conditions on nonconservative forces. A time-independent damped system does not conserve mechanical energy merely because its equations have no explicit time dependence.

For the fixed-pivot pendulum with a stationary potential and no applied work or dissipation, total mechanical energy is conserved. The pivot can exchange momentum with the support; gravity selects a direction. One must examine the actual symmetry of the specified system before asserting linear or angular momentum conservation.

If a smooth quantity $I(x)$ is conserved under the declared drift, $DI\,f_0=0$, then $\dot I=DI\,Gu$ in that model. This identifies a role for input, but does not guarantee the actuators can produce every desired change: $DI\,G$ can vanish or be constrained.

Changing coordinate labels preserves physical predictions after transformation. Changing observer is a different operation. Under a constant-velocity boost, a point’s power transforms as $F\cdot(v-U)$ and kinetic energy changes. Under accelerating or rotating observers, acceleration components acquire additional terms. Consistent mechanics is objective; individual numerical energies, powers or accelerations are not universally identical for every observer.

11 Cross-Formalism Verification and Swing Investigation

The useful comparison is concrete:

Question	Quantity to Compare	Conditions to Match
Does the inertia represent the geometry?	$\tfrac12v^{\mathsf T}Mv$ and summed Cartesian body energy	Mass locations, inertias, coordinate rates
Do the motion and loads agree?	Forward acceleration and inverse effort/body balance	Gravity signs, effort channels, external loads
Is support represented consistently?	Constraint residual and reaction forces	Copied admissible state, curvature, contact mode
Does a spatial calculation give the same endpoint motion?	Recovered $\dot p,\ddot p$ and geometric-Jacobian results	Home geometry, origins, twist ordering
What can input change immediately?	$M^{-1}B$, $B_c$, or $JB_c$	Coordinates, support, task rank, input limits
Is energy transferred between segments?	Force/moment powers and each body-energy rate	Common frame, points of application, elastic storage
Does an intervention improve impact or flight?	A declared outcome from complete alternative trajectories	Input history, event convention, uncertainty, validation data

Disagreements can reveal differing assumptions, frame conventions, load retention, integration tolerance or an implementation error. Resolve those differences before concluding that two mathematical formulations disagree physically. A check against the same formula copied twice is weaker than a Cartesian/body-balance check derived independently.

For a golf study, the connected chain is: actuator and passive forces shape body and shaft motion; support determines admissible motion and reactions; the approaching club state and contact geometry determine impact; launch velocity, orientation and spin then feed aerodynamic flight. A change that increases head speed can also change face orientation or impact location. A net-torque reduction can coexist with increasing distal speed, but that coexistence does not establish an optimal torque policy.

The next evidential step is to specify a measurable intervention and outcome, then ask which measurements could distinguish the proposed explanation from competing ones. Examples include joint reaction estimates reconciled with segment energy, measured grip motion and shaft deformation, repeated impact locations, and held-out launch outcomes. Do not equate an inverse-dynamics net moment with a unique muscle-force pattern or an instantaneous acceleration attribution with a finite causal effect.

12 Executable Verification

The following self-contained NumPy example implements the declared two-point-mass model and one stationary bilateral-constraint solve. All lengths are in meters, masses in kilograms, angles in radians and times in seconds. The selected numbers are manufactured verification parameters, not measurements of a golfer.

The constrained solver assumes an already admissible state, independent constraint rows and positive-definite inertia. It deliberately does not resolve unilateral contact, friction, impacts or inconsistent geometry. A singular solve or an inadmissible reaction requires reconsidering the model; a pseudoinverse is not a substitute for a missing contact law.

# Multibody equivalence example
import numpy as np

GRAVITY_M_S2 = 9.81


def pendulum_data(
    q: np.ndarray,
    velocity: np.ndarray,
    lengths: np.ndarray,
    masses: np.ndarray,
    gravity: float = GRAVITY_M_S2,
) -> dict[str, np.ndarray]:
    """Return the declared two-point-mass model in relative-angle coordinates."""
    q, velocity, lengths, masses = [
        np.asarray(value, dtype=float) for value in (q, velocity, lengths, masses)
    ]
    if any(value.shape != (2,) for value in (q, velocity, lengths, masses)):
        raise ValueError("Configuration, velocity, lengths and masses must be 2-vectors")
    if not all(np.all(np.isfinite(value)) for value in (q, velocity, lengths, masses)):
        raise ValueError("All state and parameter values must be finite")
    if np.any(lengths <= 0) or np.any(masses <= 0) or not np.isfinite(gravity) or gravity < 0:
        raise ValueError("Require positive lengths/masses and nonnegative finite gravity")
    l1, l2 = lengths
    m1, m2 = masses
    a1, a2 = q[0], q.sum()
    w1, w2 = velocity[0], velocity.sum()
    r1 = l1 * np.array([np.sin(a1), -np.cos(a1)])
    r2 = l2 * np.array([np.sin(a2), -np.cos(a2)])
    t1 = l1 * np.array([np.cos(a1), np.sin(a1)])
    t2 = l2 * np.array([np.cos(a2), np.sin(a2)])
    j1 = np.column_stack((t1, np.zeros(2)))
    j2 = np.column_stack((t1 + t2, t2))
    mass = np.array(
        [
            [
                (m1 + m2) * l1**2 + m2 * l2**2 + 2 * m2 * l1 * l2 * np.cos(q[1]),
                m2 * l2**2 + m2 * l1 * l2 * np.cos(q[1]),
            ],
            [m2 * l2**2 + m2 * l1 * l2 * np.cos(q[1]), m2 * l2**2],
        ]
    )
    coupling = m2 * l1 * l2 * np.sin(q[1])
    products = coupling * np.array(
        [-2 * velocity[0] * velocity[1] - velocity[1] ** 2, velocity[0] ** 2]
    )
    potential_gradient = gravity * np.array(
        [(m1 + m2) * l1 * np.sin(a1) + m2 * l2 * np.sin(a2), m2 * l2 * np.sin(a2)]
    )
    return {
        "mass": mass,
        "bias": products + potential_gradient,
        "gravity": potential_gradient,
        "points": np.array([r1, r1 + r2]),
        "jacobians": np.array([j1, j2]),
        "curvature": np.array([-(w1**2) * r1, -(w1**2) * r1 - w2**2 * r2]),
    }


def constrained_acceleration(
    mass: np.ndarray, force: np.ndarray, jacobian: np.ndarray, curvature: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
    """Solve one stationary bilateral mode; force already includes all bias loads."""
    mass, force, jacobian, curvature = [
        np.asarray(value, dtype=float) for value in (mass, force, jacobian, curvature)
    ]
    if force.ndim != 1 or curvature.ndim != 1:
        raise ValueError("Force and curvature must be vectors")
    n, k = len(force), len(curvature)
    if mass.shape != (n, n) or jacobian.shape != (k, n) or n == 0 or k == 0:
        raise ValueError("Incompatible nonempty mass, force and constraint dimensions")
    if not all(np.all(np.isfinite(value)) for value in (mass, force, jacobian, curvature)):
        raise ValueError("All solve inputs must be finite")
    if not np.allclose(mass, mass.T):
        raise ValueError("Mass matrix must be symmetric")
    np.linalg.cholesky(mass)
    block = np.block([[mass, -jacobian.T], [jacobian, np.zeros((k, k))]])
    solution = np.linalg.solve(block, np.r_[force, -curvature])
    return solution[:n], solution[n:]


lengths = np.array([0.8, 0.6])
masses = np.array([1.4, 0.5])
q = np.array([0.4, -0.7])
velocity = np.array([3.0, -2.0])
torque = np.array([2.0, -0.3])
data = pendulum_data(q, velocity, lengths, masses)
drift = np.linalg.solve(data["mass"], -data["bias"])
increment = np.linalg.solve(data["mass"], torque)
acceleration = np.linalg.solve(data["mass"], torque - data["bias"])
np.testing.assert_allclose(acceleration, drift + increment)
cartesian_acceleration = data["jacobians"] @ acceleration + data["curvature"]
loads = masses[:, None] * (cartesian_acceleration - np.array([0.0, -GRAVITY_M_S2]))
recovered = sum(j.T @ f for j, f in zip(data["jacobians"], loads, strict=True))
np.testing.assert_allclose(recovered, torque, atol=1e-12)

The sample compares a forward acceleration solve with Cartesian inertial loads projected back to joint efforts. The accompanying regression checks also perturb Cartesian geometry to verify Jacobians and curvature, test straight/folded configurations, reconcile individual body powers, and exercise the circular constraint’s input-induced reaction change.

These are mathematical and implementation checks. Numerical integration would need independent timestep and energy-error checks; a human or equipment model would need parameter identification, uncertainty analysis and validation against observations.

13 Exercises and a Working Research Standard

Derive $M$ from $m_1J_1^{\mathsf T}J_1+m_2J_2^{\mathsf T}J_2$, including its determinant. Explain why task rank loss does not make this inertia singular.
Recover the pendulum efforts from Cartesian inertial loads for a nonzero torque vector. Reconcile both individual body-power equations.
Differentiate the endpoint-circle constraint twice. Compare the full solve with one that omits $\dot J_cv$, and identify the missing physical acceleration.
Use a moving support to give an example of a constraint reaction doing nonzero work.
Derive the space twist using $\dot g g^{-1}$, recover endpoint velocity and compare it with the geometric Jacobian.
Transform from relative to absolute angles and preserve virtual work. Then use a nonlinear chart to expose the acceleration Hessian term.
Linearize a variable-inertia model about a nonzero input. Quantify the error caused by omitting $(\partial G/\partial x)\bar u$.
State a proposed golf torque intervention, retained state and loads, impact-time convention and measured outcome. Explain why a pointwise attribution alone does not answer it.
Extend the body-energy balance to a shaft mode with strain energy and damping. Identify the boundary work that couples it to the grip and head.
Identify which observations would falsify your proposed energy-transfer explanation, and which checks merely verify that its equations were implemented consistently.

The standard is to make the connections computable: declare the model, derive its motion and loads, reconcile power, propagate interventions through contact, and compare the resulting outputs with independent evidence. Each formulation then contributes something specific to understanding the swing.

Cite this page

Citation

BibTeX citation: Download .bib

@online{olson2026,
  author = {Olson, Dieter},
  title = {Multibody {Dynamics:} {Drift,} {Control,} and the {Geometry}
    of {Physical} {Law}},
  date = {2026-03-06},
  url = {https://affinedrift.com/articles/calculation-framework-comparison/multibody-drift-control-v3.html},
  langid = {en}
}

For attribution, please cite this work as:

Olson, Dieter. 2026. “Multibody Dynamics: Drift, Control, and the Geometry of Physical Law.” March 6. https://affinedrift.com/articles/calculation-framework-comparison/multibody-drift-control-v3.html.

--- title: "Multibody Dynamics: Drift, Control, and the Geometry of Physical Law" description: "Derive modeled motion through energy, body forces, twists, and task coordinates, applying those connections to investigate golf and humanoid dynamics." author: "Dieter Olson" date: "2026-03-06" date-source: "initial-publication-record" css: multibody-drift-control-v3.css format: html: toc: true toc-depth: 3 number-sections: true code-fold: false html-math-method: mathjax categories: - theory-core --- ## One Mechanism, Several Questions {#sec-introduction} A golfer's pelvis can slow while the clubhead gains speed. That observation raises several questions: what forces act between segments, where does energy go, how do current torques affect acceleration, and what would a different torque history do to impact? These questions are related, but they are not interchangeable. A segment-speed graph cannot answer all of them. Multibody dynamics gives us several ways to examine the same modeled mechanism. An energy derivation exposes the inertia and generalized efforts. Body force balances recover joint reactions. Spatial twists organize three-dimensional motion and wrenches. Task coordinates reveal what the clubhead or foot can do at its current posture. Geometry clarifies which conclusions survive a change of coordinates. The value of using several formulations is that their predictions can be compared under the **same model, state, loads, constraints and observation convention**. Agreement is a strong implementation check. It does not establish that the chosen model accurately represents a golfer. That further step requires independent measurements. ### The Shared Drift-Control Structure For a declared input vector $u$, many models have the form $$ \dot x=f_0(x,t)+G(x,t)u. $$ {#eq-affine-master} At a fixed state and time, $f_0$ is the result with those input channels set to zero; $Gu$ is the incremental state derivative. This is **control affinity**. It is linear in input at that fixed state, even when the state dependence and the resulting trajectories are nonlinear. The river analogy can help if used carefully: paddling now changes where the kayaker goes and therefore which current they encounter later. Similarly, a torque changes future configuration and velocity, hence future drift and input effectiveness. A pointwise sum of derivatives is not a sum of independently generated finite trajectories. The zero reference must be declared. Zero commanded motor voltage, zero joint torque and zero neural excitation are different interventions. Existing activation, passive tissue forces, shaft deformation and ground reactions can persist when a chosen command is zero. Nor does the split divide “physics” from “intent”: actuator forces are physical, and many mechanisms can produce the same net joint torque. ### Scope of the Comparison We first use independent local coordinates $q$, velocity $v=\dot q$, a symmetric positive-definite inertia $M(q)$, and generalized effort $\tau=B(q)u$. We then add regular stationary ideal bilateral constraints. The worked mechanism has two massless rigid rods with positive point masses at their distal ends, a fixed pivot, gravity and ideal planar revolute joints. Damping or elastic coordinates can be included without necessarily losing input affinity. Actuator dynamics require an expanded state; affinity must then be checked with respect to the actual chosen input. Contact switching, impacts, friction and hysteresis need additional rules. The [hybrid reference](../tangent-hyperplane-articles/Advanced/Hybrid_Tangent_Spaces.html) develops the event and contact distinctions. For three-dimensional floating bodies, a quaternion configuration may have more components than its generalized velocity. Then $\dot q=N(q)v$ and the equations use $\dot v$, not an unqualified identification of $v$ with $\dot q$. A humanoid also needs its base motion and a declared support mode; a joint-angle vector alone may not be a complete state. ## Configuration, Effort, and Inertia {#sec-core-concepts} ### State Is Sufficient Only for the Declared Model Configuration specifies the arrangement of bodies. Velocity specifies how that arrangement is changing. Under a well-posed second-order rigid model, $x=(q,v)$, the specified inputs and parameters determine the local evolution. A muscle activation or viscoelastic memory variable cannot be omitted if it affects that evolution independently of $q,v$. State summarizes the information a model retains; it does not uniquely reconstruct the historical causes of that state. Different torque histories can arrive at the same configuration and velocity. ### Generalized Force Is Defined by Work A generalized effort $Q$ satisfies $$ \delta W=Q^{\mathsf T}\delta q,\qquad P=Q^{\mathsf T}v. $$ An angle coordinate pairs with a moment; a translation pairs with a force. For a point load $F$ at $p(q)$, $\delta p=J\delta q$ gives $Q=J^{\mathsf T}F$. This work relation, rather than “whatever makes this coordinate accelerate,” defines generalized force. Coupling means one effort can affect several accelerations. The kinetic energy is $$ T(q,v)=\tfrac12v^{\mathsf T}M(q)v. $$ $M$ maps an acceleration increment to an inertial generalized-force increment. Its inverse maps a force increment to acceleration when other state-dependent terms and constraints are treated consistently. A surface that adds rolling resistance changes a dissipative load, not the cart's inertia merely by being rougher. Positive definiteness means every nonzero admissible independent velocity has positive kinetic energy. It requires the corresponding inertial assumptions. Redundant coordinates, massless degrees of freedom or a singular chart need separate treatment. ### The Manipulator Equation and Its Signs Write $$ M(q)\ddot q+b(q,v)=B(q)u+Q_{\rm ext},\qquad b=c(q,v)+g_V(q),\quad g_V=\frac{\partial V}{\partial q}. $$ {#eq-manipulator} Here $c=C(q,v)v$ contains the velocity-product terms of this coordinate description. The physical generalized gravitational force is $-g_V$; $g_V$ has the positive sign because it is on the left. $Q_{\rm ext}$ denotes additional loads, so gravity is not counted twice. The partition of velocity products into “Coriolis” and “centripetal” terms depends on coordinates and on the nonunique choice of $C$. Centripetal acceleration points inward. Water following a circular path needs an inward net force; an outward centrifugal term belongs to a rotating-frame description. Neither term is a new source of energy. ## A Double Pendulum Derived From Cartesian Motion {#sec-double-pendulum} ### Geometry and Velocity Let $q=(\theta_1,\theta_2)$, with $\theta_1$ measured counterclockwise from downward vertical and $\theta_2$ the second rod's angle relative to the first. Positive vertical position points upward. Define $$ r_1=l_1\begin{bmatrix}\sin\theta_1\\-\cos\theta_1\end{bmatrix}, \qquad r_2=l_2\begin{bmatrix}\sin(\theta_1+\theta_2)\\-\cos(\theta_1+\theta_2)\end{bmatrix}. $$ {#eq-dp-coords} The mass positions are $p_1=r_1$, $p_2=r_1+r_2$. The absolute angular velocities are $\omega_1=v_1$, $\omega_2=v_1+v_2$. Confusing the second relative rate with its absolute rate changes both kinetic energy and joint power. $$ J_1=\begin{bmatrix}l_1\cos\theta_1&0\\l_1\sin\theta_1&0\end{bmatrix}, \quad J_2=\begin{bmatrix} l_1\cos\theta_1+l_2\cos(\theta_1+\theta_2)&l_2\cos(\theta_1+\theta_2)\\ l_1\sin\theta_1+l_2\sin(\theta_1+\theta_2)&l_2\sin(\theta_1+\theta_2) \end{bmatrix}. $$ {#eq-dp-jacobian} Then $\dot p_i=J_i v$. Differentiating once more gives $$ \ddot p_1=J_1\ddot q-\omega_1^2r_1,\qquad \ddot p_2=J_2\ddot q-\omega_1^2r_1-\omega_2^2r_2. $$ {#eq-dp-kinematics} The last terms are the Cartesian curvature or bias accelerations. They remain even if both relative joint accelerations are zero. ### Inertia, Bias, and Potential Start from $T=\tfrac12m_1\|\dot p_1\|^2+\tfrac12m_2\|\dot p_2\|^2$. Substitution gives $M=m_1J_1^{\mathsf T}J_1+m_2J_2^{\mathsf T}J_2$, hence $$ M= \begin{bmatrix} (m_1+m_2)l_1^2+m_2l_2^2+2m_2l_1l_2\cos\theta_2& m_2l_2^2+m_2l_1l_2\cos\theta_2\\ m_2l_2^2+m_2l_1l_2\cos\theta_2&m_2l_2^2 \end{bmatrix}. $$ {#eq-dp-mass} Its determinant is $m_2l_1^2l_2^2(m_1+m_2\sin^2\theta_2)>0$. Thus straight or folded alignment does not singularize this joint inertia when all masses and lengths are positive. The off-diagonal entry can vanish or change sign for some length ratios; “weakest coupling when folded” is not a general statement about all inertial responses. The Euler–Lagrange operation $d(\partial T/\partial v)/dt-\partial T/\partial q$ yields $$ c=\beta \begin{bmatrix}-2v_1v_2-v_2^2\\v_1^2\end{bmatrix}, \qquad \beta=m_2l_1l_2\sin\theta_2. $$ {#eq-dp-coriolis} At $\sin\theta_2=0$, this vector vanishes at that instant, but $M$ is generally still coupled. It does not follow that the joints move independently. With $V=m_1g(p_1)_y+m_2g(p_2)_y$, $$ g_V=g\begin{bmatrix} (m_1+m_2)l_1\sin\theta_1+m_2l_2\sin(\theta_1+\theta_2)\\ m_2l_2\sin(\theta_1+\theta_2) \end{bmatrix}. $$ {#eq-dp-gravity} These expressions can also be recovered from Cartesian inertial loads, as the executable example below checks. [MIT's multibody notes](https://underactuated.mit.edu/multibody.html) give the standard manipulator and constraint framework; signs must be translated consistently when gravity is placed on the opposite side. ## Pointwise Attribution and Future Control {#sec-drift-control} ### The Instantaneous Split With retained additional loads independent of $u$, $$ a_0=M^{-1}(Q_{\rm ext}-b),\qquad \Delta a=M^{-1}Bu,\qquad \ddot q=a_0+\Delta a. $$ {#eq-drift-control-split} The pointwise Zero-Torque Counterfactual (ZTCF) sample uses the copied state and the declared zeroed channels. A forward ZTCF trajectory repeatedly solves the zero-input dynamics while evolving its own state. It is a different object from a series of pointwise samples evaluated along the original controlled motion. High speed does not necessarily imply large $a_0$. A free particle in Cartesian coordinates can move arbitrarily fast with zero drift acceleration. In this pendulum, aligned rods have $c=0$ regardless of the rates. Posture, direction, gravity and cancellation all matter. A norm ratio of drift and input increments also needs units, coordinate scales and a zero-denominator convention. It does not determine energy contribution, neural effort or whether a finite torque intervention improves a shot. ### Linearize the Trajectory Actually Being Used For a nominal solution $\bar x(t),\bar u(t)$, the first variation is $$ \delta\dot x=A(t)\delta x+G(\bar x,t)\delta u,\qquad A=\left.\frac{\partial f_0}{\partial x}\right|_{\bar x} +\sum_j\bar u_j\left.\frac{\partial G_j}{\partial x}\right|_{\bar x}. $$ {#eq-local-linearization} The second term is needed even though the dynamics are affine in input. If a reference path is not a solution, the error equation also contains $f_0(\bar x,t)+G(\bar x,t)\bar u-\dot{\bar x}$. Sampling only the zero-input Jacobian on a controlled path answers a different question. Frozen eigenvalues of $A(t)$ do not establish stability of the time-varying system. For example, alternating for one time unit each between $A_1=\begin{bmatrix}-1&4\\0&-1\end{bmatrix}$ and $A_2=A_1^{\mathsf T}$ gives two matrices whose eigenvalues are all $-1$, but the two-step transition $e^{A_2}e^{A_1}$ has an eigenvalue greater than two. Stability concerns the transition over time, or a suitable Lyapunov/contraction bound. ### Why a Torque Change Can Matter Elsewhere A shoulder torque can affect clubhead acceleration immediately through coupled inertia. Over time it also changes posture, velocity, shaft state and when contact occurs. These changes alter subsequent control effectiveness. To investigate a proposed “release,” specify its torque or activation history, copied state, support mode and output. Then compare the resulting trajectories at a declared clock time or at each trajectory's own impact event. ## Constraints Redirect Acceleration and Change Reactions {#sec-constrained} ### Stationary Ideal Bilateral Constraints Let $\phi(q)=0$ have $k$ independent rows locally, $J_c=\partial\phi/\partial q$, and $\kappa=\dot J_c v$. An admissible state satisfies $\phi=0$, $J_cv=0$. Its acceleration satisfies $$ M\ddot q+b=Bu+Q_{\rm ext}+J_c^{\mathsf T}\lambda,\qquad J_c\ddot q=-\kappa. $$ {#eq-cel-dynamics} The allowed accelerations form an affine subspace $\mathcal H=\{a:J_ca=-\kappa\}$ of dimension $n-k$, a hyperplane only when the rank is one. They do not generally lie in the velocity tangent space: a train following a curved track needs normal acceleration. Define $W=J_cM^{-1}J_c^{\mathsf T}$, which is positive definite under these assumptions. If $a_{\rm free}=M^{-1}(Bu+Q_{\rm ext}-b)$, substitution gives $$ \lambda=-W^{-1}(J_ca_{\rm free}+\kappa),\qquad \ddot q=a_{\rm free}-M^{-1}J_c^{\mathsf T}W^{-1}(J_ca_{\rm free}+\kappa). $$ {#eq-cel-constraint} This is equivalently the solution of the block system $$ \begin{bmatrix}M&-J_c^{\mathsf T}\\J_c&0\end{bmatrix} \begin{bmatrix}\ddot q\\\lambda\end{bmatrix} =\begin{bmatrix}Bu+Q_{\rm ext}-b\\-\kappa\end{bmatrix}. $$ {#eq-hyperplane} Use linear solves rather than explicitly forming inverses in numerical code. Redundant constraints need a rank-aware treatment; nonunique multipliers need not imply nonunique physical acceleration. ### The Constrained Input Map Let $P=I-M^{-1}J_c^{\mathsf T}W^{-1}J_c$. This projects acceleration increments onto $\ker J_c$ orthogonally in the $M$ metric: $P^2=P$, $P^{\mathsf T}M=MP$. Then $$ a_0^c=P M^{-1}(Q_{\rm ext}-b)-M^{-1}J_c^{\mathsf T}W^{-1}\kappa, \qquad B_c=P M^{-1}B, $$ {#eq-constrained-drift} $$ \ddot q=a_0^c+B_cu,\qquad \lambda(u)-\lambda(0)=-W^{-1}J_cM^{-1}Bu. $$ {#eq-constrained-ctrl} The input changes the reaction as well as the acceleration. $B_cu$ lies in the homogeneous tangent directions of the affine acceleration set; the curvature offset belongs to the anchor. It is incorrect to project the entire acceleration onto a zero-offset tangent plane. For a manufactured particle of mass $2$ kg at $r=(2,0)$ m with velocity $(0,3)$ m/s, constrained by $\phi=(r^{\mathsf T}r-4)/2=0$, $J_c=(2,0)$ and $\kappa=9$. With no other load, $a_0^c=(-4.5,0)$ m/s². An applied force $(4,6)$ N changes acceleration to $(-4.5,3)$ m/s². The reaction changes from $(-9,0)$ to $(-13,0)$ N: it cancels the new radial input while retaining the inward curvature acceleration. ### What the Zero-Work Statement Requires For these stationary ideal constraints, $$ P_c=(J_c^{\mathsf T}\lambda)^{\mathsf T}v =\lambda^{\mathsf T}J_cv=0. $$ {#eq-hyperplane-control} Zero total power does not mean zero force or zero influence on motion. For a moving constraint $\phi(q,t)=0$, $J_cv=-\phi_t$ and $P_c=-\lambda^{\mathsf T}\phi_t$, which can be nonzero. Stationary ideal Pfaffian constraints $A(q)v=0$ likewise have zero reaction power, although their velocity distribution need not integrate to a configuration submanifold. A real foot-ground constraint is unilateral and friction limited. A computed reaction requiring tension from the floor or excessive tangential force invalidates that assumed contact mode. Removing torque may cause release or slip; then the same fixed-mode affine formula no longer describes the complete intervention. These conditions matter to humanoid support and to a golfer's use of the ground. For the pendulum endpoint constrained to a circle centered at $c_0$, choose $\phi=\|p_2-c_0\|^2-R^2$. Then $J_c=2(p_2-c_0)^{\mathsf T}J_2$. Differentiating the full expression, including $J_2$, gives the required curvature term. The same block solve applies wherever the constraint row is regular and the copied state is admissible. ## Body Forces and Spatial Motion {#sec-other-frameworks} ### Newton–Euler Reveals the Loads {#sec-newton-euler} For body $i$, take all forces acting **on that body**, expressed in a common frame. Let $r_{ik}$ run from its center of mass to the point where force $F_{ik}$ acts, and $N_{ik}$ denote applied pure moments. In a body frame fixed to the center of mass, the balances are $$ m_i a_{C,i}=\sum_kF_{ik},\qquad I_{C,i}\dot\omega_i+\omega_i\times(I_{C,i}\omega_i) =\sum_k\bigl(r_{ik}\times F_{ik}+N_{ik}\bigr). $$ {#eq-ne-trans} The first equation's acceleration is the physical center-of-mass acceleration expressed in that frame, including the appropriate transport terms if obtained by differentiating body components. The moment equation includes the force lever arms. One cannot simply add joint moments about different origins. In the standard planar reduction, the out-of-plane axis is principal and the scalar rotational equation uses its moment of inertia. Parallel angular velocities alone would not prove $\omega\times I\omega=0$ for an arbitrary three-dimensional inertia tensor. The [recursive Newton–Euler inverse-dynamics algorithm](https://modernrobotics.northwestern.edu/nu-gm-book-resource/8-3-newton-euler-inverse-dynamics/) takes configuration, velocity and **prescribed acceleration**, propagates motion outward, and computes required loads and joint efforts inward. A tree recursion has linear cost in the number of links for this problem. To obtain zero-torque motion, solve a **forward-dynamics** problem. One approach obtains the bias and mass-matrix columns from inverse-dynamics evaluations, then solves $Ma=-b+Q_{\rm ext}$. An articulated-body algorithm can instead solve tree forward dynamics directly. The [forward-dynamics tutorial](https://modernrobotics.northwestern.edu/nu-gm-book-resource/8-5-forward-dynamics-of-open-chains/) explains the former construction. Setting an inverse-dynamics output torque to zero without solving for acceleration is not an algorithm. Once acceleration is known, body balance recovers reactions. Additional constraints require the corresponding coupled solve. Thus Euler–Lagrange and Newton–Euler are not restricted to disjoint physical questions; their different representations make particular calculations easier to inspect. ### Spatial Twists Need an Origin Convention {#sec-poe-screw} Let $g_{sb}=[R,p;0,1]$ map body coordinates to a fixed space frame. With angular-first ordering, the space twist satisfies $\widehat{\mathcal V}_s=\dot g_{sb}g_{sb}^{-1}$, so $$ \mathcal V_s=\begin{bmatrix}\omega_s\\v_s\end{bmatrix}, \qquad v_s=\dot p-\omega_s\times p,\qquad \dot p=v_s+\omega_s\times p. $$ {#eq-poe-fk} The linear block $v_s$ is generally **not the velocity of the body-frame origin**. The derivative of the space twist is likewise not just ordinary endpoint acceleration: $$ \ddot p=\dot v_s+\dot\omega_s\times p+\omega_s\times\dot p. $$ The [twist tutorial](https://modernrobotics.northwestern.edu/nu-gm-book-resource/3-3-2-twists-part-2-of-2/) supplies the adjoint and matrix conventions. A product of exponentials (PoE) is a kinematic representation; inertia and loads are still needed for dynamics. For our downward-vertical home posture, use $$ S_1=(0,0,1,0,0,0)^{\mathsf T},\qquad S_2=(0,0,1,-l_1,0,0)^{\mathsf T}, $$ {#eq-dp-screws} and home translation $p(0)=(0,-l_1-l_2,0)$. Then $g_{sb}(q)=e^{\widehat S_1\theta_1}e^{\widehat S_2\theta_2}g_{sb}(0)$ and $J_s=[S_1,\operatorname{Ad}_{e^{\widehat S_1\theta_1}}S_2]$. Recover $\dot p$ using the cross-position term before comparing with $J_2v$. A screw through $(l_1,0,0)$ would describe a different home geometry. For dynamics in a body center-of-mass frame, write $\mathcal V_b=(\omega_b,v_{C,b})$, $\mathcal I_b=\operatorname{diag}(I_C,mI_3)$, and $$ \operatorname{ad}_{\mathcal V_b} =\begin{bmatrix}[\omega_b]_\times&0\\{}[v_{C,b}]_\times&[\omega_b]_\times\end{bmatrix}, \qquad \mathcal W_b=\mathcal I_b\dot{\mathcal V}_b -\operatorname{ad}_{\mathcal V_b}^{\mathsf T}\mathcal I_b\mathcal V_b. $$ {#eq-ne-rot} This declares the sign rather than relying on an ambiguous coadjoint-star notation. The [rigid-body dynamics tutorial](https://modernrobotics.northwestern.edu/nu-gm-book-resource/8-2-dynamics-of-a-single-rigid-body-part-2-of-2/) uses this convention. Wrenches at other points must include their translated moments. If $\mathcal V_s=\operatorname{Ad}_g\mathcal V_b$, then $\mathcal W_s=\operatorname{Ad}_g^{-\mathsf T}\mathcal W_b$, and $\mathcal W_s^{\mathsf T}\mathcal V_s=\mathcal W_b^{\mathsf T}\mathcal V_b$. This is invariance under consistent representation of the **same** wrench and motion. It does not say that mechanical power is unchanged for an observer moving with a different velocity. ## What the Clubhead or Foot Can Do {#sec-operational-space} ### Task Dynamics as an Induced Description Let $y=\psi(q)$, $J=\partial\psi/\partial q$. Direct differentiation gives $$ \ddot y=J M^{-1}(Bu+Q_{\rm ext}+J_c^{\mathsf T}\lambda-b)+\dot Jv. $$ {#eq-os-drift} This acceleration equation remains useful when a task inertia cannot be inverted. Where $J$ has full row rank, define mobility $H=JM^{-1}J^{\mathsf T}$ and task inertia $\Lambda=H^{-1}$. Then $$ \Lambda\ddot y+\mu+p_V =F_{\rm ctrl}^{\rm eff}+F_{\rm ext}^{\rm eff}+F_c^{\rm eff}, $$ {#eq-os-dynamics} where $\mu=\Lambda(JM^{-1}c-\dot Jv)$, $p_V=\Lambda JM^{-1}g_V$, and $F^{\rm eff}=\Lambda JM^{-1}Q$ for each corresponding generalized load. These quantities depend on the full joint state; a redundant robot's task position alone does not determine them. The [task-dynamics tutorial](https://modernrobotics.northwestern.edu/nu-gm-book-resource/8-6-dynamics-in-the-task-space/) illustrates the invertible-Jacobian case. A specified generalized torque has a unique instantaneous task effect. Redundancy concerns the opposite allocation problem: many torques may produce the same task effect. In a fully actuated unconstrained model, $$ \tau=J^{\mathsf T}F^{\rm eff}+\tau_{\rm null}, \qquad JM^{-1}\tau_{\rm null}=0. $$ This dynamically weighted null space differs from an arbitrary Euclidean projection. Its instantaneous lack of task acceleration does not guarantee no future task effect: it can change posture and later dynamics. Underactuation, effort bounds and support constraints further restrict which allocations are possible. ### A Task Singularity Is Not an Inertia Singularity For the pendulum endpoint, $\det J_2=l_1l_2\sin\theta_2$. At straight or folded alignment $J_2$ loses rank, $H$ is singular, and $\Lambda=H^{-1}$ is **undefined**. Describing $\Lambda$ as a finite singular matrix is wrong. A diverging effective inertia may describe a limit approaching that posture, but it does not supply a finite force that creates an unavailable instantaneous increment. At the straight-down posture, the vertical row of $J_2$ is zero. Joint acceleration therefore supplies no vertical task-acceleration increment at that instant. Nevertheless, $(\dot J_2v)_y=l_1v_1^2+l_2(v_1+v_2)^2$ can be positive: the endpoint has inward curvature acceleration. Both joint degrees of freedom and the positive-definite $M$ remain. For a supported mechanism, use the constrained map $JB_c$ for input authority. A clubhead Jacobian can be well conditioned while input authority is restricted by contact or missing actuators. Conversely, inability to move in one direction at one instant does not settle finite-time reachability. ## Geometry Without Changing the Physical Question {#sec-geometric} ### The Kinetic Metric For unconstrained pendulum angles with full revolutions allowed, configuration space is $Q=S^1\times S^1$. Joint limits restrict that domain. The inertia defines a Riemannian kinetic metric when it is positive definite: $\langle v,w\rangle_M=v^{\mathsf T}Mw$. It measures squared speed through kinetic energy, not the energy consumed by a displacement without a duration or force history. With the Levi-Civita connection of this metric, the intrinsic equation is $$ \nabla_vv=-\operatorname{grad}_M V +(Bu+Q_{\rm ext}+Q_c)^\sharp. $$ {#eq-geo-eom} Here $\sharp$ raises a force covector using $M^{-1}$, and $\operatorname{grad}_M V=M^{-1}g_V$ in coordinates. The coordinate expansion is $\ddot q^i+\Gamma^i_{jk}v^jv^k$, with $M\Gamma(v,v)=c$. Nonzero Christoffel symbols can arise in flat space, as in polar coordinates; they are not by themselves a curvature measurement. ### The Whole State Field Transforms Intrinsically In a local chart the zero-input second-order field has components $$ f_0(q,v)= \begin{bmatrix} v\\-\Gamma(v,v)-\operatorname{grad}_M V+Q_{\rm ext}^{\sharp} \end{bmatrix}. $$ {#eq-geo-drift} This displays an intrinsic field in coordinates; the individual lower block is not an ordinary vector under arbitrary configuration changes. If $z=\chi(q)$, then $$ \dot z=D\chi\,v,\qquad \ddot z=D\chi\,\ddot q+D^2\chi[v,v]. $$ For a free Cartesian particle with coordinate $q>0$, zero acceleration becomes $\ddot z=2v^2$ in the chart $z=q^2$. The physical motion has not acquired a force. Comparing two acceleration solutions at the same $q,v$, however, cancels this Hessian term: $\Delta\ddot z=D\chi\,\Delta\ddot q$. This is why the affine acceleration-space interpretation and its tangent increments require care. The switch from relative to absolute pendulum angles is linear: $\alpha=(\theta_1,\theta_1+\theta_2)$. If $v=A^{-1}\dot\alpha$, the metric becomes $M_\alpha=A^{-\mathsf T}MA^{-1}$ and effort becomes $Q_\alpha=A^{-\mathsf T}Q$. The paired work and kinetic energy agree. Zero input describes geodesic motion only when all nonconstraint forces vanish, including potential gradients, damping and additional loads. Ideal stationary holonomic constraints then give geodesics of the induced metric on the regular constraint manifold. Nonholonomic constraints require their own constrained dynamics. ## Loads and Energy Transfer Through a Joint {#sec-hidden-loads} ### A Body Balance Is More Informative Than a Speed Sequence For a rigid body measured in an inertial frame, $T_i=\tfrac12m_i\|v_{C,i}\|^2+\tfrac12\omega_i^{\mathsf T}I_{C,i}\omega_i$. A force $F$ at point $P$ contributes $F\cdot v_P$; a pure moment $N$ contributes $N\cdot\omega_i$. Their sum is $\dot T_i$. At an ideal common joint point $P$, force $F$ on the distal body and $-F$ on the proximal body contribute $F\cdot v_P$ and $-F\cdot v_P$. They cancel globally while transferring energy between bodies. A joint actuator moment $N$ contributes combined power $N\cdot(\omega_{\rm distal}-\omega_{\rm proximal})$, which need not vanish. If contact points slip, their velocities differ and the summed contact-force power need not be zero. If a support moves, it can supply work. Forces and moments also have different units: “reaction forces are larger than torques” is not a meaningful magnitude comparison without a lever arm. Reduced coordinates eliminate ideal joint reactions from virtual work, not from the physical bodies. Recover them from body balance or a consistent redundant-coordinate formulation. Per-coordinate effort powers are coordinate-dependent bookkeeping; they are not interchangeable with segment energy balances. ### Work the Pendulum Transfer Explicitly Let $F_{12}$ be the force of rod 1 on rod 2 at $p_1$, and $\mathbf g=(0,-g)$. For the point-mass model, $$ F_{12}=m_2(\ddot p_2-\mathbf g),\qquad F_{01}=m_1(\ddot p_1-\mathbf g)+F_{12}. $$ The absolute rod rates are $\omega_1=v_1$, $\omega_2=v_1+v_2$. Since this model concentrates each rod's mass at its endpoint, its body energy balances become $$ \begin{aligned} \dot T_1&=\tau_1\omega_1-\tau_2\omega_1 -F_{12}\cdot\dot p_1+m_1\mathbf g\cdot\dot p_1,\\ \dot T_2&=\tau_2\omega_2 +F_{12}\cdot\dot p_1+m_2\mathbf g\cdot\dot p_2. \end{aligned} $$ {#eq-body-power} The force-transfer terms cancel. The two joint-2 moment terms sum to $\tau_2v_2$. Adding potential-energy rates gives $\frac{d}{dt}(T_1+T_2+V)=\tau_1v_1+\tau_2v_2$. The executable checks verify both body equations, including nonzero applied torques. A slowing first rod need not be supplying energy through $F_{12}$: gravity and actuator moments also enter its balance. For a general human segment, a decreasing angular speed does not alone establish decreasing total kinetic energy because its center-of-mass speed can change. Even a measured loss of segment energy does not identify its destination without the intervening power balances. This is the stronger version of the whip-chain argument: demonstrate the signed force and moment powers and reconcile each segment's energy change. Sequential speed peaks motivate that investigation; they do not prove it. A flexible shaft additionally stores strain energy and dissipates some energy. Its phase and boundary work must be included if the argument concerns transfer through the shaft. ### There Is No Universal Bias-Power Channel For the ideal rigid model, the kinetic-energy identity is $$ \dot T=v^{\mathsf T}(Bu+Q_{\rm ext}+Q_c-g_V),\qquad \dot E=v^{\mathsf T}(Bu+Q_{\rm ext}+Q_c),\quad E=T+V. $$ {#eq-energy-balance} The velocity products cancel the $\tfrac12v^{\mathsf T}\dot Mv$ term in differentiating $T$. Calling $(c+g_V)^{\mathsf T}v$ a positive “drift power” ignores that cancellation, the gravity sign and coordinate dependence. Gravity can convert potential to kinetic energy; inertial coupling changes how motion and energy are distributed. Neither is a new energy source. ## Symmetry, Representation, and Conservation {#sec-symmetry} Noether's theorem relates a continuous variational symmetry of the action to a conserved quantity, with appropriate conditions on nonconservative forces. A time-independent damped system does not conserve mechanical energy merely because its equations have no explicit time dependence. For the fixed-pivot pendulum with a stationary potential and no applied work or dissipation, total mechanical energy is conserved. The pivot can exchange momentum with the support; gravity selects a direction. One must examine the actual symmetry of the specified system before asserting linear or angular momentum conservation. If a smooth quantity $I(x)$ is conserved under the declared drift, $DI\,f_0=0$, then $\dot I=DI\,Gu$ in that model. This identifies a role for input, but does not guarantee the actuators can produce every desired change: $DI\,G$ can vanish or be constrained. Changing coordinate labels preserves physical predictions after transformation. Changing observer is a different operation. Under a constant-velocity boost, a point's power transforms as $F\cdot(v-U)$ and kinetic energy changes. Under accelerating or rotating observers, acceleration components acquire additional terms. Consistent mechanics is objective; individual numerical energies, powers or accelerations are not universally identical for every observer. ## Cross-Formalism Verification and Swing Investigation {#sec-cross-formalism} The useful comparison is concrete: | Question | Quantity to Compare | Conditions to Match | |:--|:--|:--| | Does the inertia represent the geometry? | $\tfrac12v^{\mathsf T}Mv$ and summed Cartesian body energy | Mass locations, inertias, coordinate rates | | Do the motion and loads agree? | Forward acceleration and inverse effort/body balance | Gravity signs, effort channels, external loads | | Is support represented consistently? | Constraint residual and reaction forces | Copied admissible state, curvature, contact mode | | Does a spatial calculation give the same endpoint motion? | Recovered $\dot p,\ddot p$ and geometric-Jacobian results | Home geometry, origins, twist ordering | | What can input change immediately? | $M^{-1}B$, $B_c$, or $JB_c$ | Coordinates, support, task rank, input limits | | Is energy transferred between segments? | Force/moment powers and each body-energy rate | Common frame, points of application, elastic storage | | Does an intervention improve impact or flight? | A declared outcome from complete alternative trajectories | Input history, event convention, uncertainty, validation data | Disagreements can reveal differing assumptions, frame conventions, load retention, integration tolerance or an implementation error. Resolve those differences before concluding that two mathematical formulations disagree physically. A check against the same formula copied twice is weaker than a Cartesian/body-balance check derived independently. For a golf study, the connected chain is: actuator and passive forces shape body and shaft motion; support determines admissible motion and reactions; the approaching club state and contact geometry determine impact; launch velocity, orientation and spin then feed aerodynamic flight. A change that increases head speed can also change face orientation or impact location. A net-torque reduction can coexist with increasing distal speed, but that coexistence does not establish an optimal torque policy. The next evidential step is to specify a measurable intervention and outcome, then ask which measurements could distinguish the proposed explanation from competing ones. Examples include joint reaction estimates reconciled with segment energy, measured grip motion and shaft deformation, repeated impact locations, and held-out launch outcomes. Do not equate an inverse-dynamics net moment with a unique muscle-force pattern or an instantaneous acceleration attribution with a finite causal effect. ## Executable Verification {#sec-practical} The following self-contained NumPy example implements the declared two-point-mass model and one stationary bilateral-constraint solve. All lengths are in meters, masses in kilograms, angles in radians and times in seconds. The selected numbers are manufactured verification parameters, not measurements of a golfer. The constrained solver assumes an already admissible state, independent constraint rows and positive-definite inertia. It deliberately does not resolve unilateral contact, friction, impacts or inconsistent geometry. A singular solve or an inadmissible reaction requires reconsidering the model; a pseudoinverse is not a substitute for a missing contact law. ```python # Multibody equivalence example import numpy as np GRAVITY_M_S2 = 9.81 def pendulum_data( q: np.ndarray, velocity: np.ndarray, lengths: np.ndarray, masses: np.ndarray, gravity: float = GRAVITY_M_S2, ) -> dict[str, np.ndarray]: """Return the declared two-point-mass model in relative-angle coordinates.""" q, velocity, lengths, masses = [ np.asarray(value, dtype=float) for value in (q, velocity, lengths, masses) ] if any(value.shape != (2,) for value in (q, velocity, lengths, masses)): raise ValueError("Configuration, velocity, lengths and masses must be 2-vectors") if not all(np.all(np.isfinite(value)) for value in (q, velocity, lengths, masses)): raise ValueError("All state and parameter values must be finite") if np.any(lengths <= 0) or np.any(masses <= 0) or not np.isfinite(gravity) or gravity < 0: raise ValueError("Require positive lengths/masses and nonnegative finite gravity") l1, l2 = lengths m1, m2 = masses a1, a2 = q[0], q.sum() w1, w2 = velocity[0], velocity.sum() r1 = l1 * np.array([np.sin(a1), -np.cos(a1)]) r2 = l2 * np.array([np.sin(a2), -np.cos(a2)]) t1 = l1 * np.array([np.cos(a1), np.sin(a1)]) t2 = l2 * np.array([np.cos(a2), np.sin(a2)]) j1 = np.column_stack((t1, np.zeros(2))) j2 = np.column_stack((t1 + t2, t2)) mass = np.array( [ [ (m1 + m2) * l1**2 + m2 * l2**2 + 2 * m2 * l1 * l2 * np.cos(q[1]), m2 * l2**2 + m2 * l1 * l2 * np.cos(q[1]), ], [m2 * l2**2 + m2 * l1 * l2 * np.cos(q[1]), m2 * l2**2], ] ) coupling = m2 * l1 * l2 * np.sin(q[1]) products = coupling * np.array( [-2 * velocity[0] * velocity[1] - velocity[1] ** 2, velocity[0] ** 2] ) potential_gradient = gravity * np.array( [(m1 + m2) * l1 * np.sin(a1) + m2 * l2 * np.sin(a2), m2 * l2 * np.sin(a2)] ) return { "mass": mass, "bias": products + potential_gradient, "gravity": potential_gradient, "points": np.array([r1, r1 + r2]), "jacobians": np.array([j1, j2]), "curvature": np.array([-(w1**2) * r1, -(w1**2) * r1 - w2**2 * r2]), } def constrained_acceleration( mass: np.ndarray, force: np.ndarray, jacobian: np.ndarray, curvature: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: """Solve one stationary bilateral mode; force already includes all bias loads.""" mass, force, jacobian, curvature = [ np.asarray(value, dtype=float) for value in (mass, force, jacobian, curvature) ] if force.ndim != 1 or curvature.ndim != 1: raise ValueError("Force and curvature must be vectors") n, k = len(force), len(curvature) if mass.shape != (n, n) or jacobian.shape != (k, n) or n == 0 or k == 0: raise ValueError("Incompatible nonempty mass, force and constraint dimensions") if not all(np.all(np.isfinite(value)) for value in (mass, force, jacobian, curvature)): raise ValueError("All solve inputs must be finite") if not np.allclose(mass, mass.T): raise ValueError("Mass matrix must be symmetric") np.linalg.cholesky(mass) block = np.block([[mass, -jacobian.T], [jacobian, np.zeros((k, k))]]) solution = np.linalg.solve(block, np.r_[force, -curvature]) return solution[:n], solution[n:] lengths = np.array([0.8, 0.6]) masses = np.array([1.4, 0.5]) q = np.array([0.4, -0.7]) velocity = np.array([3.0, -2.0]) torque = np.array([2.0, -0.3]) data = pendulum_data(q, velocity, lengths, masses) drift = np.linalg.solve(data["mass"], -data["bias"]) increment = np.linalg.solve(data["mass"], torque) acceleration = np.linalg.solve(data["mass"], torque - data["bias"]) np.testing.assert_allclose(acceleration, drift + increment) cartesian_acceleration = data["jacobians"] @ acceleration + data["curvature"] loads = masses[:, None] * (cartesian_acceleration - np.array([0.0, -GRAVITY_M_S2])) recovered = sum(j.T @ f for j, f in zip(data["jacobians"], loads, strict=True)) np.testing.assert_allclose(recovered, torque, atol=1e-12) ``` The sample compares a forward acceleration solve with Cartesian inertial loads projected back to joint efforts. The accompanying regression checks also perturb Cartesian geometry to verify Jacobians and curvature, test straight/folded configurations, reconcile individual body powers, and exercise the circular constraint's input-induced reaction change. These are mathematical and implementation checks. Numerical integration would need independent timestep and energy-error checks; a human or equipment model would need parameter identification, uncertainty analysis and validation against observations. ## Exercises and a Working Research Standard {#sec-conclusion} 1. Derive $M$ from $m_1J_1^{\mathsf T}J_1+m_2J_2^{\mathsf T}J_2$, including its determinant. Explain why task rank loss does not make this inertia singular. 2. Recover the pendulum efforts from Cartesian inertial loads for a nonzero torque vector. Reconcile both individual body-power equations. 3. Differentiate the endpoint-circle constraint twice. Compare the full solve with one that omits $\dot J_cv$, and identify the missing physical acceleration. 4. Use a moving support to give an example of a constraint reaction doing nonzero work. 5. Derive the space twist using $\dot g g^{-1}$, recover endpoint velocity and compare it with the geometric Jacobian. 6. Transform from relative to absolute angles and preserve virtual work. Then use a nonlinear chart to expose the acceleration Hessian term. 7. Linearize a variable-inertia model about a nonzero input. Quantify the error caused by omitting $(\partial G/\partial x)\bar u$. 8. State a proposed golf torque intervention, retained state and loads, impact-time convention and measured outcome. Explain why a pointwise attribution alone does not answer it. 9. Extend the body-energy balance to a shaft mode with strain energy and damping. Identify the boundary work that couples it to the grip and head. 10. Identify which observations would falsify your proposed energy-transfer explanation, and which checks merely verify that its equations were implemented consistently. The standard is to make the connections computable: declare the model, derive its motion and loads, reconcile power, propagate interventions through contact, and compare the resulting outputs with independent evidence. Each formulation then contributes something specific to understanding the swing.