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

一个方程描述一个量。但世界很少这么配合。 兔群与狼群此消彼长，RLC 网络中的电流和电压互相牵动，化学反应里的物质浓度彼此影响。只要两个未知量出现在同一组方程里，你就有了一个方程组，标量公式 $y'=ay$ 已经不够用了。

线性情形的奇迹在于：标量解 $y(t)=e^{at}y_0$ 一字不改地推广到方程组——只要你弄清楚什么是矩阵的指数 $e^{At}$。线性代数和微分方程在这里融为一体，而矩阵 $A$ 的特征结构告诉你解的长期行为、流的几何形状，以及简正模态、拍频这些物理图像背后的代数。

本章要点

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

前置知识

• 线性代数：特征值、特征向量、基的变换
• 第三章：高阶线性 ODE（任何这样的方程都可以化成二维方程组）

1. 从两个耦合方程到一个向量方程

考虑一个故意简化的生态模型：$x(t)$ 是兔群的（无量纲）数量，$y(t)$ 是狼群。兔子自身繁殖、被狼吃掉；狼只在有兔子时增长：

把未知量摞成一个向量 $\mathbf{x}=(x,y)^\top$，把系数摞成一个矩阵：

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

标量方程 $y'=ay$ 的解是 $y=e^{at}y_0$。把这个类比硬推一下：

但这只有在我们说清楚矩阵的指数是什么之后才有意义。它正是本章的中心对象。

2. 矩阵指数

2.1 用幂级数定义

模仿 $e^x = \sum x^k/k!$：

这个级数对任何方阵 $A$ 和任何 $t\in\mathbb{R}$ 都按算子范数收敛，因为 $\|(At)^k/k!\| \le (\|A\|t)^k/k!$，而后者来自一个收敛的标量级数。

2.2 真正常用的三条性质

把级数写出来之后，几乎所有论证都只用到三件事：

一个微妙的警告：$e^{A+B}=e^A e^B$ 只有在 $AB=BA$ 时才成立。这是矩阵版本的 $\sin(x+y)\ne\sin x+\sin y$——非交换性是有后果的。

2.3 实际计算的三种方式

1. 直接截断幂级数。 概念上便宜，当 $\|At\|$ 大时数值上很差。
2. 特征值分解。 若 $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)    # 同一矩阵（差一点浮点噪声）

3. 特征值方法

特征值分解的公式有一个完全等价的、不用矩阵指数的写法。

定理。 若 $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)$。线性性把剩下的活做完。

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

$A=PDP^{-1}$ 这个分解有一个非常直观的几何读法：把 $A$ 作用到一个向量上，相当于

1. $P^{-1}$：把向量改写到特征基；
2. $D$：沿每条特征方向按对应的 $\lambda_i$ 拉伸；
3. $P$：把拉伸后的向量再写回原来的基。

流 $e^{At}$ 做一模一样的事，只不过把拉伸因子换成 $e^{\lambda_i t}$。

3.2 复特征值给出旋转

实矩阵也可以有复特征值，且必然成共轭对 $\lambda=\alpha\pm\beta i$，特征向量也共轭 $\mathbf{v}=\mathbf{a}\pm i\mathbf{b}$。取 $e^{\lambda t}\mathbf{v}$ 的实部和虚部就得到两个实解：

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

4. 二维相图

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

4.1 迹–行列式速查

很多时候你不必算出特征值。$2\times 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$ 是重特征值的轨迹。

5. 重特征值与广义特征向量

上面那张分类表悄悄藏了一个麻烦。代数重数为 2 的特征值 $\lambda$ 可能有两个线性无关的特征向量（罕见——只有当 $A$ 在该子空间上等于 $\lambda I$ 时才会发生，对应星形节点），也可能只有一个（通有的亏损情形——退化节点）。

亏损时，特征值方法只给出一个解 $e^{\lambda t}\mathbf{v}$。补救办法：解

得到一个广义特征向量 $\mathbf{w}$。然后

是另一个线性无关的解。$t$ 这个因子带来的"多项式乘指数"增长，正是亏损在代数上的指纹。

这其实就是 Jordan 块现象：$A$ 在该子空间上看像 $\lambda I + N$，$N$ 幂零，因此

由于 $N$ 幂零，级数有限项就停了，留下一个 $t$ 的多项式。

6. 非齐次方程组：Duhamel 公式

加上一个外力项：

解就是参数变动法的矩阵版：

这就是 Duhamel 公式。被积函数 $e^{A(t-\tau)}\mathbf{g}(\tau)$ 应当这样读：在 $\tau$ 时刻，外力给系统一个 $\mathbf{g}(\tau)d\tau$ 的"踢"，这个踢之后在剩下的时间里按自由流 $e^{A(t-\tau)}$ 演化。整体的解就是所有这些延迟响应的叠加——一个连续时间的脉冲响应。

控制论里同一个公式以 $\mathbf{g}(t) = B\mathbf u(t)$ 的形式出现：

离可控性和可控性矩阵 $\bigl[B,\,AB,\,A^2B,\dots\bigr]$ 就只剩一步之遥了。

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

直线上的两个等质量小球，每个用劲度 $k$ 的弹簧连到墙上，相互之间用劲度 $\kappa$ 的弹簧连接：

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

牛顿定律给出

漂亮的换元是 $u=\tfrac{x_1+x_2}{\sqrt 2}$、$v=\tfrac{x_1-x_2}{\sqrt 2}$。方程解耦了：

这就是简正模态：一个同相模态，频率 $\omega_s=\sqrt{k}$（耦合弹簧从来不被拉伸），一个反相模态，频率 $\omega_a=\sqrt{k+2\kappa}$（耦合弹簧也在做功）。任何运动都是这两个模态的叠加。

戏剧性的事情发生在两个模态频率几乎相等时（$\kappa \ll k$）：只拨动 1 号小球，相当于两个模态以等幅被同时激发，而它们之间缓慢的相位差让能量整个地从 1 号转移到 2 号、再转回来。这就是拍频——两支微调失谐的音叉一同发声时听到的强弱起伏。

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

8. 稳定性：特征值判据

对线性流 $\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 矩阵的图景。

总结

练习题

基础

1. 用幂级数和特征值分解两种方式分别计算 $A=\bigl(\begin{smallmatrix}0&1\\-1&0\end{smallmatrix}\bigr)$ 的 $e^{At}$，验证两者相同。
2. 求解 $\mathbf{x}'=\bigl(\begin{smallmatrix}1&2\\0&3\end{smallmatrix}\bigr)\mathbf{x}$，$\mathbf{x}(0)=(1,1)^\top$。
3. 判断平衡点类型：(a) $A=\bigl(\begin{smallmatrix}-1&2\\-2&-1\end{smallmatrix}\bigr)$；(b) $A=\bigl(\begin{smallmatrix}1&0\\0&-2\end{smallmatrix}\bigr)$。

进阶

4. 证明 $\det e^A = e^{\mathrm{tr}\,A}$。提示：在 $\mathbb C$ 上做 Schur 上三角化。
5. 求三个等质量小球（一字排开，相邻之间和两端到墙之间都用劲度 $k$ 的弹簧连接）的简正模态频率，识别同相、反相和"呼吸"三种模态。
6. 给一个 $2\times 2$ 反例，说明 $AB\ne BA$ 时 $e^{A+B}\ne e^A e^B$。

编程

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

参考资料

• 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)

系列导航