Systems of linear ODEs and the matrix exponential

Intuition [Beginner]

Imagine two water tanks connected by a pipe. Water flows from tank $1$ to tank $2$ at one rate, and from tank $2$ back to tank $1$ at a different rate. The water level in each tank changes over time, but the two levels are coupled: what happens in tank $1$ affects tank $2$ and vice versa.

A system of ODEs captures exactly this situation. Instead of one unknown function $y(t)$, there are several — say $x_1(t)$ and $x_2(t)$ — and each derivative depends on both unknowns. Writing the unknowns as a column vector $\mathbf{x}$ and the coupling coefficients as a matrix $A$, the entire system collapses to one compact equation: ${\mathbf{x}}^{\prime }=A\mathbf{x}$.

The solution is the matrix exponential $e^{At}$, which acts like the ordinary exponential $e^{at}$ but works for matrices. It turns a coupled system into a single formula: $\mathbf{x}(t)=e^{At}\mathbf{x}(0)$. Why does this concept exist? Because most real-world systems have multiple interacting components, and the matrix exponential gives one unified way to solve them all.

Visual [Beginner]

The picture shows two curves in the $(x_1,x_2)$-plane, starting from different initial points. Both spiral inward toward the origin, because the matrix $A$ has eigenvalues with negative real parts. The spiral shape comes from the complex eigenvalues, and the inward drift comes from the negative real parts.

Each trajectory is the path traced by the tip of the vector $\mathbf{x}(t)=e^{At}\mathbf{x}(0)$ as time $t$ increases from $0$.

Worked example [Beginner]

Solve the system $x_1^{\prime }=-x_1+2x_2$, $x_2^{\prime }=-2x_1-x_2$ with initial conditions $x_1(0)=1$, $x_2(0)=0$.

Step 1. Write as a matrix equation: ${\mathbf{x}}^{\prime }=A\mathbf{x}$ where $A=\left(\begin{array}{cc}-1& 2\\ -2& -1\end{array}\right)$.

Step 2. Find the eigenvalues of $A$: $\det (A-rI)=(r+1{)}^2+4=0$, so $r=-1\pm 2i$.

Step 3. The real parts are $-1$ (decay), and the imaginary parts are $\pm 2$ (oscillation). The solution is $x_1(t)=e^{-t}\cos (2t)$, $x_2(t)=-e^{-t}\sin (2t)$.

What this tells us: the two water tanks exchange fluid in an oscillating pattern that decays over time. The matrix exponential captures both the oscillation (from the imaginary part of the eigenvalues) and the decay (from the negative real part).

Check your understanding [Beginner]

Formal definition [Intermediate+]

A first-order linear homogeneous system of ODEs with constant coefficients is

$${\mathbf{x}}^{\prime }(t)=A\mathbf{x}(t),$$

where $\mathbf{x}:\mathbb{R}\to {\mathbb{R}}^n$ (or ${\mathbb{C}}^n$) and $A$ is an $n\times n$ constant matrix.

The matrix exponential is defined by the convergent power series

$$e^{At}=\sum _{k=0}^{\infty }\frac{(At{)}^k}{k!}=I+At+\frac{A^2{t}^2}{2!}+\frac{A^3{t}^3}{3!}+\cdots ,$$

which converges for every $n\times n$ matrix $A$ and every $t\in \mathbb{R}$ (because $\Vert At{\Vert }^k/k!\le (\Vert A\Vert |t|{)}^k/k!$ and the scalar exponential series converges).

Fundamental matrix. An $n\times n$ matrix function $\Phi (t)$ is a fundamental matrix for ${\mathbf{x}}^{\prime }=A\mathbf{x}$ if each column of $\Phi$ is a solution and $\Phi (t)$ is invertible for all $t$. The matrix exponential $e^{At}$ is the unique fundamental matrix satisfying $\Phi (0)=I$.

Counterexamples to common slips

• $e^{A+B}\ne e^A{e}^B$ in general. This identity holds only when $AB=BA$. For non-commuting matrices, the correct formula is the Baker-Campbell-Hausdorff expansion.

• ${\mathbf{x}}^{\prime }=A(t)\mathbf{x}$ with variable $A(t)$ does not have solution $e^{\int A(t)dt}$. The matrix exponential formula applies only to constant-coefficient systems. Variable-coefficient systems require the Magnus expansion or other techniques.

• $e^{At}$ is not computed by exponentiating each entry of $At$. The matrix exponential is defined by the power series, not entrywise exponentiation. For diagonal matrices, $e^{At}$ does equal the diagonal matrix of exponentials, but this fails for general matrices.

