常微分方程（八）：非线性系统与相图

真实世界是非线性的： 捕食者-猎物的周期性振荡、心跳节律、神经元放电——这些现象都无法用线性方程刻画。一旦叠加原理失效，系统便展现出全新的行为：极限环、多个平衡点、双稳态、滞回效应。本章将为你提供一套几何与解析工具，让你能直接从二维相图中“读出”这些复杂动态。

总结#

• 非线性系统与线性系统存在本质差异
• Lyapunov 稳定性的直观理解：等值面、势能“碗”与吸引域
• 线性化 vs. 完整非线性图像（Hartman-Grobman 定理的实际应用）
• Lotka-Volterra 捕食-猎物模型：闭合轨道与守恒量
• 种间竞争模型：四种典型结局
• Van der Pol 振子与极限环的几何结构
• 梯度系统与 Hamilton 系统
• Poincaré-Bendixson 定理：为何二维系统不可能出现混沌

前置知识#

• 第六章：线性系统、相图分类
• 第七章：稳定性、线性化、Lyapunov 函数

从线性到非线性#

线性系统满足叠加原理：若 $\mathbf{x}_1$ 和 $\mathbf{x}_2$ 是解，则任意线性组合 $c_1\mathbf{x}_1 + c_2\mathbf{x}_2$ 也是解。这一性质是第 1–6 章所有工具（如指数试探解、特征向量、基本矩阵）的理论基石。

非线性系统打破了这一规则，代价是通常无法获得解析解；但换来的是无价之宝：

• 多个平衡点，每个都可能具有不同的稳定性类型
• 极限环——孤立且稳定的周期轨道（在线性系统中不可能存在）
• 双稳态与滞回——系统对初始条件具有“记忆”能力
• 敏感依赖性——在三维及以上系统中可导致混沌（见第九章）

物理、生物、化学、神经科学和经济学中几乎所有有趣的系统本质上都是非线性的。

Lyapunov 稳定性的直观理解#

Lyapunov 函数 $V(\mathbf{x})$ 是一个沿轨迹单调不增的标量函数（即 $\dot V \leq 0$ ）。从几何角度看，$V$ 的等值面构成一组嵌套的“碗”，而系统轨迹则从外向内穿过这些碗壁。

一旦你将 Lyapunov 稳定性理解为“轨迹沿着势能碗下滑”，所有相关定理就变得显而易见：

• 若 $\dot V \leq 0$ ，轨迹永远不会“上坡”，因此会始终停留在包含初始点的最小“碗”内（稳定）。
• 若 $\dot V < 0$ 严格成立，轨迹将持续下滑直至抵达碗底（渐近稳定）。
• 若 $\dot V > 0$ ，轨迹将向外爬升，系统不稳定。

LaSalle 不变集原理对此进行了推广：当 $\dot V$ 在某个集合上恒为零时，轨迹最终会收敛到该集合中最大的不变子集。

相图与零线#

对于系统 $x' = f(x,y),\ y' = g(x,y)$ ：

• $x$ -零线是曲线 $f(x,y) = 0$ ，在此处 $\dot x = 0$ ，轨迹垂直穿过该线。
• $y$ -零线是曲线 $g(x,y) = 0$ ，轨迹水平穿过该线。
• 平衡点位于两条零线的交点处。
• 各区域中 $f$ 和 $g$ 的符号决定了局部流向。

零线分析是手绘相图最经济高效的方法。

线性化：局部为真，全局意外#

在双曲平衡点附近，Jacobian 矩阵的特征值完全决定了局部相图结构（Hartman-Grobman 定理）。然而，一旦远离平衡点，情况可能截然不同。阻尼单摆就是一个鲜明例证：

Lotka-Volterra 捕食-猎物模型#

模型#

• $x$ ：猎物数量，$y$ ：捕食者数量
• 平凡平衡点 $(0,0)$ ：鞍点
• 共存平衡点 $(\gamma/\delta,\ \alpha/\beta)$ ：Jacobian 特征值为 $\pm i\sqrt{\alpha\gamma}$ （中心）

