常微分方程模型 — 从解析求解到定性分析

主讲：复旦大学数学科学学院 蔡志杰

基本概念

含有未知函数导数的等式

$$\frac{\mathrm{d}x}{\mathrm{d}t}=f(t,x)\text{\hspace{1em}(1)}$$

微分方程中未知函数的导数的最高阶数称为这个方程的阶。

设 $\mathbf{x}=\mathbf{x}(t)$ 是一元向量值函数，

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{f}(t,\mathbf{x})\text{\hspace{1em}(2)}$$

是一阶常微分方程组，其中 $\mathbf{f}(t,\mathbf{x})$ 是已知的向量值函数。

含有未知函数导数满足的隐式方程 (组)

$$\mathbf{F}\left(\frac{{\mathrm{d}}^k\mathbf{x}}{\mathrm{d}{t}^k},\frac{{\mathrm{d}}^{k-1}\mathbf{x}}{\mathrm{d}{t}^{k-1}},\cdots ,\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t},\mathbf{x},t\right)=0\text{\hspace{1em}(3)}$$

称为 $k$ 阶隐式常微分方程 (组)。

如果函数 $\varphi (t)$ 在 $a<t<b$ 内一阶连续可微，且

$$\frac{\mathrm{d}\varphi (t)}{\mathrm{d}t}=f(t,\varphi (t))(a<t<b)\text{\hspace{1em}(4)}$$

则称函数 $x=\varphi (t)$ 为常微分方程的解，区间 $a<t<b$ 称为解 $x=\varphi (t)$ 的定义区间。

显式解 vs. 隐式解 vs. 数值解：能写成 $x=\varphi (t)$ 的称为显式解；只能写成 $G(t,x)=0$ 形式的称为隐式解（如分离变量法常得到隐式解）；当解析式无法求出或过于复杂时，退而求其次求数值解——这正是本文后半部分 Euler 法、Runge-Kutta 法要解决的问题。

解的存在唯一性：对初值问题 $x^{\prime }=f(t,x),x(t_0)=x_0$，若 $f$ 在 $(t_0,x_0)$ 邻域内连续且关于 $x$ 满足 Lipschitz 条件，则局部存在唯一解（Picard-Lindelöf 定理）。这是后续数值方法"有意义"的前提——如果解本身不唯一，数值格式收敛到哪一支都无从谈起。

基本步骤

1. 根据实际问题建立相应的常微分方程，给出适当的初始条件；
2. 给出解析解或者数值解；
3. 对常微分方程解的形态做出分析（即本文最后的"定性分析"：不求出具体解析式，而是研究解随时间的长期行为——是否趋于平衡、是否振荡、是否发散）；
4. 用得到的结果解释实际问题，对问题的发展变化趋势做出预测。

CUMCM2024A, 2022A, 2019A, 2019B, 2016A, 2014A

模型的归结

案例1 凶杀案何时发生

某天在一住宅内发生了一起凶杀案。接到报告后，法医于 $23:00$ 赶到现场，测量到死者的体温是 $30.8^{\circ }\mathrm{C}$，再过 $1$ 小时后测得其体温为 $29.1^{\circ }\mathrm{C}$，当时室温是 $28^{\circ }\mathrm{C}$，试估计受害者的死亡时间。

牛顿冷却定律：物体的冷却速率与物体温度和环境温度的温度差成正比

$$\frac{\mathrm{d}T}{\mathrm{d}t}=-k(T-T_h)\text{\hspace{1em}(5)}$$

初始条件 $T(0)=T_0$。

求解与结果：这是变量可分离方程，解得

$$T(t)=T_h+(T_0-T_h){\mathrm{e}}^{-kt}$$

设 $23:00$ 为 $t=0$，此时 $T(0)=30.8$；$t=1$（即 $0:00$）时 $T(1)=29.1$，环境温度 $T_h=28$。代入：

$$29.1-28=(30.8-28){\mathrm{e}}^{-k}⟹{\mathrm{e}}^{-k}=\frac{1.1}{2.8}⟹k=\ln \frac{2.8}{1.1}\approx 0.933(\text{/h})$$

再利用人体正常体温 $37^{\circ }\mathrm{C}$ 作为死亡时刻 $t^{\ast }$（相对 $23:00$，$t^{\ast }<0$）的初值：

$$37-28=2.8{\mathrm{e}}^{-k{t}^{\ast }}⟹t^{\ast }=-\frac{1}{k}\ln \frac{9}{2.8}\approx -1.25\text{h}$$

即死亡时间约在 $23:00$ 之前 $1.25$ 小时，也就是约 21:45。

建模要点：该模型假设 $k$ 为常数，与体表面积/体型等因素有关，实际法医学中常需结合经验常数修正；这是"单参数辨识 + 反推初值"的典型数学建模范式，也是 CUMCM 中"根据观测点反推模型参数"的常见套路。

案例2 导弹跟踪问题

