Prev
Unit 9 · Ordinary Differential Equations: Initial Value Problems
Linear Systems and the Matrix Exponential
65m100 XP
Next
🛰️Ordinary Differential Equations: Initial Value Problems

Linear Systems and the Matrix Exponential

x′ = Ax, modes and stability, computing e^{At}, exact discretisation for digital control, transient growth, the Lyapunov equation

Nineteen dubious ways, 1978
By the 1970s, control engineers had a routine problem. A spacecraft, an aircraft or a chemical plant is modelled by linear differential equations x′=Ax+Bu\mathbf x' = A\mathbf x + B\mathbf u, but the controller is a digital computer that reads sensors and sets actuators at fixed instants. To design it, the continuous model must become a difference equation, and the matrix that does that is the exponential eATe^{AT}. Every textbook gave a formula for it. Programs implementing those formulas gave different answers for the same matrix — some of them wrong in every digit. In 1978 two numerical analysts, Cleve Moler and Charles Van Loan, published a survey in SIAM Review with a deliberately unsettling title: Nineteen Dubious Ways to Compute the Exponential of a Matrix. They examined series methods, ODE solvers, polynomial and eigenvector methods and matrix decompositions, and showed how each could fail. None was satisfactory for every matrix. A few were good enough to recommend. The paper became one of the most influential in numerical linear algebra. Twenty-five years later the journal republished it with an update, Twenty-Five Years Later. Moler, meanwhile, had written a small interactive matrix program for his students — MATLAB — and its expm function implemented the method they recommended.
A linear system with constant coefficients, x′=Ax\mathbf x' = A\mathbf x, x(0)=x0\mathbf x(0) = \mathbf x_0, has the exact solution
x(t)=eAtx0,eAt=I+At+(At)22!+(At)33!+⋯\mathbf x(t) = e^{At}\mathbf x_0, \qquad e^{At} = I + At + \frac{(At)^2}{2!} + \frac{(At)^3}{3!} + \cdots
the same series as the scalar exponential, which converges for every square matrix. Every linearised model in engineering has this form: a structure's vibration (Unit 6), small deviations of a spacecraft's attitude or an aircraft from trim, an electrical network, the error dynamics of a navigation filter. For such systems the time-stepping of lessons 9-1 to 9-4 is not needed at all: the exact solution is available.
Derivation
Why eAtx0e^{At}\mathbf x_0 solves the system, and what it looks like in modes
It solves the equation. Differentiate the series term by term; each term loses one power of tt and gains a factor AA:
ddt∑k≥0Aktkk!=∑k≥1Aktk−1(k−1)!=A∑j≥0Ajtjj!=A eAt.\frac{d}{dt}\sum_{k\ge0}\frac{A^kt^k}{k!} = \sum_{k\ge1}\frac{A^kt^{k-1}}{(k-1)!} = A\sum_{j\ge0}\frac{A^jt^j}{j!} = A\,e^{At}.
So x(t)=eAtx0\mathbf x(t) = e^{At}\mathbf x_0 satisfies x′=Ax\mathbf x' = A\mathbf x, and at t=0t = 0 the series is II, so x(0)=x0\mathbf x(0) = \mathbf x_0.
In modes. If A=VΛV−1A = V\Lambda V^{-1}, then A2=VΛV−1VΛV−1=VΛ2V−1A^2 = V\Lambda V^{-1}V\Lambda V^{-1} = V\Lambda^2V^{-1}, and in general Ak=VΛkV−1A^k = V\Lambda^kV^{-1}. Put that in every term of the series:
eAt=V(∑kΛktkk!)V−1=V eΛt V−1,eΛt=diag⁡(eλ1t,…,eλnt).e^{At} = V\Big(\sum_k\frac{\Lambda^kt^k}{k!}\Big)V^{-1} = V\,e^{\Lambda t}\,V^{-1}, \qquad e^{\Lambda t} = \operatorname{diag}(e^{\lambda_1t}, \dots, e^{\lambda_nt}).
With c=V−1x0\mathbf c = V^{-1}\mathbf x_0 this is x(t)=∑icieλitvi\mathbf x(t) = \sum_ic_ie^{\lambda_it}\mathbf v_i, the modal sum below — the same modes as lesson 6-1's guess veλt\mathbf v e^{\lambda t}, now added up to fit the initial state.
💡Modes and stability
If A=VΛV−1A = V\Lambda V^{-1} with eigenvalues λi\lambda_i and eigenvectors vi\mathbf v_i (Unit 6), then eAt=VeΛtV−1e^{At} = Ve^{\Lambda t}V^{-1}, and x(t)=∑ici eλit vi.\mathbf x(t) = \sum_i c_i\, e^{\lambda_i t}\,\mathbf v_i. The solution is a sum of modes. A real eigenvalue gives exponential growth or decay with time constant 1/∣λ∣1/|\lambda|; a complex pair σ±iω\sigma \pm i\omega gives an oscillation at frequency ω\omega with envelope eσte^{\sigma t}. The system is asymptotically stable exactly when every eigenvalue has Re⁡λ<0\operatorname{Re}\lambda < 0. For a damped oscillator x′′+2ζωnx′+ωn2x=0x'' + 2\zeta\omega_n x' + \omega_n^2x = 0 the eigenvalues are −ζωn±iωn1−ζ2-\zeta\omega_n \pm i\omega_n\sqrt{1 - \zeta^2}: damping sets the decay, the natural frequency the oscillation. Eigenvalues spread over many orders of magnitude are what made the equations of lesson 9-4 stiff.
📐Moler and Van Loan's test matrix
A=(−4924−6431)A = \begin{pmatrix} -49 & 24 \\ -64 & 31 \end{pmatrix} has eigenvalues −1-1 and −17-17, so eAe^A is small: its entries are below 1.5. Summing the Taylor series, the terms Ak/k!A^k/k! grow to 9.3×1069.3 \times 10^6 before they decay, and the answer is what remains after they cancel — catastrophic cancellation (lesson 1-2) on a grand scale. In single precision the result is wrong by 87%; in double precision it loses seven of its sixteen digits. The eigenvector formula is perfect here — but for a nearly defective matrix with two eigenvalues 10−910^{-9} apart, the eigenvector matrix has condition number 2×1092 \times 10^9 and the formula loses eight digits. scipy.linalg.expm is accurate to 10−1510^{-15} on both.
💡How expm works: scaling and squaring
The method Moler and Van Loan recommended, refined since by Nicholas Higham and Awad Al-Mohy, uses the identity eA=(eA/2s)2se^{A} = \left(e^{A/2^s}\right)^{2^s}:
  1. Scale: choose ss so that ∥A/2s∥\|A/2^s\| is below about 1 — the series of that small matrix has no large terms to cancel.
  2. Approximate eA/2se^{A/2^s} with a rational (Padé) approximation of degree up to 13, accurate to double precision for such a small argument.
  3. Square the result ss times.
