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

系列导读：这是一个为期 7 个月、共 18 章的常微分方程深度系列。我们将从零开始，用最直观的方式理解 ODE，并配以大量可视化和 2020 年最新应用案例（包括 COVID-19 疫情建模）。每个概念都从实际问题出发，用 Python 实现，让数学真正"活"起来！

引言：变化无处不在

如果你观察周围的世界，会发现一切都在变化： - 咖啡在变凉 ☕ - 人口在增长 👶 - 单摆在摆动 ⏰ - 病毒在传播 🦠 - 心脏在跳动 ❤️

微分方程就是描述"变化"的数学语言。它不仅是数学工具，更是理解自然规律的钥匙。

为什么学习微分方程？

2020 年， COVID-19 疫情席卷全球。科学家如何预测疫情发展？如何评估隔离措施的效果？答案就是微分方程！

今天，我们将： 1. 理解什么是微分方程（用最简单的例子） 2. 探索它的历史起源 3. 看 3 个真实应用案例 4. 用 Python 实现第一个微分方程

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

1.1 牛顿冷却定律：一杯咖啡的故事

想象你泡了一杯 90 ° C 的咖啡，放在 20 ° C 的房间里。几分钟后，你回来发现咖啡变凉了。

问题：咖啡的温度是如何随时间变化的？

牛顿的观察（ 1701 年）：物体的冷却速度与它和环境的温度差成正比。

用数学表达：

这就是一个微分方程！

解读每个符号（小白版）：

-：时刻的温度 - 比如（刚泡好） -（10分钟后）

-：温度的变化率 - 通俗理解：温度下降的"速度" - 单位：度/分钟 - 比如：表示每分钟降温5度

-：环境温度（ 20 ° C） - 就是房间的温度 - 咖啡最终会降到这个温度

-：冷却常数 - 与杯子材质有关（保温杯小，纸杯大） - 越大，冷得越快

公式的含义（翻译成人话）：

"咖啡的降温速度，和它比环境温度高多少成正比"

具体例子： - 如果咖啡是90°C，环境20°C，温差70°C → 降温很快 - 如果咖啡是30°C，环境20°C，温差10°C → 降温很慢 - 负号表示温度在"下降"

1.2 用 Python 模拟咖啡冷却

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义微分方程
def coffee_cooling(T, t, k, T_env):
    """
    牛顿冷却定律
    T: 当前温度
    t: 时间
    k: 冷却常数
    T_env: 环境温度
    """
    dTdt = -k * (T - T_env)
    return dTdt

# 参数设置
T0 = 90.0        # 初始温度（° C）
T_env = 20.0     # 环境温度（° C）
k = 0.1          # 冷却常数（ 1/分钟）
t = np.linspace(0, 60, 500)  # 0-60 分钟

# 求解微分方程
T_solution = odeint(coffee_cooling, T0, t, args=(k, T_env))