某部队的一导弹基地发现正北方向 $120\mathrm{km}$ 处海面上有敌艇一艘以 $90km/h$ 的速度向正东方向行驶。该基地立即发射导弹跟踪追击敌艇。导弹速度为 $450km/h$，自动导航系统使导弹在任一时刻都能对准敌艇。试问导弹在何时何处击中敌艇？

导弹在 $t$ 时刻的位置为 $P(x(t),y(t))$，取基地为原点，正东为 $x$ 轴正向，正北为 $y$ 轴正向，记 $H=120$（初始北向距离），$v_e=90$（敌艇速度），$v_w=450$（导弹速度）。敌艇 $t$ 时刻位置为 $(v_{e}t,H)$。

$${\left(\frac{\mathrm{d}x}{\mathrm{d}t}\right)}^2+{\left(\frac{\mathrm{d}y}{\mathrm{d}t}\right)}^2=v_w^2\text{\hspace{1em}(6)}$$

$$\frac{\mathrm{d}y}{\mathrm{d}x}=\frac{H-y}{v_{e}t-x}\text{\hspace{1em}(7)}$$

初始条件 $x(0)=0,y(0)=0$。

补充说明：式 (6) 表示导弹速度大小恒定为 $v_w$；式 (7) 表示导弹速度方向始终指向敌艇当前位置（自动导航条件）。这是经典的**追击曲线（pursuit curve）**问题。

求解思路：直接联立 (6)(7) 含 $t$ 较繁琐，标准处理是消去参数 $t$，把方程化为只含 $x,y$ 的二阶方程：

1. 由 (7) 变形：$v_{e}t-x=\frac{(H-y)}{y^{\prime }}$，两边对 $x$ 求导，并利用弧长关系 $\frac{\mathrm{d}s}{\mathrm{d}x}=\frac{v_w}{v_e}\sqrt{1+y^{\prime 2}}$（因为导弹走过的弧长与时间成正比、而 $v_{e}t$ 对 $x$ 的导数隐含着这个比例关系），可以消去 $t$，得到一个只含 $y^{\prime },y^{\prime \prime }$ 的一阶方程（令 $p=y^{\prime }$）：

$$\frac{v_e}{v_w}x{p}^{\prime \prime }\text{（形式）}⟹\text{可分离变量求出 }p=y^{\prime }\text{ 关于 }x\text{ 的表达式}$$

2. 再积分一次得到 $y=y(x)$（追击曲线的显式方程），是 $x$ 的幂函数与常数的组合（这是追击问题的标准结果形式）；
3. 令 $y(x)=H$（即导弹到达敌艇纬线 $y=H$，此时按照 (7) 分母为零对应"击中"）反解出击中点的 $x$ 坐标，再由弧长–时间关系 $t=\frac{s}{v_w}$ 求出击中时刻。
4. 由于 $v_w/v_e=5>1$（导弹更快），可以证明追击曲线在有限时间内与直线 $y=H$ 相交，即导弹必能击中敌艇；若 $v_w\le v_e$ 则永远追不上。

建模要点：此类问题的核心技巧是"利用速度比消去时间参数 $t$，把追击问题转化为纯粹的几何微分方程"，这是自动导引类问题（导弹制导、猎豹追羚羊等）的通用建模范式。

案例3 传染病模型

假设在传染病流行范围内只有两类人，一类是易感者，记为 S 类；另一类是感病者，记为 I 类。

假设易感者与感病者在人群中混合充分均匀，易感者感病的机会与他接触感病者的机会成正比，并且传染率为常数。不考虑出生与死亡以及人群的迁出与迁入。

$S(t)$：易感者的比例。$I(t)$：得病者的比例。

$$\frac{\mathrm{d}S}{\mathrm{d}t}=-kIS\text{\hspace{1em}(8)}$$

$$\frac{\mathrm{d}I}{\mathrm{d}t}=kIS\text{\hspace{1em}(9)}$$

这是最简单的 SI 模型。由 (8)(9) 相加可知 $S+I=\text{const}=1$（比例守恒），代入 (9) 得到关于 $I$ 的 Logistic 方程：

$$\frac{\mathrm{d}I}{\mathrm{d}t}=kI(1-I)$$

解为 $I(t)=\frac{I_0{\mathrm{e}}^{kt}}{1-I_0+I_0{\mathrm{e}}^{kt}}$，是典型的 S 形曲线，$t\to \infty$ 时 $I\to 1$，即所有人最终都会被感染（因为没有考虑康复/移出）。

现在考虑 $R(t)$：移出者（康复或隔离）的比例，$h$ 为移出率：

$$\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-kIS\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =kIS-hI\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =hI\end{array}\text{\hspace{1em}(10)}$$

这就是经典的 SIR 模型，满足 $S+I+R=1$。

基本再生数 $R_0$：由第二式 $\frac{\mathrm{d}I}{\mathrm{d}t}=I(kS-h)$，定义

