Variational Dynamics and the Moving Frame

Derives first-variation dynamics along reference motions, state transition matrices, and moving frame representations for swing sensitivity analysis.

Two swings can pass through almost the same pose and then produce different impact conditions. The difference may come from velocity, retained muscle or shaft states, contact conditions, or a changed command. Variational dynamics asks a precise local question: along a declared reference motion, how does a small change in state or input change a later outcome? It provides the derivative of a model prediction. It does not identify the golfer’s neural strategy from motion alone.

The distinction between a derivative and a finite intervention is central. The derivative obeys an exact linear equation. A finite change in a nonlinear swing generally does not. This chapter derives both the linear equation and its approximation error, develops the state transition matrix, and connects them to control authority, adjoint gradients and numerical computation. A moving frame can simplify the representation, but it does not remove physical sensitivity or guarantee stability.

First Variation Along a Feasible Motion

State, Input and Reference

Work in a regular local state chart and write \[ \dot x=F(x,u)=f(x)+G(x)u, \qquad x\in\mathbb R^n,\quad u\in\mathbb R^m. \] The input port must be specified: generalized torque, motor voltage and muscle excitation define different plants and retained states. A mechanical state commonly includes position and velocity; actuator activation, current or shaft deformation may add further states. Contact constraints must already be resolved in the chart or treated with an appropriate constrained formulation. The smooth results below hold within a fixed valid mode, not automatically through impact, slip or a coordinate singularity.

On a fixed finite interval \([t_0,T]\), let \(\bar u(t)\) be bounded and piecewise continuous, and let the feasible reference satisfy \[ \dot{\bar x}=F(\bar x,\bar u). \] If a proposed reference does not satisfy this equation, its defect is an additional forcing in the error dynamics. A planned curve and its timing are not made feasible simply by calling them nominal. Assume the state derivatives used below exist and are bounded on a tube containing the compared trajectories. State and input norms require chosen scales: adding radians, metres and activation units in an unscaled Euclidean norm has no invariant physical meaning.

A Family of Nearby Solutions

Introduce a scalar family parameter \(\varepsilon\), an initial direction \(\xi_0\), and a bounded piecewise-continuous input direction \(v(t)\): \[ x_\varepsilon(t_0)=\bar x(t_0)+\varepsilon\xi_0, \qquad u_\varepsilon(t)=\bar u(t)+\varepsilon v(t). \] Define the first variation \(\xi(t)=\partial_\varepsilon x_\varepsilon(t)|_{\varepsilon=0}\). Differentiating the differential equation and initial condition gives \[ \dot\xi=A(t)\xi+B(t)v, \qquad \xi(t_0)=\xi_0, \] where, writing the columns of \(G\) as \(g_j\), \[ A=Df(\bar x)+\sum_{j=1}^m\bar u_j Dg_j(\bar x), \qquad B=G(\bar x). \] Every term in \(A\) is an \(n\times n\) matrix. Its action is \(A\xi=Df\,\xi+(DG[\xi])\bar u\); no unspecified tensor contraction is needed. The reference command affects state sensitivity even though the original dynamics are affine in the declared input. A state-feedback comparison \(u=\pi(x,t)\) instead differentiates the closed-loop field and adds \(G D_x\pi\) to this Jacobian.

The first-variation equation is exact for the derivative \(\xi\). Under the stated local smoothness and boundedness conditions, \[ x_\varepsilon(t)=\bar x(t)+\varepsilon\xi(t)+O(\varepsilon^2) \] uniformly on a fixed finite interval for sufficiently small \(\varepsilon\). The constant depends on the interval, derivatives and amplification along the reference. This statement does not promise that one linear model remains accurate for every perturbation or indefinitely near an unstable trajectory. If the reference input lies on a bound, a two-sided family may be infeasible; then only admissible one-sided directions have a physical intervention interpretation.

The Taylor Remainder and Its Mixed Term

For a twice continuously differentiable vector field \(h\) and a state displacement \(e\), the exact first-order remainder is \[ h(\bar x+e)=h(\bar x)+Dh(\bar x)e+ \int_0^1(1-s)D^2h(\bar x+se)[e,e]\,ds. \] There is no extra factor of one-half outside this integral. The integral of \(1-s\) already equals one-half. For the scalar check \(h(x)=x^2\), its second derivative is 2 and the integral returns exactly \(e^2\).

Let \(F_0(x,t)=f(x)+G(x)\bar u(t)\). The finite error \(e=x-\bar x\) under a finite input increment \(w=u-\bar u\) satisfies \[ \dot e=Ae+Bw+\rho(t,e,w), \] with the exact forcing remainder \[ \begin{split} \rho(t,e,w)={}&\int_0^1(1-s)D_x^2F_0(\bar x+se,t)[e,e]\,ds\\ &+[G(\bar x+e)-G(\bar x)]w. \end{split} \] If \(\|D_x^2F_0\|\le H\) and \(\|DG\|\le L_G\) on the relevant segments, then \[ \|\rho\|\le\tfrac12H\|e\|^2+L_G\|e\|\,\|w\|. \] The mixed state/input term cannot be dropped merely because \(e\) is small. For example, at any fixed reference pair in \(F(x,u)=x^2+(1+x)u\), the exact remainder is \(e^2+ew\). Joint perturbations \(e,w=O(\varepsilon)\) make both terms second order. Input affinity removes a pure \(w^2\) term at frozen state; it does not remove state/input coupling along a trajectory.

The State Transition Matrix

Definition and Composition