Cost: a few dozen matrix multiplications and one linear solve, O(n3)O(n^3), with an error at the level of the condition of the problem itself. Use scipy.linalg.expm; never sum the series, and use the eigendecomposition only for symmetric (normal) matrices, where VV is orthogonal and perfectly conditioned.
Exact discretisation. A digital controller applies an input uk\mathbf u_k and holds it constant until the next sample — a zero-order hold. Over one sample period TT the solution of x′=Ax+Bu\mathbf x' = A\mathbf x + B\mathbf u is then exactly
xk+1=Φ xk+Γ uk,Φ=eAT,Γ=∫0TeAs ds  B.\mathbf x_{k+1} = \Phi\,\mathbf x_k + \Gamma\,\mathbf u_k, \qquad \Phi = e^{AT}, \quad \Gamma = \int_0^T e^{As}\,ds\;B.
Both come from one exponential of an augmented matrix: exp⁡ ⁣((AB00)T)=(ΦΓ0I)\exp\!\left(\begin{pmatrix} A & B \\ 0 & 0 \end{pmatrix}T\right) = \begin{pmatrix} \Phi & \Gamma \\ 0 & I \end{pmatrix} — a trick Van Loan published in the same year as the survey. There is no truncation error at all: the discrete model reproduces the continuous one exactly at the sampling instants, for any TT. Its eigenvalues are eλTe^{\lambda T}, inside the unit circle exactly when Re⁡λ<0\operatorname{Re}\lambda < 0 — continuous stability maps to discrete stability.
Derivation
The zero-order-hold formulas, and the augmented-matrix trick
Variation of constants. Multiply x′−Ax=Bu\mathbf x' - A\mathbf x = B\mathbf u by e−Ate^{-At}. The left side becomes an exact derivative, because ddt(e−Atx)=e−Atx′−Ae−Atx\frac{d}{dt}(e^{-At}\mathbf x) = e^{-At}\mathbf x' - Ae^{-At}\mathbf x:
ddt(e−Atx)=e−AtBu.\frac{d}{dt}\big(e^{-At}\mathbf x\big) = e^{-At}B\mathbf u.
Integrate from 00 to TT with u\mathbf u held constant, then multiply by eATe^{AT}:
x(T)=eATx(0)+∫0TeA(T−s) ds  Bu=Φ x(0)+(∫0TeAσ dσ)B u,\mathbf x(T) = e^{AT}\mathbf x(0) + \int_0^Te^{A(T - s)}\,ds\;B\mathbf u = \Phi\,\mathbf x(0) + \Big(\int_0^Te^{A\sigma}\,d\sigma\Big)B\,\mathbf u,
substituting σ=T−s\sigma = T - s in the integral. That is Φ\Phi and Γ\Gamma, with no approximation anywhere.
One exponential for both. Let M=[AB00]M = \begin{bmatrix} A & B \\ 0 & 0 \end{bmatrix}. Its powers keep the zero row: Mk=[AkAk−1B00]M^k = \begin{bmatrix} A^k & A^{k-1}B \\ 0 & 0 \end{bmatrix} for k≥1k \ge 1. Summing the exponential series block by block,
eMT=[∑kAkTkk!∑k≥1Ak−1Tkk! B0I],e^{MT} = \begin{bmatrix} \sum_k \dfrac{A^kT^k}{k!} & \sum_{k\ge1}\dfrac{A^{k-1}T^k}{k!}\,B \\ 0 & I \end{bmatrix},
and the top-right series is exactly ∫0TeAσdσ B\int_0^Te^{A\sigma}d\sigma\,B (integrate ∑jAjσj/j!\sum_j A^j\sigma^j/j! term by term). One call to expm gives Φ\Phi and Γ\Gamma together.
Euler's eigenvalues. Forward Euler uses I+ATI + AT, whose eigenvalues are 1+λT1 + \lambda T. For λ=σ+iω\lambda = \sigma + i\omega, ∣1+λT∣2=(1+σT)2+ω2T2|1 + \lambda T|^2 = (1 + \sigma T)^2 + \omega^2T^2, which exceeds 1 whenever the damping is light compared with the frequency — however small TT is.
📐Forward Euler is not a discretisation
A lightly damped structural mode (1 Hz, 2% damping) sampled at T=10T = 10 ms and pushed by 10 N for half a second. The exact map Φ\Phi has eigenvalues of magnitude 0.998744 =e−ζωnT= e^{-\zeta\omega_nT}; after 20 s its prediction agrees with a high-accuracy ODE solution to all six printed digits (−0.042333-0.042333 m). The 'obvious' model xk+1=(I+AT)xk+BTuk\mathbf x_{k+1} = (I + AT)\mathbf x_k + BT\mathbf u_k — forward Euler — has eigenvalues of magnitude 1.000717: unstable, for any step size, as lesson 9-4's stability region predicts for eigenvalues near the imaginary axis. After 20 s it predicts −2.08-2.08 m — fifty times the truth, and growing. A controller or Kalman filter designed on that model would be designed for a system that does not exist.
💡When eigenvalues mislead: transient growth
Stable eigenvalues guarantee that x(t)→0\mathbf x(t) \to 0 eventually, not that it stays small on the way. A=(−1500−2)A = \begin{pmatrix} -1 & 50 \\ 0 & -2 \end{pmatrix} has eigenvalues −1-1 and −2-2, yet ∥eAt∥\|e^{At}\| rises to 12.5 at t=0.69t = 0.69 before decaying: its eigenvectors are nearly parallel, and two decaying modes that nearly cancel at t=0t = 0 stop cancelling as one decays faster than the other. Such non-normal matrices — whose eigenvectors are far from orthogonal — appear in shear flows (the transition to turbulence), in feedback loops and in coupled structures. Lloyd N. Trefethen's pseudospectra make this precise. The engineering lesson: check ∥eAt∥\|e^{At}\|, not just the eigenvalues.
The Lyapunov equation. If a stable system is driven by white noise of intensity WW through GG, x′=Ax+Gw\mathbf x' = A\mathbf x + G\mathbf w, its covariance settles to the solution PP of
AP+PAT+GWGT=0,AP + PA^{\mathsf T} + GWG^{\mathsf T} = 0,
a linear equation for the entries of PP, solved by scipy.linalg.solve_continuous_lyapunov (the Bartels–Stewart algorithm of 1972, which first reduces AA to triangular form with the QR algorithm of lesson 6-3). For an oscillator shaken by noise it gives var⁡(x)=W/(4ζωn3)\operatorname{var}(x) = W/(4\zeta\omega_n^3) — 0.15625 for W=0.5W = 0.5, ζ=0.1\zeta = 0.1, ωn=2\omega_n = 2 — and 4000 simulated oscillators agree (0.1568). Aleksandr Lyapunov introduced the equation in 1892 to prove stability; in Unit 12 it reappears as the steady-state uncertainty of a filter, and its nonlinear relative, the Riccati equation, as the steady state of the Kalman filter.
Derivation
Where the Lyapunov equation comes from — and how the lab simulates the noise
One sample step. Over a short interval Δt\Delta t the state moves as xk+1=Φxk+wk\mathbf x_{k+1} = \Phi\mathbf x_k + \mathbf w_k, where wk\mathbf w_k is the noise collected during the interval, independent of xk\mathbf x_k, with covariance QdQ_d. Lesson 5-1's rule for a linear map, plus independence, gives the covariance one step later:
Pk+1=ΦPkΦT+Qd.P_{k+1} = \Phi P_k\Phi^{\mathsf T} + Q_d.
The continuous limit. For small Δt\Delta t, Φ≈I+AΔt\Phi \approx I + A\Delta t and Qd≈GWGTΔtQ_d \approx GWG^{\mathsf T}\Delta t (white noise adds variance in proportion to elapsed time). Then
Pk+1−Pk≈(APk+PkAT+GWGT)Δt,P_{k+1} - P_k \approx \big(AP_k + P_kA^{\mathsf T} + GWG^{\mathsf T}\big)\Delta t,
dropping terms in Δt2\Delta t^2. The covariance stops changing when the bracket is zero — the Lyapunov equation. A stable AA removes uncertainty through AP+PATAP + PA^{\mathsf T} exactly as fast as the noise adds it through GWGTGWG^{\mathsf T}.
Exact noise for any Δt\Delta t. The noise entering during one interval is ∫0ΔteA(Δt−s)Gw(s) ds\int_0^{\Delta t}e^{A(\Delta t - s)}G\mathbf w(s)\,ds (the variation-of-constants formula above with GwG\mathbf w as input), and its covariance is Qd=∫0ΔteAsGWGTeATs dsQ_d = \int_0^{\Delta t}e^{As}GWG^{\mathsf T}e^{A^{\mathsf T}s}\,ds. Van Loan showed it can be read off one more augmented exponential, exp⁡([−AGWGT0AT]Δt)\exp\Big(\begin{bmatrix} -A & GWG^{\mathsf T} \\ 0 & A^{\mathsf T} \end{bmatrix}\Delta t\Big), and Φ\Phi multiplied by its top-right block is QdQ_d. That is how the lab's 4000 oscillators are stepped without any discretisation error.
The greatest shortcoming of the human race is man's inability to understand the exponential function.
– Albert A. Bartlett, physicist, reported in Hewitt's Conceptual Physics (1981)
💡Where you will use this
Spacecraft Attitude Control and Optimal Control: every state-space design — pole placement, LQR, observers — starts from Φ\Phi and Γ\Gamma. GNSS and Navigation: the Kalman filter's state transition matrix and its discrete process noise, computed with Van Loan's augmented exponential. Structural Dynamics: exact modal time-stepping for linear structures. Electrical and Power Systems: small-signal stability of power grids is an eigenvalue problem of AA. Rocket Propulsion: linearised models of feed-system and combustion instabilities.
📚Further reading
  • Cleve Moler and Charles Van Loan, 'Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later', SIAM Review 45(1), 3–49 (2003) — the survey of the cold open, updated.
  • Nicholas Higham, Functions of Matrices: Theory and Computation (SIAM, 2008), chapter 10, 'Matrix Exponential' — scaling and squaring as expm does it.
  • Charles Van Loan, 'Computing integrals involving the matrix exponential', IEEE Transactions on Automatic Control 23(3), 395–404 (1978) — the augmented-matrix formulas for Γ and for the discrete noise covariance.
↓But there is a problem
Every problem in this unit was an initial value problem: the whole state given at one instant, then marched forward in time. A beam does not work that way. Its deflection is known at both ends — zero at two supports, or zero deflection and zero slope at a clamped end and zero moment at a free one — and the equation must be satisfied everywhere in between. A heat shield's temperature is set on both faces. The pressure in an aquifer is fixed on every boundary of the domain. These are boundary value problems, and when the unknown varies in space and time, partial differential equations. They are the next unit — and the place where most engineering simulation actually happens.