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

你身边的一切都在变化： 咖啡会凉，人口会涨，单摆会摆，病毒会传，股价会震荡，行星会绕行。这些系统都无法用“某物等于多少”来描述，而只能用“某物变化得多快”来刻画。这种描述方式正是微分方程的用武之地——学会读懂它，就等于学会了阅读物理与生物学所使用的语言。

本章将从零开始重建你的直觉。我们从一杯咖啡出发，推导出与放射性衰变、电容放电相同的方程，再逐步引入方向场、分类方法，以及那个决定常微分方程（ODE）是否具有合理解的存在唯一性定理。

总结#

• 微分方程到底是什么——超越符号，理解其本质
• ODE 与 PDE、阶数、线性性——每一种区分都有其实际意义
• 三个经典模型（冷却、衰变、振荡），你将在后续章节中反复遇见
• 方向场：在求解前就能“看见”解的行为
• 存在唯一性定理（Picard–Lindelöf），并附上一个具体反例
• 使用 Python 的 scipy 求解你的第一个数值解，并了解何时可以信任它

前置知识#

• 单变量微积分：导数、积分、链式法则
• Python 基础：NumPy 数组和 Matplotlib 绘图

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

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

1. 现在温度是多少？ —— 这是一个代数问题。
2. 此刻它冷却得多快？ —— 这是一个微积分问题。

这一行公式就是一个微分方程——它将未知函数 $T(t)$ 与其导数联系起来。未知量是函数本身，而非某个具体的数值。

负号至关重要：当 $T > T_\text{env}$ 时，右边为负，温度下降；当 $T < T_\text{env}$ 时，右边为正，温度上升。这个方程本身就蕴含了热力学第二定律。

在求解前先“读懂”模型#

我们甚至不需要写一个积分，就能从中提取出物理预测：

• 温度越接近 $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")

解析解#

有三点值得记住：

• 偏离量 $T - T_\text{env}$ 以纯指数形式衰减。
• 衰减速率完全由 $k$ 决定。若 $k$ 减半，冷却时间就加倍。
• $T$ 渐近趋于 $T_\text{env}$ ——理论上永远不会完全相等，但很快就会近到物理上不再关心。

方向场：在计算前“看见”解#

函数 $f$ 直接告诉你在 $(t, y)$ 平面上每一点处解曲线的斜率，而无需实际求解。在每个网格点画出一个具有该斜率的小箭头，就构成了方向场（也称斜率场或向量场）。真实的解曲线就是那些处处与箭头相切的轨迹。

这种可视化方法让你能分析那些根本无法写出解析解的系统。仅凭方向场，你就能回答：

• 是否存在平衡点？（箭头在何处变为水平？）
• 它是吸引还是排斥？（附近箭头是指向它还是远离它？）
• 解是否会发散、振荡，还是趋于稳定？

在整个系列中，我们会反复依赖方向场——它是理解非线性 ODE 最经济、最诚实的方式。

分类：ODE、PDE、阶数、线性性#

几个简单的标签能帮我们省下大量解释成本。

ODE 与 PDE#

本系列专注于 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$ 的直接体现：要预测一个粒子的未来，你必须知道它的初始位置和初始速度。

线性 vs. 非线性#

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

为何如此强调这一区分？因为线性是一种超能力：

• 解满足叠加原理；
• 通解可分解为齐次解加特解；
• 整套线性代数工具（特征值、变换、格林函数）均可使用。

非线性方程则没有这些便利。它们能表现出线性方程无法实现的行为——极限环、混沌、有限时间内爆破——但每个问题往往都需要专门的处理策略。

简史#

微分方程并非被“发现”的，而是被物理学“逼出来”的。微积分诞生的那个世纪，立刻就需要它。

这条发展脉络令人印象深刻：我们起初追逐解析公式，后来放弃，转而学习解流的几何。现代 ODE 理论在精神上更接近庞加莱，而非牛顿。

三个你会反复遇到的经典模型#

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

放射性衰变——宇宙的计时器#

半衰期 $t_{1/2} = (\ln 2)/\lambda$ 是一个能概括整条衰减曲线的单一数值。碳-14 的半衰期为 $5730\,\text{年}$ ，这正是考古学家能测定木炭年代的原因：测量残留的 $^{14}\mathrm{C}$ 比例，取对数即可推算出年龄。同一方程也适用于 RC 电路中的电容放电，以及药物在血液中的清除过程。

马尔萨斯人口增长#

这与衰变方程完全相同，只是将 $-\lambda$ 替换为 $r$ 。数学并不关心对象是在消亡还是在繁衍，它只关注速率常数的符号。

但它也预言了无界的指数增长——任何生物学家都会指出这是荒谬的，因为食物、空间或宿主终将耗尽。修正方案是著名的逻辑斯蒂方程（logistic equation），它如此重要，我们直接将其与原始模型并列展示。

简谐振动#

通解为 $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 章 的核心。

初值问题与解的存在性#

二者合称为初值问题（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$ 方向上的“斜率”是有界的。大多数“行为良好”的右端函数（如多项式、光滑函数）都自动满足此条件。

为何在意这一点？因为若无此保证，模型本身就存在歧义——而这种情况确实会发生。

关键结论非常明确：只有当 Picard–Lindelöf 条件成立时，以 ODE 形式写出的“物理定律”才算真正成立。 否则，未来无法由当前状态唯一确定，整个决定论模型的基础就会崩塌。

解析解 vs. 数值解#

在研究生物理课程中，你可能会手解十几个 ODE。但在专业工程与科学研究中，绝大多数 ODE 都通过数值方法求解——这毫无羞耻可言；事实上，现代建模几乎完全依赖于此。

 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，就是学会将这种语言翻译成图像、公式和预测。

练习题#

热身

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

概念

4. 证明：无论初始值 $y(0)$ 如何，$y' = -2y$ 的所有解都满足 $\lim_{t \to \infty} y(t) = 0$ 。
5. 手绘 $y' = y(1 - y)$ 的方向场，标出平衡点，并预测 $y(0) = 0.1$ 、$y(0) = 0.5$ 和 $y(0) = 1.5$ 的长期行为。（这是逻辑斯蒂方程。）
6. 为什么马尔萨斯人口模型不现实？至少提出两种应纳入模型的物理机制。
7. 证明初值问题 $y' = 3 y^{2/3}$ 、$y(0) = 0$ 存在多个解，并指出 Picard–Lindelöf 定理中哪条假设被违反。

编程

8. 实现欧拉法，在区间 $[0, 5]$ 上用步长 $h = 0.5,\,0.1,\,0.01$ 求解 $y' = -y$ 、$y(0) = 1$ 。在双对数坐标下绘制全局误差与 $h$ 的关系，观察斜率是多少？为什么？
9. 绘制 $y' = \sin t - y$ 的方向场及若干解曲线，猜测其长期行为，并用数值方法验证。
10. 复现图 6（逻辑斯蒂 vs. 指数增长）：使用 scipy.integrate.solve_ivp 同时求解两个 ODE 并叠加绘图。逻辑斯蒂解何时达到承载力 $K$ 的 99%？

下一步#

下一章我们正式开始动手解 ODE。从最简单的一阶情形入手——分离变量、积分因子、线性齐次方程的标准化套路。这些技术看起来零散，但它们的共同主题是：让方程结构暴露出可积的结构。读完下一章，你会对’怎么解’这件事不再依赖运气，而是看一眼方程就知道应该用哪一招。

参考文献#

• 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（免费视频课程）。