Linear Systems and the Matrix Exponential

With no control applied, an LTI system is just the homogeneous linear system \dot{x} = A x. In one dimension, \dot{x} = a x has the famous solution x(t) = e^{at} x_0. The miracle of state space is that the vector case looks identical — we only have to make sense of the exponential of a matrix:

\dot{x} = A x, \quad x(0) = x_0 \qquad\Longrightarrow\qquad x(t) = e^{At}\, x_0.

This one object, the matrix exponential e^{At}, propagates any initial state forward in time. The rest of the page is about what it is and how to compute it.

Defining the matrix exponential

We define e^{At} by the very same power series that defines the scalar exponential e^x = \sum_{n\ge0} x^n/n! — just feed it a matrix instead of a number:

e^{At} = I + At + \frac{(At)^2}{2!} + \frac{(At)^3}{3!} + \cdots = \sum_{n=0}^{\infty} \frac{(At)^n}{n!}.

Every term is a matrix (powers of A make sense), the series converges for every A and t, and differentiating term by term gives \frac{d}{dt}e^{At} = A e^{At} — which is exactly why x(t) = e^{At}x_0 solves \dot{x} = Ax. At t = 0 the series collapses to e^{A\cdot 0} = I, so the initial condition is met.

Computing it by diagonalization

Summing an infinite series of matrix powers by hand is hopeless — but if A is diagonalizable the series collapses beautifully.

Step 1 — diagonalize. Suppose A = V \Lambda V^{-1}, where the columns of V are eigenvectors and \Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n) holds the eigenvalues.

Step 2 — note that powers telescope. The inner factors cancel, A^k = (V\Lambda V^{-1})^k = V \Lambda^k V^{-1}, because every V^{-1}V = I in the middle collapses.

Step 3 — substitute into the series and factor. Pull the constant V and V^{-1} outside the sum:

e^{At} = \sum_{k=0}^\infty \frac{V \Lambda^k V^{-1} t^k}{k!} = V \left( \sum_{k=0}^\infty \frac{(\Lambda t)^k}{k!} \right) V^{-1} = V\, e^{\Lambda t}\, V^{-1}.

Step 4 — exponentiate the diagonal entrywise. A diagonal matrix's powers are just the entries raised to that power, so its exponential is the entrywise scalar exponential:

e^{\Lambda t} = \begin{bmatrix} e^{\lambda_1 t} & & \\ & \ddots & \\ & & e^{\lambda_n t} \end{bmatrix}.

Worked 2×2. Take A = \begin{bmatrix} -1 & 0 \\ 0 & -3 \end{bmatrix}, already diagonal with \lambda_1 = -1, \lambda_2 = -3. Then

e^{At} = \begin{bmatrix} e^{-t} & 0 \\ 0 & e^{-3t} \end{bmatrix}, \qquad x(t) = e^{At} x_0 = \begin{bmatrix} e^{-t} x_{0,1} \\ e^{-3t} x_{0,2} \end{bmatrix}.

Both coordinates decay to zero, the second three times faster. The eigenvalues are the decay rates — which is the whole story of stability.

Stability and the response to input

Because x(t) = V e^{\Lambda t} V^{-1} x_0 is built from the scalars e^{\lambda_i t}, each mode grows or shrinks according to its eigenvalue. Writing \lambda = \sigma + i\omega, the magnitude is |e^{\lambda t}| = e^{\sigma t}: it decays exactly when the real part \sigma is negative. Hence every initial state is driven to zero precisely when every eigenvalue sits in the left half-plane.

x(t) \to 0 \text{ for every } x_0 \iff \operatorname{Re}(\lambda_i) < 0 \text{ for all } i.

Switch the control back on and the same exponential solves the full equation. Treating Bu as a forcing term and applying variation of constants (an integrating factor e^{-At}) gives the complete state at any time as the free response plus a convolution of the input:

x(t) = \underbrace{e^{At} x_0}_{\text{free response}} + \underbrace{\int_0^t e^{A(t-s)} B\, u(s)\, ds}_{\text{forced response}}.

\dot{x} = Ax, x(0)=x_0 has the unique solution x(t) = e^{At}x_0.
e^{At} = \sum_{n\ge0}(At)^n/n!, with e^{A\cdot 0} = I and \frac{d}{dt}e^{At} = Ae^{At}.
If A = V\Lambda V^{-1} then e^{At} = V e^{\Lambda t} V^{-1} with e^{\Lambda t} = \operatorname{diag}(e^{\lambda_1 t}, \dots, e^{\lambda_n t}).
Stable (x(t)\to0 for all x_0) \iff \operatorname{Re}(\lambda_i) < 0 for every eigenvalue.
With input, x(t) = e^{At}x_0 + \int_0^t e^{A(t-s)}Bu(s)\,ds.

Decay or blow up

A single eigenmode evolves as the scalar x(t) = e^{\lambda t} x_0, here with x_0 = 1. Slide the eigenvalue \lambda: with \lambda < 0 the mode decays to zero — the stable case — and the more negative \lambda the faster the collapse; at \lambda = 0 it sits forever at 1; with \lambda > 0 it runs away to infinity — unstable. A whole linear system is stable exactly when every one of its modes, like this one, has a negative rate.

For all its clean theory, the matrix exponential is notoriously treacherous to compute accurately. Cleve Moler and Charles Van Loan's celebrated 1978 paper “Nineteen Dubious Ways to Compute the Exponential of a Matrix” (revisited in 2003) surveyed method after method — Taylor truncation, diagonalization, scaling and squaring, Padé approximants — and found each can fail badly: diagonalization is unreliable when eigenvectors are nearly parallel, naïve series sums lose precision for large \|At\|. The lesson is a recurring one in applied mathematics: an object can be mathematically transparent yet numerically delicate, and getting e^{At} right in software is a craft of its own.