程序设计与计算思维¶

Computer Programming and Computational Thinking

第 17 讲：时间步进法¶

2025—2026学年度春季学期

清华大学 韩文弢

问题引入¶

前面学习了离散与连续的关系：用离散的近似去逼近连续的量。

现在考虑一个新问题：一个系统随时间变化，如何用数学描述它的演化？

例如：

  1. 一批灯泡投入使用后，每天会有一定比例发生故障，如何预测未来还有多少正常工作？
  2. 一场疫情在人群中传播，感染人数先升后降，如何建模并预测走势？

这些问题都可以归结为：已知当前状态和变化规律，求系统随时间的演化——这正是常微分方程（ODE）所描述的问题。

本讲将从简单的故障模型出发，展示如何从离散递推关系出发，推导出连续的 ODE，并介绍求解 ODE 的基本数值方法——欧拉法。

准备工作¶

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from ipywidgets import interact, IntSlider

plt.rcParams['font.family'] = 'Noto Sans CJK SC'

元件故障建模¶

思考一个简单的元件故障模型，例如灯泡(实际上是从前几堂课的微观随机模型推导出这个宏观模型的)。

元件可能在任何时刻发生故障，通过每天检查一次统计当天发生故障的数量，之后再每天检查多次，这仍然是一个离散模型，最后将其转化为连续模型，在该模型中可以看到任何实际时间 $t$ 时已发生故障的数量。

每天检查一次故障（整数时间步）¶

用 $N_k$ 表示第 $k$ 天仍正常工作的灯泡（平均）数量，从初始数量 $N_0$ 开始。通过计算第 $k+1$ 天有多少个故障来求出第 $k+1$ 天结束时仍在工作数量 $N_{k+1}$ 的方程。

设 $p$ 为每个灯泡每天发生故障的概率。例如，如果 10% 的灯泡每天发生故障，则 $p = 0.1$。如果有 100 个灯泡，第一天有 10% 发生故障，那么会有 10 个灯泡故障，剩下 90 个正常。一般来说，如果 $N_k$ 中有比例为 $p$ 发生故障，预计总共有 $p \, N_k$ 个会发生故障。

因此

$$N_{k+1} = N_k - p \, N_k$$

或者

$$N_{k+1} - N_k = - p N_k.$$

在这个非常简单的模型中，实际上可以解析求解这个递推关系，找到时间 $t$ 时仍在正常工作的数量：

$$N_{k+1} = (1 - p) N_k$$

所以

$$N_k = (1 - p) N_{k-1} = (1 - p)^2 N_{k-2} = \cdots$$

因此

$$N_k = N_0 \, (1 - p)^k.$$

每天检查故障 n 次¶

现在假设想知道灯泡在半天内有多少故障。如果每天故障 10%，自然会认为半天内故障 5%。

然而，由于复利效应（如同复利利息），这并不完全正确。如果 5% 故障，然后剩余的 5% 故障，那么还剩多少比例没故障？

In [2]:
(1 - 0.05) * (1 - 0.05)
Out[2]:
0.9025

由于复利效应，总故障的比例略低于10%。与实际情况比较，结果大致正确，所以接受这个近似。

因此，假设每天检查 $n$ 次故障，每次故障的比例为 $p / n$（例如，每天 10% 且每天两次检查，则每次约 5% 故障）。

那么在第 $k$ 天第一次检查后剩余的数量为

$$N_{k + \frac{1}{n}} = \textstyle (1 - \frac{p}{n}) N_k$$

这里使用到下标是因为故障次数处于离散状态，故障次数$N$也可以写成

$$N(k + \textstyle \frac{1}{n}).$$

第二天的解是

$$N_{k+1} = N_k \, (1 - \textstyle \frac{p}{n})^n.$$

$k$ 天后剩余数量的完整解为

$$N_k = N_0 \, \textstyle (1 - \frac{p}{n})^{nk}$$

将这些公式描绘的图形绘制出得到：