$$R_0=\frac{k{S}_0}{h}$$

• 若 $R_0>1$（即 $S_0>h/k$），初期 $I$ 递增，疫情会爆发；
• 若 $R_0\le 1$，$I$ 从一开始就单调递减，疫情不会流行。

这就是流行病学中著名的"阈值现象（threshold phenomenon）"，也是 SIR 模型定性分析（而非求出解析解）能给出的最重要结论——后文"定性分析"一节会用相平面法重新推出这一结果。

案例4 单摆

在一根长为 $l$ 的细线下端，悬挂一个质量为 $m$ 的球。略微移动后，球在重力作用下来回摆动，这种装置称为单摆。

假设细线不会伸长，也没有重量，在悬点没有摩擦力，试建立单摆的运动方程。

$$\frac{{\mathrm{d}}^2\theta }{\mathrm{d}{t}^2}=-\frac{g}{l}\sin \theta \text{\hspace{1em}(11)}$$

小角度近似：当摆角 $\theta$ 很小时，$\sin \theta \approx \theta$，(11) 退化为线性方程

$${\theta }^{\prime \prime }+{\omega }^2\theta =0,\omega =\sqrt{g/l}$$

由前文二阶常系数方程的结论（$\Delta <0$ 情形），通解为 $\theta (t)=A\cos (\omega t+\phi )$，即简谐振动，周期 $T=2\pi \sqrt{l/g}$（这正是中学物理中单摆周期公式的来历）。

非线性情形：原方程 (11) 无法用初等函数求出解析解（其精确解需用椭圆积分表示），这正是"定性分析"和数值方法（如 Runge-Kutta 法）大显身手的地方：数值法可以直接对 (11) 降阶为方程组后逐步积分；定性分析则可以不求解就画出相图，判断摆是做周期摆动还是"翻转"（大能量时绕悬点整圈转动）。

方程的求解

分离变量法

$$\frac{\mathrm{d}x}{\mathrm{d}t}=f(t)g(x)$$

称为变量可分离的方程。

(1) 若 $g(x)=0$ 有实根 $x=x^{\ast }$，则

$$x(t)\equiv x^{\ast }$$

是方程的常数解，称为不动点。

(2) 若 $g(x)\ne 0$，则

$$\frac{\mathrm{d}x}{g(x)}=f(t)\mathrm{d}t$$

$${\int }_{x_0}^x\frac{\mathrm{d}\xi }{g(\xi )}={\int }_{t_0}^{t}f(\tau )\mathrm{d}\tau +C$$

例：$\frac{\mathrm{d}x}{\mathrm{d}t}=x(1-x)$（案例3中的 Logistic 方程原型）。$g(x)=x(1-x)$ 的零点 $x^{\ast }=0,1$ 是两个不动点；对 $x\ne 0,1$ 用部分分式 $\frac{1}{x(1-x)}=\frac{1}{x}+\frac{1}{1-x}$ 积分即得前文给出的 S 形解。

降阶法

$n$ 阶方程

$$F(t,x,\dot{x},\cdots ,x^{(n)})=0$$

当 $n\ge 2$ 时，统称为高阶方程。

如果方程中不显含未知函数 $x$，或者更一般地，不显含未知函数及其直到 $k-1(k\ge 1)$ 阶导数的方程

$$F(t,x^{(k)},\cdots ,x^{(n)})=0$$

令 $y=x^{(k)}$，则

$$F(t,y,\cdots ,y^{(n-k)})=0$$

是关于 $y$ 的 $n-k$ 阶方程，阶数降低了 $k$；解出 $y(t)$ 后再积分 $k$ 次即得 $x(t)$。

另一种常见降阶：若方程不显含自变量 $t$（自治方程），如 $x^{\prime \prime }=f(x,x^{\prime })$，可令 $p=x^{\prime }$，并把 $p$ 看作 $x$ 的函数，用 $x^{\prime \prime }=p\frac{\mathrm{d}p}{\mathrm{d}x}$ 降阶为关于 $p(x)$ 的一阶方程。案例2消去参数 $t$ 的处理，以及单摆方程的相图分析，用的都是这一技巧。

常数变易法

考虑一阶线性微分方程

$$\frac{\mathrm{d}x}{\mathrm{d}t}+p(t)x=q(t)$$

第一步：求对应齐次方程的通解，得 $x_h(t)=C{\mathrm{e}}^{-\int p(t)\mathrm{d}t}$。

第二步：常数变易，设 $x(t)=C(t){\mathrm{e}}^{-\int p(t)\mathrm{d}t}$，代入求得 $C^{\prime }(t)=q(t){\mathrm{e}}^{\int p(t)\mathrm{d}t}$，积分得 $C(t)$。

第三步：通解

