常微分方程（六）：线性微分方程组

一个方程描述一个量，但现实很少这么简单。 兔子和狼的数量相互制约， RLC 网络中电流电压彼此耦合，化学反应里各物质浓度互相依赖。只要两个未知量共存于一组方程，问题就变成方程组，标量公式 $$y'=ay$$ 就不够用了。

线性情况的奇妙之处在于：标量解 $y(t)=e^{at}y_0$ 可直接推广到方程组——只需理解矩阵指数 $e^{At}$ 的含义。线性代数与微分方程在此融合为一，矩阵 $$A$$ 的特征结构揭示了解的长期行为、流场几何以及简正模态和拍频等物理现象背后的本质。

总结#

写成矩阵形式 $\mathbf{x}'=A\mathbf{x}$ ，不只是记号
矩阵指数 $e^{At}$ ：定义、性质，三种计算方法
特征值法、复特征值，特征向量不足时的处理
二维相图分类与迹–行列式稳定性图
非齐次方程组和 Duhamel 公式
耦合振子：简正模态、拍频、能量交换

前置知识#

线性代数：特征值、特征向量、基变换
第三章：二阶线性 ODE （均可化为二维系统）

从两个耦合方程到一个向量方程#

x' = 2x - y, \qquad y' = x + 0.5\,y.

\mathbf{x}' = A\mathbf{x}, \qquad A = \begin{pmatrix} 2 & -1 \\ 1 & 0.5 \end{pmatrix}.

这不是简单换记号，而是换了视角。两条耦合时间序列变成平面上一条轨迹 $\mathbf{x}(t)$ ，几何代替了代数。

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

但要先定义清楚矩阵的指数才有意义。这就是本章的核心内容。

矩阵指数#

幂级数定义#

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

对任意方阵 $$A$$ 和 $t\in\mathbb{R}$ ，该级数按算子范数收敛。因为 $\|(At)^k/k!\| \le (\|A\|t)^k/k!$ ，而标量级数收敛。

矩阵指数部分和作用在向量上的轨迹，逐渐逼近真实旋转；谱范数误差随项数超指数衰减。

图 1. 左： $$A$$ 是旋转生成元时， $\sum_{k=0}^N(At)^k/k!\;\mathbf{x}_0$ 轨迹逼近单位圆。右：谱范数误差随项数超指数衰减——低阶 Padé 逼近效果好的原因。

常用三条性质#

写完级数后，论证几乎只用到三点：

性质	用途
$e^{A\cdot 0} = I$	初值条件自动满足
$\dfrac{d}{dt}e^{At} = Ae^{At}$	这是 $\mathbf{x}(t)=e^{At}\mathbf{x}_0$ 解 $\mathbf{x}'=A\mathbf{x}$ 的关键
$e^{At}e^{As}=e^{A(t+s)}$	流构成单参数群；特别地 $(e^{At})^{-1}=e^{-At}$

注意： $e^{A+B}=e^A e^B$ 仅当 $$AB=BA$$ 时成立。这是矩阵版的 $\sin(x+y)\ne\sin x+\sin y$ ——非交换性有后果。

实际计算三种方法#

截断幂级数： 概念简单，但 $\|At\|$ 大时数值差。
特征值分解： 若 $$A$$ 可对角化为 $A=PDP^{-1}$ ， $D=\mathrm{diag}(\lambda_i)$ ，则

e^{At} = P\,\mathrm{diag}(e^{\lambda_1 t},\dots,e^{\lambda_n t})\,P^{-1}.

理论分析常用此公式。 3. 缩放–平方 + Padé 逼近。 工业标准——scipy.linalg.expm、 MATLAB 的 expm 都用它。先算 $e^{At/2^s}$ 的 Padé 逼近，再平方 $$s$$ 次。对大 $\|At\|$ 稳健。

代码示例：

1
2
3
4
5
6
7
8
9
import numpy as np
from scipy.linalg import expm

A = np.array([[2.0, -1.0], [1.0, 0.5]])
print(expm(A * 0.5))        # exp(0.5·A) -- 缩放-平方 Padé

vals, P = np.linalg.eig(A)
expm_via_eig = P @ np.diag(np.exp(vals * 0.5)) @ np.linalg.inv(P)
print(expm_via_eig.real)    # 同一矩阵（差一点浮点噪声）

特征值方法#

