常微分方程（一）：微分方程的起源与直觉

你身边的一切都在变化。 咖啡在冷却，人口在增长，单摆在摆动，病毒在传播，股价在波动，行星在运行。这些系统几乎没有谁能用「某物等于多少」来描述——它们只能用「某物变化得多快」来刻画。这第二种描述方式，正是微分方程存在的理由；学会读它，就是学会读物理与生物所用的那门语言。

本章从最朴素的直觉出发重建你的认知。我们从一杯咖啡开始，推导出与放射性衰变、电容放电同一形式的方程，再向上爬到方向场、分类，以及那条决定「ODE 是否有意义」的存在唯一性定理。

本章要点

微分方程到底是什么——超越符号本身的理解
ODE 与 PDE、阶数、线性——以及为什么这些区分都值得记住
三个经典模型（冷却、衰变、振动），它们会在后续反复出现
方向场：不解方程也能「看见」解的走向
存在唯一性定理（Picard–Lindelöf），并给出一个具体的反例
用 Python 和 scipy 求解你的第一个微分方程，并知道何时该信任它

前置知识

一元微积分：导数、积分、链式法则
基本 Python：NumPy 数组与 Matplotlib 画图

1. 什么是微分方程？从一杯咖啡说起

你冲了一杯 $T_0 = 90^\circ\mathrm{C}$ 的咖啡，放在 $T_\text{env} = 20^\circ\mathrm{C}$ 的房间里。几分钟后，它凉了。两个问题同样自然：

它现在多少度？——一个代数问题。
它现在凉得多快？——一个微积分问题。

第一个问题看上去更简单，但物理学几乎总是直接回答第二个。牛顿在 1701 年给出了答案：冷却的速率与物体和环境之间的温差成正比。

$$ \frac{dT}{dt} = -k\,\bigl(T - T_\text{env}\bigr). $$

这一行就是一个微分方程——把函数 $T(t)$ 与它自己的导数关联起来的方程。要求的不是某个数，而是整个函数。

符号	含义
$T(t)$	时刻 $t$ 的温度
$dT/dt$	温度的瞬时变化率
$T_\text{env}$	环境温度（$20^\circ\mathrm{C}$）
$k > 0$	冷却常数，取决于杯子、表面积、空气流动

那个负号至关重要：当 $T > T_\text{env}$ 时右端为负，所以 $T$ 下降；当 $T < T_\text{env}$ 时右端为正，所以 $T$ 上升。这条方程已经把热力学第二定律写进自己的骨子里了。

解之前先「读」方程

不动笔，仅凭这条方程的结构，我们就能得出几条物理结论：

温度越接近 $T_\text{env}$，温差越小，冷却越慢。冷却是减速的。
取 $T = T_\text{env}$ 时右端为零——咖啡一旦达到室温，就不再变化。这是一个平衡点。
$k$ 越大（薄壁杯、通风）冷却越快；$k$ 越小（保温杯）冷却越慢。

这是整门课程的第一课：好的 ODE 在被解出之前就已经告诉你它的定性行为。

用 Python 求解

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def coffee_cooling(T, t, k, T_env):
    """牛顿冷却定律：dT/dt = -k (T - T_env)"""
    return -k * (T - T_env)

T0, T_env, k = 90.0, 20.0, 0.1
t = np.linspace(0, 60, 500)
T_sol = odeint(coffee_cooling, T0, t, args=(k, T_env))

plt.figure(figsize=(10, 5))
plt.plot(t, T_sol, color="#2563eb", linewidth=2, label="咖啡温度")
plt.axhline(T_env, color="#ef4444", linestyle="--", label=f"室温 = {T_env} °C")
plt.xlabel("时间 (分钟)"); plt.ylabel("温度 (°C)")
plt.title("牛顿冷却定律"); plt.legend(); plt.tight_layout(); plt.show()

print(f"10 分钟后: {T_sol[np.abs(t-10).argmin()][0]:.1f} °C")
print(f"30 分钟后: {T_sol[np.abs(t-30).argmin()][0]:.1f} °C")

闭合解