In [3]:
def plot_decay(n):
    p = 0.4
    N0 = 100
    T = 20
    
    # 每天检查一次 (基准)
    t_daily = np.arange(0, T + 1)
    N_daily = N0 * (1 - p) ** t_daily
    
    # 每天检查多次
    steps_per_day = 2 ** n
    t_frequent = np.arange(0, T + 1/steps_per_day, 1/steps_per_day)
    N_frequent = N0 * (1 - p/steps_per_day) ** (t_frequent * steps_per_day)
    
    plt.figure(figsize=(10, 6))
    plt.plot(t_daily, N_daily, 'o-', alpha=0.5, ms=3, label='每天一次', lw=2)
    plt.plot(t_frequent, N_frequent, 'o-', alpha=0.5, ms=2, label=f'每天{steps_per_day}次')
    
    plt.xlabel('天数 (k)')
    plt.ylabel('$N_k$')
    plt.title(f'每天{steps_per_day}次')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

# 创建 n 的滑动条，范围 0 到 8，步进 1
interact(plot_decay, n=IntSlider(min=0, max=8, step=1, value=1));

interactive(children=(IntSlider(value=1, description='n', max=8), Output()), _dom_classes=('widget-interact',)…

连续时间¶

回顾"从离散到连续"课程，越来越多的离散值被生成来跟踪，随着离散值越来越多，它们形成的曲线实际上没有明显的改变。

因此，定义一个极限对象是有意义的，也就是如果你能以正确的方式取越来越小的时间步，即取 $n \to \infty$ 的极限，会发生什么。

在该情况下，可以计算 $t$ 时刻正常工作的（平均）灯泡数量 $N(t)$，其中 $t$ 可以是任意正实数。

在微积分中，$(1 - p/n)^n$ 在 $n \to \infty$ 时收敛到 $\exp(-p)$。

另一种方法是根据差分来观察时间演化，在 $1/n$ 的时间内按比例 $p/n$ 衰减。这里将 $\delta t = 1/n$ 视为相邻连续检查之间的时间步长。如果在时间 $t$ 有 $N(t)$ 个正常工作的灯泡，并取一个长度为 $\delta t$ 的小时间步，那么应该有 $p \, \delta t$ 发生故障，因此

$$N(t + \delta t) - N(t) = -(p \, \delta t) \, N(t).$$

两边除以 $\delta t$ 得

$$\frac{N(t + \delta t) - N(t)}{\delta t} = -p \, N(t).$$

方程的左边取 $\delta t \to 0$ 的极限就是导数 $\frac{dN(t)}{dt}$ 的定义，于是得到

$$\frac{dN(t)}{dt} = - p \, N(t)$$

其中 $N(0) = N_0$ 是初始正常工作的灯泡数量。

这是一个常微分方程 (ordinary differential equations, ODE)，它是一个将 $t$ 时刻函数 $N(t)$ 的值与函数在该点的导数（斜率）联系起来的方程，这个关系必须对所有 $t$ 成立。

虽然这个方程是否有意义并不显而易见（尽管根据推导方式应该是合理的），但在微积分课程中可看出它确实有意义（在一些技术条件下），并且与初始条件一起唯一地确定了一个满足 ODE 的函数。这称为初值问题 (initial value problem, IVP)。

$N(t)$的导数是$N(t)$的倍数，因此它必须是指数函数：

$$N(t) = N_0 \exp(-p \, t)$$

这是指数函数的另一种定义方式，将其绘至图中：

In [4]:
def plot_decay1(n):
    # 物理参数
    p = 0.4
    N0 = 100
    T = 20
    
    # 每天检查一次 (基准)
    t_daily = np.arange(0, T + 1)
    N_daily = N0 * (1 - p) ** t_daily
    
    # 每天检查多次
    steps_per_day = 2 ** n
    t_frequent = np.arange(0, T + 1/steps_per_day, 1/steps_per_day)
    N_frequent = N0 * (1 - p/steps_per_day) ** (t_frequent * steps_per_day)
    
    # 连续解
    t_continuous = np.linspace(0, T, 500)
    N_continuous = N0 * np.exp(-p * t_continuous)
    
    # 开始绘图
    plt.figure(figsize=(10, 6))
    plt.plot(t_daily, N_daily, 'o-', alpha=0.5, ms=3, label='每天一次', lw=2)
    plt.plot(t_frequent, N_frequent, 'o-', alpha=0.5, ms=2, label=f'每天{steps_per_day}次')
    plt.plot(t_continuous, N_continuous, label='连续', lw=2, color='black', linestyle='--')
    
    plt.xlabel('天数(k)')
    plt.ylabel('$N_k$')
    plt.title(f'每天{steps_per_day}次')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

# 创建交互式滑动条
interact(plot_decay1, 
         n=IntSlider(min=0, max=8, step=1, value=0));

interactive(children=(IntSlider(value=0, description='n', max=8), Output()), _dom_classes=('widget-interact',)…

从图上看出，这确实是正确的极限逼近曲线。在这种情况下，$p$ 被称为速率；它是单位时间的概率，即概率与时间的比值（要获得在时间 $\delta t$ 内衰减的概率，需要将 $p$ 乘以 $\delta t$）。

总结一下发现的规律：

步的类型 时间步进 差分 解

| 整数 | $N_{k+1} = (1 - p) N_k$ | $N_{k+1} - N_k = -p N_k$ | $N_k = N_0 (1 - p)^k$ | 有理数 | $N_{k + \frac{1}{n}} = \textstyle (1 - \frac{p}{n}) N_k$ | $N_{k + \frac{1}{n}} - N_k = \textstyle (- \frac{p}{n}) N_k$ | $N_k = N_0 (1 - \frac{p}{n})^{n k}$ | 连续 | $N(t + \delta t) = (1 - p \, \delta t) N(t)$ | $\frac{dN(t)}{dt} = - p \, N(t)$ | $N(t) = N_0 \exp(-p \, t)$

SIR 模型¶

看一个更复杂的例子：用于模拟疫情或谣言在人群中传播的SIR 模型。

与作业中一样，可以建立一个完全离散的、基于代理的随机模型，在其中给出微观规则来说明单个代理如何相互交互。当用足够大的系统运行这些模型时，通常会观察到结果相当平滑。另一种方法是尝试写出描述平均值动力学的宏观离散方程。

通常，通过建立模型的连续版本来理解这类系统的行为会更容易。有些人不喜欢连续模型所以做离散模型；大型离散模型也可能在计算上很浪费。另一方面，它们也可能包含在连续框架中更难建模的效果，例如非局部效应或"不容易变成连续"的例子。

离散时间 SIR 模型¶

首先考虑 SI 模型：代理可以是易感者 (S) 和感染者 (I)。当易感者与感染者接触时，有一定概率变成感染者。

设 $S_t$ 和 $I_t$ 分别为时间 $t$ 时易感者和感染者的数量，设 $N$ 为总人数。

假设在每个时间步，每个感染者有机会与另一个人（平均）进行交互。该人从总人口中随机选择。只有当被选中的人是易感者时（这发生的概率为 $S_t / N$），并且感染尝试（概率为 $b$）成功时，才会发生新的感染。

因此，该时间步后感染者的数量变化为

$$\Delta I_t = I_{t+1} - I_t = b \, I_t \, \left(\frac{S_t}{N} \right)$$

$S_t$ 的减少量也由 $\Delta I_t$ 给出。

此外，还有恢复的可能性，一旦你被感染，每个时间步都有恒定的恢复概率 $c$。

按 $N$ 归一化很有用，因此将易感者、感染者和恢复者的比例定义为

$$s_t := \frac{S_t}{N}; \quad i_t := \frac{I_t}{N}; \quad r_t := \frac{R_t}{N}.$$

考虑概率为 $c$ 的恢复率，可以得到离散时间 SIR 模型：

$$\begin{align} s_{t+1} &= s_t - b \, s_t \, i_t \\ i_{t+1} &= i_t + b \, s_t \, i_t - c \, i_t\\ r_{t+1} &= r_t + c \, i_t \end{align}$$

下面是离散时间模型的一个示例模拟。

在编写模拟代码之前，先介绍三个新的 Python 语法知识。

Python 语法：assert 断言¶

assert 语句用于在代码中设置检查点，确保某个条件为真，否则立即报错：

assert 条件, "错误信息"

如果条件为 False，程序会抛出 AssertionError 并显示错误信息。这在验证函数输入、调试时非常有用。

Python 语法：yield 生成器¶

前面已经用过生成器表达式（如 (x**2 for x in range(10))）。现在介绍如何用 yield 关键字编写生成器函数。

普通函数用 return 返回一个值后结束。而用 yield 的函数不会结束，而是暂停，等待下一次调用时从暂停处继续执行：

In [5]:
def count_up(n):
    i = 0
    while i < n:
        yield i    # 产出 i，然后暂停
        i += 1     # 下次调用时从这里继续

调用 count_up(5) 不会立即执行函数体，而是返回一个生成器对象，可以逐个获取值：

In [6]:
count_up(5)
Out[6]:
<generator object count_up at 0x10c7a2260>
In [7]:
for x in count_up(5):   # 每次循环触发一次 yield
    print(x)            # 输出 0, 1, 2, 3, 4

0
1
2
3
4

也可以用 list() 一次性收集所有值：

In [8]:
list(count_up(5))
Out[8]:
[0, 1, 2, 3, 4]

生成器的优势在于惰性求值——它不会一次性生成所有结果，而是按需产出，节省内存。在时间步进模拟中，这与"逐步推进时间"的思想天然契合。

Python 语法：break、continue 与 for...else¶

break：立即终止整个循环，跳出循环体。

continue：跳过当前迭代的剩余语句，直接进入下一次迭代。

for...else：else 子句在循环正常结束（即没有被 break 打断）时执行。这常用于搜索场景：遍历完所有元素都没找到目标时执行 else 分支。

for x in iterable:
    if 满足退出条件:
        break          # 跳出整个循环
    if 需要跳过:
        continue       # 跳过本次，进入下一次
else:
    # 循环没有被 break，正常结束才执行
    ...

下面的例子用 for...else 在一段模拟数据中查找感染者比例的峰值时刻，用 continue 跳过异常数据点。

In [9]:
# 示例：在模拟数据中查找感染者比例峰值

data = [0.01, 0.05, 0.15, None, 0.40, 0.35, 0.20, 0.10, 0.03, 0.01]
threshold = 0.5  # 希望找到感染比例超过 threshold 的时刻

for t, val in enumerate(data):
    if val is None:          # 数据缺失
        print(f'  t={t}: 数据缺失，跳过')
        continue              # 跳过本次迭代
    if val > threshold:
        print(f'  t={t}: 感染比例 {val} 超过阈值 {threshold}，停止搜索')
        break                 # 找到了，立即终止
    print(f'  t={t}: 感染比例 {val:.2f}，继续搜索...')
else:
    # 没有被 break，说明没有找到超过阈值的时刻
    print(f'  未找到感染比例超过 {threshold} 的时刻')

print('\n--- 换一个阈值再试 ---\n')

threshold = 0.3
for t, val in enumerate(data):
    if val is None:
        continue
    if val > threshold:
        print(f'  t={t}: 感染比例 {val} 超过阈值 {threshold}，停止搜索')
        break
    print(f'  t={t}: 感染比例 {val:.2f}，继续搜索...')
else:
    print(f'  未找到感染比例超过 {threshold} 的时刻')

  t=0: 感染比例 0.01，继续搜索...
  t=1: 感染比例 0.05，继续搜索...
  t=2: 感染比例 0.15，继续搜索...
  t=3: 数据缺失，跳过
  t=4: 感染比例 0.40，继续搜索...
  t=5: 感染比例 0.35，继续搜索...
  t=6: 感染比例 0.20，继续搜索...
  t=7: 感染比例 0.10，继续搜索...
  t=8: 感染比例 0.03，继续搜索...
  t=9: 感染比例 0.01，继续搜索...
  未找到感染比例超过 0.5 的时刻

--- 换一个阈值再试 ---

  t=0: 感染比例 0.01，继续搜索...
  t=1: 感染比例 0.05，继续搜索...
  t=2: 感染比例 0.15，继续搜索...
  t=4: 感染比例 0.4 超过阈值 0.3，停止搜索
In [10]:
# 初始化
NN = 100  # 总人口
SS = NN - 1  # 初始易感者人数
II = 1  # 初始感染者人数
RR = 0  # 初始恢复者人数

# 归一化
ss, ii, rr = SS/NN, II/NN, RR/NN

# 感染和恢复概率
p_infection = 0.1
p_recovery = 0.01

# 时间步数
TT = 1000

def discrete_SIR(s0, i0, r0, T=1000, tol=1e-3):
    """离散时间 SIR 模型（生成器版本）
    
    使用 yield 逐步产出每个时间步的状态 (s, i, r)。
    当感染者比例低于 tol 时，用 break 提前终止。
    """
    assert 0 <= s0 <= 1 and 0 <= i0 <= 1 and 0 <= r0 <= 1, "比例必须在 [0, 1] 范围内"
    assert abs(s0 + i0 + r0 - 1) < 1e-10, "s + i + r 必须等于 1"
    
    s, i, r = s0, i0, r0
    yield (s, i, r)  # 产出初始状态
    
    for t in range(T):
        delta_i = p_infection * s * i
        delta_r = p_recovery * i
        
        s = s - delta_i
        i = i + delta_i - delta_r
        r = r + delta_r
        
        yield (s, i, r)  # 产出当前状态，暂停等待下一次
        
        if i < tol:  # 感染者比例极低，疫情结束
            break

# 用 list() 收集生成器的所有输出
SIR = list(discrete_SIR(ss, ii, rr))

# 提取数据
ts = range(len(SIR))
S_values = [x[0] for x in SIR]
I_values = [x[1] for x in SIR]
R_values = [x[2] for x in SIR]

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(ts, S_values, 'o-', label='S', alpha=0.2, markersize=2)
plt.plot(ts, I_values, 'o-', label='I', alpha=0.2, markersize=2)
plt.plot(ts, R_values, 'o-', label='R', alpha=0.2, markersize=2)

plt.xlabel('时间')
plt.ylabel('人口比例')
plt.title('离散时间SIR模型')
plt.legend(loc='right')
plt.grid(True, alpha=0.3)
plt.show()
No description has been provided for this image

连续时间 SIR 模型¶

现在可以像故障模型一样处理相同的过程，这次取长度为 $\delta t$ 的时间步，并将概率 $b$ 和 $c$ 替换为 $\beta$ 和 $\gamma$。取 $\delta t \to 0$ 的极限可以得到

$$\begin{align} \frac{ds(t)}{dt} &= -\beta \, s(t) \, i(t) \\ \frac{di(t)}{dt} &= +\beta \, s(t) \, i(t) &- \gamma \, i(t)\\ \frac{dr(t)}{dt} &= &+ \gamma \, i(t) \end{align}$$

可将上述方程视为具有 S、I 和 R 的化学反应模型。项 $s(t) i(t)$ 被称为质量作用形式的相互作用。

注意，这些（简单的）非线性 ODE 没有已知的解析解作为时间的函数（然而，参数解是已知的）。

Python 语法：*args 可变参数¶

在下面的代码中，solve_ivp 函数的调用方式是：

solve_ivp(sir_model, [0, 100], y0, args=(beta, gamma), t_eval=t)

其中 sir_model 是我们自己定义的导数函数，签名为 sir_model(t, y, beta, gamma)。solve_ivp 内部会自动调用它，但 solve_ivp 本身怎么知道要传 beta 和 gamma 呢？

答案在 args=(beta, gamma) 这个参数。solve_ivp 内部会把这个元组解包，等价于：

sir_model(t, y, *args)  # 等价于 sir_model(t, y, beta, gamma)

这里的 *args 就是 Python 的可变参数机制：在函数定义中，*args 会把多余的位置参数收集成一个元组；在函数调用中，* 会把一个元组（或列表）解包成独立参数。

类似地，**kwargs 用于关键字参数的收集与解包，将多余的关键字参数收集成字典。

def f(*args, **kwargs):
    print('位置参数:', args)   # 元组
    print('关键字参数:', kwargs)  # 字典

f(1, 2, 3, x=10, y=20)
# 位置参数: (1, 2, 3)
# 关键字参数: {'x': 10, 'y': 20}

下面是连续时间模型的一个示例模拟。

In [11]:
# 使用 scipy.integrate 进行更精确的 ODE 求解

# SIR 模型（连续时间）
# 注意：solve_ivp 要求导数函数的签名为 f(t, y, ...)，即 t 在前
def sir_model(t, y, beta, gamma):
    """SIR 模型的导数函数"""
    s, i, r = y
    dsdt = -beta * s * i
    didt = beta * s * i - gamma * i
    drdt = gamma * i
    return [dsdt, didt, drdt]

# 参数
beta = 0.4  # 感染率
gamma = 0.1  # 恢复率

# 初始条件
s0 = 0.99
i0 = 0.01
r0 = 0.0
y0 = [s0, i0, r0]

# 求解 ODE
t_span = (0, 100)  # 时间范围
t_eval = np.linspace(0, 100, 500)  # 需要输出解的时间点
sol = solve_ivp(sir_model, t_span, y0, args=(beta, gamma), t_eval=t_eval)

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='S')
plt.plot(sol.t, sol.y[1], label='I')
plt.plot(sol.t, sol.y[2], label='R')

plt.xlabel('时间')
plt.ylabel('人口比例')
plt.title('连续时间SIR模型 (solve_ivp)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
No description has been provided for this image

时间步进：欧拉法¶

上面展示了如何将常微分方程 (ODE) 视为在取时间步的离散时间模型中自然产生的。

如果有人给出一个 ODE ，应该如何数值求解它？

实际上，做的是相反的过程：离散化方程，将其简化为取时间步的系统！

假设微分方程是

$$\dot{x} = f(x)$$

最简单的方法是欧拉方法：使用一个小的（但不要太小）时间步 $h$（也就是上面所说的 $\delta t$）来近似导数：

$$\dot{x} \simeq \frac{x(t + h) - x(t)}{h}$$

得到

$$x(t + \delta t) \simeq x(t) + h \, f(x).$$

具有恒定时间步的欧拉方法就是以下算法：

$$x_{k+1} = x_k + h \, f(x_k)$$

如果ODE中有多个变量，可以将变量打包成向量并使用相同的方法。

对于 ODE

$$\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x})$$