特征值分解公式可以完全不用矩阵指数重写。

定理： 若 $$A$$ 的特征对 $(\lambda_i,\mathbf{v}_i)$ 满足 $\mathbf{v}_i$ 线性无关， $\mathbf{x}'=A\mathbf{x}$ 的通解为
$\mathbf{x}(t) = c_1 e^{\lambda_1 t}\mathbf{v}_1 + c_2 e^{\lambda_2 t}\mathbf{v}_2 + \cdots + c_n e^{\lambda_n t}\mathbf{v}_n.$

证明很简单：对每对特征对， $\frac{d}{dt}\bigl(e^{\lambda t}\mathbf{v}\bigr)=\lambda e^{\lambda t}\mathbf{v}=A\bigl(e^{\lambda t}\mathbf{v}\bigr)$ 。线性性完成其余部分。

几何：特征向量是天然坐标轴#

$A=PDP^{-1}$ 的分解有直观的几何解释：用 $$A$$ 作用一个向量，相当于三步操作：

$P^{-1}$ ：把向量转到特征基；
$$D$$ ：沿特征方向按 $\lambda_i$ 拉伸；
$$P$$ ：把结果转回原基。

流 $e^{At}$ 做同样的事，只是拉伸因子换成 $e^{\lambda_i t}$ 。

特征值分解的几何过程：单位圆经过"旋转到特征基—沿轴缩放—旋转回去"三步变成它在 A 下的像。

图 2. 把 $A=PDP^{-1}$ 拆成三步：特征向量让线性映射（以及它的流）退化为对角缩放的“天然坐标轴”。

复特征值给出旋转#

\mathbf{x}_1(t) = e^{\alpha t}(\mathbf{a}\cos\beta t - \mathbf{b}\sin\beta t), \qquad \mathbf{x}_2(t) = e^{\alpha t}(\mathbf{a}\sin\beta t + \mathbf{b}\cos\beta t).

几何上看： $\beta$ 是 $(\mathbf{a},\mathbf{b})$ -平面内的旋转角频率， $\alpha$ 决定轨道是向外螺旋（ $\alpha>0$ ）、向内螺旋（ $\alpha<0$ ）、还是闭合曲线（ $\alpha=0$ ）。

二维相图#

平面上的 $\mathbf{x}'=A\mathbf{x}$ ，特征值 $\lambda_1,\lambda_2$ 决定了原点附近的几何特性。分类表值得记住。

特征值	类型	稳定性
$\lambda_1<\lambda_2<0$ （实）	稳定节点——轨道全指向原点	渐近稳定
$0<\lambda_1<\lambda_2$ （实）	不稳定节点——轨道全离开原点	不稳定
$\lambda_1<0<\lambda_2$ （实）	鞍点——两稳定两不稳定方向	不稳定
$\alpha\pm\beta i$ ， $\alpha<0$	稳定螺旋	渐近稳定
$\alpha\pm\beta i$ ， $\alpha>0$	不稳定螺旋	不稳定
$\pm\beta i$ （纯虚）	中心——闭轨	Lyapunov 稳定但不渐近
$\lambda_1=\lambda_2$ ，两个特征向量	星形节点	$\lambda$ 符号决定
$\lambda_1=\lambda_2$ ，一个特征向量	退化节点	$\lambda$ 符号决定

图 3. 相图动物园。绿色是稳定方向，红色是不稳定方向。特征值直接决定图形类型，其余只是旋转和缩放。

迹–行列式速查#

\lambda^2 - \tau\lambda + \delta = 0, \qquad \tau = \mathrm{tr}\,A = \lambda_1+\lambda_2, \quad \delta = \det A = \lambda_1\lambda_2,

判别式 $\Delta = \tau^2 - 4\delta$ 。 $(\tau,\delta)$ 的位置直接判断平衡点类型：

$\delta < 0$ → 实特征值异号 → 鞍点。
$\delta > 0,\ \Delta > 0$ → 实特征值同号 → 节点（ $\tau$ 符号决定稳定性）。
$\delta > 0,\ \Delta < 0$ → 复特征值 → 螺旋（ $\tau$ 符号决定稳定性）。
$\delta > 0,\ \tau = 0$ → 纯虚特征值 → 中心。
抛物线 $\Delta = 0$ 是重特征值轨迹。

