常微分方程（七）：稳定性理论

轻轻推一下系统，它是会恢复原状、逐渐漂离，还是会彻底崩溃？ 这个问题的答案，决定了桥梁能否扛过风暴、生态系统能否从干旱中恢复，以及经济能否走出危机。稳定性理论正是为此而生——而且它根本不需要解出微分方程。我们将学会直接从相平面的几何结构中，读出系统的命运。

总结#

• 三个精确定义：Lyapunov 稳定、渐近稳定、不稳定
• 利用 Jacobian 矩阵进行线性化，以及 Hartman-Grobman 定理
• Lyapunov 直接方法：通过构造类能量函数来证明稳定性
• LaSalle 不变集原理：处理边界情况的利器
• 迹-行列式平面：对所有二维线性系统进行统一分类
• 四种典型分岔：鞍-结点、跨临界、叉型、Hopf
• 实际应用：单摆、捕食者-猎物模型、倒立摆控制

前置知识#

• 第六章：线性系统、特征值、相图
• 多元微积分：偏导数、Jacobian 矩阵

先看图：六张图讲完线性系统#

稳定性本质上是一个几何问题，描述的是轨迹在相空间中如何运动。仅用六张图，就能完整刻画二维线性系统的所有行为。

对于非线性系统，这六种图形依然会出现——只不过仅在每个平衡点的局部邻域内成立。阻尼单摆和 Lotka-Volterra 捕食者-猎物模型就是绝佳例证：

稳定性的精确定义#

考虑系统 $\mathbf{x}' = \mathbf{f}(\mathbf{x})$ ，其平衡点为 $\mathbf{x}^*$ （即 $\mathbf{f}(\mathbf{x}^*) = \mathbf{0}$ ）。

Lyapunov 稳定： 对任意 $\varepsilon > 0$ ，都存在 $\delta > 0$ ，使得只要初始状态满足 $\|\mathbf{x}(0) - \mathbf{x}^*\| < \delta$ ，就有 $\|\mathbf{x}(t) - \mathbf{x}^*\| < \varepsilon$ 对所有 $t > 0$ 成立。
直觉：附近的轨迹永远不会跑远。

渐近稳定： 不仅是 Lyapunov 稳定，还要求 $\mathbf{x}(t) \to \mathbf{x}^*$ （当 $t \to \infty$ ）。
直觉：附近的轨迹不仅不跑远，最终还会回到平衡点。

不稳定： 不满足 Lyapunov 稳定的条件。

吸引域是指所有能收敛到 $\mathbf{x}^*$ 的初始条件构成的集合。渐近稳定性是一个局部性质，而吸引域的大小则告诉你这个“局部”究竟有多广。

为什么要区分两种稳定？ 中心（即闭合轨道）是 Lyapunov 稳定的，但不是渐近稳定的——轨迹虽然始终在附近打转，却永远不会真正停下来。Lotka-Volterra 模型就是最经典的例子。

线性化：Jacobian 方法#

Hartman-Grobman 定理#

如果 Jacobian 矩阵 $J$ 的所有特征值的实部都不为零（这种平衡点称为双曲的），那么非线性系统在 $\mathbf{x}^*$ 附近与它的线性化系统 $\mathbf{u}' = J\mathbf{u}$ 是拓扑等价的。

• 所有特征值实部 $< 0$ ：渐近稳定
• 存在特征值实部 $> 0$ ：不稳定
• 特征值为纯虚数：线性化方法失效，必须使用 Lyapunov 方法

实例：阻尼单摆#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

gamma, omega0 = 0.3, 1.0
def pendulum(X, t):
    x, y = X
    return [y, -gamma*y - omega0**2*np.sin(x)]

fig, ax = plt.subplots(figsize=(10, 6))
t = np.linspace(0, 50, 2000)
for x0 in np.linspace(-3*np.pi, 3*np.pi, 15):
    for y0 in np.linspace(-3, 3, 5):
        sol = odeint(pendulum, [x0, y0], t)
        ax.plot(sol[:,0], sol[:,1], 'b-', linewidth=0.3, alpha=0.5)
ax.set_title('阻尼单摆相图')
ax.grid(True, alpha=0.3); plt.tight_layout(); plt.show()

Lyapunov 直接方法#

核心思想#

无需解出 ODE，也能判断稳定性。只需构造一个类能量的标量函数 $V(\mathbf{x})$ ，并证明它沿着系统的轨迹是递减的。