欧拉方法变为

$$\mathbf{x}_{k+1} = \mathbf{x}_k + h \, \mathbf{f}(\mathbf{x}_k),$$

其中 $\mathbf{f}$ 表示将变量向量映射到 ODE 右侧向量的向量值函数。

然而，总的来说，欧拉方法并不是模拟 ODE 动态的好算法！从本节开头的图表可以看到原因：这样取时间步实际上对近似连续曲线做得很差。数值分析课程展示了如何设计更好的数值方法来更准确地近似 ODE 的真实解。

Python 语法：zip() 并行迭代¶

在下面的代码中，需要同时遍历两个子图 axes 和两个不同的步长 [1, 0.1]。Python 提供了 zip() 函数来实现并行迭代：

for ax, h in zip(axes, [1, 0.1]):
    ...

zip() 将多个可迭代对象"拉链式"地配对：第一个元素配第一个元素，第二个配第二个，以此类推。当最短的那个迭代完时停止。

list(zip(['a', 'b', 'c'], [1, 2, 3]))
# [('a', 1), ('b', 2), ('c', 3)]

之前学过的 enumerate(iterable) 等价于 zip(range(len(iterable)), iterable)。zip 更灵活，可以同时遍历任意多个序列。

In [12]:
# 欧拉方法示例：求解 dN/dt = -p * N
def euler_method(f, y0, t):
    """欧拉方法求解 ODE
    f: 导数函数 dy/dt = f(y, t)
    y0: 初始条件
    t: 时间数组"""
    n = len(t)
    y = np.zeros(n)
    y[0] = y0
    
    for i in range(n - 1):
        dt = t[i+1] - t[i]
        y[i+1] = y[i] + dt * f(y[i], t[i])
    
    return y