Key theorem with proof [Intermediate+]

Theorem (Matrix exponential solution). The unique solution of the initial value problem ${\mathbf{x}}^{\prime }=A\mathbf{x}$, $\mathbf{x}(0)={\mathbf{x}}_0$ is

$$\mathbf{x}(t)=e^{At}{\mathbf{x}}_0.$$

Moreover, $e^{At}$ satisfies:

1. $\frac{d}{dt}{e}^{At}=A{e}^{At}=e^{At}A$.
2. $e^{A\cdot 0}=I$.
3. $e^{A(t+s)}=e^{At}{e}^{As}$ for all $t,s\in \mathbb{R}$.

Proof of (1). Differentiate the power series term by term:

$$\frac{d}{dt}{e}^{At}=\frac{d}{dt}\sum _{k=0}^{\infty }\frac{A^k{t}^k}{k!}=\sum _{k=1}^{\infty }\frac{A^k\cdot k{t}^{k-1}}{k!}=\sum _{k=1}^{\infty }\frac{A^k{t}^{k-1}}{(k-1)!}.$$

Re-indexing $j=k-1$:

$$\frac{d}{dt}{e}^{At}=A\sum _{j=0}^{\infty }\frac{A^j{t}^j}{j!}=A{e}^{At}.$$

Since $A^k=A\cdot A^{k-1}=A^{k-1}\cdot A$, the same computation gives $e^{At}A$. Term-by-term differentiation is justified because the series converges uniformly on every bounded interval (the derivative series has the same radius of convergence).

Proof of (2). At $t=0$: $e^{A\cdot 0}={\sum }_{k=0}^{\infty }{A}^k\cdot 0^k/k!=I+0+0+\cdots =I$.

Proof of (3). Fix $s$ and let $F(t)=e^{A(t+s)}$. Then $F^{\prime }(t)=A{e}^{A(t+s)}=AF(t)$ by property (1), and $F(0)=e^{As}$. The function $G(t)=e^{At}{e}^{As}$ satisfies $G^{\prime }(t)=A{e}^{At}{e}^{As}=AG(t)$ and $G(0)=e^{As}$. By uniqueness of solutions to ${\mathbf{y}}^{\prime }=A\mathbf{y}$ with given initial condition, $F(t)=G(t)$.

Proof of the solution formula. The function $\mathbf{x}(t)=e^{At}{\mathbf{x}}_0$ satisfies ${\mathbf{x}}^{\prime }(t)=A{e}^{At}{\mathbf{x}}_0=A\mathbf{x}(t)$ by property (1), and $\mathbf{x}(0)=e^{A\cdot 0}{\mathbf{x}}_0=I{\mathbf{x}}_0={\mathbf{x}}_0$ by property (2). Uniqueness follows from the existence-uniqueness theorem for linear ODE systems. $\square$

Bridge. The foundational reason the matrix exponential solves the system is that the power series $\sum (At{)}^k/k!$ differentiates to $A$ times itself, and this is exactly the matrix analogue of the scalar identity $(d/dt)e^{at}=a{e}^{at}$. The central insight is that $e^{At}$ plays the role of the scalar exponential $e^{at}$, but in $n$ dimensions simultaneously. This result builds toward Lyapunov stability theory where the eigenvalues of $A$ determine whether $e^{At}$ decays, and appears again in 02.06.02 via the companion-matrix reduction where a scalar nth-order ODE becomes a system whose solution is $e^{At}$. Putting these together, the diagonalisation identity $e^{At}=P{e}^{Dt}{P}^{-1}$ (where $D$ is diagonal with the eigenvalues) is the bridge that connects the abstract power series to concrete computations: the eigenvalues of $A$ control the behaviour of every solution.

Exercises [Intermediate+]

Advanced results [Master]

Theorem 1 (Diagonalisation formula). If $A=PD{P}^{-1}$ with $D=\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_n)$, then $e^{At}=P\mathrm{diag}(e^{{\lambda }_{1}t},\dots ,e^{{\lambda }_{n}t})P^{-1}$. The columns of $P$ are eigenvectors of $A$, and each eigenvalue ${\lambda }_j={\alpha }_j+i{\beta }_j$ contributes a factor $e^{{\alpha }_{j}t}(\cos {\beta }_{j}t+i\sin {\beta }_{j}t)$ to the solution.