这条方程友善到可以手算。令 $u = T - T_\text{env}$，则 $du/dt = -ku$，分离变量并代回得：

$$ \boxed{\,T(t) = T_\text{env} + \bigl(T_0 - T_\text{env}\bigr)\,e^{-kt}.\,} $$

有三件事值得记住：

偏离量 $T - T_\text{env}$ 按纯指数衰减。
衰减速度完全由 $k$ 决定。$k$ 减半，冷却时间加倍。
$T \to T_\text{env}$ 是渐近的——理论上永远到不了，但实际上很快物理就不再关心这点差距。

图 1. 牛顿冷却律给出一个解族，每条曲线对应一个初始温度，但所有曲线最终都汇聚到那条虚线（室温）。热咖啡、冷咖啡、冻酸奶——同一物理，同一终点。

2. 方向场：解方程之前就先「看见」解

这里有一个深刻的小技巧，数学界用了 200 年才完全意识到它的价值。对任意一阶 ODE

$$ \frac{dy}{dt} = f(t, y), $$

函数 $f$ 在 $(t, y)$ 平面的每一点都告诉你解曲线在该点的斜率，根本不需要先解出 $y$。在格点上画出对应斜率的小箭头，就得到方向场（也叫斜率场或向量场）。解曲线就是处处与箭头相切的轨迹。

图 2. 衰减原型方程 $dy/dt = -y$ 的方向场。箭头给出局部斜率；紫色曲线是真实的解，沿场线穿行；红色虚线是平衡点 $y = 0$。注意所有解都被「漏斗般」推向零——这个平衡点是一个全局吸引子。

这种可视化让我们能够推理那些根本写不出闭合解的系统。仅凭一张图，你已经可以回答：

有没有平衡点？（哪里的箭头变水平？）
它是吸引还是排斥？（附近的箭头指向它还是远离它？）
解会爆炸、振荡，还是稳定下来？

整个系列里我们都会反复使用方向场——这是理解非线性 ODE 最便宜、最诚实的方式。

3. 分类：ODE、PDE、阶数、线性

几个标签可以为后面省下大量笔墨。

ODE 与 PDE

类型	自变量个数	示例
常微分方程 (ODE)	一个（通常是 $t$）	$\dfrac{dy}{dt} = -ky$
偏微分方程 (PDE)	多个（$x, y, t, \dots$）	$\dfrac{\partial u}{\partial t} = k\,\dfrac{\partial^2 u}{\partial x^2}$

本系列只讨论 ODE。PDE（热方程、波动方程、薛定谔、Navier–Stokes）会在第 13 章另起炉灶。

阶数 = 方程中出现的最高阶导数

一阶：$y' = f(t, y)$——放射性衰变、RC 电路、单种群增长。
二阶：$y'' + p(t)y' + q(t)y = g(t)$——弹簧、单摆、RLC 电路、行星轨道。

阶数为何重要？因为要确定一个唯一解，你必须给出恰好这么多个初始条件。一阶 ODE 需要 $y(0)$；二阶 ODE 同时需要 $y(0)$ 和 $y'(0)$——位置和速度。这不是数学的怪癖，而是 $F = ma$ 在告诉你：要预测一个粒子的未来，必须同时知道它的初始位置和初始速度。

图 3. 左：$y' = -y$ 给出单参数解族——选定起点，曲线就唯一确定。右：$x'' + x = 0$ 给出双参数解族——必须同时给出初始位置与初始速度。阶数 = 宇宙向你索取多少初始信息。

线性与非线性

如果 $y$ 及其各阶导数只以一次幂出现，且彼此不相乘，这个方程就是线性的。$n$ 阶线性方程的一般形式为

$$ a_n(t)\,y^{(n)} + a_{n-1}(t)\,y^{(n-1)} + \dots + a_0(t)\,y = g(t). $$

只要出现 $y^2$、$\sin y$、$y\,y'$、$\sqrt{y}$ 等任意「不老实」的项，方程就是非线性的。

为什么这么在意这个区分？因为线性是一种超能力：

解可以叠加（叠加原理）。
完整解可以分解为齐次解 + 特解。
整套线性代数机器（特征值、积分变换、Green 函数）都能直接套用。

