常微分方程（二）：一阶微分方程的求解方法

银行存款的复利、肝脏代谢一片药物、水箱里的盐慢慢稀释、电容在电源驱动下逐步充满——这些看似毫不相关的现象，背后都是同一类方程：一阶常微分方程。本章的目标只有一个：让你看见一个一阶方程的瞬间，就能判断它属于四种典型形态中的哪一种，并立刻知道该使用哪种求解技巧。这四套方法看似各自独立，其实背后是同一个思想——找到一个变量替换或乘子，把方程化成"一眼就能积出来"的形式。

本章要点

如何识别可分离方程并直接积分得到通解
用积分因子把线性方程整理成"完美导数"的形式
恰当方程的判定条件，以及它的几何意义：解就是某个势函数的等高线
伯努利替换 $v = y^{1-n}$ 如何把一整族非线性方程线性化
五个真实应用：药物半衰期、RC 充电、盐水混合、Logistic 增长、复利增长——每一个都是可以直接复用的小型案例

前置知识

第一章的概念：什么是 ODE、初值问题、方向场长什么样
标准积分技巧：换元法、分部积分、部分分式分解

1. 可分离方程：最自然的形式

1.1 形态识别

一阶 ODE 是可分离的，当且仅当右端可以分解为"$x$ 的函数"乘以"$y$ 的函数"：

$$ \frac{dy}{dx} \;=\; g(x)\,h(y). $$

正是这种因式分解，让我们能把所有 $y$ 移到左边、所有 $x$ 移到右边：

$$ \frac{dy}{h(y)} \;=\; g(x)\,dx \qquad\Longrightarrow\qquad \int \frac{dy}{h(y)} \;=\; \int g(x)\,dx + C. $$

常数 $C$ 是族参数——每取一个值，就挑出族中的一条解曲线。下图中，$\frac{dy}{dx} = -x/y$ 经过分离后得到 $y\,dy + x\,dx = 0$，积分即为 $x^2 + y^2 = C$，对应一族同心圆。方向场上的小段处处与这些圆相切——这就是"解"在几何上的全部含义。

分离变量把一个二维方向场坍缩成一族一维曲线 $x^2+y^2=C$。

1.2 例题：指数增长与衰减

求解 $\dfrac{dy}{dx} = ky$，其中 $k$ 是常数。

分离： $\;\dfrac{dy}{y} = k\,dx$
两边积分： $\;\ln|y| = kx + C_1$
取指数： $\;y = A e^{kx}$，其中 $A = \pm e^{C_1}$ 把符号和常数一起吸收

更深的一点：这个解里 $k$ 控制的是时间尺度，不是形状。两种典型解读：

$k > 0$：每隔 $\ln 2 / k$ 翻一倍——无资源限制下的种群、连续复利、失控的链式反应
$k < 0$：半衰期 $\ln 2 / |k|$——放射性衰变、电容放电、肝脏代谢药物

1.3 应用：药物代谢与半衰期法则

绝大多数药物遵循一阶清除：从血液中被清除的速率，正比于此刻血液中的含量：

$$ \frac{dC}{dt} = -kC, \qquad C(t) = C_0 e^{-kt}. $$

半衰期 $t_{1/2}$ 是浓度降为初值一半所需的时间。代入 $C_0/2 = C_0 e^{-k t_{1/2}}$ 解出来就是那个万能公式：

$$ t_{1/2} = \frac{\ln 2}{k} \approx \frac{0.693}{k}. $$

布洛芬的 $t_{1/2} \approx 2$ 小时。从一颗 400 mg 的剂量出发：

时间	残余	含义
0 h	400 mg	峰浓度
2 h	200 mg	一个半衰期
4 h	100 mg	两个半衰期
6 h	50 mg	已落到治疗窗口以下

这就解释了为什么布洛芬建议每 4–6 小时服用一次：再长，浓度跌出有效区间；再短，则会蓄积。临床上常说"五个半衰期清除完毕"（约 99% 已代谢），也是这个公式的直接推论。

1.4 应用：Logistic 方程

纯指数增长 $P' = rP$ 是一种幻觉：现实中没有任何东西可以无限增长。Verhulst 在 1838 年加了一个刹车——环境容纳量 $K$，得到了至今仍是数学生态学基石的方程：