For piecewise continuous \(A(t)\), define the state transition matrix by \[ \partial_t\Phi(t,s)=A(t)\Phi(t,s),\qquad\Phi(s,s)=I. \] Each column is the homogeneous solution from one coordinate basis vector, so \(\xi(t)=\Phi(t,s)\xi(s)\) when \(v=0\). The matrix is itself a solution of a matrix initial-value problem. Its coordinate representation changes if the basis or state scaling changes.

On an interval where the linear equation exists, uniqueness gives \[ \Phi(t,r)\Phi(r,s)=\Phi(t,s), \qquad \Phi(t,s)^{-1}=\Phi(s,t). \] For the first identity, both products propagate the same arbitrary initial vector and solve the same initial-value problem. For the second, set the endpoint equal to the start in the composition identity. These are identities for the exact fundamental solution, including backward propagation on the common interval. A numerically tiny singular value can make inversion inaccurate even though the exact matrix remains nonsingular.

For constant \(A\), \(\Phi(t,s)=\exp(A(t-s))\). In a time-varying system, \(\exp(\int_s^t A(r)dr)\) is not generally the answer. A sufficient condition for that simplification is pairwise commutation of the matrices at different times.

Liouville’s Identity

Jacobi’s determinant derivative, for invertible \(\Phi\), reads \[ \frac{d}{dt}\det\Phi=\det\Phi\,\operatorname{tr}(\Phi^{-1}\dot\Phi) =\det\Phi\,\operatorname{tr}A. \] The trace is unchanged by similarity. Equivalently, the adjugate is \(\operatorname{adj}\Phi=(\det\Phi)\Phi^{-1}\), without an inverse transpose. This derivative identity is distinct from the matrix determinant lemma for a rank-one update. Integrating the scalar equation gives \[ \det\Phi(t,s)=\exp\left(\int_s^t\operatorname{tr}A(r)\,dr\right)>0. \] It is a valuable independent check on a computed STM, but a correct determinant alone does not establish that its individual entries or singular directions are correct.

Time Ordering and the Peano–Baker Series

The matrix initial-value problem is equivalent to the Volterra equation \[ \Phi(t,s)=I+\int_s^t A(r)\Phi(r,s)\,dr. \] Repeated substitution produces a sum of ordered products. With \(\Phi_0=I\) and \(\Phi_{k+1}(t,s)=\int_s^t A(r)\Phi_k(r,s)dr\), \[ \Phi(t,s)=I+\int_s^t A(r_1)dr_1+ \int_s^t\int_s^{r_1}A(r_1)A(r_2)\,dr_2dr_1+\cdots. \] Later times occur on the left. This is a time-ordered sum, not an unrestricted average over all matrix permutations. If \(\|A(r)\|\le L\) in an induced norm and \(0\le t-s\le\Delta\), integration over the ordered simplex gives \[ \|\Phi_k(t,s)\|\le\frac{(L\Delta)^k}{k!}. \] The exponential majorant proves absolute and uniform convergence on this bounded interval. Passing to the integral equation identifies the sum with its unique solution. More generally, local integrability of \(\|A\|\) suffices; the short treatment by Baake and Schlaegel supplies that version.

A direct check makes order tangible. Apply \(A_1=\begin{pmatrix}0&1\\0&0\end{pmatrix}\) for one unit of time and then \(A_2=\begin{pmatrix}0&0\\1&0\end{pmatrix}\) for one unit. Each squares to zero, so \[ \Phi(2,0)=(I+A_2)(I+A_1)=\begin{pmatrix}1&1\\1&2\end{pmatrix}. \] Reversing the intervals gives \(\begin{pmatrix}2&1\\1&1\end{pmatrix}\). Neither ordering is recovered by pretending the two generators commute. The Magnus representation organizes a local logarithm \(\Phi=\exp\Omega\) using integrals and commutators; a truncated Magnus series needs its own convergence and error analysis. An exact infinite representation and a finite numerical approximation are different claims.

Variation of Constants and Outcome Sensitivity

Derivation of the Forced Solution

Write \(\xi(t)=\Phi(t,t_0)c(t)\). Differentiation cancels the homogeneous terms and yields \(\dot c=\Phi(t_0,t)B(t)v(t)\). Integration and composition give \[ \xi(t)=\Phi(t,t_0)\xi_0+ \int_{t_0}^t\Phi(t,s)B(s)v(s)\,ds. \] This is exact for the variational equation. Numerical integration and quadrature introduce additional errors. The first term propagates the initial variation; the second propagates each infinitesimal input contribution from its application time to the observation time. Initial-state effects and input effects superpose within this one linear problem.

For a differentiable fixed-time output \(y=h(x(T))\), define \(C_T=Dh(\bar x(T))\). Then \[ \delta y=C_T\Phi(T,t_0)\xi_0+ \int_{t_0}^T K(T,s)v(s)\,ds, \qquad K(T,s)=C_T\Phi(T,s)B(s). \] The output kernel, not \(B\) or \(\Phi\) alone, determines the first-order effect on the chosen task. Clubface orientation, normal clubhead velocity and strike location select different output rows. Their coordinates and units must be declared, and orientation differences require a valid local orientation chart.

Earlier inputs are not always more influential. In a scalar stable model \(\dot\xi=-a\xi+v\), \(a>0\), the terminal kernel is \(e^{-a(T-s)}\), so equal short input increments applied later have greater terminal-state influence. For a double integrator, the position kernel is \(T-s\) while the velocity kernel is 1. The same actuator and horizon therefore yield different timing sensitivities for different tasks. Coupling, constraints and actuator dynamics can rotate and filter these effects further.

Impact Time Is an Additional Variable