$$x(t)=C{\mathrm{e}}^{-\int p(t)\mathrm{d}t}+{\mathrm{e}}^{-\int p(t)\mathrm{d}t}\int q(t){\mathrm{e}}^{\int p(t)\mathrm{d}t}\mathrm{d}t$$

例：$x^{\prime }+2tx=t$。这里 $p(t)=2t,q(t)=t$，$\int p\mathrm{d}t=t^2$。齐次解 $x_h=C{\mathrm{e}}^{-t^2}$；$C(t)=\int t{\mathrm{e}}^{t^2}\mathrm{d}t=\frac{1}{2}{\mathrm{e}}^{t^2}+C$，故通解 $x(t)=\frac{1}{2}+C{\mathrm{e}}^{-t^2}$（可直接验证特解 $x_p=1/2$ 满足原方程）。

二阶齐次常系数线性方程

$$\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^2}+p\frac{\mathrm{d}x}{\mathrm{d}t}+qx=0,{\lambda }^2+p\lambda +q=0,\Delta =p^2-4q$$

$\Delta$	根的情况	通解	力学含义（如阻尼振动）
$\Delta >0$	两不等实根 ${\lambda }_1,{\lambda }_2$	$C_1{\mathrm{e}}^{{\lambda }_{1}t}+C_2{\mathrm{e}}^{{\lambda }_{2}t}$	过阻尼，单调衰减无振荡
$\Delta =0$	二重根 ${\lambda }_1={\lambda }_2=-p/2$	$(C_1+C_{2}t){\mathrm{e}}^{{\lambda }_{1}t}$	临界阻尼
$\Delta <0$	共轭复根 $\alpha \pm \mathrm{i}\beta$	${\mathrm{e}}^{\alpha t}(C_1\cos \beta t+C_2\sin \beta t)$	欠阻尼，振荡衰减（$p>0$）或增幅（$p<0$）

例：单摆小角度方程 ${\theta }^{\prime \prime }+{\omega }^2\theta =0$ 对应 $p=0,q={\omega }^2>0$，$\Delta =-4{\omega }^2<0$，$\alpha =0,\beta =\omega$，故 $\theta (t)=C_1\cos \omega t+C_2\sin \omega t$——纯振荡不衰减（因为没有阻尼项 $p{\theta }^{\prime }$），与前面的物理直觉一致。

数值计算方法

基本概念

$$\frac{\mathrm{d}y}{\mathrm{d}x}=f(x,y),y(x_0)=y_0$$

数值解法：寻求解 $y(x)$ 在一系列离散节点 $x_1<x_2<\cdots <x_n<x_{n+1}<\cdots$ 上的近似值 $y_1,y_2,\cdots$。相邻两个节点的间距 $h_i=x_i-x_{i-1}$ 称为步长（通常取定步长 $h$）。

误差的两个层面：

• 局部截断误差：假设 $y_n$ 精确，只看走一步 $y_n\to y_{n+1}$ 引入的误差，记为 $O(h^{p+1})$；
• 整体（全局）截断误差：累积 $N=(x_n-x_0)/h$ 步之后的总误差，记为 $O(h^p)$——因为累积了 $N\propto 1/h$ 步局部误差。

一个方法若局部截断误差是 $O(h^{p+1})$，则称其为 $p$ 阶方法。

Euler 法

从初始点 $P_0(x_0,y_0)$ 出发沿切线方向推进：

$$y_{n+1}=y_n+hf(x_n,y_n)$$

向量形式：${\mathbf{y}}_{n+1}={\mathbf{y}}_n+h\mathbf{f}(x_n,{\mathbf{y}}_n)$。

误差阶：对精确解做泰勒展开 $y(x_n+h)=y_n+hf(x_n,y_n)+\frac{h^2}{2}{y}^{\prime \prime }(\xi )$，可见 Euler 法每步丢弃了 $O(h^2)$ 的项，故局部截断误差 $O(h^2)$，是一阶方法，整体误差 $O(h)$——精度较低，步长减半误差只减半，实际计算中很少直接用于精细问题，但因形式简单、稳定性分析直观，常作为教学起点与其他方法的比较基准。

隐式改进：梯形法/后退 Euler 法：把 $f(x_{n+1},y_{n+1})$ 也考虑进来（如 $y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_{n+1})]$）可以提高精度到二阶，并显著改善稳定性，但每步需要解一个（通常非线性的）方程，计算量更大——这也是为什么显式 Runge-Kutta 法（不需要解方程）在实践中更受欢迎的原因之一。

Runge-Kutta 法

基本思路

$$\mathbf{y}(x_{n+1})=\mathbf{y}(x_n)+{\int }_{x_n}^{x_{n+1}}\mathbf{f}(x,\mathbf{y}(x))\mathrm{d}x$$

$${\int }_{x_n}^{x_{n+1}}\mathbf{f}(x,\mathbf{y}(x))\mathrm{d}x\approx h\sum _{i=1}^r{c}_i\mathbf{f}(x_n+{\lambda }_{i}h,\mathbf{y}(x_n+{\lambda }_{i}h))$$