迹-行列式平面上的稳定性图：抛物线 Delta=0 分隔节点与螺旋，det=0 轴是鞍点，tr=0 且 det>0 是中心。

图 4. 稳定性图。单点 $(\tau,\delta)$ 就能判断平衡类型，无需计算特征值。中心（ $\tau=0,\ \delta>0$ ）是非通有情况：任意小扰动都会变成螺旋。

重特征值与广义特征向量#

分类表里藏了个麻烦。代数重数为 2 的 $\lambda$ ，可能有两个线性无关特征向量（罕见，只有 $A=\lambda I$ 时成立，对应星形节点），也可能只有一个（通例，对应退化节点）。

(A - \lambda I)\,\mathbf{w} = \mathbf{v},

\mathbf{x}_2(t) \;=\; e^{\lambda t}\bigl(t\,\mathbf{v} + \mathbf{w}\bigr)

是另一个独立解。 $$t$$ 带来的多项式乘指数增长，正是亏损的代数特征。

e^{(\lambda I + N)t} = e^{\lambda t}\bigl(I + tN + \tfrac12 t^2 N^2 + \cdots\bigr),

$$N$$ 幂零，级数有限项终止，留下 $$t$$ 的多项式。

亏损 2x2 系统：流向唯一的特征方向剪切，第二个独立解多了一项 t·e^{λt}。

图 5. 左：亏损系统 $\dot{\mathbf{x}}=\bigl(\begin{smallmatrix}-1&1\\0&-1\end{smallmatrix}\bigr)\mathbf{x}$ 。唯一特征方向（绿）是轨道贴着进来的方向；广义方向（红虚线）是 $t\mathbf v$ 栖身的地方。右：两个基本解分量随时间变化——注意先升后降的 $te^{\lambda t}$ 形状。

非齐次方程组： Duhamel 公式#

\mathbf{x}' = A\mathbf{x} + \mathbf{g}(t), \qquad \mathbf{x}(0)=\mathbf{x}_0.

\boxed{\;\mathbf{x}(t) \;=\; e^{At}\mathbf{x}_0 \;+\; \int_0^t e^{A(t-\tau)}\,\mathbf{g}(\tau)\,d\tau.\;}

这就是 Duhamel 公式。看被积函数 $e^{A(t-\tau)}\mathbf{g}(\tau)$ ： $\tau$ 时刻，外力给系统一个 $\mathbf{g}(\tau)d\tau$ 的“踢”，随后按自由流 $e^{A(t-\tau)}$ 演化。整体解是所有延迟响应的叠加，即连续时间脉冲响应。

\mathbf{x}(t) = e^{At}\mathbf{x}_0 + \int_0^t e^{A(t-\tau)} B\,\mathbf u(\tau)\,d\tau,

离可控性与矩阵 $\bigl[B,\,AB,\,A^2B,\dots\bigr]$ 只差一步。

应用：耦合振子、简正模态与拍频#

直线上的两个等质量小球，分别用劲度 $$k$$ 的弹簧连到墙上，中间用劲度 $\kappa$ 的弹簧连接：

1
2
墙 —/\/\/— [m] —/\/\/— [m] —/\/\/— 墙
       k        κ        k

\ddot x_1 = -k\,x_1 + \kappa(x_2 - x_1), \qquad \ddot x_2 = -k\,x_2 + \kappa(x_1 - x_2).

\ddot u = -k\,u, \qquad \ddot v = -(k+2\kappa)\,v.

这就是简正模态。同相模态频率 $\omega_s=\sqrt{k}$ ，耦合弹簧不拉伸；反相模态频率 $\omega_a=\sqrt{k+2\kappa}$ ，耦合弹簧做功。任何运动都是两者的叠加。

当 $\kappa \ll k$ 时，模态频率接近。只拨动 1 号小球，两模态等幅激发，缓慢相位差让能量从 1 号转到 2 号再转回。这就是拍频，类似两支微调失谐音叉发声时的强弱起伏。

三条堆叠的时间序列：纯对称模态、纯反对称模态、以及只拨动一个小球时出现的拍频图样。

图 6. 上：对称模态，两小球同相，频率 $\omega_s$ 。中：反对称模态，反相，频率更高 $\omega_a$ 。下：通有初值，两模态轻微失谐，干涉产生拍频，能量周期性转移，包络为虚线 $\cos\bigl(\tfrac{\omega_a-\omega_s}{2}t\bigr)$ 。