# 参数
p = 0.4
N0 = 100
T = 20

# 定义导数函数：dN/dt = -p * N
def dN_dt(N, t):
    return -p * N

# 使用不同的时间步长
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for ax, h in zip(axes, [1, 0.1]):
    t = np.arange(0, T + h, h)
    N_euler = euler_method(dN_dt, N0, t)
    t_continuous = np.linspace(0, T, 500)
    N_exact = N0 * np.exp(-p * t_continuous)
    
    ax.plot(t, N_euler, 'o-', label=f'欧拉方法 (h={h})', alpha=0.7)
    ax.plot(t_continuous, N_exact, label='实际解', lw=2)
    ax.set_xlabel('天数')
    ax.set_ylabel('$N(t)$')
    ax.set_title(f'欧拉方法时间步 h = {h}')
    ax.legend()
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
No description has been provided for this image

可以看到，时间步越小，欧拉方法的近似就越准确。

使用 scipy.integrate.solve_ivp 求解 ODE¶

手动实现欧拉法有助于理解原理，但在实际应用中，应使用经过充分验证的库函数。

SciPy 提供了 solve_ivp（Solve Initial Value Problem），它是求解常微分方程初值问题的推荐接口（替代了旧版的 odeint）。其基本用法为：

from scipy.integrate import solve_ivp