使得所有轨道均为闭合曲线。其时间序列与相平面如下所示：

循环过程的文字描述#

1. 猎物丰富 → 捕食者食物充足 → 捕食者数量上升。
2. 捕食者众多 → 猎物被大量捕食 → 猎物数量骤降。
3. 猎物稀少 → 捕食者食物短缺 → 捕食者数量下降。
4. 捕食者减少 → 猎物压力减轻 → 猎物数量回升，回到第 1 步。

模型局限#

• 结构不稳定：哪怕加入微小扰动项，闭合轨道也会被破坏。
• 缺乏环境承载力：若无捕食者，猎物将无限增长。
• 忽略空间分布与时滞效应。

这些缺陷推动了下一节更现实模型的发展。

竞争模型：四种可能结局#

乘积 $\alpha_{12}\,\alpha_{21}$ （即种间相互抑制强度）决定了系统最终走向哪一种情形。

这正是竞争排斥原理的数学表达。

Van der Pol 振子：非线性阻尼催生极限环#

其精妙之处在于阻尼系数 $-\mu(1 - x^2)$ ：

• 当 $|x| < 1$ 时，阻尼为负——系统向自身注入能量。
• 当 $|x| > 1$ 时，阻尼为正——能量逐渐耗散。

内部轨迹向外扩张，外部轨迹向内收缩，最终所有轨迹都收敛到同一个稳定极限环——这是一个孤立的周期轨道，吸引其吸引域内的所有初始点。

实际应用场景包括：心脏起搏细胞、神经元动作电位、电子振荡电路、间歇泉喷发、商业周期等。

梯度系统与 Hamilton 系统：两种特殊世界#

梯度系统：$\mathbf{x}' = -\nabla V$ #

• 轨迹沿势能 $V$ 的最陡下降方向运动。
• $\dot V = -|\nabla V|^2 \leq 0$ ，势能始终不增。
• 不可能存在周期轨道（绕行一圈后势能不可能更低）。
• 机器学习中的梯度下降法正是其离散版本。

Hamilton 系统：$x' = \dfrac{\partial H}{\partial y},\ y' = -\dfrac{\partial H}{\partial x}$ #

• $H$ 沿任意轨迹守恒（$\dot H = 0$ ）。
• 相空间体积守恒（Liouville 定理）。
• 轨道即为 $H$ 的等值线。
• 无阻尼单摆是典型实例。

这两类系统分别代表了耗散谱的两个极端。

Poincaré-Bendixson 定理：二维系统无法产生混沌#

定理（Poincaré-Bendixson）：对于光滑的二维连续动力系统，若有界轨迹不趋向任何平衡点，则必收敛于一条闭合轨道。

换句话说，在二维系统中，长期行为只有两种可能：趋于平衡点或趋于周期轨道，混沌无处容身。

其背后的关键是 Jordan 曲线定理——闭合轨道将平面分割为内外两部分，从而“困住”轨迹。一旦引入第三维，轨迹便可从轨道“上方”逃逸，混沌的大门由此开启（见第九章）。

Bendixson 判据（排除闭合轨道）：若在单连通区域内，$\partial f/\partial x + \partial g/\partial y$ 恒为非零且符号不变，则该区域内不存在闭合轨道。

数值方法速查#

1
2
3
4
5
6

def rk4_step(f, x, t, h):
    k1 = h * np.array(f(x, t))
    k2 = h * np.array(f(x + k1/2, t + h/2))
    k3 = h * np.array(f(x + k2/2, t + h/2))
    k4 = h * np.array(f(x + k3,    t + h))
    return x + (k1 + 2*k2 + 2*k3 + k4) / 6

练习题#

概念题

1. 为什么 $y' = y^2$ 不满足叠加原理？请给出具体反例。
2. 为何二维连续系统不可能出现混沌？
3. 手绘 Duffing 方程 $x' = y,\ y' = -x + x^3 - 0.2y$ 的相图。