这与量子两能级系统的 Rabi 振荡、耦合 LC 电路的能量交换、弱耦合摆钟同步现象，数学本质相同。

稳定性：特征值判据#

线性流 $\mathbf{x}'=A\mathbf{x}$ 的长时间行为由 $$A$$ 的谱完全决定：

定理。
全部 $\mathrm{Re}(\lambda) < 0\ \Longrightarrow\ e^{At}\to 0$ ，原点渐近稳定。
任一 $\mathrm{Re}(\lambda) > 0\ \Longrightarrow$ 原点不稳定。
$\mathrm{Re}(\lambda) = 0$ 需进一步分析（中心、 Jordan 块）。

Liouville 公式更细致： $\det \Phi(t) = \det \Phi(0)\,\exp\bigl(\int_0^t \mathrm{tr}\,A\,d\tau\bigr)$ 。
$\mathrm{tr}\,A < 0$ 表示相空间体积收缩，是耗散系统。
$\mathrm{tr}\,A = 0$ 表示体积守恒，对应 Hamilton 系统。
$\mathrm{tr}\,A > 0$ 表示体积膨胀。

第七章通过线性化将这套理论推广到非线性系统。
Hartman–Grobman 定理指出，远离边界情形时，非线性系统在不动点附近的局部行为等同于其 Jacobi 矩阵的行为。

总结#

概念	核心
向量形式 $\mathbf{x}'=A\mathbf{x}$	解是 $e^{At}\mathbf{x}_0$
矩阵指数	幂级数定义，特征值分解或 Padé 缩放-平方法计算
特征值方法	$\sum c_k e^{\lambda_k t}\mathbf{v}_k$ ，特征模态线性组合
复特征值	实解为 $e^{\alpha t}\bigl(\cos\beta t,\ \sin\beta t\bigr)$ ，螺旋形式
重特征值（亏损）	广义特征向量引入 $te^{\lambda t}$ 解
相图	特征值对应节点 / 螺旋 / 鞍点 / 中心，迹-行列式图概括
非齐次	Duhamel： $\mathbf{x}=e^{At}\mathbf{x}_0+\int_0^t e^{A(t-\tau)}\mathbf{g}(\tau)d\tau$
稳定性	全部 $\mathrm{Re}(\lambda)<0\ \Leftrightarrow$ 渐近稳定

练习题#

基础

用幂级数和特征值分解算 $A=\bigl(\begin{smallmatrix}0&1\\-1&0\end{smallmatrix}\bigr)$ 的 $e^{At}$ ，验证结果一致。
解 $\mathbf{x}'=\bigl(\begin{smallmatrix}1&2\\0&3\end{smallmatrix}\bigr)\mathbf{x}$ ，初值 $\mathbf{x}(0)=(1,1)^\top$ 。
判断平衡点类型：(a) $A=\bigl(\begin{smallmatrix}-1&2\\-2&-1\end{smallmatrix}\bigr)$ ；(b) $A=\bigl(\begin{smallmatrix}1&0\\0&-2\end{smallmatrix}\bigr)$ 。

进阶

证明 $\det e^A = e^{\mathrm{tr}\,A}$ 。提示：在 $\mathbb C$ 上对 $$A$$ 做 Schur 上三角化。
求三个等质量小球（一字排开，弹簧劲度 $$k$$ ）的简正模态频率，区分同相、反相和呼吸模态。
找一个 $2\times 2$ 反例，说明 $AB\ne BA$ 时 $e^{A+B}\ne e^A e^B$ 。

编程

写个相图绘制器，输入任意 $2\times 2$ 矩阵，自动标注平衡点类型（迹-行列式判别法）。
模拟耦合振子系统，画拍频周期 $T_{\text{beat}} = 2\pi/(\omega_a - \omega_s)$ 随 $\kappa$ 变化的曲线，对比小 $\kappa$ 近似。

一个具体算例：上三角矩阵的特征值与模态分解#

A = \begin{pmatrix} 3 & 1 \\ 0 & 2 \end{pmatrix}.

特征多项式 $\det(A - \lambda I) = (3 - \lambda)(2 - \lambda) - 0 = 0$ ，所以 $\lambda_1 = 3$ 、 $\lambda_2 = 2$ 。