Golf impact is often an event rather than a fixed clock time. Suppose a smooth scalar guard \(g(x,t)=0\) is crossed transversely at \(T\), so \(g_xF+g_t\ne0\). A fixed-time state variation \(\xi(T)\) shifts the event by \[ \delta T=-\frac{g_x\xi(T)}{g_xF+g_t}. \] For a pre-impact output \(h(x(T),T)\), its event variation is \[ \delta y=h_x\xi(T)+(h_xF+h_t)\delta T. \] All derivatives are evaluated on the reference event. Parameter-dependent guards add their direct parameter terms. A reset or collision map requires its derivative and the appropriate event-time correction; grazing contact can invalidate this smooth sensitivity. A calculation at a fixed nominal impact time must not silently stand in for sensitivity of the actual collision outcome.

Finite Error and the Limits of Linear Prediction

Let \(z\) solve \(\dot z=Az+Bw\) with the same initial displacement as the finite nonlinear error \(e\). Their difference \(d=e-z\) obeys \[ \dot d=Ad+\rho(t,e,w),\qquad d(t_0)=0, \qquad d(t)=\int_{t_0}^t\Phi(t,s)\rho(s,e(s),w(s))\,ds. \] The nonlinear forcing \(\rho\) is not the entire derivative difference \(\dot e-\dot z\): that difference also contains \(Ad\). Combining the exact identity with the remainder bound gives \[ \|d(t)\|\le\int_{t_0}^t\|\Phi(t,s)\| \left[\tfrac12H\|e(s)\|^2+L_G\|e(s)\|\|w(s)\|\right]ds. \] For \(\|\Phi(t,s)\|\le K\) on the comparison interval, and suprema \(E=\sup\|e\|\), \(U=\sup\|w\|\), this is bounded by \(K(t-t_0)(HE^2/2+L_GEU)\). Stability by itself does not imply \(K\le1\) in an arbitrary norm; stable nonnormal systems can amplify errors transiently.

Twice continuous differentiability on a compact tube suffices for this quadratic bound. A cubic bound after subtracting the quadratic Taylor polynomial requires stronger control of the Hessian, such as a Lipschitz Hessian. The finite-error estimate concerns nonlinear derivatives in the chosen chart; it is not necessarily intrinsic manifold curvature. Even a scalar system on a flat line can have a nonzero remainder.

Re-linearizing can improve a local prediction, but it does not alone prove quadratic convergence of an optimizer. Newton and suitable SQP/DDP convergence results require the relevant second-order information, regularity, nondegeneracy and sufficiently accurate steps near an appropriate solution. Constraints, regularization, line searches and omitted second derivatives affect the rate. iLQR and DDP should not inherit a universal \(2^{-2^k}\) error law from the Taylor remainder. Nonlinear rollout and gradient verification remain necessary.

Discrete Flow and Its Derivative

Zero-Order-Held Inputs

For a fixed interval \(h>0\), let \(\varphi_h(x_k,u_k)\) denote the nonlinear solution at its endpoint with initial state \(x_k\) and held input \(u_k\). It satisfies \[ \varphi_h(x_k,u_k)=x_k+ \int_{t_k}^{t_k+h}F(x(s),u_k)\,ds. \] The trajectory inside the integral is the corresponding nonlinear solution. Drift is retained. Continuous input affinity does not imply that \(\varphi_h\) is affine in the held input: \(\dot x=xu\) gives \(\varphi_h(x,u)=xe^{uh}\).

Along the held-input reference, integrate \[ \dot X=AX,\quad X(t_k)=I, \qquad \dot U=AU+B,\quad U(t_k)=0. \] Then \(A_k=D_x\varphi_h=X(t_k+h)\) and \(B_k=D_u\varphi_h=U(t_k+h)\), or equivalently \[ A_k=\Phi(t_k+h,t_k),\qquad B_k=\int_{t_k}^{t_k+h}\Phi(t_k+h,s)B(s)\,ds. \] These are exact derivatives of the exact flow. Freezing \(A,B\) at the beginning of a nonlinear interval is an additional approximation. Differentiating a numerical step gives derivatives of that numerical map; for a fixed Runge–Kutta grid, differentiating its stages consistently agrees with applying those stages to the augmented sensitivity equations.

A Singular-Safe Linear Formula

For constant linear \(A,B\), the block exponential gives \[ \exp\left(h\begin{pmatrix}A&B\\0&0\end{pmatrix}\right) =\begin{pmatrix}A_d&B_d\\0&I\end{pmatrix}, \qquad A_d=e^{Ah}. \] Here \(B_d=\int_0^h e^{As}B\,ds\). Only when \(A\) is invertible may it also be written \(A^{-1}(e^{Ah}-I)B\); even then, explicitly forming an inverse or subtracting nearly equal matrices may be unattractive numerically. The block formula handles an integrator. For \(A=\begin{pmatrix}0&1\\0&0\end{pmatrix}\) and \(B=(0,1)^T\), \[ A_d=\begin{pmatrix}1&h\\0&1\end{pmatrix},\qquad B_d=\begin{pmatrix}h^2/2\\h\end{pmatrix}. \] The discrete transition from \(k\) to \(\ell\) is \(A_{\ell-1}\cdots A_k\), with the earliest interval on the right. Exact smooth-flow Jacobians are invertible; arbitrary discrete maps and coarse numerical approximations need not share that property.

A Checked Pendulum Example

Model and Analytic Transition