非线性方程没有这些便利。它们能做线性方程做不到的事——极限环、混沌、有限时间爆破——但每一个非线性问题往往都得量身定制。

图 4. 同样大小的两幅图，完全不同的世界。左：线性方程 $y' = -y + \sin t$ 的所有解都被拽到同一条红色周期吸引子上——未来彻底「忘记」初始条件。右：非线性方程 $y' = y^2 - t$ 中，$y(0)$ 的微小变化导致命运迥异，部分解甚至在有限时间内逃逸到 $+\infty$。欢迎来到非线性的世界。

4. 微分方程简史

微分方程并不是「被发现」的，而是被物理「逼出来」的。微积分诞生的那个世纪也正好需要它。

年代	人物	里程碑
1666	牛顿	微积分 + 运动定律。$F = ma$ 本身就是二阶 ODE
1690s	伯努利家族	悬链线、最速降线——把物理问题翻译成 ODE
1739	欧拉	第一个系统化的 ODE 理论；欧拉方法（仍是现代积分器的种子）
1820s	柯西	严格的存在性证明——解到底存不存在？
1890	皮卡、林德勒夫	唯一性定理的现代形式
1880s+	庞加莱	定性理论：放弃追求公式，转而研究解流的几何形状
1963	洛伦兹	混沌从一个三方程的天气玩具模型中诞生

这条弧线令人印象深刻：我们从追逐公式开始，然后放弃，转而学着研究解流的几何。现代 ODE 理论在精神上更接近庞加莱，而不是牛顿。

5. 三个反复出现的经典模型

ODE 之所以无所不在，是因为少数几个方程会以不同的物理面貌反复出现。下面是三个最经典的。

图 5. 全科学中被复用最多的三条方程。每条都只有一行，却共同描述了热交换、考古测年和机械工程中所有弹簧–质量系统。

5.1 放射性衰变——宇宙的秒表

每个放射性原子在单位时间内衰变的概率相同，因此总衰变率与剩余原子数成正比：

$$ \frac{dN}{dt} = -\lambda N \qquad\Longrightarrow\qquad N(t) = N_0\,e^{-\lambda t}. $$

半衰期 $t_{1/2} = (\ln 2)/\lambda$ 是一个能概括整条曲线的单一数字。碳-14 的 $t_{1/2} = 5730\,\text{年}$，正是因此考古学家可以为一块木炭测年：测出残留的 $^{14}\mathrm{C}$ 比例，取对数，就读出年龄。同一条方程也描述了 RC 电路中电容的放电，以及药物在血液中的清除。

5.2 马尔萨斯人口增长

若一个数量为 $P$ 的种群有恒定的人均出生率减死亡率 $r$，则

$$ \frac{dP}{dt} = rP \qquad\Longrightarrow\qquad P(t) = P_0\,e^{rt}. $$

这与衰变方程完全相同，只是把 $-\lambda$ 换成 $r$。数学不在乎一个东西是死还是生，它在乎的是速率常数的符号。

它同时预测无界的指数增长——任何生物学家都会告诉你这是胡说，因为食物、空间或宿主迟早会耗尽。修正方案是 Logistic 方程，它太重要，我们直接把它和它的天真表亲画在一起。

图 6. 左：马尔萨斯的指数模型一路冲向无穷；Logistic 模型则弯下腰，饱和于环境容量 $K$。右：增长率的「相视图」。Logistic 曲线有两个平衡点（$P = 0$ 与 $P = K$），中间是增长区，$K$ 之上是衰减区。我们将在第二章细致剖析它。

5.3 简谐振动

弹簧上的质点 $m$，受 Hooke 力 $F = -kx$ 与牛顿第二定律 $F = ma$ 共同支配：

$$ m\,\frac{d^2 x}{dt^2} = -kx \qquad\Longleftrightarrow\qquad x'' + \omega_0^2\,x = 0, \quad \omega_0 = \sqrt{k/m}. $$

通解为 $x(t) = A\cos(\omega_0 t + \varphi)$——一条纯正弦波，自然角频率为 $\omega_0$。这条方程会在小角度单摆、LC 电路、分子振动，乃至几乎所有线性系统的本征模中反复出现。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