思路解读：这一步的本质是用数值积分（求积公式）去近似右端的定积分——$c_i$ 相当于求积公式的权重，${\lambda }_i$ 相当于取样点在区间内的相对位置（$0\le {\lambda }_i\le 1$）。唯一的困难是被积函数中出现的 $\mathbf{y}(x_n+{\lambda }_{i}h)$ 是未知的（因为我们还没求出解），Runge-Kutta 法用 ${\mathbf{y}}_n$ 加上前面已算出的 ${\mathbf{K}}_j$ 的组合去预测这个值，这就是 ${\mathbf{K}}_i$ 的递推定义的来源。

2阶 Runge-Kutta 法

$$\begin{array}{rl}{\mathbf{y}}_{n+1}& ={\mathbf{y}}_n+h(c_1{\mathbf{K}}_1+c_2{\mathbf{K}}_2)\\ {\mathbf{K}}_1& =\mathbf{f}(x_n,{\mathbf{y}}_n)\\ {\mathbf{K}}_2& =\mathbf{f}(x_n+{\lambda }_{2}h,{\mathbf{y}}_n+{\mu }_{21}h{\mathbf{K}}_1)\end{array}$$

把 ${\mathbf{K}}_2$ 按 $h$ 展开，与精确解泰勒展开比较到 $h^2$ 项，得阶条件：

$$1-c_1-c_2=0,\frac{1}{2}-c_2{\lambda }_2=0,\frac{1}{2}-c_2{\mu }_{21}=0$$

由后两式知 ${\lambda }_2={\mu }_{21}$（这一相等关系保证了 ${\mathbf{K}}_2$ 是在同一个"预测点"上求值，是显式 RK 法的一般规律：$c_i={\sum }_j{a}_{ij}$），三个方程、四个未知数 $(c_1,c_2,{\lambda }_2,{\mu }_{21})$，有一个自由度，取 $c_2=a\ne 0$：

$$c_2=a,c_1=1-a,{\lambda }_2={\mu }_{21}=\frac{1}{2a}$$

常用取值：

$a$	方法	公式
$a=1$	改进 Euler 法（梯形法的显式版，Heun 法）	${\mathbf{y}}_{n+1}={\mathbf{y}}_n+\frac{h}{2}({\mathbf{K}}_1+{\mathbf{K}}_2)$，${\mathbf{K}}_2=\mathbf{f}(x_n+h,{\mathbf{y}}_n+h{\mathbf{K}}_1)$
$a=1/2$	中点法	${\mathbf{y}}_{n+1}={\mathbf{y}}_n+h{\mathbf{K}}_2$，${\mathbf{K}}_2=\mathbf{f}(x_n+\frac{h}{2},{\mathbf{y}}_n+\frac{h}{2}{\mathbf{K}}_1)$

四阶 Runge-Kutta 法

$$\begin{array}{rl}{\mathbf{y}}_{n+1}& ={\mathbf{y}}_n+\frac{h}{6}({\mathbf{K}}_1+2{\mathbf{K}}_2+2{\mathbf{K}}_3+{\mathbf{K}}_4)\\ {\mathbf{K}}_1& =\mathbf{f}(x_n,{\mathbf{y}}_n)\\ {\mathbf{K}}_2& =\mathbf{f}\left(x_n+\frac{h}{2},{\mathbf{y}}_n+\frac{h}{2}{\mathbf{K}}_1\right)\\ {\mathbf{K}}_3& =\mathbf{f}\left(x_n+\frac{h}{2},{\mathbf{y}}_n+\frac{h}{2}{\mathbf{K}}_2\right)\\ {\mathbf{K}}_4& =\mathbf{f}(x_n+h,{\mathbf{y}}_n+h{\mathbf{K}}_3)\end{array}$$

局部截断误差 $O(h^5)$，整体误差 $O(h^4)$，是四阶方法，每步需计算 4 次 $\mathbf{f}$。这是工程与建模竞赛中最常用的数值积分格式。

方法对比小结

方法	每步计算 $f$ 次数	局部误差	整体误差	特点
Euler 法	1	$O(h^2)$	$O(h)$	最简单，精度低
2阶 RK	2	$O(h^3)$	$O(h^2)$	精度与工作量折中
4阶 RK	4	$O(h^5)$	$O(h^4)$	最常用，精度高

数值方法的选用建议：对光滑、非刚性（non-stiff）方程，经典 RK4 定步长足够；若方程刚性较强（不同分量变化速率差异巨大，如化学反应网络），显式 RK 法可能需要极小步长才能稳定，此时应改用隐式方法（如隐式 Euler、Radau 方法）或专门的刚性求解器（如 MATLAB 的 ode15s、ode23s）。

定性分析