计算题

4. 对系统 $x' = x - x^3,\ y' = -y$ ，找出所有平衡点并分类。
5. 证明 $V = x^2 + y^2$ 是系统 $x' = -x + y^2,\ y' = -y$ 的 Lyapunov 函数。
6. 验证 $H = \delta x - \gamma\ln x + \beta y - \alpha\ln y$ 在 Lotka-Volterra 模型中守恒。

编程题

7. 复现图 5 中的四种竞争情形，并为每个吸引域着色。
8. 数值估算 Van der Pol 振子在 $\mu \in \{0.1, 0.5, 1, 3, 10\}$ 时的周期 $T(\mu)$ 。
9. 比较 Euler 与 RK4 在 Van der Pol 方程上的精度，找出 Euler 方法失效的 $\mu$ 值。

极限环稳定性：Floquet 乘子#

平衡点的稳定性由 Jacobian 特征值判断，而极限环 $\gamma(t)$ 的稳定性则由 Floquet 乘子 决定——它回答的问题是：“若将轨迹横向扰动一小段距离，绕行一周后该扰动会被放大还是缩小？”

对上述方程在一个周期内积分，得到单值矩阵 $\Phi(T)$ 。其特征值 $\mu_i$ 即为 Floquet 乘子。其中必有一个 $\mu_1 = 1$ （对应沿轨道的相位平移），其余乘子决定横向稳定性：

• 若其余所有 $|\mu_i| < 1$ ，则为稳定极限环（吸引）。
• 若存在任一 $|\mu_i| > 1$ ，则为不稳定。

数值上，可同时积分极限环及其扰动的基本矩阵：

1
2
3
4
5
6
7
8
9

def monodromy(f, jac, gamma0, T, n=2000):
    from scipy.integrate import solve_ivp
    d = len(gamma0)
    def rhs(t, y):
        x, Phi = y[:d], y[d:].reshape(d, d)
        return np.concatenate([f(x), (jac(x) @ Phi).ravel()])
    y0 = np.concatenate([gamma0, np.eye(d).ravel()])
    sol = solve_ivp(rhs, (0, T), y0, dense_output=False, rtol=1e-10)
    return sol.y[d:, -1].reshape(d, d)

以 $\mu = 1$ 的 Van der Pol 系统为例，其非平凡 Floquet 乘子约为 $0.04$ ，表明极限环具有极强的吸引力——轨迹通常只需绕行几圈便会迅速锁定到环上。

一个具体算例：兔子-狐狸的 30 个时间单位#

把 Lotka-Volterra 写成具体数字，远比抽象推导更能让人记住"中心轨道"的味道。我取这样一组参数：兔子初值 $x_0 = 50$ ，狐狸初值 $y_0 = 25$ ，繁殖率 $\alpha = 0.6$ ，被捕食率 $\beta = 0.02$ ，狐狸自然死亡率 $\gamma = 0.4$ ，狐狸增殖率 $\delta = 0.01$ 。代入共存平衡点公式 $(\gamma/\delta,\ \alpha/\beta) = (40,\ 30)$ ，初值正好落在右下方一点点。

我用 RK4、步长 $h = 0.01$ 跑到 $t = 30$ ，关键节点的数值是这样的：

到 $t = 30$ 几乎回到出发点——周期约 $T \approx 2\pi/\sqrt{\alpha\gamma} = 2\pi/\sqrt{0.24} \approx 12.83$ 时间单位，30 大约对应 2.34 个周期，残差完全来自 RK4 的截断误差，而不是动力学本身。守恒量 $H$ 浮动在第四位小数，说明步长选得够细。

这一组数字最直观的一点：兔子峰值 71.8（$t \approx 15$ ）远高于狐狸的稳态 30，狐狸峰值 41.9（$t \approx 20$ ）也明显高于兔子的稳态 40——也就是说，振幅由初始的偏离量决定，不是方程参数本身。如果我换初值 $(45, 28)$ ，轨迹会缩成一个更小的环，$H$ 取另一个常数。这正是"中心"和"螺旋"的本质区别：螺旋有渐近半径，中心则任何半径都能存在。