omega0 = 2 * np.pi  # 自然频率：1 Hz
t = np.linspace(0, 3, 500)

def spring(state, t):
    x, v = state
    return [v, -omega0**2 * x]

sol = odeint(spring, [1.0, 0.0], t)  # x(0) = 1, v(0) = 0

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(t, sol[:, 0], color="#10b981", linewidth=2)
ax1.set_xlabel("时间 (s)"); ax1.set_ylabel("位移 x")
ax1.set_title("简谐振动")

ax2.plot(sol[:, 0], sol[:, 1], color="#7c3aed", linewidth=2)
ax2.set_xlabel("x"); ax2.set_ylabel("v = dx/dt")
ax2.set_title("相图——能量守恒给出一个闭合椭圆")
ax2.set_aspect("equal")
plt.tight_layout(); plt.show()

右图是一张相图，第一次让我们看到二阶系统天然生活在 $(x, v)$ 二维空间里。能量守恒迫使轨迹必须是一条闭合椭圆；这个几何事实将成为第 7 章的核心。

6. 初值问题与「解到底存不存在」

一条裸的 ODE 有无穷多个解——一整族（一参数或 $n$ 参数）曲线。要从中挑出唯一的一条，就需要一个初始条件：

$$ \frac{dy}{dt} = f(t, y), \qquad y(t_0) = y_0. $$

二者合在一起就是初值问题（IVP）。随之而来的问题，数学界用了两百年才认真对待：

这个 IVP 真的有解吗？如果有，唯一吗？

现代答案就是 Picard–Lindelöf 定理。

定理（Picard–Lindelöf, 1890）. 若 $f(t, y)$ 在 $(t_0, y_0)$ 附近的一个矩形区域内关于 $t$ 连续、关于 $y$ 满足 Lipschitz 条件，则 IVP 在包含 $t_0$ 的某个开区间上存在唯一解。

「关于 $y$ Lipschitz」是指存在常数 $L$ 使得 $|f(t, y_1) - f(t, y_2)| \le L\,|y_1 - y_2|$——通俗地说，$f$ 沿 $y$ 方向的「斜率」是有界的。大多数「老实」的右端项（多项式、光滑函数）都自动满足。

为什么要在意？因为没有这条保证，模型本身就是有歧义的——这种事确实会发生。

图 7. 左：对 $y' = -y$，$f$ 处处光滑；穿过任意一点恰好有一条解。红色加粗曲线是从 $(0, 1)$ 出发的唯一解。右：对 $y' = 3\,y^{2/3}$，$f$ 连续但在 $y = 0$ 处不 Lipschitz（斜率无穷大）。从原点出发，无穷多条解都合法：$y(t) = 0$ 是一条，$y(t) = (t-c)^3$（任意 $c \ge 0$）也是一条，方程自己根本无法决定。物理上：这意味着模型坏了。

结论很尖锐：一条以 ODE 形式写下的「物理定律」，只有当 Picard–Lindelöf 成立时才真的是一条定律。 否则，未来无法由现在唯一决定，决定论模型的全部意义就崩塌了。

7. 解析解 vs 数值解

	解析解	数值解
得到的是什么	公式：$y(t) = C e^{-kt}$	离散表 $(t_n, y_n)$
从哪里来	笔、纸，闭合形式的积分	时间步进：$y_{n+1} = y_n + h\,f(t_n, y_n) + \dots$
优点	精确；揭示结构（谁依赖谁）	适用于任何写得出来的 ODE
缺点	大多数 ODE 没有闭合解	只是近似；需要细心控制步长

在物理研究生课程上你或许会手算十几条 ODE。但在工业工程与科研里，绝大多数 ODE 都是数值求解的——这并不丢人，事实上现代建模几乎都依赖数值方法。

最简单的数值方法欧拉法，本质上就是把导数定义离散化：

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np
import matplotlib.pyplot as plt

def euler(f, y0, t):
    """y' = f(y, t) 的前向欧拉法"""
    y = np.zeros(len(t)); y[0] = y0
    for i in range(len(t) - 1):
        h = t[i+1] - t[i]
        y[i+1] = y[i] + h * f(y[i], t[i])
    return y