Use a point mass \(m\) on a massless rod of length \(\ell\), a fixed pivot, no damping, and angle \(\theta\) measured from downward vertical. Positive pivot torque \(\tau\) increases \(\theta\): \[ \ddot\theta=-\frac g\ell\sin\theta+\frac{\tau}{m\ell^2}. \] At the downward zero-torque equilibrium, with state \((\theta,\dot\theta)\), \[ A=\begin{pmatrix}0&1\\-\omega_0^2&0\end{pmatrix}, \quad B=\begin{pmatrix}0\\1/(m\ell^2)\end{pmatrix}, \quad \omega_0=\sqrt{g/\ell}. \] Since \(A^2=-\omega_0^2I\), separating even and odd exponential terms gives \[ \Phi(t,0)=\begin{pmatrix} \cos(\omega_0t)&\sin(\omega_0t)/\omega_0\\ -\omega_0\sin(\omega_0t)&\cos(\omega_0t) \end{pmatrix}. \] Differentiation verifies \(\dot\Phi=A\Phi\), and its determinant is 1. The eigenvalues are \(\pm i\omega_0\). This undamped linear oscillator is stable, not asymptotically stable. Stability of the nonlinear downward equilibrium can be established separately with its local energy minimum. Purely imaginary linear eigenvalues alone do not settle a general nonlinear system’s stability.

Numbers, Units and Impulse Interpretation

For this illustrative calculation set \(m=1\) kg, \(\ell=1\) m and rounded \(g=10\) m/s\(^2\), distinct from the 9.81 convention in other chapters. The small-amplitude period is \(2\pi/\sqrt{10}=1.986918\) s. The following values are generated by the program below and independently checked against a matrix-exponential implementation.

\(t\) (s)	\(\Phi_{11}\)	\(\Phi_{12}\) (s)	\(\Phi_{21}\) (1/s)	\(\Phi_{22}\)
0.25	0.703441	0.224760	-2.247601	0.703441
0.50	-0.010342	0.316211	-3.162109	-0.010342

The upper-right STM entry converts initial angular velocity to angle and has units of seconds; the lower-left has reciprocal-second units. A matrix norm that mixes these coordinates needs an explicit scale. For initial displacement \((0.1\,\mathrm{rad},0)^T\), the first-order response at 0.5 s is given below. It is not an exact finite-amplitude nonlinear pendulum trajectory.

Variation	\(\delta\theta\) (rad)	\(\delta\dot\theta\) (rad/s)
Initial Angle	-0.001034	-0.316211
Torque Impulse	0.224760	0.703441

The second row assumes zero initial variation and an ideal torque impulse at 0.25 s with angular impulse \(J_\tau=1\) N m s. The jump is \(\Delta\dot\theta=J_\tau/(m\ell^2)\), and the linear response at 0.5 s is \(\Phi(0.25,0)B J_\tau\). The Dirac distribution has units of 1/s, so its coefficient is torque times time, not torque. The application time 0.25 s is not a quarter of the stated period; that quarter is approximately 0.49673 s. A finite actuator pulse must also satisfy amplitude, rate and retained actuator-state limits. A unit impulse illustrates a linear response calculation and is not evidence for an executable muscle command.

Parameter Sensitivity and Adjoint Gradients

Forward Sensitivities

Let \(\dot x=F(x,u(t,p),p,t)\) with \(p\in\mathbb R^r\). All partial derivatives below hold the other arguments fixed. With fixed initial time, \[ S=\frac{\partial x}{\partial p},\qquad \dot S=F_x S+F_p+F_u u_p,\qquad S(t_0)=\frac{\partial x_0}{\partial p}. \] The direct-control term vanishes only when the prescribed input is independent of \(p\). If the control is state feedback, first form the closed-loop field and differentiate its state and parameter dependence. Variation of constants adds both the propagated initial sensitivity and the integral of the direct forcing. Assuming a zero initial sensitivity without declaring it can reverse an optimization conclusion.

For a differentiable observable \(h(x(T),p)\) at fixed \(T\), the derivative is \(h_xS(T)+h_p\). A maximum over time or a minimum distance can be nonsmooth when the active extremum changes or ties. Event-defined quantities require the event-time calculation above. Sensitivity is a local derivative, not a replacement for nonlinear simulation across an arbitrary parameter range.

Derivation of the Adjoint

For a scalar objective at fixed \(t_0,T\), \[ J(p)=\psi(x(T),p)+\int_{t_0}^T L(x,u,p,t)\,dt, \] define the column adjoint as the derivative of remaining cost with respect to the current state under the specified future input policy. With an open-loop prescribed input, it satisfies \[ -\dot\lambda=F_x^T\lambda+L_x, \qquad \lambda(T)=\psi_x. \] Here \(L_x\) and \(\psi_x\) are column gradients. To verify the signs, differentiate \(\lambda^TS\): \[ \frac{d}{dt}(\lambda^TS)=-L_x^TS+\lambda^T(F_p+F_uu_p). \] Substitute this identity into the direct derivative of \(J\) and integrate. In row-gradient form, \[ \begin{split} \frac{dJ}{dp}={}&\psi_p+\lambda(t_0)^TS(t_0)\\ &+\int_{t_0}^T\left[L_p+\lambda^TF_p+ (L_u^T+\lambda^TF_u)u_p\right]dt. \end{split} \] The terminal plus sign follows from this convention. An alternative Hamiltonian convention can change signs consistently; changing only the terminal sign breaks the identity. With \(T=T(p)\), add \((\psi_x^TF+L)T_p\) at the endpoint, and add direct \(\psi_tT_p\) if the terminal cost explicitly depends on time. Variable initial times and event resets require their corresponding boundary terms.

The adjoint is advantageous when there are many parameters but few scalar objectives. Its state dimension is \(n\), while forward sensitivity uses \(nr\) entries; parameter-gradient accumulation still has a cost. Backward evaluation also needs the forward trajectory, through storage, dense output or checkpointing. Reconstructing an unstable reverse-time trajectory without stored information can be inaccurate.