直觉与陷阱：为什么"中心"在数值实验里几乎活不下来#

刚学这套理论的同学最容易踩的坑，是把 Lotka-Volterra 的中心当成稳健现象。它结构不稳定：只要给方程加一个微小扰动项（比如让兔子自身有承载力 $\dot x = \alpha x(1-x/K) - \beta xy$ ，$K = 1000$ ），原本闭合的轨道立刻变成稳定螺旋——所有轨迹都会缠绕到唯一的内部不动点。

我自己第一次用 Euler 方法跑这个系统时也被坑过：步长 $h = 0.1$ 看起来挺合理，但 100 个时间单位后轨道半径增长了 30%，看上去像是不稳定螺旋。换成 RK4 才发现，那个"螺旋"完全是 Euler 方法的数值耗散反着加了能量——一个非保结构积分器在保守系统上必然失败。教训是：遇到 Hamilton 结构或者中心型平衡点，不要用 Euler / RK4，要用 leapfrog 或辛 Runge-Kutta。这条原则在第十一章和第十七章会反复出现。

更深的陷阱是：纯虚特征值 $\pm i\sqrt{\alpha\gamma}$ 在线性化层面给出的"中心"，对完整非线性方程并不一定成立。Hartman-Grobman 定理在这种非双曲情形是失效的。Lotka-Volterra 之所以仍然是闭轨，靠的是守恒量 $H$ 的存在；一旦扰动破坏了 $H$ ，闭轨立刻消失。

反例与应用：极限环不是"非线性的中心"#

很多人把极限环和中心搞混。它们的区别非常关键：

• 中心：相图里有连续一族闭合轨道（每个初值对应一条），能量层级是连续参数；任意小扰动可以破坏这种结构。
• 极限环：只有一条孤立的闭合轨道，邻近所有轨迹都被它吸引（或排斥）；扰动后可能形变，但不会消失。

Van der Pol 振子是极限环的标准例子：心脏起搏细胞的去极化-复极化节律就是这种结构——细胞被噪声、温度、激素扰动，但搏动周期始终被锁在同一条极限环上，这是生命系统鲁棒性的物理基础。如果心脏靠的是 Lotka-Volterra 那种中心，任何一次外界扰动都会改变心率，人就死定了。

工业上的反例是激光腔的弛豫振荡、电网的次同步振荡、化工反应器的温度震荡——这些都是稳定极限环。设计者关心的不是"会不会振"，而是"振幅落在哪条环上"。Andronov-Hopf 分岔是从稳定不动点诞生极限环的标准机制，第七章的 Hopf 分岔正是它在 ODE 层面的一面镜子。

数值常见陷阱：刚性、辛结构与相位漂移#

陷阱一：刚性系统使用显式方法。Lotka-Volterra 模型尚可应付，但若修改为 $\dot x = -1000(x - \cos t) - \sin t$ ，显式 RK4 方法为保证稳定性需满足步长 $h < 2/1000$ ，即每秒模拟需 500 步以上。此时应改用 LSODA 或 BDF 等隐式方法。诊断依据：Jacobian 特征值跨越多个数量级。

陷阱二：Hamilton 系统使用非辛积分器。以单摆为例，其守恒量为 $H = p^2/2 - \cos q$ 。使用 RK4 积分 $10^4$ 步后，能量会出现明显漂移。而 leapfrog 或 symplectic Euler 方法虽有震荡，但能量始终保持有界。这也正是 PDE-ML 系列中辛网络章节的动机所在。

陷阱三：长期积分中的相位漂移。在极限环或准周期轨道上，传统数值方法会累积线性增长的相位误差。若希望在 $10^4$ 个周期后仍能绘制清晰相图，标准 RK4 无法胜任。此时应采用几何积分器（geometric integrator），并配合较大的步长。

接入机器学习：Neural ODE 的相图即网络行为#