sol = solve_ivp(fun, t_span, y0, ...)

其中：

参数 含义
fun 导数函数，签名为 fun(t, y) 或 fun(t, y, *args)，注意 $t$ 在前
t_span 时间范围元组，如 (t_start, t_end)
y0 初始条件（标量或数组）
args 传给 fun 的额外参数元组
t_eval 希望输出解的时间点数组（可选，不指定则由求解器自动选择）
method 求解方法，默认 'RK45'（Runge-Kutta 4/5 阶），也可选 'Euler'、'RK23' 等

返回值 sol 是一个 OdeResult 对象：

属性 含义
sol.t 求解器实际使用的时间点数组
sol.y 解的数组，形状为 (n_vars, n_times)
sol.success 是否成功求解
sol.message 求解状态的描述信息

对比：欧拉法 vs solve_ivp¶

用同一个故障衰减问题 $\frac{dN}{dt} = -0.4 N$，对比手动欧拉法和 solve_ivp 的结果。

In [13]:
p = 0.4
N0 = 100
T = 20

def dN_dt(t, N):          # solve_ivp 要求 t 在前
    return -p * N

# 手动欧拉法（大步长）
h_euler = 2.0
t_euler = np.arange(0, T + h_euler, h_euler)
N_euler = np.zeros(len(t_euler))
N_euler[0] = N0
for i in range(len(t_euler) - 1):
    N_euler[i+1] = N_euler[i] + h_euler * dN_dt(t_euler[i], N_euler[i])