Theorem 2 (Jordan canonical form and the matrix exponential). For any matrix $A$, there exists an invertible $P$ such that $A=PJ{P}^{-1}$ where $J$ is the Jordan form: a block-diagonal matrix of Jordan blocks $J_i(\lambda )=\lambda I+N_i$ where $N_i$ is nilpotent. Then $e^{At}=P{e}^{Jt}{P}^{-1}$ and each Jordan block contributes $e^{\lambda t}{\sum }_{k=0}^{m-1}{N}^k{t}^k/k!$ where $m$ is the block size. This handles non-diagonalisable matrices.

Theorem 3 (Lyapunov stability theorem). The equilibrium $\mathbf{x}=0$ of ${\mathbf{x}}^{\prime }=A\mathbf{x}$ is:

• Asymptotically stable if and only if every eigenvalue of $A$ has negative real part.
• Stable (but not asymptotically) if and only if every eigenvalue has non-positive real part and every eigenvalue with zero real part is semisimple (its algebraic and geometric multiplicities coincide).
• Unstable otherwise.

The Lyapunov equation $A^{T}P+PA=-Q$ (for any positive-definite $Q$) has a unique positive-definite solution $P$ if and only if $A$ is asymptotically stable.

Theorem 4 (Controllability and the Kalman rank condition). The system ${\mathbf{x}}^{\prime }=A\mathbf{x}+B\mathbf{u}$ is controllable (any initial state can be driven to any target state in finite time by a suitable input $\mathbf{u}$) if and only if the controllability matrix $\mathcal{C}=[B\mid AB\mid A^{2}B\mid \cdots \mid A^{n-1}B]$ has full row rank $n$. This is the Kalman rank condition.

Theorem 5 (Putzer's algorithm). The matrix exponential can be computed without finding eigenvectors, using only the eigenvalues ${\lambda }_1,\dots ,{\lambda }_n$ and the powers of $A$. Define $P_0=I$, $P_k=(A-{\lambda }_{k}I)P_{k-1}$, and solve the cascade of first-order linear ODEs $r_1^{\prime }={\lambda }_1{r}_1$, $r_k^{\prime }={\lambda }_k{r}_k+r_{k-1}$ for $k\ge 2$, with $r_1(0)=1$, $r_k(0)=0$ for $k\ge 2$. Then $e^{At}={\sum }_{k=1}^n{r}_k(t)P_{k-1}$.

Theorem 6 (Non-homogeneous systems and variation of constants). The solution of ${\mathbf{x}}^{\prime }=A\mathbf{x}+\mathbf{f}(t)$ with $\mathbf{x}(t_0)={\mathbf{x}}_0$ is

$$\mathbf{x}(t)=e^{A(t-t_0)}{\mathbf{x}}_0+{\int }_{t_0}^t{e}^{A(t-s)}\mathbf{f}(s)ds.$$

This is the variation of constants (or Duhamel) formula.

Synthesis. The foundational reason the matrix exponential works is that the power series $\sum (At{)}^k/k!$ converges for every matrix and differentiates to $A$ times itself, and this is exactly the content that makes $\mathbf{x}(t)=e^{At}\mathbf{x}(0)$ solve the system. The central insight is that the eigenvalues of $A$ control the qualitative behaviour: decay, growth, oscillation, or a mix. Putting these together with the Jordan form, even non-diagonalisable matrices are handled by the nilpotent correction $e^{\lambda t}\sum N^k{t}^k/k!$. The bridge is between the abstract power series and the concrete eigenvalue decomposition $A=PD{P}^{-1}$ that reduces $e^{At}$ to scalar exponentials. This identification appears again in the companion matrix reduction of 02.06.02 where scalar nth-order equations become systems, the pattern generalises to the non-homogeneous Duhamel formula, and builds toward Lyapunov stability and controllability where the eigenvalue locations determine whether the system can be stabilised or controlled.

Full proof set [Master]

Proposition 1 (The power series for $e^{At}$ converges). For any $n\times n$ matrix $A$ and any $t\in \mathbb{R}$, the series ${\sum }_{k=0}^{\infty }(At{)}^k/k!$ converges.

Proof. Let $\Vert \cdot \Vert$ be any submultiplicative matrix norm (e.g., the operator norm). Then $\Vert (At{)}^k/k!\Vert \le (\Vert A\Vert \cdot |t|{)}^k/k!$. The series $\sum (\Vert A\Vert |t|{)}^k/k!=e^{\Vert A\Vert |t|}$ converges. By the Weierstrass M-test, the matrix series converges absolutely, hence converges. $\square$

Proposition 2 ($e^{A+B}=e^A{e}^B$ when $AB=BA$). If $AB=BA$, then $e^{A+B}=e^A{e}^B$.

Proof. When $AB=BA$, the binomial theorem applies to matrices: $(A+B{)}^k={\sum }_{j=0}^k(\genfrac{}{}{0ex}{}{k}{j})A^j{B}^{k-j}$. So

$$e^{A+B}=\sum _{k=0}^{\infty }\frac{1}{k!}\sum _{j=0}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right)A^j{B}^{k-j}=\sum _{k=0}^{\infty }\sum _{j=0}^k\frac{A^j}{j!}\frac{B^{k-j}}{(k-j)!}.$$