A Boundary-Term Check

In nondimensional variables, take \(\dot x=px\), \(x(0)=p\), and \(J=\tfrac12x(T)^2+\int_0^T\tfrac12x(t)^2dt\), with fixed \(T\). The exact solution is \(x(t)=pe^{pt}\) and \(S(t)=e^{pt}(1+pt)\). The adjoint obeys \(-\dot\lambda=p\lambda+x\), \(\lambda(T)=x(T)\). Its gradient is \[ \frac{dJ}{dp}=\lambda(0)+\int_0^T\lambda(t)x(t)\,dt. \] The \(\lambda(0)\) term is required because \(dx(0)/dp=1\). The accompanying program checks this against finite differences of the independently integrated objective at \(p=0.4\), \(T=0.7\). Agreement over a useful finite-difference step range checks the implementation; arbitrarily tiny steps can instead expose subtractive cancellation.

Discrete Adjoint and Backpropagation

For \(x_{k+1}=f_k(x_k,p)\) and \(J=\psi(x_N,p)+\sum_k L_k(x_k,p)\), with input dependencies included in \(f_k,L_k\), \[ \lambda_N=\psi_x,\qquad \lambda_k=(D_xf_k)^T\lambda_{k+1}+(L_k)_x. \] The parameter gradient adds \(\psi_p\), the initial-state term and \(\sum_k[(L_k)_p+\lambda_{k+1}^T D_pf_k]\). This is reverse-mode differentiation of the implemented discrete graph. It need not equal a separately discretized continuous adjoint at finite step size. Automatic differentiation differentiates the operations and branches actually executed; it does not certify the physical model, smoothness of contact events or an adaptive solver’s gradient accuracy.

Volume, Distance and Stability Are Different

Which Divergence?

For the same prescribed reference input applied to neighboring initial states, Liouville’s formula describes infinitesimal coordinate volume. Its divergence is that of \(F_0=f+G\bar u\), so \[ \operatorname{tr}A=\operatorname{div}f+ \sum_j\bar u_j\operatorname{div}g_j. \] Feedback changes the field and its divergence. For a finite region, integrate the local Jacobian determinant over the initial region; one reference determinant is not generally the exact factor for the whole region. A nonconstant volume density also changes the geometric volume formula. Hamiltonian Liouville preservation refers to the canonical symplectic volume, not every unweighted position/velocity chart.

Negative trace along an interval makes its infinitesimal coordinate-volume ratio less than 1. A bound \(\operatorname{tr}A\le-c<0\) gives the uniform rate \(e^{-c(t-s)}\); strict negativity alone need not give a uniform exponential rate. Volume decrease does not imply distance contraction or stability: \(A=\operatorname{diag}(1,-2)\) has determinant factor \(e^{-t}\) while the first direction grows as \(e^t\).

For a damped pendulum \(\ddot\theta+c\dot\theta+(g/\ell)\sin\theta=\tau/(m\ell^2)\) with \(c\) in 1/s, the trajectory Jacobian has lower-left entry \(-(g/\ell)\cos\bar\theta(t)\) and trace \(-c\). The volume factor is \(e^{-c(t-s)}\), including near the unstable upright equilibrium. Damping and shrinking phase volume do not make that equilibrium stable.

Growth in a Declared Metric

For an unforced variation in a fixed scaled Euclidean chart, \[ \frac{d}{dt}\|\xi\|^2=\xi^T(A^T+A)\xi. \] The skew part cancels from this instantaneous quadratic derivative, although its rotation of directions influences later exposure to stretching. With input, add \(2\xi^TBv\). The norm is a perturbation measure, not automatically mechanical energy. For a positive metric \(P(t)\), \[ \frac{d}{dt}(\xi^TP\xi)=\xi^T(\dot P+A^TP+PA)\xi \] in the homogeneous problem. A uniformly positive and uniformly bounded metric and an inequality \(\dot P+A^TP+PA\preceq-2\alpha P\), \(\alpha>0\), establish exponential decay for this variational system. Extending a differential inequality to a nonlinear contraction claim requires a suitable domain, metric regularity and comparison between trajectories there.

Frozen eigenvalues are insufficient for general LTV stability. To verify this explicitly, let \[ D=\begin{pmatrix}-1&4\\0&-1\end{pmatrix},\quad \Omega=\begin{pmatrix}0&1\\-1&0\end{pmatrix},\quad R(t)=e^{\Omega t},\quad A(t)=R(t)DR(t)^T. \] Every frozen \(A(t)\) has eigenvalues \(-1,-1\). But \(\eta=R(t)^T\xi\) satisfies \(\dot\eta=(D-\Omega)\eta\), whose eigenvalues are \(-1\pm\sqrt3\). One direction grows exponentially. This smooth example separates instantaneous spectra from the propagated dynamics. Even a constant Hurwitz matrix can produce transient growth in a chosen Euclidean norm if it is nonnormal.

Control Authority Over a Horizon

At one instant, \(Bv\) is an increment in the state derivative. It is not an instantaneous finite state displacement. Its column space describes possible derivative directions only with unrestricted signed input directions; actual actuator bounds restrict that set. A fully actuated mechanical system with \(d\) configuration coordinates often has \(2d\) states and a rank-\(d\) input matrix. Thus \(\operatorname{rank}B<n\) by itself does not define mechanical underactuation. Underactuation concerns independent generalized-force authority relative to the mechanical degrees of freedom.