基本要求：

1. $V(\mathbf{x}^*) = 0$ ，且在其他地方 $V(\mathbf{x}) > 0$ （正定）
2. 沿轨迹有 $\dot V = \nabla V \cdot \mathbf{f}(\mathbf{x}) \leq 0$

稳定性判定定理：

几何直观#

轨迹会从外向内穿越 $V$ 的等高线。由于 $V$ 在 $\mathbf{x}^*$ 处取得最小值，轨迹无法逃出足够小的等高线区域；如果 $V$ 严格递减，轨迹就会一路滑到底部。

如何寻找 $V$ #

• 物理能量：对力学系统，直接用动能加势能
• 二次型：设 $V = \mathbf{x}^T P \mathbf{x}$ ，其中 $P$ 满足 Lyapunov 方程 $A^T P + PA = -Q$
• 试探法：从简单的 $V = x^2 + y^2$ 开始，计算 $\dot V$ ，再调整系数

实例：单摆的能量#

这说明下垂位置是稳定的。而 LaSalle 原理（见下节）可以将这一结论进一步提升为渐近稳定。

LaSalle 不变集原理#

有时 $\dot V \leq 0$ ，但它在一个完整的集合上恒为零，而不仅仅在平衡点处为零。此时标准的 Lyapunov 方法只能保证稳定性，无法得出收敛性。

定理： 令 $E = \{\mathbf{x} : \dot V = 0\}$ ，并设 $M$ 是 $E$ 中最大的不变子集。那么，每一条有界的轨迹最终都会趋近于 $M$ 。

特别地，如果 $M = \{\mathbf{x}^*\}$ ，那么 $\mathbf{x}^*$ 就是渐近稳定的。

以阻尼单摆为例，$\dot V = 0$ 要求 $\omega = 0$ 。但在直线 $\omega = 0$ 上，动力学方程给出 $\dot\omega = -\omega_0^2 \sin\theta$ ，除非同时有 $\theta = 0$ ，否则该项不为零。因此，唯一的不变集就是 $M = \{(0,0)\}$ ，于是我们立刻得到渐近稳定性。

迹-行列式分类#

对于二维线性系统 $\mathbf{x}' = A\mathbf{x}$ ，记 $\tau = \operatorname{tr}(A)$ ，$\Delta = \det(A)$ 。其特征值为 $\lambda_{1,2} = \tfrac12(\tau \pm \sqrt{\tau^2 - 4\Delta})$ ，由此可得如下分类：

一张图就足以概括整个表格：

分岔：当图形本身发生变化#

缓慢调节一个参数 $r$ ，系统的平衡点可能会凭空产生、相互湮灭，或者交换稳定性。四种“标准形式”足以局部刻画所有余维为 1 的分岔现象。

Hopf 分岔是自然界中所有自持振荡现象的核心机制——从心跳节律到变星的脉动，皆源于此。

应用 1：Lotka-Volterra 捕食者-猎物模型#

它保证了所有轨道都是闭合的。因此，该模型呈现出周期性的种群数量波动（见图 2 右侧）。

应用 2：倒立摆控制#

倒立平衡点本身是一个鞍点。通过施加线性反馈控制 $u = -K\mathbf{x}$ ，可以将闭环系统的特征值全部移到复平面的左半部分，从而将原来的鞍点转变为一个稳定焦点。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

import numpy as np
from scipy.linalg import solve_continuous_are
from scipy.integrate import odeint

g, l, m, M_cart = 9.81, 0.5, 0.1, 1.0
A = np.array([[0,1,0,0],[0,0,-m*g/M_cart,0],
              [0,0,0,1],[0,0,(M_cart+m)*g/(M_cart*l),0]])
B = np.array([[0],[1/M_cart],[0],[-1/(M_cart*l)]])
Q = np.diag([10, 1, 100, 1]); R = np.array([[0.1]])
P = solve_continuous_are(A, B, Q, R)
K = np.linalg.inv(R) @ B.T @ P
print("闭环特征值：", np.linalg.eigvals(A - B @ K))

练习题#

基础题。

1. 判断下列系统的稳定性：(a) $x' = -x,\ y' = -2y$ ；(b) $x' = y,\ y' = -\sin x - 0.5y$ 。
2. 使用 $V = x^2 + y^2$ 分析系统 $x' = -x + y^2,\ y' = -y$ 。
3. 找出方程 $\dot x = rx - x^3$ 的所有分岔点。