This double sum equals $({\sum }_j{A}^j/j!)({\sum }_m{B}^m/m!)=e^A{e}^B$ by the Cauchy product formula, since both series converge absolutely. $\square$

Connections [Master]

• n-th-order linear ODE with constant coefficients 02.06.02. The companion matrix reduction transforms every scalar nth-order ODE into a first-order system ${\mathbf{x}}^{\prime }=A\mathbf{x}$, and the eigenvalues of the companion matrix are exactly the roots of the characteristic polynomial. The matrix exponential $e^{At}$ generalises the scalar exponential solutions $e^{rt}$ to the coupled setting, recovering the scalar theory when $n=1$.

• First-order linear and separable ODEs 02.08.01. The existence-uniqueness theorem for first-order linear ODE systems provides the theoretical foundation for the solution formula $\mathbf{x}(t)=e^{At}\mathbf{x}(0)$. The scalar integrating-factor method from that unit is the $1\times 1$ case of the matrix exponential formula.

• Second-order linear ODE with constant coefficients 02.08.02. A second-order equation $y^{\prime \prime }+b{y}^{\prime }+cy=0$ becomes a $2\times 2$ system via the companion matrix $A=\left(\begin{array}{cc}0& 1\\ -c& -b\end{array}\right)$. The eigenvalues of $A$ are the roots of $r^2+br+c=0$, and $e^{At}$ produces the same oscillatory, exponential, or critically-damped solutions as the scalar characteristic-polynomial method.

Historical & philosophical context [Master]

Sylvester 1883, in On the Equation to the Secular Inequalities in the Planetary Theory ^{[Sylvester1883]}, introduced the matrix exponential in the context of solving systems of linear differential equations arising from celestial mechanics. His formulation recognised that the power series $\sum (At{)}^k/k!$ generalises the scalar exponential to matrices and satisfies the key differential-equation property.

Peano 1888, in Integration par series des equations differentielles lineaires ^[Peano1888], gave the rigorous power-series definition of the matrix exponential and proved the fundamental properties: term-by-term differentiation yields $d/dt{e}^{At}=A{e}^{At}$, and the composition $e^{A(t+s)}=e^{At}{e}^{As}$. Peano's treatment was the first to establish the matrix exponential as the canonical solution operator for linear ODE systems with constant coefficients. The modern synthesis — diagonalisation, Jordan form, stability theory, and the Kalman controllability criterion — was developed through the mid-twentieth century, with the Jordan-form computation of $e^{At}$ appearing in Gantmacher's Theory of Matrices (1959) and the controllability theory in Kalman's foundational papers of 1960-1963.

Bibliography [Master]

@article{Sylvester1883,
  author = {Sylvester, James Joseph},
  title = {Sur les quantites formant un groupe de nonnes analogues aux quaternions de Hamilton},
  journal = {Comptes Rendus de l'Academie des Sciences, Paris},
  volume = {94},
  year = {1883},
  pages = {1336--1340},
}

@article{Peano1888,
  author = {Peano, Giuseppe},
  title = {Integration par series des equations differentielles lineaires},
  journal = {Mathematische Annalen},
  volume = {32},
  year = {1888},
  pages = {450--456},
}

@book{BoyceDiPrima2012,
  author = {Boyce, William E. and DiPrima, Richard C.},
  title = {Elementary Differential Equations and Boundary Value Problems},
  publisher = {Wiley},
  year = {2012},
  edition = {10th},
}

@book{CoddingtonLevinson1955,
  author = {Coddington, Earl A. and Levinson, Norman},
  title = {Theory of Ordinary Differential Equations},
  publisher = {McGraw-Hill},
  year = {1955},
}

@book{Gantmacher1959,
  author = {Gantmacher, Felix R.},
  title = {The Theory of Matrices},
  publisher = {Chelsea},
  year = {1959},
}