t = np.linspace(0, 5, 100)
y_exact = np.exp(-t)
y_euler = euler(lambda y, t: -y, 1.0, t)

plt.figure(figsize=(8, 5))
plt.plot(t, y_exact, color="#2563eb", linewidth=2, label=r"精确解：$y = e^{-t}$")
plt.plot(t, y_euler, color="#ef4444", linewidth=2, linestyle="--", label="欧拉近似")
plt.xlabel("t"); plt.ylabel("y(t)")
plt.title("解析解 vs 数值解")
plt.legend(); plt.tight_layout(); plt.show()

欧拉法是所有更高级方法的种子（Runge–Kutta、Adams–Bashforth、用于轨道传播的辛积分器）。第 11 章我们会拆解 scipy.integrate.solve_ivp 里那些现代求解器到底在做什么——但底层逻辑永远一样：估计导数、走一步、重复、控制误差。

总结

概念	核心要点
微分方程	把一个未知函数与它的导数关联起来的方程
ODE vs PDE	一个自变量 vs 多个自变量
阶数	出现的最高阶导数——也等于你要给的初始条件个数
线性 vs 非线性	线性可叠加、有完整理论；非线性可以混沌
方向场	解之前先「看见」流的方向
IVP + Picard–Lindelöf	初始条件 + Lipschitz 的 $f$ ⇒ 唯一解
解析解 vs 数值解	能解析则解析；不能则数值（绝大多数情况）

本章一句话：物理学之所以把定律写成微分方程，是因为大自然说的就是「变化率」这门方言；学 ODE，就是学着把这门方言翻回成图像、公式与预测。

练习题

热身

验证 $y(t) = C e^{3t}$（任意常数 $C$）是 $y' = 3y$ 的解。
在 $20^\circ\mathrm{C}$ 的房间里，咖啡 10 分钟内从 $90^\circ\mathrm{C}$ 冷却到 $60^\circ\mathrm{C}$。求冷却常数 $k$，并预测 $t = 30$ min 时的温度。
某放射性样品半衰期 10 年。100 g 样品 25 年后剩多少？

概念

证明：无论 $y(0)$ 为何，$y' = -2y$ 的所有解都满足 $y \to 0$（当 $t \to \infty$）。
手画 $y' = y(1 - y)$ 的方向场。找出所有平衡点；预测 $y(0) = 0.1$、$y(0) = 0.5$、$y(0) = 1.5$ 时的长期行为。（这就是 Logistic 方程。）
为什么马尔萨斯人口模型不现实？至少提出两种应该被加进模型的物理机制。
证明 IVP $y' = 3 y^{2/3}$、$y(0) = 0$ 有不止一个解，并指出违反了 Picard–Lindelöf 的哪条假设。

编程

实现欧拉法，用步长 $h = 0.5,\,0.1,\,0.01$ 求解 $y' = -y$、$y(0) = 1$ 在 $[0, 5]$ 上的解。在 log–log 坐标上画出全局误差与 $h$ 的关系。斜率是多少？为什么？
画出 $y' = \sin t - y$ 的方向场和若干解曲线。猜测它的长期行为，并用数值结果验证。
重现图 6（Logistic vs 指数）：用 scipy.integrate.solve_ivp 同时求解两条 ODE 并叠图。Logistic 解在何时达到 $K$ 的 99%？

参考资料

Boyce, DiPrima & Meade, Elementary Differential Equations and Boundary Value Problems, Wiley（第 11 版, 2017）。
Strogatz, Nonlinear Dynamics and Chaos, CRC Press（第 2 版, 2015）——第 1–2 章是本章绝佳的伴读。
Tenenbaum & Pollard, Ordinary Differential Equations, Dover（1985）——经典的详尽参考。
Hirsch, Smale & Devaney, Differential Equations, Dynamical Systems, and an Introduction to Chaos, Academic Press（第 3 版, 2012）。
MIT OCW 18.03SC——Differential Equations（免费视频公开课）。

系列导航


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