For unrestricted square-integrable inputs in an LTV model, the terminal controllability Gramian is \[ W(T,t_0)=\int_{t_0}^T\Phi(T,s)B(s)B(s)^T\Phi(T,s)^T\,ds. \] For any \(a\), \(a^TWa=\int\|B^T\Phi^Ta\|^2ds\ge0\). Positive definiteness means no terminal direction is orthogonal to all propagated channels and is equivalent to full finite-horizon controllability of this linear problem. It does not guarantee bounded-input feasibility or global nonlinear controllability. For a positive input-cost matrix \(R_u\), replace \(BB^T\) by \(BR_u^{-1}B^T\); this makes the chosen effort weighting explicit.

For the unit double integrator over length \(T>0\), \[ W=\begin{pmatrix}T^3/3&T^2/2\\T^2/2&T\end{pmatrix}, \qquad \det W=T^4/12>0. \] Its rank-one instantaneous channel can influence both terminal position and velocity through propagation. Yet as the horizon shortens, some corrections become expensive. With \(|u|\le u_{\max}\) and zero initial error, \(|\delta q(T)|\le u_{\max}T^2/2\) and \(|\delta v(T)|\le u_{\max}T\). Controllability rank alone hides these practical limits.

Singular values require state and input scaling. For a full-column-rank \(B\), \(\kappa_2(B^TB)=\kappa_2(B)^2=(\sigma_{\max}(B)/\sigma_{\min}(B))^2\). Redundant columns make \(B^TB\) singular even when useful output authority remains. Neither a large condition number nor a low rank alone proves saturation for a specified task. Evaluate the desired output change, feasible channel set, horizon, uncertainties and retained actuator dynamics together. The finite-horizon kernel supplies a disciplined way to investigate earlier preparation versus later correction in the swing; it supplies no universal timing prescription without a validated model and measured boundary conditions.

Moving Frames and Geometric Meaning

Let \(\xi=E(t)\eta\), with a differentiable invertible basis matrix \(E\). Then \[ \dot\eta=\left(E^{-1}AE-E^{-1}\dot E\right)\eta+E^{-1}Bv. \] The basis-rate term is essential. Setting \(E=\Phi(t,t_0)\) eliminates the homogeneous coefficient, but the basis itself then carries all the original stretching and can become ill-conditioned. This coordinate choice is not a physical stabilization mechanism.

The STM is the derivative, or tangent pushforward, of the specified flow with respect to initial state. It need not be parallel transport for the kinetic metric or any Levi-Civita connection. On the Euclidean line, \(\dot x=ax\) gives a flow derivative \(e^{a(t-s)}\), whereas Euclidean parallel transport preserves a vector’s length. A connection must be specified before making a parallel-transport claim. Under a smooth coordinate change \(y=\chi(x)\), the same flow derivative has matrix \(D\chi(\bar x(t))\Phi(t,s)D\chi(\bar x(s))^{-1}\). Finite subtraction on a manifold is chart-dependent; tangent vectors and finite displacements should not be interchanged without that declaration.

Reliable Numerical Computation

Integrating What the Task Requires

Integrating the full STM adds \(n^2\) state variables. A dense product \(A\Phi\) costs \(O(n^3)\) per right-hand-side evaluation, plus the cost of obtaining \(A\). If only a few initial directions are needed, propagate an \(n\times r\) sensitivity block instead. If only a few scalar output gradients are needed, an adjoint may be preferable. Sparse products, Jacobian-vector products and matrix-exponential actions can avoid forming dense matrices. No single method is uniformly cheapest across dimensions, structures, output requests and tolerances.

Analytic derivatives and automatic differentiation should be checked against independent directional differences where feasible. A derivative can faithfully differentiate the wrong force law. Integrate the reference state and sensitivities consistently, use scales appropriate to each block, refine tolerances or the time step, and compare the final task output and derivative. The MIT trajectory-linearization treatment connects these derivatives to local trajectory feedback while keeping the reference feasible.

Matrix Exponentials and Conditioning

For constant coefficients, scaling and squaring approximates \(e^Z\) through \[ e^Z\approx\left[r_m(2^{-s}Z)\right]^{2^s}, \qquad r_m=N_m/D_m. \] Modern algorithms choose a diagonal \([m/m]\) Padé degree and scaling according to error estimates, then solve the matrix-polynomial system rather than forming its denominator inverse. The Al-Mohy–Higham analysis explains how unnecessary scaling can worsen rounding error. A fixed \([7/8]\) rule and threshold are not a general guarantee of machine accuracy. Backward error, conditioning of the matrix-exponential problem and roundoff together determine forward accuracy; consult the numerical library’s documented algorithm and verify the application.

For a dense fixed-degree approximant, polynomial evaluation and a solve have cubic scaling, with additional cubic work for each squaring. Constants, selected degree, structure and reuse matter. An adaptive ODE solver’s cost depends on accepted and rejected steps, not simply the number of requested output samples. For time-varying coefficients, repeatedly applying frozen exponentials also introduces a time-discretization error.

Stiffness Is Not the Same as Loss of Information

Stiffness is a numerical stability restriction relative to the solution and accuracy sought. Widely separated decaying modes can create it; an eigenvalue ratio is a diagnostic in some problems, not a universal definition. An explicit solver may need a short step for stability after the fast transient is already small. Suitable implicit methods address that restriction, but do not restore information lost through extreme expansion or contraction.

For \(A=\operatorname{diag}(-1,-1000)\), \(\Phi=\operatorname{diag}(e^{-t},e^{-1000t})\) has condition number \(e^{999t}\) in these scaled coordinates. In contrast, \(\Phi=e^{-100t}I\) has condition number 1 even when every entry is very small. The smallest singular value falling below an absolute machine-epsilon threshold is not by itself a singularity criterion. Relative singular-value separation, absolute error scales, underflow and the intended operation all matter. Inversion or recovery of suppressed directions can be unreliable long before a forward output of interest becomes unusable.