进阶题。

4. 证明总能量是阻尼单摆的一个 Lyapunov 函数，并应用 LaSalle 原理推导其渐近稳定性。
5. 分析 Van der Pol 振子：证明原点不稳定，但存在一个稳定的极限环。
6. 对 Lotka-Volterra 竞争模型，推导出物种共存与竞争排斥的条件。

编程题。

7. 编写程序，根据迹和行列式自动分类二维线性系统的平衡点，并复现图 3。
8. 制作 Hopf 分岔的动画：扫描参数 $\mu$ ，观察极限环如何从无到有并逐渐增大。

一个具体算例：六种相图的迹-行列式判别表#

文中给出了迹-行列式分类原则，但对刚学的人来说，"$\Delta > 0,\ \tau < 0,\ \tau^2 > 4\Delta$ " 这种条件太抽象。我把六种典型相图配上具体矩阵、$\tau$ 、$\Delta$ 、特征值，做成一张可以背的表：

读这张表的方法：

• 第一关看 $\Delta$ 的符号。$\Delta < 0$ 立刻锁定鞍点，不必看其他。
• 第二关看 $\tau^2 - 4\Delta$ 。正号且 $\Delta > 0$ 是结点（实根同号），负号是螺旋或中心（复根），零是退化或星形。
• 第三关看 $\tau$ 符号决定稳定性（负稳定，正不稳定，零是中心）。

手算验证一例。取退化结点 $A = \begin{pmatrix} -2 & 1 \\ 0 & -2 \end{pmatrix}$ 。$\det(A - \lambda I) = (-2 - \lambda)^2 = 0$ ，$\lambda = -2$ （重根）。解 $(A + 2I)\mathbf{v} = 0$ ：$\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}\mathbf{v} = 0$ ，只得 $\mathbf{v} = (1, 0)^\top$ 。亏损一个特征向量，需要广义特征向量 $\mathbf{w}$ ：解 $(A + 2I)\mathbf{w} = \mathbf{v}$ ，即 $\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}\mathbf{w} = \binom{1}{0}$ ，得 $\mathbf{w} = (0, 1)^\top$ 。通解：$\mathbf{x}(t) = c_1 e^{-2t}(1, 0)^\top + c_2 e^{-2t}\bigl(t (1, 0)^\top + (0, 1)^\top\bigr)$ 。

注意 $t e^{-2t}$ 这一项：它先增长（$t < 1/2$ ）再衰减（$t > 1/2$ ）。所以退化节点的轨迹会先远离原点再回来，这是教科书图里画的"弯钩"形状的来源。这种"先冲一下再回收"是亏损 Jordan 块的标志，跟稳定结点的"直线进圆点"是肉眼可辨的区别。

直觉与陷阱：中心是"骗人"的#

迹-行列式平面里，中心对应一条线 $\tau = 0,\ \Delta > 0$ ——是测度为零的集合。在物理实验里你不会真的见到一个完美的中心：地球大气、电路板上的电感、肌肉里的兴奋纤维，没有任何一个能让 $\tau$ 精确等于 0。这意味着——

陷阱一：理论中心实验上不存在。任何阻尼、任何能量耗散、任何模型不确定性都会把 $\tau$ 推离 0，瞬间把中心变成稳定螺旋（$\tau < 0$ ）或不稳定螺旋（$\tau > 0$ ）。Lotka-Volterra 模型给的中心，加任何小扰动就破坏闭合轨道——这正是它在生态学里被批评"过于简化"的根本原因。

陷阱二：纯虚特征值时 Hartman-Grobman 失效。线性化告诉你"中心"，但完整非线性方程可能给出稳定螺旋、不稳定螺旋、甚至极限环。教科书会用 $\dot r = a r^3$ （在极坐标下）这种例子说明：线性部分是 $\pm i$ ，但 $a > 0$ 的三次项让原点不稳定，$a < 0$ 让它稳定。只看特征值你会判错——这是 Lyapunov 方法不可替代的根本理由。

陷阱三：数值方法在中心上"伪稳定"或"伪不稳定"。Euler 方法每一步注入正能量，让中心看起来不稳定；隐式 Euler 注入负能量，看起来稳定。连两步 Adams-Bashforth 和 RK4 都有微小的能量漂移。要正确数值模拟一个真正的中心（比如 Hamilton 系统），必须用辛积分器（symplectic integrator），第十一章和第十七章会反复强调这件事。