Neural ODE $\dot z = f_\theta(z, t)$ 将分类任务转化为“将每个输入点通过流形映射至正确区域”。训练完成后，$f_\theta$ 的相图本身就直接体现了网络的行为：

• 拓扑障碍真实存在。Dupont et al. (2019) 指出，Neural ODE 无法学习需要两条不相交曲线交叉的映射——因为 ODE 流是同胚映射。解决方案是将状态空间扩充至 $\mathbb{R}^{d+k}$ （即 ANODE），这正是相平面理论的直接推论。
• 吸引子可视化提供可解释性。训练一个二维 Neural ODE 分类器后，绘制 $f_\theta$ 的向量场，即可看到决策边界由不动点及其稳定流形自然构成，无需借助 SHAP 等事后解释工具。
• 稳定性正则化。在损失函数中加入 $\|J_\theta\|_F^2$ 项，可限制流的局部拉伸程度。这是一种源于本章相图理论的隐式对抗鲁棒性机制，而非机器学习领域的经验之谈。

总结#

下一步#

下一章遇到 ODE 最迷人也最反常识的一面——混沌。Lorenz 在 1961 年发现：一个三维确定性系统可以产生看似随机的行为，对初值的依赖敏感到无法长期预测天气。我们会从单变量分岔（鞍-结分岔、Hopf 分岔）讲到 Lorenz 方程的奇异吸引子。

工程实践片段：竞争模型在产品市场里的具象化#

把竞争方程从生态搬到商业，参数立刻有了实际含义。设两家产品的用户数 $x$ 、$y$ ，自然增长率 $r_1 = r_2 = 0.3$ （月度），各自市场容量 $K_1 = 100$ 万、$K_2 = 80$ 万。互相挤占系数 $\alpha_{12}$ （产品二每万用户挤占产品一多少潜在市场）和 $\alpha_{21}$ 决定结局。当 $\alpha_{12} \alpha_{21} < 1$ 时双方共存，这是大多数细分市场的常态——比如 iOS 与 Android、淘宝与京东、微信与 QQ。当 $\alpha_{12} \alpha_{21} > 1$ 时则进入双稳态，谁先抢到用户谁赢——这正是社交网络早期"网络效应"主导的真实写照。

我读 Fortune 100 公司的财报建模时一个反复出现的教训：用线性外推预测市场份额几乎一定错。用户增长是非线性的，Bass 扩散模型本质就是 Lotka-Volterra 的特殊情形。把"竞争系数"接到广告投入、价格弹性、用户切换成本这些可观测变量上，才能把这套模型从黑板搬进董事会。

零线分析的工程价值在这里也很直接：产品经理盯着的"用户数 vs 留存率"二维相图，本质上就是 $x$ -零线和 $y$ -零线的交点。两条曲线的几何关系决定了哪个产品最终主导，跟生态学里的物种竞争是同一份数学。这是我特别喜欢这套理论的地方——抽象到极致反而能跨学科迁移。

一句话提醒：相图比方程更重要#

学非线性动力学最大的认知跃迁，是从"找解析解"切换到"读相图"。Poincaré 在十九世纪末就强调过这一点：天体力学里的三体问题没有解析解，但相空间的几何结构清楚地告诉我们哪些初值会发散、哪些会准周期、哪些会被某条不变流形捕获。把同样的思路搬到生态、化学、神经科学，零线、不动点、极限环、不变流形就是工程师真正需要的语言。一旦学会把方程在脑子里翻译成相图，所有的"数值实验失败"基本都能在画图阶段提前预见——这是我自己学了十几年最有用的一句话。

参考文献#

• Strogatz, Nonlinear Dynamics and Chaos, Westview Press (2015)
• Murray, Mathematical Biology I, Springer (2002)
• Hirsch, Smale, & Devaney, Differential Equations, Dynamical Systems, and Chaos, Academic Press (2012)
• Verhulst, Nonlinear Differential Equations and Dynamical Systems, Springer (1996)