# solve_ivp（RK45 自适应步长）
t_eval = np.linspace(0, T, 200)
sol = solve_ivp(dN_dt, (0, T), [N0], t_eval=t_eval)

# 精确解
t_exact = np.linspace(0, T, 500)
N_exact = N0 * np.exp(-p * t_exact)

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(t_euler, N_euler, 'rs--', ms=8, label=f'欧拉法 (h={h_euler})')
plt.plot(sol.t, sol.y[0], 'b.-', label='solve_ivp (RK45)')
plt.plot(t_exact, N_exact, 'k-', lw=2, label='精确解')

plt.xlabel('天数')
plt.ylabel('$N(t)$')
plt.title('欧拉法 vs solve_ivp 对比')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
No description has been provided for this image

可以看到，solve_ivp 的默认求解器（RK45，即 4 阶 Runge-Kutta 方法）即使在很少的输出点上也几乎与精确解重合，而欧拉法在大步长下偏差显著。

  • 欧拉法：最简单的 ODE 数值方法，适合教学理解原理，但精度低、稳定性差。
  • solve_ivp：工程实践中应优先使用，内置多种高阶方法（默认 RK45），支持自适应步长，精度和稳定性远优于手动欧拉法。
  • 两者的核心思想相同：都是时间步进——从当前状态出发，根据导数函数推算下一时刻的状态。区别在于 solve_ivp 使用了更精巧的步进公式。

本讲小结¶

  • 从离散到连续：元件故障模型展示了如何从离散递推关系 $N_{k+1} = (1-p)N_k$ 出发，通过不断缩小时间步长，取极限得到连续 ODE $\frac{dN}{dt} = -pN$，其解析解为指数衰减 $N(t) = N_0 e^{-pt}$。
  • SIR 模型：用于描述疫情传播的经典模型，包含易感者 (S)、感染者 (I)、恢复者 (R) 三个群体，可建立离散和连续两个版本。连续 SIR 模型是非线性 ODE 组，一般没有已知的解析解。
  • 欧拉法与 solve_ivp：欧拉法 $x_{k+1} = x_k + h \, f(x_k)$ 是求解 ODE 最基本的数值方法，直观但精度低。实际应用中应使用 scipy.integrate.solve_ivp，它默认采用 RK45 高阶方法，精度和稳定性远优于手动欧拉法。
  • Python 语法：本讲引入了 yield 生成器（惰性产出时间步结果）、break（提前终止循环）、assert（断言验证输入）、zip()（并行迭代）和 *args/**kwargs（可变参数）。