$$ \frac{dP}{dt} = r P\left(1 - \frac{P}{K}\right). $$

这是可分离方程。左边用部分分式：

$$ \int \!\left(\frac{1}{P} + \frac{1/K}{1 - P/K}\right)\,dP \;=\; \int r\,dt, $$

取指数后即得 Logistic 曲线：

$$ P(t) \;=\; \frac{K}{1 + \left(\dfrac{K}{P_0} - 1\right)e^{-rt}}. $$

不必真的解出来，仅看方程本身就能读出三件事：

$P \ll K$：括号近似为 $1$，近似指数增长
$P = K/2$：右端等于 $rK/4$，是整条曲线的最大增速
$P \to K$：括号趋于 $0$，种群饱和

下图右侧把 $\dot P$ 作为 $P$ 的函数画出来——这是一条向下的抛物线，零点在 $P=0$ 与 $P=K$，峰值正好在 $K/2$。

不同初值下的 Logistic 曲线，以及增长率 dP/dt 关于 P 的抛物线。

左：所有解最终都收敛到 $K$。右：增长率 $\dot P = rP(1-P/K)$ 在 $P=K/2$ 处达到最大值。

2. 一阶线性方程：积分因子法

2.1 标准形式

$$ \frac{dy}{dx} + P(x)\,y \;=\; Q(x). $$

这是应用数学的"主力工种"。任何"某个量的变化率，正比于自身加上一个外部驱动"的过程，最终都长这个样子。

2.2 关键技巧

两边同乘 积分因子

$$ \mu(x) \;=\; e^{\int P(x)\,dx}. $$

这样设计 $\mu$，是为了让乘积法则把左边收成一个完整的导数：

$$ \frac{d}{dx}\bigl[\mu(x)\,y\bigr] \;=\; \mu(x)\,Q(x). $$

直接积一次再除以 $\mu$：

$$ \boxed{\,y(x) \;=\; \frac{1}{\mu(x)} \left[\int \mu(x)\,Q(x)\,dx + C\right].\,} $$

整套方法就这一行。它为什么管用？因为 $\mu$ 满足 $\mu' = \mu P$，正好契合乘积法则的需求。

2.3 例题：$y' + 2xy = x$

步骤	计算	结果
1. 识别	$P = 2x$，$\;Q = x$	—
2. 积分因子	$\mu = \exp\!\int 2x\,dx$	$\mu = e^{x^2}$
3. 乘 $\mu$	左端化为 $\frac{d}{dx}[e^{x^2}y]$	右端 $= x e^{x^2}$
4. 积分	$\int x e^{x^2} dx = \tfrac12 e^{x^2}$	$e^{x^2}y = \tfrac12 e^{x^2} + C$
5. 解出	两边除以 $e^{x^2}$	$y = \tfrac12 + C e^{-x^2}$

当 $x \to \pm\infty$ 时 $C e^{-x^2} \to 0$，所有解最终都汇入平衡值 $y = 1/2$；积分常数 $C$ 只在原点附近的一段瞬态里起作用。

2.4 积分因子在做什么：三幅图说清楚

取 $y' + 2y = 4$（于是 $P = 2$，$\mu = e^{2x}$，平衡点 $y=2$）：

y’ + 2y = 4 的方向场与解族、积分因子 e^{2x}、以及左端坍缩为完整导数。

左：每条解 $y = 2 + Ce^{-2x}$ 都迅速被吸引到平衡 $y=2$。中：积分因子 $e^{2x}$ 不过是一条干净的指数。右：乘上 $\mu$ 之后，$\mu(x)y(x)$ 与 $\int 4\mu\,dx$ 完全吻合——左端真的成了一个原函数。

2.5 应用：RC 电路充电

电阻 $R$ 与电容 $C$ 串联，外加电压源 $V_s$，按基尔霍夫电压定律：

$$ RC\,\frac{dV_c}{dt} + V_c \;=\; V_s. $$

对照标准形式：$P(t) = 1/RC$，$Q(t) = V_s/RC$。积分因子 $\mu(t) = e^{t/RC}$。当 $V_c(0)=0$ 时，最终结果是

$$ V_c(t) \;=\; V_s\bigl(1 - e^{-t/\tau}\bigr), \qquad \tau = RC. $$

