51. Systems of linear ODEs via eigen-decomposition (bridge to calculus course)
This closing lesson returns to where three Calculus IV lessons pointed here in the first place. §13.8 solved \vec X'=A\vec X using eigenvalues and eigenvectors as a self-contained, from-scratch tool. Now that Modules 19 and 21 have built the full theory around them, this lesson revisits the same equation and shows every step of that calculus-course solution was really Module 19's diagonalization, applied to a differential equation instead of a static computation.
The equation, and the diagonalization shortcut
For \vec X'=A\vec X with A diagonalizable, A=PDP^{-1} (§19.1). Substituting \vec Y=P^{-1}\vec X (a change of basis, §17.3, into eigenvector coordinates):
\vec X'=A\vec X\ \Longrightarrow\ P\vec Y'=AP\vec Y\ \Longrightarrow\ \vec Y'=P^{-1}AP\vec Y=D\vec Y
D\vec Y=D diagonal decouples the system entirely: each y_i'=\lambda_iy_i is an independent scalar equation (Calculus I's y'=ky), with solution y_i(t)=c_ie^{\lambda_it}. Transforming back, \vec X=P\vec Y:
\vec X(t)=c_1e^{\lambda_1t}\vec v_1+c_2e^{\lambda_2t}\vec v_2+\cdots+c_ne^{\lambda_nt}\vec v_n
— exactly §13.8's solution formula, now derived from diagonalization directly rather than guessed via the ansatz \vec X=e^{\lambda t}\vec v. The two routes agree because they're the same fact: diagonalizing A is finding the directions along which the system decouples into independent scalar exponentials.
The matrix exponential
This solution has a compact closed form using the matrix exponential, e^{At}=\sum_{k=0}^\infty\frac{(At)^k}{k!} (the ordinary Taylor series for e^x, with x replaced by the matrix At — well defined because matrix powers and scalar multiplication are both already in hand, §16.2). The general solution is simply \vec X(t)=e^{At}\vec X(0), and diagonalization computes it in closed form exactly as §19.1's power formula did for A^k:
e^{At}=Pe^{Dt}P^{-1},\qquad e^{Dt}=\operatorname{diag}(e^{\lambda_1t},\dots,e^{\lambda_nt})
(substituting A=PDP^{-1} into the series and cancelling every interior P^{-1}P, identically to §19.1's A^k=PD^kP^{-1} derivation — the only difference is summing an infinite series of powers instead of taking one finite power.) This single formula reproduces §13.8's node/saddle/spiral/center classification directly: the eigenvalues' signs (real, same/opposite sign; complex) determine how e^{\lambda_it} behaves as t\to\infty, exactly as that lesson described geometrically.
Doing it in Python
import numpy as np
A = np.array([[1., 1.], [4., 1.]]) # from Calculus IV Section 13.8
eigvals, eigvecs = np.linalg.eig(A)
print("eigenvalues:", [round(v.real, 4) for v in eigvals])
print("classification:", "saddle" if eigvals[0]*eigvals[1] < 0 else "node/spiral/center")
# X(0) = (1, 0); solve for c1, c2 using X(0) = c1*v1 + c2*v2
c = np.linalg.solve(eigvecs, np.array([1., 0.]))
print("coefficients:", [round(v, 4) for v in c.tolist()])
# X(t) at t=0.5, reconstructed from the eigen-solution
t = 0.5
X_t = sum(c[i] * np.exp(eigvals[i]*t) * eigvecs[:, i] for i in range(2))
print("X(0.5) via eigen-solution:", [round(v.real, 4) for v in X_t])
eigenvalues: [3.0, -1.0]
classification: saddle
coefficients: [1.118, -1.118]
X(0.5) via eigen-solution: [2.5441, 3.8752]
Confirming the matrix exponential formula matches directly:
import numpy as np
import scipy.linalg as la
A = np.array([[1., 1.], [4., 1.]])
X0 = np.array([1., 0.])
t = 0.5
X_expm = la.expm(A * t) @ X0
print("X(0.5) via matrix exponential:", [round(v, 4) for v in X_expm.tolist()])
X(0.5) via matrix exponential: [2.5441, 3.8752]
Worked example
Solve \vec X'=A\vec X, \vec X(0)=(2,1), for A=\begin{pmatrix}-1&0\\0&-3\end{pmatrix} (already diagonal).
A is diagonal, so P=I, D=A — eigenvalues -1,-3 read directly off the diagonal (§19.1's shortcut), eigenvectors \vec e_1,\vec e_2.
\vec X(t)=c_1e^{-t}(1,0)+c_2e^{-3t}(0,1)
\vec X(0)=(c_1,c_2)=(2,1).
\boxed{\vec X(t)=(2e^{-t},\,e^{-3t})}
Sanity check. Both eigenvalues negative (same sign) \Rightarrow stable node (§13.8's classification): every component decays to \vec0 as t\to\infty, confirmed directly by inspection since e^{-t}\to0 and e^{-3t}\to0. The second component decays faster (-3 vs -1), so the trajectory bends toward the x-axis as t grows — consistent with §13.8's "trajectories align with the less negative (slower-decaying) eigendirection" behavior for a stable node.
Your turn
1. For A diagonalizable with a repeated positive eigenvalue \lambda=2 (algebraic multiplicity 2, and — crucially — geometric multiplicity 2 also, so A is not defective), what does every solution do as t\to\infty?
2. Why does §19.4's Jordan form matter for solving \vec X'=A\vec X when A is defective?
3. True or false: e^{At} is always invertible, for any square matrix A and any t.
Solutions
1. Every solution is a combination of e^{2t}\vec v_1 and e^{2t}\vec v_2 (both eigenvectors sharing the same eigenvalue) — so every component of every solution grows without bound as t\to\infty (an unstable node, the higher-dimensional version of §13.8's classification, since both directions grow at the identical rate e^{2t}).
2. Diagonalization requires a full eigenbasis (§19.1), which a defective matrix doesn't have — the derivation at the top of this lesson breaks down at "P^{-1}AP=D" specifically. §19.4's Jordan form A=PJP^{-1} still applies unconditionally, and solving \vec Y'=J\vec Y for a Jordan block produces solutions with an extra t\cdot e^{\lambda t} term (from the block's off-diagonal 1, coupling consecutive equations instead of fully decoupling them) — the same te^{\lambda t} term that appears in a calculus course's "repeated root" case for a single second-order ODE, now understood as a Jordan block of size 2 in disguise.
3. True. \det(e^{At})=e^{\operatorname{tr}(At)} (a fact about the matrix exponential, not shown here) — and e^x\neq0 for every real (or complex) x, so this determinant is never zero, regardless of A or t. Geometrically, the flow of a linear ODE system never collapses space to a lower dimension; it only stretches, rotates, or shrinks it, always reversibly.
Check yourself in code
Solve \vec X'=A\vec X, \vec X(0)=(1,1) for A=\begin{pmatrix}2&0\\0&-1\end{pmatrix}, and evaluate \vec X(1).
Print exactly this:
eigenvalues: [-1.0, 2.0]
X(1) = [7.3891, 0.3679]
import numpy as np
A = np.array([[2., 0.], [0., -1.]])
eigvals = np.linalg.eigvalsh(A)
print("eigenvalues:", eigvals.tolist())
# X(0) = (1, 1) = 1*e1 + 1*e2, so c1=1, c2=1 directly
t = 1.0
X_t = np.array([1*np.exp(2*t), 1*np.exp(-1*t)])
# print "X(1) = [...]" with X_t rounded to 4 decimal places
import numpy as np
A = np.array([[2., 0.], [0., -1.]])
eigvals = np.linalg.eigvalsh(A)
print("eigenvalues:", eigvals.tolist())
t = 1.0
X_t = np.array([1*np.exp(2*t), 1*np.exp(-1*t)])
print("X(1) =", [round(v, 4) for v in X_t.tolist()])
\vec X'=A\vec X solves by diagonalizing A=PDP^{-1}, decoupling into independent scalar exponentials in eigenvector coordinates — exactly Calculus IV §13.8's formula, now derived from Module 19's machinery rather than guessed. The matrix exponential e^{At}=Pe^{Dt}P^{-1} packages the whole solution into one formula, and Jordan form (§19.4) extends it, with an extra te^{\lambda t} term, to defective matrices diagonalization alone cannot handle.
This closes the Linear Algebra course. Every module built toward this point: Module 16's elimination and determinants, Module 17's abstract vector spaces, Module 18's transformations, Module 19's eigenvalues, Module 20's geometry, Module 21's spectral theory, Module 22's numerical methods, and Module 23's bilinear structures all converge in this closing module's five applications — Markov chains, PCA, graph theory, regression, and this lesson's differential equations — each one a direct, traceable payoff of a tool built somewhere in the eight modules before it.