Variational Dynamics and the Moving Frame
Even though a nonlinear system curves and twists through state space, tiny errors around a reference trajectory evolve by a linear rule that we can compute and exploit. This chapter derives that rule—the variational equation—from first principles, constructs its exact solution via the Peano-Baker series, and shows how it provides a “moving local coordinate system” that travels along with the nominal motion. We examine rigorous proofs of its fundamental properties, connect it to sensitivity analysis and optimal control, and confront the numerical challenges in computing it accurately.
Here is the key idea of this chapter in everyday language.
Imagine you are a pilot following a planned flight path from New York to London. Before takeoff, you computed the exact route, the exact fuel burn, the exact timing. That is your nominal trajectory—the plan.
Now, the real flight will not follow the plan exactly. Wind gusts push you off course. The question is: if you start slightly off-plan (maybe 100 meters north of where you expected to be), how does that error evolve over time? Does it grow? Shrink? Rotate?
The variational equation answers this question precisely. It is a linear rule that tells you, at each instant, how small deviations from the plan evolve. The rule changes as you move along the flight path (different winds at different altitudes), but at each instant it is a simple, linear relationship: twice the error leads to twice the deviation, and errors in different directions combine by addition.
This is enormously powerful. It means you can analyze the sensitivity of any trajectory without having to re-simulate the entire system for every possible perturbation. One matrix—the state transition matrix—captures all possible small deviations at once.
From Global Nonlinearity to Local Linearity
Setup: Nominal and Perturbed Trajectories
Given the control-affine dynamics \[\begin{equation} \dot{\x} = f(\x) + G(\x)\,\uvec, \quad \x\in\R^n,\;\uvec\in\R^m, \label{eq:ch2:affine} \end{equation}\] where \(f:\R^n\to\R^n\) and \(G:\R^n\to\R^{n\times m}\) are smooth, suppose we have a nominal trajectory \(\bar{\x}(t)\) generated by a nominal input \(\bar{\uvec}(t)\): \[\begin{equation} \dot{\bar{\x}}(t) = f\!\bigl(\bar{\x}(t)\bigr) + G\!\bigl(\bar{\x}(t)\bigr)\,\bar{\uvec}(t). \label{eq:ch2:nominal} \end{equation}\] Consider a nearby trajectory \(\x(t) = \bar{\x}(t) + \dx(t)\) under a perturbed input \(\uvec(t) = \bar{\uvec}(t) + \du(t)\) with \(\|\dx\|\) small. Then \(\x\) satisfies \[\begin{equation} \dot{\bar{\x}} + \dot{\dx} = f(\bar{\x}+\dx) + G(\bar{\x}+\dx)(\bar{\uvec}+\du). \label{eq:ch2:perturbed} \end{equation}\]
Derivation of the Variational Equation From First Principles
We expand the right-hand side of~\(\eqref{eq:ch2:perturbed}\) in a Taylor series about \((\bar{\x}, \bar{\uvec})\).
Expansion of \(f(\bar{\x}+\dx)\)
Using the multivariate Taylor theorem, \[\begin{equation} f(\bar{\x}+\dx) = f(\bar{\x}) + \frac{\partial f}{\partial \x}\bigg|_{\bar{\x}}\,\dx + \frac{1}{2}\int_0^1(1-\tau)\frac{\partial^2 f}{\partial\x^2}\bigg|_{\bar{\x}+\tau\dx}(\dx,\dx)\,\dd\tau. \label{eq:ch2:f-taylor} \end{equation}\]
The first derivative is the Jacobian \(J_f(\bar{\x}) := \partial f/\partial\x|_{\bar{\x}}\), which is an \(n\times n\) matrix. The second derivative term represents the Hessian applied to the perturbation: for each state component \(i\), we have \[ \left[\frac{\partial^2 f_i}{\partial\x^2}(\dx,\dx)\right]_i = \sum_{j,k}\frac{\partial^2 f_i}{\partial x_j\partial x_k}\dx_j\,\dx_k. \] The integral bound on the second derivative is \(O(\|\dx\|^2)\). Let us denote \[ R_f(\dx) := \frac{1}{2}\int_0^1(1-\tau)\frac{\partial^2 f}{\partial\x^2}\bigg|_{\bar{\x}+\tau\dx}(\dx,\dx)\,\dd\tau. \] Then \[\begin{equation} f(\bar{\x}+\dx) = f(\bar{\x}) + J_f(\bar{\x}))\dx + R_f(\dx), \quad \|R_f(\dx)\| = O(\|\dx\|^2). \label{eq:ch2:f-expanded} \end{equation}\]
Expansion of \(G(\bar{\x}+\dx)(\bar{\uvec}+\du)\)
Similarly, expand \(G(\bar{\x}+\dx)\) as \[\begin{equation} G(\bar{\x}+\dx) = G(\bar{\x}) + \frac{\partial G}{\partial \x}\bigg|_{\bar{\x}}\,\dx + O(\|\dx\|^2). \label{eq:ch2:G-taylor} \end{equation}\] Thus \[\begin{align} G(\bar{\x}+\dx)(\bar{\uvec}+\du) &= \left[G(\bar{\x}) + \frac{\partial G}{\partial \x}\bigg|_{\bar{\x}}\,\dx + O(\|\dx\|^2)\right](\bar{\uvec}+\du)\\ &= G(\bar{\x})\bar{\uvec} + G(\bar{\x})\du + \frac{\partial G}{\partial \x}\bigg|_{\bar{\x}}(\dx\otimes\bar{\uvec}) + O(\|\dx\|^2) + O(\|\dx\|\|\du\|). \label{eq:ch2:G-expanded} \end{align}\]
Here the term \(\frac{\partial G}{\partial \x}(\dx\otimes\bar{\uvec})\) involves the directional derivative of \(G\) with respect to \(\x\) applied to \(\dx\), scaled by \(\bar{\uvec}\). In index form, the \(i\)-th component is \(\sum_{j,k}\frac{\partial G_{ij}}{\partial x_k}\dx_k\,\bar u_j\).
Combining Terms
Substituting~\(\eqref{eq:ch2:f-expanded}\) and~\(\eqref{eq:ch2:G-expanded}\) into~\(\eqref{eq:ch2:perturbed}\): \[\begin{align} \dot{\bar{\x}} + \dot{\dx} &= f(\bar{\x}) + J_f(\bar{\x})\dx + R_f(\dx)\\ &\quad+ G(\bar{\x})\bar{\uvec} + G(\bar{\x})\du + \frac{\partial G}{\partial \x}\bigg|_{\bar{\x}}(\dx\otimes\bar{\uvec}) + O(\|\dx\|^2). \end{align}\] Subtracting the nominal equation~\(\eqref{eq:ch2:nominal}\): \[\begin{equation} \dot{\dx} = J_f(\bar{\x})\dx + \frac{\partial G}{\partial \x}\bigg|_{\bar{\x}}(\dx\otimes\bar{\uvec}) + G(\bar{\x})\du + R_f(\dx) + O(\|\dx\|^2). \label{eq:ch2:diff-before-drop} \end{equation}\]
Dropping Higher-Order Terms
To obtain the linearized (first-order) equation, we retain only terms linear in \(\dx\) and \(\du\) and drop all \(O(\|\dx\|^2)\) terms: \[\begin{equation} \dot{\dx} = J_f(\bar{\x})\dx + \frac{\partial G}{\partial \x}\bigg|_{\bar{\x}}(\dx\otimes\bar{\uvec}) + G(\bar{\x})\du. \label{eq:ch2:variational-unsimplified} \end{equation}\]
The \(R_f(\dx)\) and \(O(\|\dx\|^2)\) terms represent curvature of the vector field. They vanish in the limit \(\|\dx\|\to 0\), making them negligible for infinitesimal perturbations. Note that the term \(\frac{\partial G}{\partial \x}(\dx\otimes\bar{\uvec})\) is linear in \(\dx\) (when \(\bar{\uvec}\) is fixed along the nominal trajectory), so it is retained.
Compact Matrix Form
Define \[\begin{equation} \mat{A}(t) := \frac{\partial}{\partial \x}\bigl[f(\x)+G(\x)\bar{\uvec}(t)\bigr]\bigg|_{\x=\bar{\x}(t)}, \qquad \mat{B}(t) := G\!\bigl(\bar{\x}(t)\bigr). \label{eq:ch2:AB-def} \end{equation}\]
The matrix \(\mat{A}(t)\) is the Jacobian of the vector field \(f + G\bar{\uvec}\) along the nominal trajectory. Explicitly, \[\begin{equation} \mat{A}(t) = J_f(\bar{\x}(t)) + \frac{\partial G}{\partial\x}\bigg|_{\bar{\x}(t)} \, (\mathbf{1} \otimes \bar{\uvec}(t))^T, \end{equation}\] where the second term arises from the control input contribution to the state sensitivity. Then~\(\eqref{eq:ch2:variational-unsimplified}\) becomes:
\[\begin{equation} \boxed{ \dot{\dx}(t) = \mat{A}(t)\,\dx(t) + \mat{B}(t)\,\du(t), } \label{eq:ch2:variational} \end{equation}\]
This is the variational equation: a linear, time-varying ODE that governs the evolution of infinitesimal perturbations.
The error incurred by dropping the \(O(\|\dx\|^2)\) terms in the derivation is the residual. For initial conditions with \(\|\dx(t_0)\| \ll 1\), this residual remains small (growing only quadratically with \(\|\dx\|\)) over a finite time interval. This is the precise justification for using linearization to study local stability and design feedback control.
The variational equation~\(\eqref{eq:ch2:variational}\) is a linear, time-varying system that lives in the tangent space along \(\bar{\x}(t)\). It describes how infinitesimal perturbations propagate. The matrices \(\mat{A}(t)\) and \(\mat{B}(t)\) change as the nominal trajectory evolves, reflecting the fact that the tangent space itself moves with the base point. Notice that \(\mat{A}\) depends on the nominal input \(\bar{\uvec}(t)\): different reference inputs generate different tangent-space dynamics!
The State Transition Matrix
Definition and Fundamental ODE
When \(\du = 0\), the variational equation reduces to the homogeneous system \[\begin{equation} \dot{\dx} = \mat{A}(t)\,\dx. \label{eq:ch2:homogeneous} \end{equation}\] The solution is characterized by the state transition matrix (STM) \(\Phi(t,t_0)\), defined as the unique \(n\times n\) matrix satisfying \[\begin{equation} \boxed{ \dot\Phi(t,t_0) = \mat{A}(t)\,\Phi(t,t_0), \quad \Phi(t_0,t_0) = \mat{I}_n. } \label{eq:ch2:stm-ode} \end{equation}\]
Once \(\Phi(t,t_0)\) is known, the solution to~\(\eqref{eq:ch2:homogeneous}\) with initial condition \(\dx(t_0)\) is \[\begin{equation} \dx(t) = \Phi(t,t_0)\,\dx(t_0). \label{eq:ch2:stm-solution} \end{equation}\]
The STM is not a solution in the usual sense; it is a fundamental solution matrix. Each column of \(\Phi(t,t_0)\) is a solution to the homogeneous equation with a unit initial condition. Specifically, if \(\dx_i(t_0) = \mathbf{e}_i\) (the \(i\)-th unit vector), then \(\dx_i(t) = \Phi(t,t_0)\mathbf{e}_i\) is the \(i\)-th column of \(\Phi(t,t_0)\).
Properties and Rigorous Proofs
The state transition matrix inherits fundamental properties from the structure of linear ODEs:
Semigroup Property
For any \(t_0 \leq t_1 \leq t_2\), \[\begin{equation} \Phi(t_2,t_0) = \Phi(t_2,t_1)\,\Phi(t_1,t_0). \label{eq:ch2:semigroup} \end{equation}\]
Proof. Let \(\Psi(t) := \Phi(t,t_1)\,\Phi(t_1,t_0)\). We compute \[\begin{align} \dot{\Psi}(t) &= \frac{\dd}{\dd t}\left[\Phi(t,t_1)\,\Phi(t_1,t_0)\right]\\ &= \mat{A}(t)\,\Phi(t,t_1)\,\Phi(t_1,t_0) = \mat{A}(t)\,\Psi(t). \end{align}\]
At \(t=t_1\), we have \(\Psi(t_1) = \Phi(t_1,t_1)\,\Phi(t_1,t_0) = \mat{I}_n\,\Phi(t_1,t_0) = \Phi(t_1,t_0)\).
Both \(\Psi(t)\) and \(\Phi(t,t_0)\) satisfy the same linear ODE \(\dot{\mat{X}} = \mat{A}(t)\mat{X}\). By uniqueness of solutions to linear ODEs, if they are equal at \(t_1\), they are equal everywhere. Thus \(\Psi(t) = \Phi(t,t_0)\). Evaluating at \(t = t_2\) yields \(\Phi(t_2,t_0) = \Phi(t_2,t_1)\,\Phi(t_1,t_0)\), completing the proof.
Invertibility
The state transition matrix is invertible for all \(t, t_0\), and \[\begin{equation} \Phi(t,t_0)^{-1} = \Phi(t_0,t). \label{eq:ch2:inverse} \end{equation}\]
Proof. From the semigroup property, \[ \Phi(t,t_0)\,\Phi(t_0,t) = \Phi(t,t) = \mat{I}_n. \] Thus \(\Phi(t_0,t)\) is indeed the right inverse. Similarly, from \(\Phi(t_0,t)\,\Phi(t,t_0) = \Phi(t_0,t_0) = \mat{I}_n\), it is also the left inverse. Hence \(\Phi(t,t_0)^{-1} = \Phi(t_0,t)\).
Liouville’s Formula
The determinant of the state transition matrix satisfies \[\begin{equation} \det\Phi(t,t_0) = \exp\!\left(\int_{t_0}^{t}\trace\mat{A}(s)\,\dd s\right). \label{eq:ch2:liouville} \end{equation}\]
Proof. We use the matrix determinant lemma (also called Jacobi’s formula): \[\begin{equation} \frac{\dd}{\dd t}\det(\Phi(t,t_0)) = \trace\left(\text{adj}(\Phi)\,\frac{\dd\Phi}{\dd t}\right), \label{eq:ch2:jacobi} \end{equation}\] where \(\text{adj}(\Phi)\) is the adjugate matrix. For a nonsingular matrix, \(\text{adj}(\Phi) = \det(\Phi)\,\Phi^{-T}\), so \[\begin{align} \frac{\dd}{\dd t}\det(\Phi) &= \det(\Phi)\trace(\Phi^{-T}\,\dot{\Phi})\\ &= \det(\Phi)\trace(\Phi^{-T}\,\mat{A}\,\Phi)\\ &= \det(\Phi)\trace(\mat{A}\,\Phi\,\Phi^{-T})\\ &= \det(\Phi)\trace(\mat{A}). \end{align}\]
Here we used the cyclic property of the trace: \(\trace(ABC) = \trace(BCA) = \trace(CAB)\).
Integrating from \(t_0\) to \(t\), with initial condition \(\det(\Phi(t_0,t_0)) = \det(\mat{I}_n) = 1\): \[\begin{equation} \det(\Phi(t,t_0)) = \exp\!\left(\int_{t_0}^{t}\trace(\mat{A}(s))\,\dd s\right). \end{equation}\]
The state transition matrix \(\Phi(t,t_0)\) accumulates exact infinitesimal linear contributions along the trajectory. Each infinitesimal step \(\dd(\dx)=\mat{A}(t)\dx\,\dd t\) is exact in its tangent space. The integral \[ \dx(t_1) = \Phi(t_1,t_0)\dx(t_0) + \int_{t_0}^{t_1}\Phi(t_1,s)\mat{B}(s)\du(s)\,\dd s \] is the exact accumulation of these exact local effects. Residuals arise from curvature (when we compare the true finite-displacement trajectory to the linear prediction), not from integration itself.
The Peano-Baker Series and Exact Solution
While the ODE~\(\eqref{eq:ch2:stm-ode}\) defines the STM implicitly, we can construct an explicit series solution called the Peano-Baker series (also known as the Dyson series in physics).
Derivation of the Series
Suppose we have a formal integral representation of~\(\eqref{eq:ch2:stm-ode}\): \[\begin{equation} \Phi(t,t_0) = \mat{I} + \int_{t_0}^{t}\mat{A}(s_1)\,\Phi(s_1,t_0)\,\dd s_1. \label{eq:ch2:stm-integral} \end{equation}\]
Verify: differentiating both sides with respect to \(t\) gives \(\dot{\Phi} = \mat{A}(t)\Phi(t,t_0)\), and at \(t=t_0\), \(\Phi(t_0,t_0) = \mat{I}\). This integral equation is equivalent to the ODE.
Now, substitute the integral equation into itself iteratively: \[\begin{align} \Phi(t,t_0) &= \mat{I} + \int_{t_0}^{t}\mat{A}(s_1)\left[\mat{I} + \int_{t_0}^{s_1}\mat{A}(s_2)\Phi(s_2,t_0)\,\dd s_2\right]\dd s_1\\ &= \mat{I} + \int_{t_0}^{t}\mat{A}(s_1)\,\dd s_1\\ &\quad+ \int_{t_0}^{t}\mat{A}(s_1)\int_{t_0}^{s_1}\mat{A}(s_2)\Phi(s_2,t_0)\,\dd s_2\,\dd s_1. \end{align}\]
Continuing this recursion, we obtain the formal series:
\[\begin{equation} \boxed{ \Phi(t,t_0) = \sum_{k=0}^{\infty}\Phi_k(t,t_0), } \label{eq:ch2:peano-baker} \end{equation}\]
where \[\begin{align} \Phi_0(t,t_0) &:= \mat{I},\\ \Phi_1(t,t_0) &:= \int_{t_0}^{t}\mat{A}(s_1)\,\dd s_1,\\ \Phi_2(t,t_0) &:= \int_{t_0}^{t}\int_{t_0}^{s_1}\mat{A}(s_1)\mat{A}(s_2)\,\dd s_2\,\dd s_1,\\ \Phi_k(t,t_0) &:= \int_{t_0}^{t}\int_{t_0}^{s_1}\cdots\int_{t_0}^{s_{k-1}}\mat{A}(s_1)\mat{A}(s_2)\cdots\mat{A}(s_k)\,\dd s_k\cdots\dd s_1. \label{eq:ch2:peano-terms} \end{align}\]
This is the Peano-Baker series: an infinite product series where \(\Phi_k\) represents the cumulative effect of \(k\) applications of \(\mat{A}\), with all possible orderings integrated over the time interval.
Convergence
For time-varying linear systems with bounded coefficients, the series converges uniformly:
Peano-Baker Convergence
Let \(\mat{A}(t)\) be piecewise continuous with \(\|\mat{A}(t)\| \leq M\) for all \(t\in[t_0,T]\). Then the Peano-Baker series converges absolutely and uniformly to the unique solution of~\(\eqref{eq:ch2:stm-ode}\).
Proof. We bound the series term by term. For the \(k\)-th term: \[\begin{align} \|\Phi_k(t,t_0)\| &\leq \int_{t_0}^{t}\int_{t_0}^{s_1}\cdots\int_{t_0}^{s_{k-1}}\|\mat{A}(s_1)\|\cdots\|\mat{A}(s_k)\|\,\dd s_k\cdots\dd s_1\\ &\leq M^k \int_{t_0}^{t}\int_{t_0}^{s_1}\cdots\int_{t_0}^{s_{k-1}}\,\dd s_k\cdots\dd s_1\\ &= M^k \frac{(t-t_0)^k}{k!}. \end{align}\]
Thus the series is bounded by the exponential series: \[ \sum_{k=0}^{\infty}\|\Phi_k(t,t_0)\| \leq \sum_{k=0}^{\infty}\frac{M^k(t-t_0)^k}{k!} = e^{M(t-t_0)}, \] which is finite and independent of the ordering in \(\mat{A}\). Hence the series converges uniformly on any finite interval \([t_0,T]\), and the sum satisfies the ODE~\(\eqref{eq:ch2:stm-ode}\).
The Peano-Baker series highlights a crucial feature: even though \([\mat{A}(t_1),\mat{A}(t_2)] \neq 0\) in general (the matrices do not commute), the series still provides the exact solution. The summation over all orderings effectively accounts for non-commutativity. This is in stark contrast to the Magnus expansion from perturbation theory, which organizes the series differently to emphasize commutators.
Numerical Computation of the State Transition Matrix
For practical implementation, we rarely sum the Peano-Baker series explicitly. Instead, we integrate the matrix ODE~\(\eqref{eq:ch2:stm-ode}\) numerically.
Direct Integration of the Matrix ODE
The most straightforward approach is to view \(\Phi(t,t_0)\) as an \(n\times n\) matrix and integrate \[\begin{equation} \dot{\Phi} = \mat{A}(t)\,\Phi, \quad \Phi(t_0) = \mat{I}_n, \label{eq:ch2:stm-direct} \end{equation}\] using a standard ODE solver. This requires integrating \(n^2\) coupled scalar ODEs (the entries of \(\Phi\)).
At each time step, we compute the Jacobian \(\mat{A}(t)\) by evaluating the partial derivatives of the nominal dynamics at \((\bar{\x}(t),\bar{\uvec}(t))\). This requires either automatic differentiation, manual Jacobian formulas, or finite-difference approximations of the derivatives. In production code, automatic differentiation is often most reliable and accurate.
The Matrix Exponential and Padé Approximation
For time-invariant systems (constant \(\mat{A}\)), the STM has a closed form: \[\begin{equation} \Phi(t,t_0) = \exp\!\bigl[\mat{A}(t-t_0)\bigr]. \label{eq:ch2:stm-exponential} \end{equation}\]
The matrix exponential \(\exp(\mat{A}\tau) = \sum_{k=0}^{\infty}\frac{(\mat{A}\tau)^k}{k!}\) is a fundamental object in control theory. Computing it accurately is essential but non-trivial.
Padé Approximation
The Padé rational approximation \([\ell/m]\) of \(\exp(x)\) is the ratio of two polynomials that matches the Taylor series up to order \(\ell+m\): \[\begin{equation} [\ell/m]_{\exp}(x) = \frac{N_\ell(x)}{D_m(x)}, \quad N_\ell, D_m \text{ are polynomials of degree } \ell, m. \label{eq:ch2:pade} \end{equation}\]
For the matrix exponential, we compute \(N_\ell(\mat{A}\tau)\) and \(D_m(\mat{A}\tau)\) as matrix polynomials and then solve the linear system \[\begin{equation} D_m(\mat{A}\tau)\,\Phi(t,t_0) = N_\ell(\mat{A}\tau), \label{eq:ch2:pade-solve} \end{equation}\] which gives \(\Phi\) with spectral accuracy.
A common choice is the \([7/8]\) Padé approximation with scaling and squaring: 1. Scale: Find the smallest integer \(q\) such that \(\|\mat{A}\tau / 2^q\| < \rho\) for some threshold \(\rho\) (typically \(0.5\)). 2. Compute: Apply the \([7/8]\) Padé approximant to \(\mat{A}\tau/2^q\). 3. Square: Repeatedly square the result \(q\) times to recover \(\exp(\mat{A}\tau)\).
This approach achieves unit roundoff accuracy in double precision with minimal computational cost.
Computational Cost
For an \(n\times n\) STM: - Direct integration: Each step of an RK4 integrator performs \(\sim 4\) matrix multiplications (one per stage), each costing \(O(n^3)\). For \(N\) steps: \(O(4Nn^3)\) flops. - Matrix exponential (scaling & squaring): Computing \(N_\ell(\mat{A})\) costs \(O(n^3)\); solving the linear system costs \(O(n^3)\) (via LU factorization). Squaring \(q\) times adds \(qO(n^3)\). Total per call: \(O((q+2)n^3)\) flops. - Evaluation cost per time step: Regardless of method, if we evaluate \(\Phi\) at \(N\) time points, the dominant cost is \(O(Nn^3)\).
For large \(n\) (e.g., \(n=100\)), the \(O(n^3)\) factor is significant. Specialized structure in \(\mat{A}\) (sparse, low-rank, block-diagonal) can be exploited to reduce cost. For general dense problems, direct integration with adaptive step-size control is often most efficient in practice.
Computing \(\Phi\) for stiff systems (those with widely separated time scales) is notoriously ill-conditioned. We discuss remedies in Section~\(\ref{sec:ch2:stability}\).
Variation of Constants Formula
When both initial perturbations \(\dx(t_0) \neq 0\) and control perturbations \(\du(t) \neq 0\) are present, the complete solution of the variational equation is given by the variation of constants (also called Duhamel formula):
\[\begin{equation} \boxed{ \dx(t) = \Phi(t,t_0)\,\dx(t_0) + \int_{t_0}^{t} \Phi(t,s)\,\mat{B}(s)\,\du(s)\,\dd s. } \label{eq:ch2:duhamel} \end{equation}\]
Derivation
Consider the non-homogeneous equation \[\begin{equation} \dot{\dx} = \mat{A}(t)\dx + \mat{B}(t)\du(t). \end{equation}\]
We seek a particular solution using the method of variation of parameters. Assume \[\begin{equation} \dx(t) = \Phi(t,t_0)\,\mathbf{c}(t), \end{equation}\] where \(\mathbf{c}(t)\) is a time-varying coefficient vector to be determined. Then \[\begin{equation} \dot{\dx} = \dot{\Phi}(t,t_0)\,\mathbf{c} + \Phi(t,t_0)\,\dot{\mathbf{c}} = \mat{A}(t)\Phi(t,t_0)\,\mathbf{c} + \Phi(t,t_0)\,\dot{\mathbf{c}}. \end{equation}\]
Matching with the desired dynamics: \[\begin{equation} \mat{A}(t)\Phi(t,t_0)\,\mathbf{c} + \Phi(t,t_0)\,\dot{\mathbf{c}} = \mat{A}(t)\Phi(t,t_0)\,\mathbf{c} + \mat{B}(t)\du. \end{equation}\]
This simplifies to \[\begin{equation} \Phi(t,t_0)\,\dot{\mathbf{c}} = \mat{B}(t)\du, \quad \Rightarrow \quad \dot{\mathbf{c}} = \Phi(t_0,t)\mat{B}(t)\du, \end{equation}\] using \(\Phi(t_0,t) = \Phi(t,t_0)^{-1}\). Integrating from \(t_0\) to \(t\): \[\begin{equation} \mathbf{c}(t) = \mathbf{c}(t_0) + \int_{t_0}^{t}\Phi(t_0,s)\mat{B}(s)\du(s)\,\dd s. \end{equation}\]
With the initial condition \(\dx(t_0) = \Phi(t_0,t_0)\mathbf{c}(t_0) = \mathbf{c}(t_0)\), we have \[\begin{equation} \dx(t) = \Phi(t,t_0)\left[\dx(t_0) + \int_{t_0}^{t}\Phi(t_0,s)\mat{B}(s)\du(s)\,\dd s\right]. \end{equation}\]
Using the semigroup property \(\Phi(t,t_0)\Phi(t_0,s) = \Phi(t,s)\): \[\begin{equation} \dx(t) = \Phi(t,t_0)\dx(t_0) + \int_{t_0}^{t}\Phi(t,s)\mat{B}(s)\du(s)\,\dd s, \end{equation}\] which is the Duhamel formula.
Physical Interpretation
The Duhamel formula decomposes the state perturbation into two contributions:
- Homogeneous term \(\Phi(t,t_0)\dx(t_0)\): The initial perturbation \(\dx(t_0)\) is propagated forward via the STM. This is the response to initial conditions alone.
- Forced term \(\int_{t_0}^{t}\Phi(t,s)\mat{B}(s)\du(s)\,\dd s\): At each time \(s\), a control perturbation \(\du(s)\) is applied. It enters through the actuator channels \(\mat{B}(s)\), and its effect is then propagated forward from time \(s\) to time \(t\) by the STM \(\Phi(t,s)\). The integral sums the contributions of all past control inputs, weighted by their propagation matrices.
The Duhamel formula is the engine behind sensitivity analysis, adjoint methods, and trajectory optimization. It tells us how each input perturbation at time \(s\) contributes to the final state perturbation at time \(t\), weighted by the transition matrix \(\Phi(t,s)\). This is the precise sense in which input effects “add up”—they superpose linearly in the tangent space. For control design, this formula shows which times and which channels have the most leverage: early times (with large \(\Phi(t,s)\)) and actuator-rich directions (large column spaces of \(\mat{B}(s)\)) are most influential.
Discrete-Time Variational Dynamics
For computational implementation, we discretize the continuous-time dynamics. Given a sampling period \(h\) and zero-order-hold (ZOH) input, the discrete-time flow map is \[\begin{equation} \x_{k+1} = f_d(\x_k, \uvec_k), \label{eq:ch2:discrete-flow} \end{equation}\] where \(f_d\) is the time-\(h\) solution operator (flow map) of the continuous-time system~\(\eqref{eq:ch2:affine}\). Formally, \[\begin{equation} f_d(\x_k,\uvec_k) := \Phi_c\!\bigl(t_k+h, t_k; \uvec_k\bigr), \end{equation}\] where \(\Phi_c\) denotes the flow of the continuous-time system with zero-order-held input \(\uvec(t) = \uvec_k\) for \(t\in[t_k,t_k+h)\).
Linearization and Discrete Jacobians
Linearizing the discrete flow map around the nominal sequence \((\bar{\x}_k, \bar{\uvec}_k)\): \[\begin{equation} \dx_{k+1} = \mat{A}_k\,\dx_k + \mat{B}_k\,\du_k, \label{eq:ch2:discrete-variational} \end{equation}\] where \[\begin{equation} \mat{A}_k = \frac{\partial f_d}{\partial \x}\bigg|_{(\bar{\x}_k,\bar{\uvec}_k)}, \qquad \mat{B}_k = \frac{\partial f_d}{\partial \uvec}\bigg|_{(\bar{\x}_k,\bar{\uvec}_k)}. \label{eq:ch2:discrete-jacobians} \end{equation}\]
Computing Discrete Jacobians via Matrix Exponential
For a control-affine continuous system, the discrete flow under ZOH control satisfies \[\begin{equation} \x_{k+1} = \x_k + \int_{t_k}^{t_k+h}\Phi_c(t_k+h,\tau)\,G(\bar{\x}(\tau))\,\bar{\uvec}_k\,\dd\tau, \end{equation}\] where \(\Phi_c(t_k+h,\tau)\) is the state transition matrix of the linearized system along the nominal trajectory.
At the nominal trajectory point \(\bar{\x}_k\), assuming the input is held constant at \(\bar{\uvec}_k\) over the interval: \[\begin{equation} \bar{\x}_{k+1} = \bar{\x}_k + \int_{t_k}^{t_k+h}\Phi_c\!\bigl(t_k+h,\tau;\bar{\uvec}_k\bigr)G(\bar{\x}(\tau))\bar{\uvec}_k\,\dd\tau. \end{equation}\]
For the Jacobians: \[\begin{equation} \mat{A}_k = \Phi_c(t_k+h,t_k;\bar{\uvec}_k)\bigg|_{\bar{\x}_k}, \label{eq:ch2:Ak-stm} \end{equation}\] which is the STM over one sampling interval.
\[\begin{equation} \mat{B}_k = \int_{t_k}^{t_k+h}\Phi_c(t_k+h,\tau;\bar{\uvec}_k)\,G(\bar{\x}(\tau))\,\dd\tau, \label{eq:ch2:Bk-integral} \end{equation}\] which is an integral of the STM weighted by the input matrix.
For time-invariant systems with constant \(A\) and \(B\): \[\begin{equation} \mat{A}_k = \exp(A h), \quad \mat{B}_k = \left(\int_0^h \exp(A s)\,\dd s\right)B = A^{-1}\left[\exp(Ah)-\mat{I}\right]B. \label{eq:ch2:AB-discrete-formula} \end{equation}\]
Discrete State Transition Matrix
The discrete state transition matrix from stage \(k\) to stage \(\ell > k\) is the product of individual Jacobians: \[\begin{equation} \Phi_{d,\ell,k} = \mat{A}_{\ell-1}\,\mat{A}_{\ell-2}\cdots\mat{A}_k, \quad \ell > k, \qquad \Phi_{d,k,k} = \mat{I}_n. \label{eq:ch2:discrete-stm} \end{equation}\]
This satisfies a discrete semigroup property: \(\Phi_{d,\ell,k} = \Phi_{d,\ell,j}\Phi_{d,j,k}\) for \(k \leq j \leq \ell\).
Unlike the continuous STM, the discrete STM is not a fundamental solution to a matrix ODE; it is a pure product of matrices. For numerical work, computing \(\Phi_{d,\ell,k}\) by successive multiplication is standard, but for large \(\ell-k\), this accumulates rounding errors and can become numerically unstable (Section~\(\ref{sec:ch2:stability}\)).
A Worked Example: Linearized Pendulum
To concretize the theory, we compute the variational equation and STM explicitly for a simple pendulum linearized about the downward equilibrium.
Pendulum Dynamics
Consider a simple pendulum with mass \(m\), length \(\ell\), and a motor torque \(\tau\) applied at the pivot: \[\begin{equation} \ddot\theta = \frac{g}{\ell}\sin\theta + \frac{\tau}{m\ell^2}. \label{eq:ch2:pend-ode} \end{equation}\]
Writing this as a first-order system with \(x_1 = \theta\) and \(x_2 = \dot\theta\): \[\begin{equation} \begin{pmatrix}\dot x_1\\\dot x_2\end{pmatrix} = \begin{pmatrix}x_2\\\frac{g}{\ell}\sin x_1 + \frac{u}{m\ell^2}\end{pmatrix}, \label{eq:ch2:pend-first} \end{equation}\] where \(u = \tau\) is the control input.
Nominal Trajectory
Suppose the nominal trajectory is the equilibrium at the downward position with zero control: \[\begin{equation} \bar{\x}(t) = \begin{pmatrix}0\\0\end{pmatrix}, \quad \bar u(t) = 0. \label{eq:ch2:pend-nominal} \end{equation}\]
Verify: \(\dot{\bar x}_1 = 0\), \(\dot{\bar x}_2 = \frac{g}{\ell}\sin(0) + 0 = 0\). This is indeed a fixed point.
Jacobians
At the point \((\bar{\x}, \bar u) = (0,0)\): \[\begin{equation} \mat{A} = \frac{\partial}{\partial \x}\begin{pmatrix}x_2\\\frac{g}{\ell}\sin x_1\end{pmatrix}\bigg|_{(0,0)} = \begin{pmatrix}0 & 1\\\frac{g}{\ell}\cos(0) & 0\end{pmatrix} = \begin{pmatrix}0 & 1\\\frac{g}{\ell} & 0\end{pmatrix}. \label{eq:ch2:pend-A} \end{equation}\]
\[\begin{equation} \mat{B} = \frac{\partial}{\partial u}\begin{pmatrix}x_2\\\frac{g}{\ell}\sin x_1 + \frac{u}{m\ell^2}\end{pmatrix}\bigg|_{(0,0)} = \begin{pmatrix}0\\\frac{1}{m\ell^2}\end{pmatrix}. \label{eq:ch2:pend-B} \end{equation}\]
The variational equation is thus \[\begin{equation} \begin{pmatrix}\delta\dot x_1\\\delta\dot x_2\end{pmatrix} = \begin{pmatrix}0 & 1\\\frac{g}{\ell} & 0\end{pmatrix}\begin{pmatrix}\delta x_1\\\delta x_2\end{pmatrix} + \begin{pmatrix}0\\\frac{1}{m\ell^2}\end{pmatrix}\delta u. \label{eq:ch2:pend-variational} \end{equation}\]
Eigenvalues and State Transition Matrix
For the time-invariant case, we compute the eigenvalues of \(\mat{A}\): \[\begin{equation} \det(\lambda \mat{I} - \mat{A}) = \det\begin{pmatrix}\lambda & -1\\-\frac{g}{\ell} & \lambda\end{pmatrix} = \lambda^2 + \frac{g}{\ell} = 0. \end{equation}\]
Thus \(\lambda = \pm i\sqrt{g/\ell}\), which are purely imaginary. Let \(\omega_0 := \sqrt{g/\ell}\) (the natural frequency).
The matrix exponential can be computed using the formula for a 2D system: \[\begin{equation} \exp(\mat{A}t) = \cos(\omega_0 t)\mat{I} + \frac{\sin(\omega_0 t)}{\omega_0}\mat{A} = \begin{pmatrix}\cos(\omega_0 t) & \frac{\sin(\omega_0 t)}{\omega_0}\\ -\omega_0\sin(\omega_0 t) & \cos(\omega_0 t)\end{pmatrix}. \label{eq:ch2:pend-stm} \end{equation}\]
Verify: \(\frac{\dd}{\dd t}\Phi(t,0) = \begin{pmatrix}-\omega_0\sin(\omega_0 t) & \cos(\omega_0 t)\\\ -\omega_0^2\cos(\omega_0 t) & -\omega_0\sin(\omega_0 t)\end{pmatrix}\). Comparing with \(\mat{A}\Phi(t,0)\): \[\begin{equation} \begin{pmatrix}0 & 1\\\frac{g}{\ell} & 0\end{pmatrix}\begin{pmatrix}\cos(\omega_0 t) & \frac{\sin(\omega_0 t)}{\omega_0}\\ -\omega_0\sin(\omega_0 t) & \cos(\omega_0 t)\end{pmatrix} = \begin{pmatrix}-\omega_0\sin(\omega_0 t) & \cos(\omega_0 t)\\ \omega_0^2\cos(\omega_0 t) & \omega_0\sin(\omega_0 t)\end{pmatrix}. \end{equation}\]
Good (noting \(\omega_0^2 = g/\ell\)).
Numerical Example
For concreteness, let \(\ell = 1\) m and \(g = 10\) m/s\(^2\), so \(\omega_0 = \sqrt{10} \approx 3.162\) rad/s. The natural period is \(T = 2\pi/\omega_0 \approx 1.987\) s.
At time \(t = 0.5\) s: \[\begin{align} \omega_0 t &= \sqrt{10} \times 0.5 \approx 1.581 \text{ rad},\\ \cos(\omega_0 t) &\approx \cos(1.581) \approx -0.0108,\\ \sin(\omega_0 t) &\approx \sin(1.581) \approx 1.0000,\\ \frac{\sin(\omega_0 t)}{\omega_0} &\approx \frac{1.0000}{3.162} \approx 0.3162. \end{align}\]
Thus \[\begin{equation} \Phi(0.5, 0) \approx \begin{pmatrix}-0.0108 & 0.3162\\-3.162 & -0.0108\end{pmatrix}. \label{eq:ch2:pend-stm-numeric} \end{equation}\]
An initial perturbation \(\delta\x_0 = (0.1\text{ rad}, 0)^\T\) (small angular displacement) would evolve as \[\begin{equation} \delta\x(0.5) = \Phi(0.5,0)\delta\x_0 \approx \begin{pmatrix}-0.0108\\-0.3162\end{pmatrix}, \end{equation}\] meaning the angle perturbation has nearly oscillated back through zero, and a velocity perturbation of about \(-0.32\) rad/s has developed.
Duhamel Integral With Impulse Control
Suppose we apply an impulse of torque \(\delta u(t) = u_0\,\delta(t-t^*)\) at time \(t^* = 0.25\) s. The response at time \(t = 0.5\) s is \[\begin{equation} \delta\x(0.5) = \Phi(0.5, 0.25)\mat{B}\,u_0, \end{equation}\] where \[\begin{equation} \Phi(0.5,0.25) = \Phi(0.25, 0) = \begin{pmatrix}\cos(\sqrt{10} \times 0.25) & \frac{\sin(\sqrt{10} \times 0.25)}{\sqrt{10}}\\ -\sqrt{10}\sin(\sqrt{10} \times 0.25) & \cos(\sqrt{10} \times 0.25)\end{pmatrix}. \end{equation}\]
With \(\sqrt{10} \times 0.25 \approx 0.7906\) rad: \[\begin{equation} \Phi(0.5,0.25) \approx \begin{pmatrix}0.7038 & 0.2499\\-0.7905 & 0.7038\end{pmatrix}. \end{equation}\]
If \(u_0 = 1\) N\(\cdot\)m and \(m = 1\) kg, then \(\mat{B}u_0 = (0, 1)^\T\), and \[\begin{equation} \delta\x(0.5) \approx \begin{pmatrix}0.2499\\0.7038\end{pmatrix}, \end{equation}\] meaning a torque impulse at quarter-period induces a displacement of about 0.25 rad and velocity of 0.70 rad/s.
The numerical values are illustrative. In practice, one would solve this on a computer using ODE solvers (RK4, Runge-Kutta-Fehlberg) or leverage the analytic matrix exponential for this simple 2D system.
Sensitivity Analysis and Parameter Dependence
Parameter Sensitivity via the Variational Equation
Often we wish to understand how the state depends on a parameter \(p \in \mathbb{R}\), where \[\begin{equation} \dot{\x} = f(\x, p), \quad x\in\R^n, p \in \R^{n_p}. \label{eq:ch2:param-system} \end{equation}\]
Taking the partial derivative with respect to \(p\): \[\begin{equation} \frac{\partial}{\partial p}\dot{\x} = \frac{\partial}{\partial p}f(\x,p) = \frac{\partial f}{\partial \x}\frac{\partial \x}{\partial p} + \frac{\partial f}{\partial p}. \end{equation}\]
Let \(S(t) := \frac{\partial \x(t)}{\partial p}\) be the sensitivity (a matrix of shape \(n \times n_p\)). Then \[\begin{equation} \dot S = \frac{\partial f}{\partial \x}\bigg|_{\x(t),p} S + \frac{\partial f}{\partial p}\bigg|_{\x(t),p}, \quad S(t_0) = \frac{\partial \x(t_0)}{\partial p} = 0 \quad \text{(if } x(t_0) \text{ is independent of } p\text{)}. \label{eq:ch2:sensitivity-ode} \end{equation}\]
This is a linear, non-homogeneous ODE in \(S\), with the solution \[\begin{equation} S(t) = \int_{t_0}^{t}\Phi(t,s)\frac{\partial f}{\partial p}\bigg|_{\x(s),p}\,\dd s, \label{eq:ch2:sensitivity-duhamel} \end{equation}\] which is directly an application of the Duhamel formula with \(\mat{A}(s) = \frac{\partial f}{\partial\x}|_{\x(s),p}\), \(\mat{B}(s) = \frac{\partial f}{\partial p}|_{\x(s),p}\), and input \(\du = \dd p\).
Computing \(S(t)\) allows us to quickly evaluate the sensitivity of any observable (e.g., final position, maximum altitude, minimum distance) to parameter changes without re-simulating the full nonlinear system.
Adjoint Method for Gradient Computation
Now suppose we have a cost functional \[\begin{equation} J(p) = \ell(x(T), p) + \int_{t_0}^{T}\ell_t(\x(t),\uvec(t),p)\,\dd t, \label{eq:ch2:cost} \end{equation}\] and we wish to compute \(\frac{\partial J}{\partial p}\).
Rather than propagating sensitivities forward (as in~\(\eqref{eq:ch2:sensitivity-ode}\)), the adjoint method propagates a costate vector \(\lambda\) backward in time. Define \[\begin{equation} \lambda(t) := \frac{\partial J}{\partial \x}\bigg|_{\x(t)}. \label{eq:ch2:costate} \end{equation}\]
The adjoint equation is \[\begin{equation} \boxed{ \dot\lambda = -\frac{\partial f}{\partial \x}^T\lambda - \frac{\partial \ell_t}{\partial \x}, \quad \lambda(T) = -\frac{\partial \ell}{\partial \x}\bigg|_{x(T)}. } \label{eq:ch2:adjoint} \end{equation}\]
This equation is integrated backward from \(T\) to \(t_0\). Once \(\lambda(t)\) is known, the gradient is \[\begin{equation} \frac{\partial J}{\partial p} = \int_{t_0}^{T}\left[\lambda^T\frac{\partial f}{\partial p} + \frac{\partial \ell_t}{\partial p}\right]\dd t + \frac{\partial \ell}{\partial p}\bigg|_{x(T)}. \label{eq:ch2:grad-adjoint} \end{equation}\]
The adjoint method is far more efficient than forward sensitivity propagation when the cost depends on many parameters but only a few states. Forward sensitivity would require \(n_p\) parallel integrations (one per parameter), each of dimension \(n\). The adjoint requires only one backward integration of dimension \(n\), plus quadrature along the trajectory. This is the foundation of gradient-based optimization in optimal control, and it connects directly to the state transition matrix: the propagation of \(\lambda\) is governed by \(-\mat{A}(t)^T\), the transpose of the linearization.
Connection to Backpropagation
In machine learning, backpropagation through a neural network is exactly the adjoint method applied to the loss function. If the network is viewed as a dynamical system (e.g., via Neural ODEs), then the backward pass computes the adjoint equation \(\dot\lambda = -\mat{A}(t)^T\lambda\), and the gradient is accumulated via quadrature. This is not a coincidence: both variational dynamics and neural network training are governed by the same differential geometry and calculus of variations.
Modern automatic differentiation frameworks (PyTorch, JAX, TensorFlow) implement the adjoint method (or its discrete analogue) under the hood when computing gradients of dynamical systems with respect to parameters. This is a striking example of how geometric insights from control theory have been incorporated into machine learning infrastructure.
Liouville’s Formula and Phase Space Volume
Liouville’s formula (Theorem~\(\ref{thm:ch2:liouville}\)) has a beautiful geometric interpretation related to the evolution of phase-space volumes.
Interpretation via Divergence
Consider a small volume element \(\delta V\) in state space centered at \(\bar{\x}(t)\). Under the flow of the nonlinear system~\(\eqref{eq:ch2:affine}\), this volume element evolves as \[\begin{equation} \frac{\delta V(t)}{\delta V(t_0)} = \det(\Phi(t,t_0)) = \exp\left(\int_{t_0}^{t}\trace\mat{A}(s)\,\dd s\right) = \exp\left(\int_{t_0}^{t}\nabla \cdot f\bigg|_{\x(s)}\,\dd s\right). \label{eq:ch2:volume-evolution} \end{equation}\]
Here \(\nabla \cdot f = \sum_i \frac{\partial f_i}{\partial x_i} = \trace(\frac{\partial f}{\partial \x})\) is the divergence of the vector field \(f\).
Dissipative Systems
If \(\trace(\mat{A}(t)) < 0\) everywhere along the trajectory, then \[\begin{equation} \det(\Phi(t,t_0)) = \exp\left(\int_{t_0}^{t}\trace\mat{A}(s)\,\dd s\right) < 1, \end{equation}\] meaning phase-space volumes shrink exponentially. Such systems are called dissipative or contracting. Most nonlinear systems with friction or damping fall into this category.
Example: Pendulum With Friction
For the pendulum with damping coefficient \(c\): \[\begin{equation} \ddot\theta + c\dot\theta + \frac{g}{\ell}\sin\theta = \tau/(m\ell^2). \end{equation}\]
In first-order form, \(\mat{A} = \begin{pmatrix}0 & 1\\-(g/\ell) & -c\end{pmatrix}\). Thus \(\trace(\mat{A}) = -c < 0\), so volumes shrink at rate \(e^{-ct}\).
For numerical integration, tracking \(\det(\Phi)\) via Liouville’s formula provides a diagnostic: if \(\det(\Phi)\) grows unexpectedly, it may indicate numerical instability or an error in the Jacobian computation.
Numerical Stability and Ill-Conditioning
Computing the STM accurately is challenging, especially for stiff systems (those with widely separated time scales). We discuss the challenges and remedies.
Ill-Conditioning in Stiff Systems
A system is stiff if the Jacobian \(\mat{A}(t)\) has eigenvalues with vastly different magnitudes or real parts. For example, if \(\lambda_1 = -1\) and \(\lambda_2 = -1000\), the stiffness ratio is \(1000:1\).
For such systems, naive numerical integration with explicit methods (RK4) becomes extremely slow: the step size must be tiny to resolve the fastest mode (with \(\lambda_2\)), even though the solution is dominated by the slower mode (with \(\lambda_1\)). Implicit methods (BDF, Rosenbrock) handle stiffness better but can still suffer from conditioning issues when computing the STM.
Consider the STM: if a column \(\Phi_{\cdot,j}(t,t_0)\) corresponds to an initial condition along an eigenvector of a stable eigenvalue, that column will decay as \(e^{\text{Re}(\lambda)t}\). For large times and stable systems, \(\Phi\) becomes numerically singular (entries underflow to zero), making it difficult to compute \(\Phi^{-1}\) accurately.
QR-Based Propagation
A robust remedy is to propagate the STM using a QR decomposition: \[\begin{equation} \Phi(t,t_0) = Q(t,t_0)\,R(t,t_0), \label{eq:ch2:qr-stm} \end{equation}\] where \(Q\) is orthonormal and \(R\) is upper triangular. The QR factorization is updated at each step without explicitly computing \(\Phi\).
- Propagate: Integrate the ODE \(\dot Q = \mat{A}(t)Q\) with \(Q(t_0,t_0) = \mat{I}\).
- Orthonormalize: At each step, compute \([Q_{\text{new}}, R_{\text{step}}] = \text{qr}(\tilde{Q})\), where \(\tilde{Q}\) is the current approximate matrix from the integrator. Store \(R = R \cdot R_{\text{step}}\).
- Recover: At any time, \(\Phi = Q \cdot R\) can be reconstructed, but more importantly, \(R\) remains well-conditioned (bounded) throughout the integration.
This approach prevents underflow and overflow and maintains numerical stability for arbitrarily long integration times.
Singular Value Monitoring
Another diagnostic is to monitor the singular values of \(\Phi(t,t_0)\) as the integration progresses: \[\begin{equation} \Phi = U\,\Sigma\,V^T, \quad \Sigma = \diag(\sigma_1, \ldots, \sigma_n), \quad \sigma_1 \geq \cdots \geq \sigma_n \geq 0. \label{eq:ch2:svd-stm} \end{equation}\]
For a stable system, the smallest singular value \(\sigma_n\) decays exponentially, reflecting the contraction in stable directions. If \(\sigma_n\) becomes smaller than machine epsilon (\(\sim 10^{-16}\) in double precision), further numerical integration loses information, and the computed STM becomes unreliable. At that point, it is best to restart the integration from a recent checkpoint or use higher precision arithmetic.
Best practices for stable STM computation: 1. Use automatic differentiation to compute \(\mat{A}(t)\) accurately. 2. For stiff systems, use implicit integrators (BDF or Rosenbrock) rather than explicit ones. 3. For very long integration times, use QR-based propagation to maintain stability. 4. Monitor singular values of \(\Phi\) and switch to higher precision if needed. 5. When \(\det(\Phi)\) becomes very small or very large, verify against Liouville’s formula.
Qualitative Interpretation of the Jacobians
\(\mat{A}(t)\): Local State Sensitivity
The eigenvalues and singular values of \(\mat{A}(t)\) reveal the local perturbation landscape:
- Eigenvalues with positive real part: indicate locally unstable directions—perturbations that amplify exponentially. If \(\lambda = \alpha + i\beta\) with \(\alpha > 0\), perturbations grow as \(e^{\alpha t}\) while oscillating with frequency \(\beta\).
- Eigenvalues with negative real part: indicate locally contracting directions. Perturbations decay as \(e^{\alpha t}\) with \(\alpha < 0\).
- The symmetric part \(\sym(\mat{A}) = \frac{1}{2}(\mat{A}+\mat{A}^T)\) governs instantaneous growth/decay of perturbation energy. The quadratic form \(\|\dx\|^2\) evolves as \[\begin{equation} \frac{\dd}{\dd t}\|\dx\|^2 = 2\dx^T\dot{\dx} = 2\dx^T\mat{A}\dx = 2\dx^T[\sym(\mat{A})]\dx + 2\dx^T[\skewpart(\mat{A})]\dx. \end{equation}\] The skew part \(\skewpart(\mat{A}) = \frac{1}{2}(\mat{A}-\mat{A}^T)\) does not affect energy directly; it causes rotation.
For control design, the eigenvalues of \(\mat{A}(t)\) tell us which directions are naturally unstable (requiring control action to stabilize) and which are naturally stable (requiring control only for disturbance rejection).
Residuals Under Finite Perturbations
When we use the variational equation to predict the evolution of a finite (not infinitesimal) perturbation \(\dx_0\), the true trajectory \(\x(t) = \bar{\x}(t)+\dx(t)\) differs from the linear prediction by a residual.
Residual Integral Form
The true state satisfies \[\begin{equation} \dot{\x} = f(\x) + G(\x)\uvec = f(\bar{\x}+\dx) + G(\bar{\x}+\dx)(\bar{\uvec}+\du). \end{equation}\]
Expanding as before but keeping second-order terms: \[\begin{equation} \dot{\x} = f(\bar{\x}) + J_f(\bar{\x})\dx + \frac{1}{2}\frac{\partial^2 f}{\partial\x^2}(\dx,\dx) + G(\bar{\x})\du + O(\|\dx\|^2) + \cdots \end{equation}\]
The true state perturbation \(\dx_{\text{true}}\) satisfies \[\begin{equation} \dot{\dx}_{\text{true}} = J_f(\bar{\x})\dx_{\text{true}} + \frac{1}{2}\frac{\partial^2 f}{\partial\x^2}(\dx_{\text{true}}, \dx_{\text{true}}) + O(\|\dx_{\text{true}}\|^2) + G(\bar{\x})\du. \end{equation}\]
Let \(\dx_{\text{lin}}\) be the solution from the linearized equation: \[\begin{equation} \dot{\dx}_{\text{lin}} = J_f(\bar{\x})\dx_{\text{lin}} + G(\bar{\x})\du. \end{equation}\]
The residual is \(r(\dx) := \dot{\dx}_{\text{true}} - \dot{\dx}_{\text{lin}} = \frac{1}{2}\frac{\partial^2 f}{\partial\x^2}(\dx_{\text{true}}, \dx_{\text{true}}) + O(\|\dx_{\text{true}}\|^3)\).
Integrating the error equation: \[\begin{equation} \dx_{\text{true}}(t) - \dx_{\text{lin}}(t) = \int_{t_0}^{t}\Phi(t,s)\,r(\dx_{\text{true}}(s))\,\dd s = \frac{1}{2}\int_{t_0}^{t}\Phi(t,s)\frac{\partial^2 f}{\partial\x^2}|_{\x(s)}(\dx_{\text{true}}(s), \dx_{\text{true}}(s))\,\dd s + O(\|\dx_{\text{true}}\|^3). \label{eq:ch2:residual-integral} \end{equation}\]
Residual Bound
Under \(C^2\) smoothness, the Hessian is bounded: \[\begin{equation} \left\|\frac{\partial^2 f}{\partial\x^2}\right\| \leq C_H \end{equation}\] in a neighborhood of the nominal trajectory. Thus \[\begin{equation} \|r(\dx)\| \leq \tfrac{1}{2}\,C_H\,\|\dx\|^2. \label{eq:ch2:residual-hessian} \end{equation}\]
For small \(\|\dx\|\), the residual is much smaller than the main term, justifying linearization. As \(\|\dx\|\) grows, the residual becomes significant. Assuming \(\|\Phi(t,s)\| \lesssim 1\) (e.g., stable system), we have \[\begin{equation} \|\Delta x(t)\| \lesssim \int_{t_0}^{t}C_H\,\|\dx_{\text{true}}(s)\|^2\,\dd s. \label{eq:ch2:residual-growth} \end{equation}\]
If we denote \(M := \max_{s\in[t_0,t]}\|\dx_{\text{true}}(s)\|\), then \[\begin{equation} \|\Delta x(t)\| \lesssim C_H\,M^2\,(t-t_0). \label{eq:ch2:residual-final} \end{equation}\]
This is not a defect of linearization but a geometric feature: it measures the curvature of the vector field along the trajectory.
Residual Scaling and Iterative Refinement
The \(O(\|\dx\|^2)\) scaling is why iterative linearization-based methods (DDP, SQP) converge quadratically near a solution: each re-linearization eliminates the dominant residual term, leaving only higher-order curvature. After one iteration, \(M\) shrinks quadratically, so after \(k\) iterations, the residual is smaller by a factor of \(2^{-2^k}\), leading to superlinear (quadratic) convergence.
This quadratic convergence is one of the key advantages of locally optimal methods like DDP and SQP over open-loop trajectory optimization methods that do not re-linearize: they can refine a solution very quickly once they are in a neighborhood of the optimum.
Structure-Preserving Numerical Integration
The numerical methods discussed in Section~\(\ref{sec:ch2:numerical-stm}\) (Runge–Kutta methods, Padé approximants) are general-purpose. For specific classes of systems—particularly Hamiltonian mechanical systems and systems evolving on Lie groups—specialized integrators that preserve geometric structure can dramatically improve accuracy and long-term behavior.
Why Structure Preservation Matters
Consider a Hamiltonian system with \(H(\q, \vec{p}) = T(\vec{p}) + V(\q)\). The continuous dynamics preserve the Hamiltonian (\(\dot{H} = 0\)), the symplectic form, and phase-space volume (Liouville’s theorem, Section~\(\ref{sec:ch2:liouville}\)). A generic Runge–Kutta integrator preserves none of these. Over long simulations, this manifests as: - Energy drift: the numerically computed energy \(H(\q_k, \vec{p}_k)\) drifts monotonically upward or downward, even though the true energy is constant. - Phase-space expansion: trajectories spiral outward in phase space, destroying the qualitative behavior of the system. - Broken symmetries: conservation laws (angular momentum, for instance) are violated, leading to unphysical trajectories.
For short simulations (a single golf downswing at 0.5 s), these effects are typically negligible. For repeated simulations, parameter sweeps, or optimization loops that iterate many times, accumulation of drift can corrupt the results.
Symplectic Integrators
A symplectic integrator is a numerical scheme that preserves the symplectic 2-form \(\dd\q \wedge \dd\vec{p}\). The simplest example is the symplectic Euler method: \[\begin{align} \vec{p}_{k+1} &= \vec{p}_k - h \frac{\partial V}{\partial \q}(\q_k), \label{eq:ch2:symplectic-euler-p} \\ \q_{k+1} &= \q_k + h \frac{\partial T}{\partial \vec{p}}(\vec{p}_{k+1}). \label{eq:ch2:symplectic-euler-q} \end{align}\]
The key property: the map \((\q_k, \vec{p}_k) \mapsto (\q_{k+1}, \vec{p}_{k+1})\) has Jacobian with determinant exactly 1 (for all step sizes \(h\)), preserving phase-space volume exactly. Symplectic integrators do not conserve energy exactly, but the energy error remains bounded for all time (it oscillates rather than drifting).
Higher-order symplectic methods (Störmer–Verlet, Yoshida compositions) achieve \(O(h^p)\) accuracy while maintaining the symplectic property.
Lie Group Integrators
When the state space includes rotation groups (SO(3) for rigid-body attitude, SE(3) for pose), standard integrators that parameterize rotations via Euler angles or quaternions face singularities or constraint drift. Lie group integrators work directly on the group:
- Compute the velocity in the Lie algebra \(\mathfrak{g}\) (tangent space at the identity).
- Integrate in \(\mathfrak{g}\) using standard methods.
- Map back to the group via the exponential map: \(R_{k+1} = R_k \cdot \exp(\omega_k \, h)\).
This guarantees that \(R_{k+1} \in \text{SO}(3)\) exactly (no normalization or re-projection needed). The Crouch-Grossman and Munthe-Kaas methods provide systematic frameworks for constructing Lie group integrators of arbitrary order.
Variational Integrators
Variational integrators derive the discrete equations of motion from a discrete variational principle (discrete Hamilton’s principle) rather than discretizing the continuous equations. The resulting integrator automatically preserves the symplectic structure and momentum maps.
For the discrete Lagrangian \(L_d(\q_k, \q_{k+1}) \approx \int_{t_k}^{t_{k+1}} L(\q, \dot{\q}) \, \dd t\), the discrete Euler-Lagrange equations yield: \[\begin{equation} D_2 L_d(\q_{k-1}, \q_k) + D_1 L_d(\q_k, \q_{k+1}) = 0, \label{eq:ch2:discrete-euler-lagrange} \end{equation}\] which implicitly defines \(\q_{k+1}\) given \(\q_{k-1}\) and \(\q_k\).
For trajectory optimization and DDP (Chapter~\(\ref{ch:optimal-control}\)), variational integrators have a specific advantage: the STM \(\Phi\) computed from the variational integrator inherits the symplectic structure. This means the sensitivities used in the Riccati recursion are geometrically consistent, improving convergence of the optimization. For systems on Lie groups (robotics, spacecraft), the combination of Lie group integration with variational principles provides a fully structure-preserving computational framework.
Chapter Summary
- Variational Equation from First Principles: The variational equation \(\dot{\dx} = \mat{A}(t)\dx + \mat{B}(t)\du\) emerges from a careful Taylor expansion of the nonlinear dynamics, keeping first-order terms in \(\dx\) and \(\du\) and dropping \(O(\|\dx\|^2)\) residuals. This is exact for infinitesimal perturbations and an excellent first-order approximation for small but finite perturbations.
- State Transition Matrix: The STM \(\Phi(t,t_0)\) is the unique solution to \(\dot\Phi = \mat{A}(t)\Phi\) with \(\Phi(t_0,t_0) = \mat{I}\). It characterizes the evolution of any perturbation via \(\dx(t) = \Phi(t,t_0)\dx(t_0)\).
- Fundamental Properties: The STM satisfies three key properties: semigroup (\(\Phi(t_2,t_0) = \Phi(t_2,t_1)\Phi(t_1,t_0)\)), invertibility (\(\Phi^{-1}(t,t_0) = \Phi(t_0,t)\)), and Liouville’s formula for the determinant.
- Peano-Baker Series: The STM admits an explicit infinite series solution (Peano-Baker series), which converges uniformly for bounded \(\mat{A}(t)\). This series elegantly handles non-commutativity in the matrix product.
- Numerical Integration: Computing \(\Phi\) typically requires integrating the \(n \times n\) matrix ODE. For time-invariant systems, the matrix exponential can be computed via Padé approximation with scaling and squaring. For stiff systems, QR-based propagation ensures numerical stability.
- Variation of Constants (Duhamel) Formula: When both initial and input perturbations are present, the solution is \(\dx(t) = \Phi(t,t_0)\dx(t_0) + \int_{t_0}^{t}\Phi(t,s)\mat{B}(s)\du(s)\dd s\). This shows how effects superpose linearly in the tangent space and is fundamental to sensitivity analysis and optimal control.
- Discrete-Time Variational Dynamics: For computational implementation, the continuous-time variational equation is discretized. The discrete STM is a product of Jacobians: \(\Phi_{d,\ell,k} = \prod_{j=k}^{\ell-1}\mat{A}_j\).
- Sensitivity and Adjoint Methods: Parameter sensitivities \(\frac{\partial\x}{\partial p}\) satisfy a linearized ODE whose solution is given by the Duhamel formula. The adjoint method computes cost gradients efficiently by propagating a costate backward in time, governed by \(\dot\lambda = -\mat{A}(t)^T\lambda\).
- Residuals and Convergence: Finite perturbations incur residuals of order \(O(\|\dx\|^2)\), reflecting manifold curvature. Iterative re-linearization (as in DDP) eliminates the dominant residual, leading to quadratic convergence.
- Numerical Challenges: Stiff systems and long integration times can render \(\Phi\) numerically singular. QR-based propagation and singular value monitoring provide robust diagnostics and solutions.
- Matrix Interpretation: \(\mat{A}(t)\) encodes local sensitivity (eigenvalues reveal stable/unstable directions); \(\mat{B}(t)\) encodes local control authority (column space reveals reachable directions). Together they are the raw material for all subsequent stability and optimality analyses.
The variational equation is not merely a linear approximation: it is an exact first-order differential-geometric object describing the behavior of tangent vectors along a trajectory in a curved (nonlinear) vector field. The state transition matrix is the parallel transport operator in this geometry, evolving vectors along the trajectory while respecting the local geometry. This perspective connects classical control theory to modern differential geometry and is essential for understanding advanced topics like Riemannian optimization, geometric tracking control, and neural network training dynamics.
Exercises
- STM of a linear system. For the LTI system \(\mat{A} = \begin{bmatrix} 0 & 1 \\ -\omega^2 & -2\zeta\omega \end{bmatrix}\) (damped harmonic oscillator), compute \(\Phi(t, 0) = e^{\mat{A}t}\) analytically using eigenvalue decomposition. Verify the semigroup property \(\Phi(t_2, 0) = \Phi(t_2, t_1)\Phi(t_1, 0)\).
- Peano-Baker truncation error. For the system in Exercise 1, compute the first three terms of the Peano-Baker series~\(\eqref{eq:ch2:peano-baker}\). Estimate the truncation error at \(t = 1\) and compare with the exact exponential.
- Liouville’s formula. Verify Liouville’s formula (Theorem~\(\ref{thm:ch2:liouville}\)) for the 2D system \(\dot{x}_1 = x_2\), \(\dot{x}_2 = -x_1 - 0.1 x_2\) by computing both \(\det\Phi(t, 0)\) directly and \(\exp(\int_0^t \trace\mat{A}(s)\,\dd s)\).
- Sensitivity of a nonlinear system. For the Van der Pol oscillator \(\ddot{x} - \mu(1 - x^2)\dot{x} + x = 0\) with \(\mu = 1\), compute the variational equation along the trajectory starting from \((x_0, \dot{x}_0) = (2, 0)\). Integrate the STM numerically for one period. What do the singular values of \(\Phi\) reveal?
- Input sensitivity. For the control-affine system \(\dot{\x} = \mat{A}\x + \mat{B}\uvec\) with \(\mat{A}\), \(\mat{B}\) from the pendulum example in Section~\(\ref{sec:ch2:pendulum-example}\), use the variation of constants formula to compute the response to a unit pulse input \(\uvec(t) = \delta(t - t_0)\).
- Symplectic integrator comparison. Implement both RK4 and the symplectic Euler method for the undamped pendulum \(\ddot\theta + \sin\theta = 0\). Integrate for 1000 periods. Plot the energy \(H(\theta, \dot\theta)\) as a function of time for both methods. Explain the qualitative difference.
- Discrete-time STM. For the discrete system \(\x_{k+1} = \mat{A}_k \x_k\) with \(\mat{A}_k\) periodic (\(\mat{A}_{k+N} = \mat{A}_k\)), express the one-period STM as a product. What is the relationship between the eigenvalues of this product and the stability of the system?
- Magnus expansion. For a time-varying \(\mat{A}(t) = \mat{A}_0 + \varepsilon \mat{A}_1 \sin(\omega t)\), compute the first-order Magnus expansion and compare with the Peano-Baker series at \(t = 2\pi/\omega\). Under what conditions does the Magnus expansion converge faster?