应用与反例：分岔图在工业控制里的真实价值#

迹-行列式平面看似抽象，对工程师却是一张"可达性地图"——它告诉你通过调节参数可以把系统推向哪种行为。具体例子：

• 核反应堆点火。中子通量 $n$ 与延迟中子前驱体浓度 $C$ 满足 $\dot n = (\rho - \beta) n / \Lambda + \lambda C$ ，$\dot C = \beta n / \Lambda - \lambda C$ ，其中 $\rho$ 是反应性。$\rho < 0$ 是稳定结点（反应缓慢衰减），$\rho = 0$ 是退化稳定（临界），$\rho > 0$ 是不稳定结点（功率指数发散——切尔诺贝利的根源）。控制棒就是在 $\rho$ 上手动调节。
• 激光腔阈值。腔内光子数 $N$ 与粒子反转 $\Delta n$ 的方程组在泵浦功率小于阈值时是稳定结点（光子衰减为零），高于阈值后通过跨临界分岔变成稳定螺旋（激光振荡建立）。这是激光器设计中"阈值"概念的数学本质。
• 生态系统崩塌。三种群食物链模型在某个环境参数（资源、温度、捕捞率）变化下，可以从稳定结点经过 Hopf 分岔进入极限环，再经过倍周期级联进入混沌——这正是渔业管理里"过度捕捞导致系统失稳"的数学描述。

反例：迹-行列式分类只对二维线性系统普适。三维及以上：你需要看完整谱（多个实部交叠），鞍点不再只是"两稳两不稳"，可以是"鞍-焦点"（Shilnikov 同宿轨道是混沌的发源地）。非线性系统：分类只在双曲平衡点附近局部成立，全局可能有多个不动点、极限环、奇异吸引子并存——这就是为什么 Lorenz 系统三个不动点都不稳定，但轨迹却被局限在一个有限的混沌吸引子上。这种全局-局部分离正是第八章和第九章的主题。

数值实现：用 SciPy 解 Lyapunov 方程#

在纸上，方程 $A^\top P + PA = -Q$ 看起来像是一个可以直接求解的矩阵方程。但手算最多只能应付二维情形。维度更高时，就得求助于 scipy.linalg.solve_continuous_lyapunov：

1
2
3
4
5
6
7
8

import numpy as np
from scipy.linalg import solve_continuous_lyapunov, eigvals

A = np.array([[0, 1], [-2, -3]])      # 阻尼振子
Q = np.eye(2)
P = solve_continuous_lyapunov(A, -Q)   # 符号约定：A P + P A^T = -Q
print("P =", P)
print("P 正定？", np.all(eigvals(P).real > 0))

其底层采用的是 Bartels-Stewart 算法：先对 $A$ 进行 Schur 分解（$A = U T U^\top$ ），再将未知量变换为 $\tilde P = U^\top P U$ ，从而将原问题转化为一个上三角线性方程组，最后按列回代求解。该算法的时间复杂度为 $O(n^3)$ ，比直接将方程展开为 $n^2 \times n^2$ 的 Kronecker 形式快了整整 $n^3$ 倍。

常见陷阱： 如果矩阵 $A$ 有特征值落在虚轴上，Lyapunov 方程就会变成奇异的，此时 SciPy 会抛出 LinAlgError 异常。这个数值上的“报错”，恰恰对应着理论上的非双曲平衡点——在那里，Hartman-Grobman 定理失效。同一个现象，只是用了两种不同的语言来描述。

当线性化失效：中心流形与高阶项#

Hartman-Grobman 定理要求 Jacobian 的所有特征值实部非零。一旦有特征值的实部为零，线性部分就只能告诉你中性方向的信息，而真正的稳定性则由非线性项决定。

一个简单例子：$\dot x = -x^3$ 。在 $x=0$ 处，Jacobian 为 0，线性分析完全沉默。但若取 $V(x) = x^2/2$ ，则 $\dot V = -x^4 \le 0$ ，由此可证系统是稳定的。这正是 Lyapunov 方法能看见而特征值分析所遗漏的情形。

中心流形定理： 假设系统在原点有 $n_s$ 个稳定方向和 $n_c$ 个中心（中性）方向，那么存在一个 $n_c$ 维的不变流形 $W^c$ （称为中心流形），使得原点附近的轨迹都会指数级快速地坍缩到 $W^c$ 上。于是，整个系统的稳定性问题就简化为在 $W^c$ 上分析降维后的动力学。下一章要讲的 Hopf 分岔标准型，正是通过这种方法推导出来的。