时间常数 $\tau$ 是每个电气工程师都背得滚瓜烂熟的：

时刻	$V_c/V_s$	工程上的读法
$\tau$	63.2%	“一个时间常数”
$3\tau$	95.0%	在大多数逻辑电路里可视为充满
$5\tau$	99.3%	工程意义上"完全充满"

把符号翻一下、令 $V_s = 0$，同一个公式就描述了放电过程——也解释了相机闪光灯为什么需要等一两秒才能再次拍照。

3. 恰当方程

3.1 几何思想

把方程写成微分形式：

$$ M(x,y)\,dx + N(x,y)\,dy = 0. $$

假设存在一个"势函数" $F(x,y)$ 满足

$$ \frac{\partial F}{\partial x} = M, \qquad \frac{\partial F}{\partial y} = N. $$

那么 $M\,dx + N\,dy$ 就是 $F$ 的全微分 $dF$，方程 $dF = 0$ 等价于"$F$ 沿着解曲线保持常数"，于是整个解族就是

$$ F(x,y) \;=\; C. $$

解曲线就是 $F$ 的等高线。这个观点其实非常深刻：求解 ODE 的工作被替换成了"找出一个标量函数 $F$，让它的等高线图本身就是解族"。

3.2 恰当性判定

势函数存在当且仅当混合偏导相等（即 $F_{xy} = F_{yx}$）：

$$ \boxed{\,\dfrac{\partial M}{\partial y} \;=\; \dfrac{\partial N}{\partial x}.\,} $$

如果不满足，方程就不是恰当的——但还有挽救的办法（见 §3.4）。

3.3 例题：恰当方程的求解步骤

考虑 $(2x + y)\,dx + (x + 2y)\,dy = 0$。先验证：$M_y = 1 = N_x$，恰当。然后还原 $F$：

将 $M$ 关于 $x$ 积分：$\;F = x^2 + xy + g(y)$
关于 $y$ 求导，与 $N$ 比对：$\;x + g'(y) = x + 2y \Rightarrow g'(y) = 2y \Rightarrow g(y) = y^2$
因此 $F(x,y) = x^2 + xy + y^2$，解族为 $x^2 + xy + y^2 = C$

几何上这是一族倾斜的椭圆。下图把梯度场 $\nabla F = (M, N)$ 与等高线一并画出；梯度处处与等高线垂直——这就是"$dF = 0$ 沿解成立"在几何上的精确含义：梯度告诉你"上山方向"，而解曲线始终待在同一高度。

F(x,y) = x^2 + xy + y^2 的等高线与梯度场，二者处处垂直。

每一条等高线就是一条解。灰色箭头是 $\nabla F$，与等高线垂直——这正是"沿解 $dF = 0$“的几何表达。

3.4 救场：把非恰当方程变恰当

若 $M_y \neq N_x$，则乘以 $\mu(x,y)$ 修正它。两种最常见、有显式公式的情况：

条件	取
$(M_y - N_x)/N$ 仅依赖 $x$	$\mu(x) = \exp\!\int \dfrac{M_y - N_x}{N}\,dx$
$(N_x - M_y)/M$ 仅依赖 $y$	$\mu(y) = \exp\!\int \dfrac{N_x - M_y}{M}\,dy$

这与线性方程的积分因子是同一个构造——事实上，线性方程的方法只是这条更普遍方法的特例。

4. 伯努利方程

4.1 形态

$$ \frac{dy}{dx} + P(x)\,y \;=\; Q(x)\,y^{n}, \qquad n \neq 0, 1. $$

$n = 0$ 时它已经是线性的；$n = 1$ 时退化为可分离的。其余情形里，$y^n$ 这一项破坏了线性。

4.2 替换方法

令 $v = y^{1-n}$，则 $\dfrac{dv}{dx} = (1-n)\,y^{-n}\,\dfrac{dy}{dx}$。把原方程乘以 $(1-n) y^{-n}$ 重新整理：

$$ \frac{dv}{dx} + (1-n)P(x)\,v \;=\; (1-n)Q(x). $$

这是关于 $v$ 的线性方程。用第 2 节的积分因子求出 $v$，再通过 $y = v^{1/(1-n)}$ 回代即可。

4.3 替换的几何含义