前面几节要么求出解析解，要么用数值方法逐点计算近似解——但很多实际问题（尤其是非线性方程，如单摆的大角度摆动、SIR 模型、捕食者-被捕食者模型）既求不出解析解，逐点数值计算也无法直接给出"长期行为是怎样的"这类结构性结论。定性分析的目标就是：不求解方程，直接从方程本身的结构判断解的长期趋势、是否稳定、是否振荡。

相空间与自治系统

对于自治系统（右端不显含自变量 $t$）

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{f}(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^n$$

$\mathbf{x}$ 所在的空间 ${\mathbb{R}}^n$ 称为相空间（phase space）；解 $\mathbf{x}(t)$ 在相空间中画出的曲线称为轨线（trajectory）；所有轨线构成的图像称为相图（phase portrait）。$n=2$ 时称为相平面，是最直观、也最常分析的情形（如单摆的 $(\theta ,{\theta }^{\prime })$ 平面、SIR 模型的 $(S,I)$ 平面）。

自治系统的一个关键性质：相空间中轨线互不相交（除非在同一条轨线上），因为过每一点的解唯一（存在唯一性定理）。这保证了相图是"结构清晰"的——一旦画出，就能预判从任意初值出发解的走向，无需重新求解。

平衡点与线性化

满足 $\mathbf{f}({\mathbf{x}}^{\ast })=0$ 的点 ${\mathbf{x}}^{\ast }$ 称为系统的平衡点（equilibrium），对应常数解 $\mathbf{x}(t)\equiv {\mathbf{x}}^{\ast }$（与分离变量法中"不动点"概念一致，是其向量形式的推广）。

判断平衡点附近解的行为，标准做法是线性化：令 $\mathbf{u}=\mathbf{x}-{\mathbf{x}}^{\ast }$（小扰动），泰勒展开到一阶：

$$\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}\approx J({\mathbf{x}}^{\ast })\mathbf{u},J({\mathbf{x}}^{\ast })={\frac{\partial \mathbf{f}}{\partial \mathbf{x}}|}_{{\mathbf{x}}^{\ast }}(\text{雅可比矩阵})$$

平衡点附近的行为由 $J({\mathbf{x}}^{\ast })$ 的特征值 $\lambda$ 决定（这是 Hartman–Grobman 定理的内容：只要 $\lambda$ 都不含零实部，线性化系统与原系统在局部拓扑等价）。二维情形的分类：

特征值 ${\lambda }_1,{\lambda }_2$	平衡点类型	稳定性
均为负实数	稳定结点（node）	渐近稳定
均为正实数	不稳定结点	不稳定
一正一负	鞍点（saddle）	不稳定
共轭复数，实部 $<0$	稳定焦点（螺旋收敛）	渐近稳定
共轭复数，实部 $>0$	不稳定焦点	不稳定
纯虚数	中心（center），闭轨线	临界（李雅普诺夫稳定，非渐近）

李雅普诺夫稳定性

对一般非线性系统，即使线性化失效（特征值含零实部）或需要判断全局（而非局部）稳定性，可用 Lyapunov 第二方法：若能构造一个函数 $V(\mathbf{x})$（通常理解为"能量"）满足

$$V({\mathbf{x}}^{\ast })=0,V(\mathbf{x})>0(\mathbf{x}\ne {\mathbf{x}}^{\ast }),\frac{\mathrm{d}V}{\mathrm{d}t}=\nabla V\cdot \mathbf{f}(\mathbf{x})\le 0$$

则 ${\mathbf{x}}^{\ast }$ 是（李雅普诺夫）稳定的；若 $\dot{V}<0$ 严格成立（$\mathbf{x}\ne {\mathbf{x}}^{\ast }$），则是渐近稳定的。这一方法的优点是不需要求解方程，只需要找到一个合适的 $V$（好比证明系统的"能量"始终耗散）。

极限环

对于非线性自治系统，还可能出现既非平衡点、又不是简单振荡的孤立闭合轨线——极限环（limit cycle）：附近的轨线都会螺旋逼近这条闭轨（稳定极限环）或远离它（不稳定极限环）。这是非线性系统特有的现象（线性系统只能有"中心"型的一族闭轨，不可能孤立）。经典例子是范德波尔振荡子 $x^{\prime \prime }-\mu (1-x^2)x^{\prime }+x=0$，描述自激振荡电路。

案例串讲：用相平面法重新分析前面的模型

（1）单摆方程 (11) 的相图

令 $x_1=\theta ,x_2={\theta }^{\prime }$，方程化为一阶方程组

$$x_1^{\prime }=x_2,x_2^{\prime }=-\frac{g}{l}\sin x_1$$

平衡点：$x_2=0,\sin x_1=0\Rightarrow x_1=k\pi$。