# 绘图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 左图：温度随时间变化
ax1.plot(t, T_solution, 'b-', linewidth=2.5, label='Coffee Temperature')
ax1.axhline(y=T_env, color='r', linestyle='--', linewidth=2, label=f'Room Temp = {T_env}° C')
ax1.fill_between(t, T_solution.flatten(), T_env, alpha=0.3, color='orange')
ax1.set_xlabel('Time (minutes)', fontsize=12)
ax1.set_ylabel('Temperature (° C)', fontsize=12)
ax1.set_title('Newton\'s Law of Cooling', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 标注关键点
ax1.plot(0, T0, 'go', markersize=10, label='Start')
ax1.text(2, T0-2, f'Initial: {T0}° C', fontsize=10, color='green')

# 右图：冷却速率
dT_dt = -k * (T_solution.flatten() - T_env)
ax2.plot(t, dT_dt, 'r-', linewidth=2.5)
ax2.set_xlabel('Time (minutes)', fontsize=12)
ax2.set_ylabel('Cooling Rate (° C/min)', fontsize=12)
ax2.set_title('Rate of Temperature Change', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)

plt.tight_layout()
plt.savefig('/Users/kchen/Desktop/Blog/source/_posts/常微分方程系列/第 1 章/图 1_咖啡冷却.png', 
            dpi=300, bbox_inches='tight')
plt.close()

print("✓ 图 1 已生成：咖啡冷却")
print(f"初始温度：{T0}° C")
print(f"10 分钟后：{T_solution[np.abs(t-10).argmin()][0]:.2f}° C")
print(f"30 分钟后：{T_solution[np.abs(t-30).argmin()][0]:.2f}° C")

关键洞察： 1. 温度指数衰减地接近环境温度 2. 冷却速率越来越慢（温度差越小，冷却越慢） 3. 理论上永远不会完全达到环境温度（但会非常接近）

1.3 微分方程的精确解

对于这个简单方程，我们能求出精确解：

验证： - 时：✓ - 时：✓

二、微分方程的分类与基本概念

2.1 什么是"常"微分方程？

定义：只含有一个自变量（通常是时间）的微分方程。

对比： - 常微分方程（ ODE）：（一个变量） - 偏微分方程（ PDE）：（多个变量）

2.2 阶数：导数的最高次数

一阶 ODE：最高导数是一阶

例子：人口增长、放射性衰变、 RC 电路

二阶 ODE：最高导数是二阶

例子：弹簧振动、行星运动、 RLC 电路

2.3 线性 vs 非线性

线性 ODE：及其导数都是一次的

非线性 ODE：含有、、等

为什么重要？ - 线性方程：有成熟理论，容易求解 - 非线性方程：更真实，但往往只能数值求解

三、历史之旅：微分方程的诞生

3.1 牛顿与莱布尼茨（ 1600s）

1666 年：牛顿发明微积分，立即用于描述运动。

牛顿第二定律本身就是微分方程：

这是物理学中最重要的微分方程！

3.2 伯努利家族（ 1700s）

雅各布·伯努利和约翰·伯努利兄弟研究了许多经典 ODE： - 悬链线问题（悬挂链条的形状） - 最速降线问题（最快下降路径）

丹尼尔·伯努利（他们的侄子/儿子）： - 研究弹性体振动 - 提出叠加原理

3.3 欧拉（ 1700s）

1739 年：欧拉系统化了 ODE 理论，引入： - 常系数线性方程的通解 - 欧拉方法（第一个数值算法！）

3.4 拉普拉斯（ 1800s）

1785 年：拉普拉斯变换 - 将微分方程转化为代数方程 - 工程师的最爱！

3.5 庞加莱（ 1880s）

革命性转变：不再追求精确解，而是研究解的性质 - 稳定性理论 - 相空间（ Phase Space） - 为现代动力系统奠基

# 可视化：历史时间线
fig, ax = plt.subplots(1, 1, figsize=(14, 6))

# 历史事件
events = [
    (1666, "Newton\nInvents Calculus", 'blue'),
    (1690, "Bernoulli Brothers\nBrachistochrone", 'green'),
    (1739, "Euler\nSystematizes ODEs", 'red'),
    (1785, "Laplace\nTransform", 'purple'),
    (1881, "Poincar é\nQualitative Theory", 'orange'),
    (1963, "Lorenz\nChaos Theory", 'brown'),
    (2020, "COVID-19\nEpidemic Models", 'crimson')
]

years = [e[0] for e in events]
names = [e[1] for e in events]
colors = [e[2] for e in events]

# 绘制时间线
ax.plot([1650, 2030], [0, 0], 'k-', linewidth=2)

for i, (year, name, color) in enumerate(events):
    # 交替上下放置
    y_pos = 0.3 if i % 2 == 0 else -0.3
    
    # 画点
    ax.plot(year, 0, 'o', markersize=15, color=color, zorder=3)
    
    # 画连线
    ax.plot([year, year], [0, y_pos*0.8], color=color, linewidth=2, linestyle='--')
    
    # 标注
    ax.text(year, y_pos, name, ha='center', va='bottom' if y_pos > 0 else 'top',
            fontsize=10, fontweight='bold', color=color,
            bbox=dict(boxstyle='round,pad=0.5', facecolor='white', edgecolor=color, linewidth=2))

ax.set_xlim(1640, 2040)
ax.set_ylim(-0.5, 0.5)
ax.set_xlabel('Year', fontsize=12, fontweight='bold')
ax.set_title('History of Differential Equations: Key Milestones', fontsize=14, fontweight='bold')
ax.axis('off')

plt.tight_layout()
plt.savefig('/Users/kchen/Desktop/Blog/source/_posts/常微分方程系列/第 1 章/图 2_历史时间线.png', 
            dpi=300, bbox_inches='tight')
plt.close()

print("✓ 图 2 已生成：历史时间线")

四、三个经典应用案例

4.1 案例 1：放射性衰变

背景：放射性元素以固定的"半衰期"衰减。

模型：-：时刻的原子数 -：衰变常数

解：

应用：碳-14 定年法（考古学）

# 碳-14 衰变模拟
fig, ax = plt.subplots(1, 1, figsize=(10, 6))

t_half = 5730  # 碳-14 半衰期（年）
lambda_c14 = np.log(2) / t_half
t = np.linspace(0, 30000, 1000)
N = np.exp(-lambda_c14 * t)

ax.plot(t, N, 'b-', linewidth=2.5, label='C-14 Remaining')
ax.axhline(y=0.5, color='r', linestyle='--', linewidth=2, label='50% (1 half-life)')
ax.axhline(y=0.25, color='orange', linestyle='--', linewidth=2, label='25% (2 half-lives)')
ax.axvline(x=t_half, color='r', linestyle=':', alpha=0.5)
ax.axvline(x=2*t_half, color='orange', linestyle=':', alpha=0.5)

ax.set_xlabel('Time (years)', fontsize=12)
ax.set_ylabel('Fraction Remaining', fontsize=12)
ax.set_title('Radioactive Decay: Carbon-14', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

# 标注半衰期
ax.text(t_half, 0.5, f'  {t_half} years', fontsize=10, va='bottom', color='red')

plt.tight_layout()
plt.savefig('/Users/kchen/Desktop/Blog/source/_posts/常微分方程系列/第 1 章/图 3_放射性衰变.png', 
            dpi=300, bbox_inches='tight')
plt.close()

print("✓ 图 3 已生成：放射性衰变")

4.2 案例 2：指数人口增长（马尔萨斯模型）

1798 年：托马斯·马尔萨斯提出最简单的人口模型。

假设：出生率恒定，死亡率恒定-：人口数量 -：净增长率（出生率 - 死亡率）

解：

问题：指数增长是不可持续的！（资源有限）

改进： Logistic 模型（下一章）

# 人口增长对比：不同增长率
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

t = np.linspace(0, 100, 500)
P0 = 1000  # 初始人口

# 左图：不同增长率
growth_rates = [0.01, 0.02, 0.03, 0.04]
colors = ['blue', 'green', 'orange', 'red']

for r, color in zip(growth_rates, colors):
    P = P0 * np.exp(r * t)
    ax1.plot(t, P, color=color, linewidth=2.5, label=f'r = {r} (doubling: {np.log(2)/r:.1f} yrs)')

ax1.set_xlabel('Time (years)', fontsize=12)
ax1.set_ylabel('Population', fontsize=12)
ax1.set_title('Exponential Growth: Effect of Growth Rate', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_yscale('log')

# 右图：半对数坐标
r = 0.03
P = P0 * np.exp(r * t)
ax2.semilogy(t, P, 'b-', linewidth=2.5)
ax2.set_xlabel('Time (years)', fontsize=12)
ax2.set_ylabel('Population (log scale)', fontsize=12)
ax2.set_title('Exponential Growth on Log Scale\n(Becomes a Straight Line!)', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('/Users/kchen/Desktop/Blog/source/_posts/常微分方程系列/第 1 章/图 4_人口增长.png', 
            dpi=300, bbox_inches='tight')
plt.close()

print("✓ 图 4 已生成：人口增长")

4.3 案例 3：简谐振动（弹簧-质量系统）

胡克定律：弹簧的恢复力与位移成正比

牛顿第二定律：

整理得：

其中是角频率。

解：-：振幅 -：初相位

这是简谐运动！

# 简谐振动模拟
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

# 参数
m = 1.0      # 质量（ kg）
k = 10.0     # 弹簧常数（ N/m）
omega = np.sqrt(k / m)
A = 1.0      # 振幅（ m）
phi = 0      # 初相位

t = np.linspace(0, 10, 1000)
x = A * np.cos(omega * t + phi)
v = -A * omega * np.sin(omega * t + phi)

# 上图：位移-时间
ax1.plot(t, x, 'b-', linewidth=2.5, label='Displacement x(t)')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.axhline(y=A, color='r', linestyle='--', linewidth=1.5, alpha=0.5, label=f'Amplitude = {A}m')
ax1.axhline(y=-A, color='r', linestyle='--', linewidth=1.5, alpha=0.5)
ax1.set_xlabel('Time (seconds)', fontsize=12)
ax1.set_ylabel('Displacement (m)', fontsize=12)
ax1.set_title(f'Simple Harmonic Motion: ω = {omega:.2f} rad/s, Period = {2*np.pi/omega:.2f}s', 
             fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 下图：相空间（位移-速度）
ax2.plot(x, v, 'g-', linewidth=2.5)
ax2.plot(x[0], v[0], 'ro', markersize=10, label='Start')
ax2.arrow(x[50], v[50], x[51]-x[50], v[51]-v[50], 
          head_width=0.3, head_length=0.2, fc='blue', ec='blue')
ax2.set_xlabel('Displacement x (m)', fontsize=12)
ax2.set_ylabel('Velocity v (m/s)', fontsize=12)
ax2.set_title('Phase Space Portrait (Circular Orbit)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_aspect('equal')

plt.tight_layout()
plt.savefig('/Users/kchen/Desktop/Blog/source/_posts/常微分方程系列/第 1 章/图 5_简谐振动.png', 
            dpi=300, bbox_inches='tight')
plt.close()

print("✓ 图 5 已生成：简谐振动")

五、微分方程的几何意义：方向场

5.1 什么是方向场？

对于一阶 ODE：

在平面的每一点，方程告诉我们解曲线在该点的斜率。

方向场：在网格点上画出这些斜率（小箭头）。

5.2 例子：

def direction_field_example():
    fig, ax = plt.subplots(1, 1, figsize=(10, 8))
    
    # 定义 ODE
    def f(t, y):
        return y - t
    
    # 创建网格
    t_vals = np.linspace(-2, 4, 20)
    y_vals = np.linspace(-2, 4, 20)
    T, Y = np.meshgrid(t_vals, y_vals)
    
    # 计算每点的斜率
    dY = f(T, Y)
    dT = np.ones_like(dY)
    
    # 归一化箭头长度
    M = np.sqrt(dT**2 + dY**2)
    dT_norm = dT / M
    dY_norm = dY / M
    
    # 绘制方向场
    ax.quiver(T, Y, dT_norm, dY_norm, M, cmap='viridis', alpha=0.6)
    
    # 绘制几条解曲线
    from scipy.integrate import odeint
    
    def ode_func(y, t):
        return y - t
    
    t_span = np.linspace(-2, 4, 500)
    initial_conditions = [-1, 0, 1, 2, 3]
    
    for y0 in initial_conditions:
        sol = odeint(ode_func, y0, t_span)
        ax.plot(t_span, sol, 'r-', linewidth=2, alpha=0.8)
    
    ax.set_xlabel('t', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title('Direction Field: dy/dt = y - t', fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-2, 4)
    ax.set_ylim(-2, 4)
    
    plt.tight_layout()
    plt.savefig('/Users/kchen/Desktop/Blog/source/_posts/常微分方程系列/第 1 章/图 6_方向场.png', 
                dpi=300, bbox_inches='tight')
    plt.close()
    
    print("✓ 图 6 已生成：方向场")

direction_field_example()

关键洞察： - 红色曲线是积分曲线（解） - 它们顺着箭头方向流动 - 不同初值→不同解曲线

六、初值问题（ IVP）

6.1 定义

初值问题：给定 ODE 和初始条件，求解

存在唯一性定理（ Picard-Lindel ö f）：如果满足某些条件（连续且 Lipschitz 连续），则 IVP 有唯一解。

6.2 为什么初值重要？

没有初值，微分方程有无穷多个解（通解）。有了初值，确定唯一一个解（特解）。

物理意义： - 初始位置和初始速度→确定物体的整个运动轨迹 - 初始人口→确定未来人口变化

七、数值解 vs 解析解

7.1 解析解

用公式精确表示： $某个公式$

优点：精确、优美
缺点：大多数 ODE 求不出解析解！

7.2 数值解

用计算机逐步逼近：

优点：适用于任何 ODE
缺点：只是近似、需要计算资源

现实：工程和科学中 99%的 ODE 都是数值求解！

# 对比：解析解 vs 数值解
fig, ax = plt.subplots(1, 1, figsize=(10, 6))

# ODE: dy/dt = -y, y(0) = 1
# 解析解: y(t) = e^(-t)

t = np.linspace(0, 5, 100)
y_analytical = np.exp(-t)

# 数值解：欧拉方法
def euler_method(f, y0, 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

def f(y, t):
    return -y

y_numerical = euler_method(f, 1.0, t)

# 绘图
ax.plot(t, y_analytical, 'b-', linewidth=3, label='Analytical Solution:$y = e^{-t}$')
ax.plot(t, y_numerical, 'r--', linewidth=2, label='Numerical Solution (Euler)', alpha=0.7)
ax.plot(t, np.abs(y_analytical - y_numerical), 'g:', linewidth=2, label='Error', alpha=0.7)

ax.set_xlabel('Time t', fontsize=12)
ax.set_ylabel('y(t)', fontsize=12)
ax.set_title('Analytical vs Numerical Solution', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/Users/kchen/Desktop/Blog/source/_posts/常微分方程系列/第 1 章/图 7_数值 vs 解析.png', 
            dpi=300, bbox_inches='tight')
plt.close()

print("✓ 图 7 已生成：数值 vs 解析解")

八、展望 2020：从牛顿到 COVID-19

8.1 疫情建模预告

2020 年初， COVID-19 爆发。全世界都在问： - 疫情会持续多久？ - 隔离有用吗？ - 何时达到峰值？

答案在SIR 模型（第 13 章重点）：

这是一个微分方程组！

8.2 神经 ODE（ Neural ODEs）

2018 年突破：用神经网络学习微分方程！

传统：设计，求解神经 ODE：让神经网络学习！

应用：时间序列、生成模型、物理建模

九、总结与展望

本章关键要点

微分方程描述变化：温度、人口、振动...
ODE vs PDE：单变量 vs 多变量
阶数：最高导数的次数
线性 vs 非线性：决定求解难度
方向场：可视化解的流动
初值问题：确定唯一解
数值解：现实中的主要工具

为什么重要？

微分方程是： - 物理学的语言（牛顿定律） - 生物学的工具（种群模型） - 工程学的基础（控制系统） - 经济学的模型（增长理论） - 医学的武器（疫情预测）

下一章预告

《一阶微分方程》

我们将学习： - 可分离方程： - 线性方程： - 恰当方程： - 混合问题：两种液体的混合 - RC 电路：电容器充放电

每种方程都有特定技巧！

练习题

基础题

验证是的解（是任意常数）。
咖啡从 90 ° C 冷却到 60 ° C 用了 10 分钟，环境温度 20 ° C 。求冷却常数。
放射性元素的半衰期是 10 年。 100 克样品， 20 年后剩多少？

进阶题

证明：所有满足的解都趋向于 0（当）。
对于，画出方向场并猜测解的行为。
为什么指数人口模型不现实？提出一个改进方案。

编程题

用 Python 实现欧拉方法，求解，，。
制作一个动画 GIF，展示不同初值的解曲线在方向场中的演化。
模拟单摆运动（不考虑阻尼），绘制相空间轨迹。

参考资料

Boyce, W. E., & DiPrima, R. C. (2012). Elementary Differential Equations. Wiley.
Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos. CRC Press.
MIT OpenCourseWare 18.03SC - Differential Equations
Tenenbaum, M., & Pollard, H. (1985). Ordinary Differential Equations. Dover.
Imperial College COVID-19 Response Team (March 2020). Impact of NPIs.

代码资源

本章所有代码和图表： - GitHub: [ode-series]/chapter-1 - Jupyter Notebook: 交互式版本 - 动画 GIF: 方向场演化、简谐振动 ---

本文是《常微分方程的世界》系列的第 1 章，共 18 章。 作者：[Your Name] | 发布日期： 2020-01-01 2020 年特别版：包含 COVID-19 疫情建模应用