实用建议： 在编程时，如果 np.linalg.eigvals(J) 返回的某个特征值满足 abs(real) < 1e-8，就不要仅凭特征值谱来断言稳定性。此时，要么尝试构造一个 Lyapunov 函数，要么将泰勒展开推进到更高阶。

接入 ML：Lyapunov 神经网络与稳定性约束训练#

我在实践中总结出两条有效技巧：

1. 将正定性直接嵌入网络架构： 令 $V_\theta(x) = \|\phi_\theta(x)\|^2 + \epsilon \|x\|^2$ ，其中 $\phi_\theta$ 是一个标准的 MLP。这样，$V_\theta \ge \epsilon\|x\|^2$ 就自然成立，无需额外添加惩罚项。
2. 用 SMT 求解器验证下降条件。 基于采样的损失函数不能构成严格证明。训练完成后，可以把 $V_\theta$ 和 $f$ 交给 dReal 或 Z3 这样的求解器，让它寻找反例，再据此重新训练——这就是所谓的 CEGIS（Counterexample-Guided Inductive Synthesis）循环。

值得关注的论文包括：Neural Lyapunov Control (Chang et al., 2019) 和 Learning Stability Certificates from Data (Boffi et al., 2020)。这些方法已经成功部署在腿式机器人的控制器中。

总结#

下一步#

下一章正面进入非线性世界。线性 ODE 是干净的（叠加原理、特征值分析），但真实世界很少线性。捕食者-猎物、神经元放电、化学振荡——它们的关键特征（极限环、混沌、分岔）都不会出现在线性方程里。下一章用相图和定性分析处理这些。

工程实践片段：稳定性分析的几个真实场景#

把稳定性理论从教科书搬到工程里，能讲很多有教益的故事。

故事一：直立摩托车的动态稳定性。骑过摩托车的人都知道：低速时车很难直立（鞍点），高速时反而稳如泰山（稳定螺旋）。这不是错觉。陀螺效应让车体绕垂直轴的角动量在转向时产生回正力矩，是速度的函数。把摩托车-车手系统线性化，特征值的实部在某个临界速度（约 20 km/h）从正变负——这是工程上"自稳定速度"的概念。设计师调节前轮拖距和叉杆角度，本质上就是把这个临界速度往下推，让低速更安全。

故事二：化工反应器的多稳态。放热化学反应器在某些参数下有三个稳态：低温（反应慢、产物少）、高温（反应快、危险）、中温（不稳定）。线性化告诉你前两个是稳定结点、中间是鞍点。但工厂启动时怎么从低温到高温？这就是 Lyapunov 函数派上用场的地方——计算吸引域的边界（鞍点的稳定流形），知道从哪个初始温度出发能"翻越鞍点"达到目标稳态。八十年代美国 Bhopal 化工事故部分就是因为操作员没意识到反应器已被推过吸引域边界，进入失控区。

故事三：电网的暂态稳定性。短路故障后线路重合闸时，发电机功角可能已经偏离平衡点很远。能量函数法（Lyapunov 直接法的电力系统版本）给出"临界清除时间"——故障必须在多少毫秒内切除才能保证系统重新进入吸引域。继电保护装置的整定就是按这条边界设置的。

这三个故事共享同一个抽象：稳定性是局部概念，吸引域才是全局答案。Hartman-Grobman 告诉你不动点附近的行为，但 Lyapunov 函数能告诉你从多远的地方还能回得来。这是我反复强调给学生的核心区别——会算特征值不等于懂稳定性。

故事四：神经网络训练的隐含稳定性。深度学习里的"梯度爆炸"和"梯度消失"，本质上是反向传播图的雅可比矩阵特征值的累积乘积。LSTM 的精巧设计就是把这个雅可比的最大特征值固定在 1 附近——保持中性稳定。Transformer 的 LayerNorm 在做类似的事。这是稳定性理论被搬进 ML 的一个具体例子。

参考文献#

• Strogatz, Nonlinear Dynamics and Chaos, CRC Press (2015)
• Khalil, Nonlinear Systems, Prentice Hall (2002)
• Guckenheimer & Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations, Springer (1983)
• Perko, Differential Equations and Dynamical Systems, Springer (2001)