• $x_1=0$（摆锤竖直向下）：雅可比矩阵特征值为 $\pm \mathrm{i}\sqrt{g/l}$（纯虚数）——中心，对应摆锤附近的小幅周期摆动，与前面小角度近似 ${\theta }^{\prime \prime }+{\omega }^2\theta =0$ 的振荡解完全一致；
• $x_1=\pi$（摆锤竖直向上，倒立）：特征值为 $\pm \sqrt{g/l}$（一正一负）——鞍点，对应不稳定的倒立平衡。

相图上可以看到：能量较小时轨线是绕 $x_1=0$ 的闭合曲线（周期摆动）；能量足够大时轨线不再闭合，而是穿过 $x_1=\pi$ 附近，对应摆锤"翻过头顶"整圈旋转——这正是数值解法（RK4）无法直接"看出"、但相平面定性分析一眼就能给出的结构性结论。

（2）SIR 模型 (10) 的相平面分析

由于 $S+I+R=1$，系统本质上是二维的，取 $(S,I)$ 平面，方程为

$$\frac{\mathrm{d}S}{\mathrm{d}t}=-kSI,\frac{\mathrm{d}I}{\mathrm{d}t}=I(kS-h)$$

两式相除消去 $t$：$\frac{\mathrm{d}I}{\mathrm{d}S}=-1+\frac{h}{kS}$，可直接积分得到轨线方程 $I+S-\frac{h}{k}\ln S=\text{const}$，从而不用求出 $S(t),I(t)$ 的具体时间历程，就能画出 $(S,I)$ 平面上的轨线族。

由 $\frac{\mathrm{d}I}{\mathrm{d}t}=I(kS-h)$ 可知：

• 当 $S>h/k$ 时 $I$ 递增（疫情上升）；
• 当 $S<h/k$ 时 $I$ 递减（疫情消退）；
• $S=h/k$ 处 $I$ 取极大值。

这就从相平面角度重新导出了前文的阈值现象（$R_0=k{S}_0/h>1$ 则爆发）——且不需要求出任何一个显式的时间解，是定性分析方法威力的典型体现。

（3）Lotka–Volterra 捕食者-被捕食者模型（补充案例）

$$\frac{\mathrm{d}x}{\mathrm{d}t}=\alpha x-\beta xy,\frac{\mathrm{d}y}{\mathrm{d}t}=\delta xy-\gamma y$$

（$x$：猎物种群，$y$：捕食者种群）。非零平衡点 $(x^{\ast },y^{\ast })=(\gamma /\delta ,\alpha /\beta )$ 处雅可比矩阵特征值为纯虚数 $\pm \mathrm{i}\sqrt{\alpha \gamma }$——中心，对应种群数量在平衡点附近做周期性振荡（这是生态学中"猎物-捕食者数量周期性此消彼长"现象的数学解释），可用 Lyapunov 函数 $V(x,y)=\delta x-\gamma \ln x+\beta y-\alpha \ln y$ 证明这些闭轨在全局意义上也是稳定的（$\dot{V}\equiv 0$，说明 $V$ 是守恒量）。

---

小结：定性分析与数值方法的关系

	定性分析	数值方法（Euler / RK）
目标	长期行为结构（稳定性、周期性、阈值）	具体某条轨线的逐点近似值
是否需要求解方程	否	是（逐步递推）
典型工具	平衡点分类、Lyapunov 函数、相图	差分递推公式
适用场景	判断"会不会发生""趋势如何"	需要具体数值、绘制精确轨迹

在实际建模（如 CUMCM 竞赛）中，两者通常配合使用：先用定性分析快速判断模型的合理性与解的大致形态（例如 SIR 模型判断是否爆发），再用 Runge-Kutta 等数值方法算出具体的数量-时间曲线，供预测和决策使用。

参数辨识

前面几节建立的模型 $\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{f}(t,\mathbf{x};\mathbf{\theta })$ 中，往往含有若干未知参数 $\mathbf{\theta }$（如案例1中的冷却系数 $k$，案例3 SIR 模型中的传染率 $k$、移出率 $h$）。这些参数通常无法直接测量，只能借助观测数据反推——这就是**参数辨识（parameter identification）**问题，是前文"基本步骤"第 1、2 步之间承上启下的关键环节：模型结构定了，但要真正"能用"，还需要把参数的具体数值确定下来。

问题描述

已知 $0<t_1<t_2<\cdots <t_n$ 上 $\mathbf{x}$ 的观测值

$$\mathbf{x}(t_i)={\mathbf{x}}_i(i=1,2,\cdots ,n)$$

希望根据这些信息来确定未知参数 $\mathbf{\theta }$ 的值，这就是参数辨识问题。

参数辨识一般采用最小二乘法，将参数辨识问题转化为优化问题来处理：

$$\underset{\mathbf{\theta }}{\min }E=\sum _{i=1}^n\Vert \mathbf{x}(t_i,\mathbf{\theta })-{\mathbf{x}}_i{\Vert }^2$$