QR Propagation With Correct Product Order

Suppose \(\Phi_k=Q_kR_k\) and a computed interval propagator is \(P_k\). Form \(Y=P_kQ_k\) and factor \(Y=Q_{k+1}S_k\). Then \[ \Phi_{k+1}=P_k\Phi_k=Q_{k+1}S_kR_k, \qquad R_{k+1}=S_kR_k. \] The new triangular factor multiplies on the left. During the interval, one propagates \(Y\) with \(\dot Y=AY\) from an orthonormal starting basis; \(Y\) is not orthonormal between QR steps. The accumulated \(R\) carries the stretching and is not promised to remain bounded or well-conditioned. QR can stabilize basis propagation and support logarithmic growth diagnostics, but reconstructing the full product can still overflow, underflow or lose relative information. Keep scaled factors or checkpointed short transitions when appropriate. Restarting a local transition does not magically recover information already lost from the original long-horizon map.

Structure-Preserving Integration

Symplectic Structure and Energy

For an autonomous unforced canonical Hamiltonian system, the exact flow preserves its Hamiltonian and symplectic form. A symplectic numerical map satisfies \((D\Psi_h)^T\mathbb J D\Psi_h=\mathbb J\) and hence preserves canonical volume. Volume preservation alone is weaker in dimensions above two. A generic explicit RK4 method is not symplectic, but implicit midpoint and Gauss Runge–Kutta methods are symplectic under their defining assumptions; Hairer’s symplectic-integration lectures derive these examples.

For separable \(H(q,p)=T(p)+V(q)\), kick-drift symplectic Euler gives \[ p_{k+1}=p_k-h\nabla V(q_k),\qquad q_{k+1}=q_k+h\nabla T(p_{k+1}). \] It preserves symplectic structure in exact arithmetic but does not conserve the original energy exactly or remain stable for every step. For the unit oscillator with state order \((q,p)\), its matrix is \[ S_h=\begin{pmatrix}1-h^2&h\\-h&1\end{pmatrix}, \qquad \det S_h=1,\quad\operatorname{tr}S_h=2-h^2. \] When \(0<h<2\), its eigenvalues lie on the unit circle and it preserves the positive quadratic \(\tfrac12(q^2+p^2-hqp)\). At \(h>2\) one eigenvalue has magnitude greater than 1, so energy can grow without bound despite symplecticity. At \(h=2\) the repeated eigenvalue needs its Jordan structure checked. Long-time near-energy conservation requires suitable small-step, regularity and bounded-trajectory assumptions; it is not an all-time, all-step guarantee.

Numerical energy error need not drift monotonically or always produce an outward spiral. A short downswing is not an automatic accuracy certificate, and independently restarted parameter sweeps do not accumulate integration time merely because many runs are executed. Check force balance, convergence and task outputs at the actual time scales. Forced, dissipative or impact dynamics need their own work and jump balances rather than an inappropriate energy-conservation test.

Rotations and Variational Integrators

For a body-to-space rotation \(R\) and body-resolved angular velocity \(\omega_b\), \(\dot R=R\widehat\omega_b\), where \(\widehat\omega_b z=\omega_b\times z\). A frozen-velocity update \(R_{k+1}=R_k\exp(h\widehat\omega_{b,k})\) stays on \(SO(3)\) in exact arithmetic and is first order for general time-varying velocity. Space-resolved velocity uses left multiplication. Higher order needs appropriate stage and noncommutativity treatment; simply integrating three angular components and exponentiating does not automatically supply it. Floating-point matrix products can still require monitoring of orthogonality. Euler-angle chart singularities and quaternion normalization drift are different issues.

A discrete action \(\sum_k L_d(q_k,q_{k+1})\), varied with fixed endpoints, gives \[ D_2L_d(q_{k-1},q_k)+D_1L_d(q_k,q_{k+1})=0. \] For a regular discrete Lagrangian these equations define a local map. Its discrete Legendre transforms \(p_k=-D_1L_d(q_k,q_{k+1})\) and \(p_{k+1}=D_2L_d(q_k,q_{k+1})\) place it in canonical phase space, where the unforced map is symplectic. Momentum preservation follows for symmetries actually retained by the discrete Lagrangian; see West’s variational-integrator account. Forced, constrained and nonsmooth extensions require their corresponding formulations. Differentiating a symplectic map preserves its tangent symplectic identity in canonical coordinates, but does not guarantee faster optimization or make a generic position/velocity Jacobian canonically symplectic.

Reproducible Checks

The standalone Python program uses NumPy and SciPy. Its functions expose the pendulum transition, singular-safe ZOH sensitivity, ordered QR product, scalar remainder, boundary-sensitive adjoint and symplectic-Euler matrix. The published table values are generated from it. Independent tests compare its pendulum result with matrix-exponential and ODE calculations, its QR result with direct noncommuting products, its ZOH result with a polynomial solution and its adjoint with finite differences of a separately integrated objective.

# Variational chapter example
import numpy as np
from scipy.integrate import quad
from scipy.linalg import expm

GRAVITY_M_S2 = 10.0
LENGTH_M = 1.0
MASS_KG = 1.0


def pendulum_stm(time: float) -> np.ndarray:
    """Downward point-mass pendulum, state (angle, angular velocity)."""
    frequency = np.sqrt(GRAVITY_M_S2 / LENGTH_M)
    cosine, sine = np.cos(frequency * time), np.sin(frequency * time)
    return np.array([[cosine, sine / frequency], [-frequency * sine, cosine]])