取 $y' + y = y^2$（$n = 2$，$v = 1/y$）。替换后变成线性方程 $v' - v = -1$，解是绕平衡点 $v = 1$ 的指数族。在 $v$-平面里，整族是规整、笔直的；而在原 $y$-平面里，它们弯曲、甚至会在有限时间内炸到无穷。

y’ + y = y^2 在原 (x, y) 平面与替换后 (x, v) 平面（v = 1/y）的对比。

同样五条轨迹。左 $(x, y)$：非线性、$y=1$ 是不稳定平衡。右 $(x, v)$：线性、$v=1$ 是稳定平衡（从下方对应 $y=1$）。

同样的替换还能"还原"出 Logistic 方程：它是 $n = 2$、$P(t) = -r$、$Q(t) = -r/K$ 的伯努利方程。也就是说，§1.4 里我们费力做的部分分式积分，本质上只是一次伯努利替换 + 线性求解。

5. 综合应用：水箱里的盐

这是经典的"房室模型”——同样的骨架在流行病学、药代动力学、化学工程里都会再次出现。

问题设定。 1000 L 水箱初装纯水。浓度 $2$ g/L 的盐水以 $5$ L/min 流入；水箱充分搅拌，以同样 $5$ L/min 流出，体积保持不变。求盐量 $Q(t)$。

建模。 流入减流出：

$$ \underbrace{\frac{dQ}{dt}}{\text{变化率}} = \underbrace{5 \cdot 2}{\text{流入: g/min}}

\underbrace{5 \cdot \tfrac{Q}{1000}}_{\text{流出: g/min}} = 10 - \tfrac{Q}{200}, \qquad Q(0) = 0. $$

求解。 这是线性方程：$Q' + \tfrac{1}{200}Q = 10$。积分因子 $\mu = e^{t/200}$，结果为

$$ Q(t) \;=\; 2000\bigl(1 - e^{-t/200}\bigr). $$

读懂答案。 平衡值 $Q^\ast = 2000$ g，对应浓度 $2$ g/L——正好等于流入浓度。时间常数 $\tau = 200$ min 控制收敛快慢。把同一方程的初值改为 $Q(0) > Q^\ast$，描述的就是一个偏咸的水箱被慢慢冲淡——下图把两种情形画在一起，它们殊途同归地走向同一个平衡。

两种情形，同一个平衡。从 $Q(0)=0$ 灌满与从 $Q(0)=3000$ 冲淡，最终都到达流入浓度 $2$ g/L；只是偏离的方向不同。

6. 数值方法预览

只要 $f(x,y)$ 稍微复杂一点，闭式技巧就会失效——而真实世界里的方程绝大多数都属于"复杂"那一类。这时候的办法是：沿着方向场一步一步走。

6.1 欧拉方法

最朴素的想法：每一步用切线代替曲线。

$$ y_{n+1} \;=\; y_n + h\,f(t_n, y_n). $$

单步局部截断误差：$O(h^2)$
走 $1/h$ 步后的全局误差：$O(h)$，故称"一阶方法"

便宜、好写，但精度通常不够。

6.2 四阶 Runge-Kutta（RK4）

数值 ODE 的标准主力。每步评估四次斜率再加权平均：

$$ k_1 = f(t_n, y_n), \qquad k_2 = f\!\left(t_n + \tfrac{h}{2},\, y_n + \tfrac{h}{2}k_1\right), $$$$ k_3 = f\!\left(t_n + \tfrac{h}{2},\, y_n + \tfrac{h}{2}k_2\right), \qquad k_4 = f(t_n + h,\, y_n + h k_3), $$$$ y_{n+1} \;=\; y_n + \frac{h}{6}\bigl(k_1 + 2k_2 + 2k_3 + k_4\bigr). $$

全局误差 $O(h^4)$——同样步长下比欧拉精度高出四个数量级，代价只是多三次函数求值。

6.3 看方向场，能预判数值解的形态

在敲代码之前，先把右端画出来。两条表面相似的方程，长期行为可能截然不同：

左（线性）：唯一的直线解 $y=x+1$，其余指数发散。右（伯努利）：三个平衡点 $y=-1, 0, +1$，所有轨迹被压向 $\pm 1$，与初值符号无关。同属一阶方程，几何形态完全不同。

6.4 SciPy 三行就能搞定