特征向量。对 $\lambda_1 = 3$ ：解 $(A - 3I)\mathbf{v} = 0$ ，即 $\begin{pmatrix} 0 & 1 \\ 0 & -1 \end{pmatrix}\mathbf{v} = 0$ ，得 $\mathbf{v}_1 = (1, 0)^\top$ 。对 $\lambda_2 = 2$ ： $\begin{pmatrix} 1 & 1 \\ 0 & 0 \end{pmatrix}\mathbf{v} = 0$ ，得 $\mathbf{v}_2 = (1, -1)^\top$ 。

x_1(t) = 3 e^{3t} + 2 e^{2t}, \qquad x_2(t) = -2 e^{2t}.

e^{At} = P \begin{pmatrix} e^{3t} & 0 \\ 0 & e^{2t} \end{pmatrix} P^{-1} = \begin{pmatrix} e^{3t} & e^{3t} - e^{2t} \\ 0 & e^{2t} \end{pmatrix}.

取 $$t = 1$$ ： $e^{3} \approx 20.09$ ， $e^{2} \approx 7.39$ ，所以 $e^{A} \approx \begin{pmatrix} 20.09 & 12.70 \\ 0 & 7.39 \end{pmatrix}$ 。代入 $\mathbf{x}(0) = (5, -2)^\top$ 得 $\mathbf{x}(1) \approx (5 \cdot 20.09 - 2 \cdot 12.70,\ -2 \cdot 7.39) = (74.95,\ -14.78)$ ，与 $3 e^{3} + 2 e^{2} = 60.27 + 14.78 = 75.05$ 略有偏差——是我手算 $$e^3, e^2$$ 时只保留两位小数引入的误差，scipy 算出来是 $$74.97$$ ，两边匹配到三位。

\mathbf{x}(t) = \underbrace{3 e^{3t}\binom{1}{0}}_{\text{模态 1（增长更快）}} + \underbrace{2 e^{2t}\binom{1}{-1}}_{\text{模态 2}}.

$t \to \infty$ 时模态 1 完全主导， $\mathbf{x}(t)/\|\mathbf{x}(t)\|$ 收敛到 $\mathbf{v}_1 = (1, 0)^\top$ 的方向——这就是幂法（power method）的本质：任何随机初值在 $$A$$ 反复作用下都会被最大模特征值的特征方向"吸"过去。Google PageRank 的早期算法核心就是这个事实。

直觉与陷阱：上三角不是"特殊"，特征值也不"等于对角元"#

很多人记规则记成"上三角矩阵的特征值就是对角元"。这话对实矩阵成立（任何上三角实矩阵 $\det(A - \lambda I)$ 展开就是对角元乘积），但容易让人误以为这是某种偶然的简单情形。其实这反映了一个深层事实：任何复方阵都酉相似于上三角矩阵——这就是 Schur 分解 $$A = U T U^*$$ ， $$T$$ 是上三角， $$U$$ 是酉的。Schur 分解保证了求特征值在数值上等价于求上三角矩阵的对角元，这是 LAPACK dgeev 内部 QR 算法收敛后做的事。

陷阱一：亏损情形。如果上面那个矩阵改成 $\begin{pmatrix} 3 & 1 \\ 0 & 3 \end{pmatrix}$ ，特征值都是 3，但 $(A - 3I) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}$ 只有一个特征向量 $(1, 0)^\top$ 。这时通解必须包含 $t e^{3t}$ 项——多项式增长是 Jordan 块的签名。许多偏微分方程数值离散后正好得到这种亏损矩阵（对流方程的迎风差分），如果你不加广义特征向量，初值在亏损方向上会被忽略。

陷阱二：条件数。如果 $$P$$ 接近奇异（特征向量近乎平行）， $P^{-1}$ 范数巨大， $e^{At} = P D P^{-1}$ 的浮点误差会爆炸。Moler 和 Van Loan 在 Nineteen Dubious Ways 那篇综述里专门讨论过：高度非正常矩阵（normal 矩阵指 $$A A^* = A^* A$$ ，对称、酉都属于此类）下，特征值法不可信，必须用 scaling-and-squaring。

应用与反例：模态分解在工业中的真实分量#

模态分解不是教科书装饰，而是实打实的工程语言：