def zoh(matrix: np.ndarray, channel: np.ndarray, interval: float) -> tuple[np.ndarray, np.ndarray]:
    """Exact linear held-input blocks, including singular state matrices."""
    states, inputs = channel.shape
    augmented = np.zeros((states + inputs, states + inputs))
    augmented[:states, :states] = matrix
    augmented[:states, states:] = channel
    result = expm(interval * augmented)
    return result[:states, :states], result[:states, states:]


def qr_product(intervals: list[np.ndarray]) -> tuple[np.ndarray, np.ndarray]:
    """Finite illustrative product; accumulated factors can be ill-conditioned."""
    if not intervals:
        raise ValueError("Supply at least one square interval propagator")
    orthogonal = np.eye(intervals[0].shape[0])
    triangular = orthogonal.copy()
    for interval in intervals:
        orthogonal, factor = np.linalg.qr(interval @ orthogonal)
        triangular = factor @ triangular
    return orthogonal, triangular


def joint_remainder(displacement: float, increment: float) -> float:
    """Exact remainder for F(x,u)=x*x+(1+x)*u at a frozen pair."""
    return displacement**2 + displacement * increment


def scalar_adjoint_gradient(parameter: float, terminal: float) -> float:
    """x'=p*x, x(0)=p; half squared terminal and running state cost."""

    def state(time: float) -> float:
        return float(parameter * np.exp(parameter * time))

    def adjoint(time: float) -> float:
        forced = quad(
            lambda moment: np.exp(parameter * (moment - time)) * state(moment),
            time,
            terminal,
        )[0]
        return float(np.exp(parameter * (terminal - time)) * state(terminal) + forced)

    return adjoint(0.0) + quad(lambda time: adjoint(time) * state(time), 0.0, terminal)[0]


def symplectic_euler_matrix(step: float) -> np.ndarray:
    """Unit oscillator, kick then drift, state (q,p)."""
    return np.array([[1 - step**2, step], [-step, 1]])


initial_response = pendulum_stm(0.5) @ np.array([0.1, 0.0])
impulse_response = pendulum_stm(0.25) @ np.array([0.0, 1.0 / (MASS_KG * LENGTH_M**2)])
if __name__ == "__main__":
    for time in [0.25, 0.5]:
        print(time, pendulum_stm(time))
    print("Initial response:", initial_response)
    print("Impulse response:", impulse_response)
    print("Adjoint gradient:", scalar_adjoint_gradient(0.4, 0.7))

What This Establishes for the Swing

Variational analysis connects preparation and correction through a specified trajectory, retained state, input port and output. Earlier effort may establish velocities or elastic and activation states that shape later motion. A later input can still affect some impact variables while having little authority over others. Constraint reactions can redirect acceleration without supplying total mechanical work; actuator dynamics can delay or filter a command. These mechanisms enter the reference dynamics and its derivatives, rather than being inferred from a single speed peak.

To investigate a proposed release mechanism, first declare the model and event, then identify measured and inferred states, compute the appropriate finite-horizon output sensitivities, impose feasible inputs and uncertainties, and check representative finite nonlinear interventions. Compare those predictions with data that distinguish competing mechanisms. A locally sensitive direction does not prove human causation, and a locally insensitive output does not prove that earlier muscular work was irrelevant. The mathematical tools sharpen the research question while keeping model evidence separate from observations of a golfer.

Exercises

1. A Joint Taylor Check. For \(F(x,u)=x^2+(1+x)u\), derive the exact finite-error equation about a feasible reference. Compare equal scaling of state/input perturbations with a fixed finite input change as the state displacement tends to zero.
2. Ordered Propagation. Compute both orders of the two nilpotent interval generators above. Compare with \(\exp(A_1+A_2)\) and identify the commutator term responsible for the discrepancy.
3. A Moving Stable Spectrum. Differentiate \(\Phi(t,0)=R(t)\exp((D-\Omega)t)\) for the rotating example and verify its initial condition. Show that one solution grows although every frozen eigenvalue is negative.
4. Pendulum Limits. Using \(m=\ell=1\) and \(g=10\), compare the nonlinear pendulum and first variation from angles 0.1, 0.05 and 0.025 rad over 0.5 s. State the norm and integration tolerance; investigate the special odd-symmetry error order at this equilibrium.
5. Bounded Pulses. Replace the ideal torque impulse by a rectangular pulse with the same angular impulse. Vary duration and enforce a declared torque limit. Separate the short-pulse approximation from actuator feasibility.
6. The Initial-State Gradient. Reproduce the scalar adjoint check with \(x(0)=p\), then repeat with fixed \(x(0)\). Identify the changed boundary term and verify gradients across several finite-difference steps.
7. A Short-Horizon Task. Derive the double-integrator Gramian and amplitude bounds. Compare the achievable terminal position and velocity corrections as \(T\) decreases; explain why full rank does not imply a uniformly easy correction.
8. Volume Versus Distance. For \(A=\operatorname{diag}(1,-2)\), compute the determinant, singular values and response from each coordinate direction. Repeat for a damped pendulum at its upright equilibrium.
9. A QR Product. Implement the interval factorization and compare the correct left accumulation with the reversed order using noncommuting propagators. Monitor reconstruction error and factor magnitudes as the horizon grows.
10. Symplecticity Versus Stability. Verify the unit-oscillator symplectic identity and modified quadratic for symplectic Euler. Compare \(h=0.5\), 2 and 3, including the repeated-eigenvalue case. Do not infer stability from determinant 1.
11. A Golf Output Contract. Choose one pre-impact outcome, declare its units, frame, event and retained actuator states, and write its first-order input kernel. Describe the measurements needed to test a timing claim and a finite intervention that could falsify it.

Cite this page