牛顿冷却定律：$T' = -k(T - T_{\text{env}})$，$T(0) = 90$°C，$T_{\text{env}} = 20$°C，$k = 0.1$ /min。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
from scipy.integrate import solve_ivp
import numpy as np

def cooling(t, T, k=0.1, T_env=20):
    return -k * (T - T_env)

sol = solve_ivp(cooling, [0, 60], [90.0], method="RK45", dense_output=True)

t = np.linspace(0, 60, 300)
T = sol.sol(t)[0]
print(f"t=10 min 时: {T[50]:.1f} C")
print(f"t=30 min 时: {T[150]:.1f} C")

solve_ivp 会自适应步长。绝大多数问题不需要你手写 RK4；你需要知道的是：什么时候标准求解器会吃力（刚性方程——见第 11 章）。

7. 总结：四问筛选法

碰到一个一阶方程，按下表顺序自问自答，第一个"是"就告诉你方法。

提问	若是	方法
能写成 $y' = g(x)\,h(y)$ 吗？	可分离	移项后两边积分
能写成 $y' + P(x)y = Q(x)$ 吗？	线性	积分因子 $\mu = e^{\int P\,dx}$
能写成 $M\,dx + N\,dy = 0$ 且 $M_y = N_x$ 吗？	恰当	还原势函数 $F$；解为 $F = C$
能写成 $y' + Py = Qy^n$（$n \notin \{0,1\}$）吗？	伯努利	替换 $v = y^{1-n}$，再走第二行
都不行？	—	数值方法（第 11 章），或寻找一个聪明的换元

这四套方法并不彼此独立。线性方程的积分因子正是"修正非恰当方程的积分因子"在特定情形下的退化；伯努利替换则把非线性问题归约为线性。所以学好这四种方法，本质上是在反复练习同一个观念——找一个变量替换或乘子，把方程化成可以一眼积分的形式——只是在四种典型场景下落地了四次。

8. 练习题

基础题。

求解 $y' = y^2 \sin x$，$y(0) = 1$。解会在何处爆破？
求解 $y' + y = e^{-x}$。指出哪一部分是瞬态、哪一部分是稳态。
验证 $(2xy + 3)\,dx + (x^2 + 1)\,dy = 0$ 是恰当方程并求解。
1000 L 水箱中初有 50 kg 盐。纯水以 10 L/min 流入，混合液以同样速率流出。求 $Q(t)$，并求盐量减半所需的时间。

进阶题。

求解伯努利方程 $y' + y = y^3$。画出相线、指出所有平衡点。
对 $y' = 3y^{2/3}$，$y(0) = 0$，给出至少两条不同的解（定义在 $[0, \infty)$ 上）。这为什么不与存在唯一性定理矛盾？
细菌种群满足 $P' = 0.5\,P(1 - P/1000)$，$P(0) = 50$。显式求出 $P(t)$，并求增长率达到最大值的时刻。

编程题。

从零实现欧拉方法和 RK4，比较它们在 $y' = -y$，$y(0) = 1$，$[0, 5]$ 区间上、步长 $h = 0.5, 0.1, 0.01$ 时的全局误差。在 log-log 坐标下绘制误差关于 $h$ 的曲线，验证斜率分别为 $1$ 和 $4$。
自己复现 §6.3 的方向场对比：用 numpy.meshgrid 与 matplotlib.pyplot.quiver 画方向场，再用 scipy.integrate.solve_ivp 计算轨迹。试一试 $y' = y - y^5$，并指出它的所有平衡点。

参考资料

Boyce & DiPrima，Elementary Differential Equations and Boundary Value Problems，Wiley (2012)。本科段的标准参考；第 2–3 章按经典严谨度覆盖了上面的全部内容。
Zill，A First Course in Differential Equations with Modeling Applications，Cengage (2017)。偏应用风格；混合槽与电路例题处理得格外干净。
Strogatz，Nonlinear Dynamics and Chaos，Westview (2014)。第 2 章对 Logistic 方程的相线分析，是这一章后必读的延伸。
SciPy 文档：scipy.integrate 。solve_ivp 的方法选择（RK45、Radau、BDF）与事件检测都在此查阅。

系列导航


上一章	第一章：微分方程的起源与直觉
当前	第二章：一阶微分方程的求解方法
下一章	第三章：高阶线性微分方程