结构动力学：高层建筑的有限元模型动辄 $$10^5$$ 自由度，但前 20 个模态贡献了 90% 的位移能量。台北 101 防风方案就是先做模态分析，找到第一阶横向模态频率约 0.15 Hz，调谐质量阻尼器精确匹配这个频率把振幅压低 40%。
量子化学：分子振动光谱本质上是简正模态——把势能在平衡构型展开到二阶，对力常数矩阵做特征分解，得到的 $\omega_i$ 就是红外吸收峰的频率。
机器学习：PCA 是协方差矩阵的特征分解；图卷积网络（GCN）用图拉普拉斯的特征向量做谱滤波；DMD（动态模态分解）用线性算子的特征值预测流体流场。

但模态分解不是万能。有两类系统它会出问题：

第一，强非线性。线性模态在 Duffing 振子 $\ddot x + x + \alpha x^3 = 0$ 大振幅下完全失效——非线性把模态混合（modal coupling），频率随振幅变化。这时要用非线性正规模（Shaw-Pierre nonlinear normal modes）或 Koopman 算子谱。

第二，非正常算子。 $$A$$ 的特征值都在左半平面不代表 $\|e^{At}\|$ 单调下降。瞬态增长（transient growth）现象在剪切流稳定性、神经动力学、生态网络里都很关键——Trefethen 的伪谱理论（pseudospectra）就是为这类问题设计的。Couette 流的特征值全负，但小扰动可以放大 $$10^4$$ 倍再衰减，这是流体湍流亚临界转捩的核心机制。光看特征值会错过这一切。

下一步#

下一章引入稳定性理论——研究’扰动一下，系统会回到平衡点还是飞走’。这是工程师最关心的问题（控制系统、机械振动、化工反应器都需要稳定性分析）。我们会从线性化局部分析开始，然后用 Lyapunov 函数处理非线性情况。

工程实践片段：耦合振子在工程中的真实身影#

线性系统组的理论看似抽象，落到具体工程里非常具象。

例一：建筑抗震设计。一栋十层楼简化成十个质量块用刚度矩阵 $$K$$ 连接，方程组就是 $M\ddot{\mathbf{x}} + C\dot{\mathbf{x}} + K\mathbf{x} = 0$ 。求 $M^{-1}K$ 的特征值得到固有频率，最低几阶通常对应整体摆动模态。设计师必须避开 $0.5\sim 5$ Hz 这个地震主能量带——日本东京大学的研究表明，台风和地震的能量谱在低频区差异巨大，所以南方抗台风楼和北方抗震楼的最优固有频率完全不同。这是用矩阵特征值"算"出来的硬约束。

例二：电子滤波器级联。多级 RLC 滤波器的状态变量耦合形式正好是线性方程组 $\dot{\mathbf{x}} = A\mathbf{x} + B u$ 。设计带通滤波器时，工程师手里就拿着 $$A$$ 矩阵的特征值在复平面上拖动——这就是模拟电路圈子里说的"极点位置决定一切"。Sallen-Key 拓扑、状态变量滤波器、双二阶节，每一种都是把六阶传递函数分解成多个二阶节的特征对，理论上对应分块对角化。

例三：电力系统稳定性。北美东部联网的电网有上千台同步发电机，每台用 $\delta_i$ （功角）和 $\omega_i$ （转速偏差）描述，整体动力学是一个数千维线性系统。慢相干模态（0.1~1 Hz）对应区域间的电力波动，工程师通过把 $$A$$ 矩阵的特征值往左半平面推（加阻尼）来保稳定。1996 年北美西部大停电的根因就是某个区域间模态阻尼变成负值——纯数学上的"特征值跨过虚轴"，物理上就是几千万人停电。

这三个场景共享同一份数学：线性 ODE 组的特征值和特征向量决定了能不能稳、稳得多快、振哪些频率。理解这一点之后，看似杂乱的工程领域立刻被一根线串起来。这也是我反复推荐学生认真学习 $e^{At}$ 的原因——它不是一个数学符号，是物理世界的"流"。

参考文献#

Hirsch, Smale, & Devaney, Differential Equations, Dynamical Systems, and an Introduction to Chaos, Academic Press (2012)
Strang, Linear Algebra and Learning from Data, Wellesley-Cambridge Press (2019)
Moler & Van Loan, “Nineteen Dubious Ways to Compute the Exponential of a Matrix,” SIAM Review 45(1) (2003)
Perko, Differential Equations and Dynamical Systems, Springer (2001)