即寻找一组参数 $\mathbf{\theta }$，使模型在各观测时刻给出的预测值 $\mathbf{x}(t_i,\mathbf{\theta })$ 与实际观测值 ${\mathbf{x}}_i$ 的误差平方和最小。

与前文方法的衔接

$\mathbf{x}(t_i,\mathbf{\theta })$ 从哪里来？ 这正是本文前面几节的作用所在：

• 若方程能求出解析解（如案例1的 $T(t)=T_h+(T_0-T_h){\mathrm{e}}^{-kt}$），$\mathbf{x}(t_i,\mathbf{\theta })$ 直接代入解析式计算即可，$E(\mathbf{\theta })$ 是 $\mathbf{\theta }$ 的显式函数；
• 若方程无解析解（如非线性的 SIR 模型、单摆大角度方程），则对每一组试探参数 $\mathbf{\theta }$，都需要用 Runge-Kutta 等数值方法先把方程从 $t_0$ 积分到各观测时刻 $t_i$，得到 $\mathbf{x}(t_i,\mathbf{\theta })$，再代入 $E$ 中评估——这意味着参数辨识实际上是一个**"外层优化 + 内层数值积分"的嵌套过程**：外层不断调整 $\mathbf{\theta }$，内层每次都要重新做一次 RK4 积分。

案例1 重新表述为参数辨识问题：未知参数 $\mathbf{\theta }=k$，观测数据为 $(t_1,T_1)=(0,30.8)$、$(t_2,T_2)=(1,29.1)$（$n=2$ 个数据点，$1$ 个未知参数，恰好能唯一确定，无需优化，直接解方程即可，这是最小二乘法的退化情形）。若有更多测温数据点（$n>1$ 个未知参数个数），则 $E(k)={\sum }_i(T_h+(T_0-T_h){\mathrm{e}}^{-k{t}_i}-T_i{)}^2$ 一般无法使 $E=0$，只能求使 $E$ 最小的 $k$——这才是最小二乘法真正发挥作用（消除测量误差影响）的典型场景。

求解方法

（1）线性最小二乘：若 $\mathbf{x}(t,\mathbf{\theta })$ 关于 $\mathbf{\theta }$ 是线性的（如某些可通过换元转化为线性关系的模型，例如对 $T(t)$ 两边取对数后 $\ln (T-T_h)=\ln (T_0-T_h)-kt$ 关于 $k$ 就是线性的），$E(\mathbf{\theta })$ 是 $\mathbf{\theta }$ 的二次函数，令 ${\nabla }_{\mathbf{\theta }}E=0$ 可直接解出正规方程组（线性代数问题），一步得到解析解，这也是案例1中"先取对数、再线性拟合"这类做法的理论依据。

（2）非线性最小二乘（一般情形）：$E(\mathbf{\theta })$ 关于 $\mathbf{\theta }$ 非线性时，通常没有解析解，需要迭代求解，常用算法包括：

• 梯度下降法：沿 $-{\nabla }_{\mathbf{\theta }}E$ 方向迭代更新 $\mathbf{\theta }$，简单但收敛慢；
• Gauss–Newton 法：利用残差 $r_i=\mathbf{x}(t_i,\mathbf{\theta })-{\mathbf{x}}_i$ 对 $\mathbf{\theta }$ 的一阶灵敏度（雅可比矩阵 $\partial r_i/\partial \mathbf{\theta }$）构造线性近似，收敛速度快于梯度下降；
• Levenberg–Marquardt 法：在 Gauss–Newton 基础上加入阻尼项，兼顾梯度下降的稳健性和 Gauss–Newton 的收敛速度，是实践中最常用的非线性最小二乘算法（MATLAB 的 lsqcurvefit、lsqnonlin，Python 的 scipy.optimize.curve_fit、least_squares 底层都基于此）。

灵敏度方程：若用 Gauss–Newton / L-M 法，需要 $\partial \mathbf{x}(t_i,\mathbf{\theta })/\partial \mathbf{\theta }$。当 $\mathbf{x}(t,\mathbf{\theta })$ 本身由 ODE 数值积分得到而无解析式时，可以对原方程关于 $\mathbf{\theta }$ 求偏导，得到灵敏度方程（sensitivity equations）：

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \mathbf{x}}{\partial \mathbf{\theta }}\right)=\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\cdot \frac{\partial \mathbf{x}}{\partial \mathbf{\theta }}+\frac{\partial \mathbf{f}}{\partial \mathbf{\theta }}$$

这是关于 $\partial \mathbf{x}/\partial \mathbf{\theta }$ 的一组新的（线性）常微分方程，可以和原方程 ${\mathbf{x}}^{\prime }=\mathbf{f}$ 联立后一起用 RK4 数值积分求出，从而避免每次调整 $\mathbf{\theta }$ 都用差分近似梯度（既省计算量又更精确）。实践中也常直接用有限差分近似梯度（简单但需多算几次 ODE，精度和效率都不如灵敏度方程），初学建模时更常用。

完整讲义