前言

0.1 为什么用微分方程建模

微分方程是描述变化率的数学语言。在现实世界中，绝大多数现象的本质不是静态的"是什么"，而是动态的"怎么变"——人口的增长速率、疾病的传播速度、弹簧的振动频率、化学反应的消耗速率……这些都可以用微分方程精确刻画。

数学建模竞赛中，微分方程模型是最核心、最高频的题型之一。掌握微分方程的数值求解能力，意味着你可以将复杂的现实问题转化为可计算的模型。

上图展示了三种典型的微分方程解曲线：

• 左图：$\frac{dy}{dt} = ky,\;k>0$，指数增长模型。增长率与当前值成正比，解为 $y(t) = y_0 e^{kt}$。适用于无限制条件下的种群增长、复利计算等场景。

• 中图：$\frac{dy}{dt} = ky,\;k<0$，指数衰减模型。适用于放射性衰变、冷却过程、药物代谢等。

• 右图：$\frac{d^2y}{dt^2} = -y$，简谐振动方程。解为正弦函数，描述了弹簧振子、单摆小角度摆动等周期运动。

这三种方程虽然简单，但它们是构建复杂模型的基石。后续章节中，你会看到更复杂的模型本质上都是这些基本模式的组合与扩展。

0.2 常用微分方程模型预览

在数学建模中，以下两类模型极为常见：

左图 — Lotka-Volterra 捕食者-猎物模型：

​其中 x 为猎物数量，y 为捕食者数量。猎物有自然增长率 α，被捕食率为 βxy；捕食者因捕食而增长 δxy，自然死亡率为 γy。

右图 — SIR 传染病模型：

​S 为易感者，I 为感染者，R 为康复者。β 为感染率，γ 为康复率。这是新冠疫情期间被广泛讨论的经典模型。

这两类模型都是耦合的常微分方程组，通常没有解析解，必须借助数值方法求解——这正是本教程要解决的问题。

0.3 教程结构

本教程以 SciPy 的 solve_ivp 为核心工具，系统讲解常微分方程的数值求解方法：

0.4 环境配置

依赖安装

pip install numpy scipy matplotlib

版本要求

solve_ivp 在 SciPy 1.0+ 中引入，是官方推荐的 ODE 求解接口（旧的 odeint 已被标记为 legacy）。

可视化中文适配

本教程所有示例均配置了 matplotlib 中文支持：

import matplotlib.pyplot as plt

plt.rcParams.update({
    'font.sans-serif': ['SimHei', 'Microsoft YaHei', 'DejaVu Sans'],
    'axes.unicode_minus': False,
})

• SimHei（黑体）和 Microsoft YaHei（微软雅黑）是 Windows 系统自带的中文字体

• axes.unicode_minus 确保负号正常显示

验证安装

import numpy as np
import scipy
import matplotlib

print(f"NumPy: {np.__version__}")
print(f"SciPy: {scipy.__version__}")
print(f"Matplotlib: {matplotlib.__version__}")

确保运行后版本号满足上述要求，即可开始后续章节的学习。

0.5 前置知识

阅读本教程需要具备：

1. 高等数学基础：导数、积分、泰勒展开的基本概念

2. 线性代数基础：向量、矩阵、特征值

3. Python 编程：熟悉 NumPy 数组操作、基本的数据可视化

4. ODE 基本概念：了解什么是微分方程、初值问题、解析解与数值解的区别

不需要深入的数学分析背景——教程中的所有数学内容都会配合代码和可视化进行直观解释。

0.6 解析解 vs 数值解

现实情况：绝大多数建模问题中的微分方程没有解析解，或者解析解形式过于复杂而不实用。数值解法是唯一可行的途径。

本教程的重点就是：给定一个微分方程（组），如何用 SciPy 高效、准确地求出数值解，并用直观的可视化呈现结果。

第1章 初识 solve_ivp

1.1 为什么选择 solve_ivp

SciPy 1.0 引入了 solve_ivp 作为求解常微分方程初值问题的统一接口。在 scipy 1.0 之前，最常用的函数是 odeint，但它已被标记为 legacy，不再推荐使用。

solve_ivp 的优势：

• 统一的函数签名：所有求解方法共享同一套 API

• 多种求解器：支持显式 Runge-Kutta、BDF、Radau 等不同算法

• 事件检测：可以在求解过程中检测特定事件（如到达某个状态）

• 密集输出：提供任意时间点的连续插值

from scipy.integrate import solve_ivp
import numpy as np

# solve_ivp 的最简调用
# 求解 dy/dt = 0.5*y, y(0) = 1
def exp_growth(t, y, k):
    return k * y

sol = solve_ivp(exp_growth, [0, 5], [1.0], args=(0.5,), dense_output=True)

print(f"时间步: {sol.t}")
print(f"解: {sol.y[0]}")
print(f"状态: {sol.message}")

输出：

时间步: [0.         0.03138145 0.13825591 0.32678526 0.63876289 1.05884731
         1.60909647 2.29634592 3.19885228 4.46787954 5.        ]
解: [ 1.          1.01594513  1.07141356  1.17396838  1.37110218  1.69380744
      2.39031659  3.71374384  6.39638404 12.18494287 14.84131591]
状态: The solver successfully reached the end of the integration interval.

1.2 函数签名详解

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45',
                           t_eval=None, dense_output=False,
                           events=None, args=None,
                           rtol=1e-3, atol=1e-6,
                           max_step=0.0)

返回的 Bunch 对象包含以下重要字段：

sol.t        # 时间点数组（求解器实际选择的时间步）
sol.y        # 解数组，sol.y[i] 对应第 i 个方程在所有时间点的解
sol.t_events # 各事件触发的时间点列表
sol.success  # 布尔值，求解是否成功
sol.message  # 求解状态描述
sol.nfev     # 右端函数的调用次数（性能指标）
sol.sol      # 密集输出插值函数（仅当 dense_output=True）

1.3 指数增长模型

1.3.1 问题描述

最简单的 ODE——指数增长模型：

​其解析解为 $y(t) = y_0 e^{kt}$。当 k>0 时呈指数增长，k<0 时呈指数衰减。

from scipy.integrate import solve_ivp
import numpy as np

def exp_growth(t, y, k):
    return k * y

# 求解
t_span = [0, 5]
y0 = [1.0]
sol = solve_ivp(exp_growth, t_span, y0, args=(0.5,), dense_output=True)

# 与解析解对比
t_fine = np.linspace(0, 5, 200)
y_numerical = sol.sol(t_fine)
y_analytical = np.exp(0.5 * t_fine)

1.3.2 数值精度验证

左图展示了求解器的输出行为：

• 蓝色插值曲线：通过 dense_output=True 获得的连续解，与解析解（黑色虚线）几乎完全重合

• 红色圆点：求解器自适应选择的实际计算点——注意这些点不是均匀分布的，solve_ivp 会根据解的变化率自动调整步长。在增长较快的区域（tt 较大时）步长更大，在增长较慢的区域（tt 较小时）步长更密集

• 总共只用了 11 个计算点就达到了高精度，这就是自适应步长的威力

右图展示了数值解与解析解的绝对误差（对数刻度）：

• 误差普遍在 10⁻⁶量级以下，远低于默认的相对容差 rtol=1e-3

• 误差曲线上的周期性尖峰反映了自适应步长算法的步长调整过程——每当求解器估计误差超过阈值时，就会减小步长重新计算，导致误差骤降

• 误差在 t 较大时略有增长，这是因为指数增长放大了累积误差

1.3.3 放射性衰变——指数衰减的应用

左图展示了放射性衰变模型 dN/dt=−λN,N(0)=100，其中 λ=0.3：

• 蓝色曲线（数值解）与黑色虚线（解析解）完全重合，肉眼无法区分

• 红色虚线标记了半衰期 t1/2​=ln2/λ≈2.3——此时原子核数量恰好减少一半

• 经过约 4 个半衰期（t≈9.2），剩余量已不足初始值的 7%；经过 10 个半衰期后，几乎完全衰变殆尽

右图是相线（Phase line）——一种一维的相图，横坐标为 dy/dt，纵坐标为 y：

• 直线斜率为负，表示 y 越大，下降越快

• 当 y→0 时，dy/dt→0——零点是系统的平衡点（不动点），一旦到达就不会再变化

• 相线分析是理解微分方程定性行为的有力工具，后续章节中会扩展到二维相图

1.4 不同求解器方法的对比

solve_ivp 内置了五种求解方法：

1.4.1 Logistic 增长模型

用 logistic 增长模型来对比各方法的精度：

其中 r=2.0 为增长率，K=10.0 为环境承载力。其解析解为：

from scipy.integrate import solve_ivp
import numpy as np

def logistic(t, y, r, K):
    return r * y * (1 - y / K)

r, K = 2.0, 10.0
y0 = [0.5]

# 对比三种显式方法
for method in ['RK45', 'RK23', 'DOP853']:
    sol = solve_ivp(logistic, [0, 5], y0, args=(r, K),
                    method=method, rtol=1e-6, atol=1e-9, dense_output=True)
    print(f"{method}: {sol.nfev} 次函数调用")

输出：

RK45: 158 次函数调用
RK23: 371 次函数调用
DOP853: 173 次函数调用

1.4.2 精度与效率分析

上图对比了三种显式方法在同一问题上的表现（均设置 rtol=1e-6）：

• 左图 RK45：误差在 10⁻⁶∼10⁻⁸之间波动，函数调用 158 次，是默认方法的最佳平衡

• 中图 RK23：误差约在 10⁻⁵量级，但函数调用高达 371 次——低阶方法在相同精度下需要更小的步长，反而更慢

• 右图 DOP853：误差约在 10⁻⁷∼10⁻⁸之间，函数调用 173 次——高阶方法每步开销大，但步长可以更大，综合效率优秀

选择建议：

• 默认用 RK45——大多数非刚性问题的最佳选择

• 需要更高精度时尝试 DOP853——虽然每步开销大，但总步数少

• 避免用 RK23——在相同精度下通常不如 RK45 高效

• 遇到求解器报错或步长过小的警告时，可能是刚性问题，改用 BDF 或 LSODA（第4章详细讲解）

1.5 t_eval 参数的使用

solve_ivp 默认只在求解器自适应选择的时间点输出解，但有时你需要在特定时间点获得解值：

# 方法1：用 t_eval 指定输出时间点
sol = solve_ivp(exp_growth, [0, 5], [1.0], args=(0.5,),
                t_eval=np.linspace(0, 5, 50))

print(f"输出 {len(sol.t)} 个点")  # 50 个点

# 方法2：用 dense_output 进行任意点插值
sol = solve_ivp(exp_growth, [0, 5], [1.0], args=(0.5,), dense_output=True)

# 在任意时间点求值（甚至超出积分区间，但外推不可靠）
t_query = [0.5, 1.5, 2.5, 3.5, 4.5]
y_query = sol.sol(t_query)

推荐做法：通常 dense_output=True 更灵活，且插值精度足够高（对于 RK45 等方法是 4 阶 Hermite 插值）。

1.6 实用技巧速查

1.7 本节小结

• solve_ivp 是 SciPy 1.0+ 推荐的 ODE 求解接口，统一了所有求解方法的 API

• 基本调用格式：solve_ivp(fun, [t0, tf], y0, args=(), method='RK45')

• 函数签名必须是 f(t, y, *args) → dydt，其中 y 是数组，即使只有一个方程

• 返回值 sol.t 是求解器自适应选择的时间步，不等于 t_span 的等分点

• dense_output=True 返回连续插值函数 sol.sol(t)，可在任意时间点查询

• rtol 和 atol 控制精度——减小它们可以获得更高精度但计算更慢

• nfev（函数调用次数）是衡量求解效率的关键指标

• 默认 RK45 适用于大多数非刚性问题，遇到报错时考虑 BDF 或 LSODA

第2章 高阶方程降阶

2.1 核心思想：降阶为方程组

solve_ivp 只能求解一阶常微分方程或一阶方程组。对于 n 阶 ODE：

标准降阶技巧：引入新变量 $y_1 = y,\ y_2 = y',\ y_3 = y'',\ \dots,\ y_n = y^{(n-1)}$，得到等价的一阶方程组：

初始条件也相应转换为：$y_1(0) = y(0),\ y_2(0) = y'(0),\ \dots,\ y_n(0) = y^{(n-1)}(0)$。

from scipy.integrate import solve_ivp
import numpy as np

# 示例：二阶方程 y'' + y = 0, y(0)=1, y'(0)=0
# 降阶：令 y1=y, y2=y'
# 得方程组：y1' = y2,  y2' = -y1

def harmonic(t, y):
    y1, y2 = y
    return [y2, -y1]

sol = solve_ivp(harmonic, [0, 10], [1.0, 0.0], dense_output=True)

2.2 阻尼谐振子

2.2.1 物理模型

质量为 m 的物体通过劲度系数为 k 的弹簧连接到墙壁，并受到与速度成正比的阻尼力 −cv：

令 $v = \frac{dx}{dt}$，降阶为一阶方程组：

def harmonic_oscillator(t, state, m, c, k):
    x, v = state
    return [v, (-c * v - k * x) / m]

m, c, k = 1.0, 2.0, 10.0
y0 = [1.0, 0.0]  # x(0)=1, v(0)=0
sol = solve_ivp(harmonic_oscillator, [0, 10], y0, args=(m, c, k))

2.2.2 三种阻尼状态

左图比较了三种阻尼情形下的位移响应 x(t)：

• 无阻尼（蓝色）：c=0，系统以固有频率 $\omega_n = \sqrt{k/m} = \sqrt{10} \approx 3.16$ 做等幅振荡，永不衰减

• 欠阻尼（红色）：$c=2 < c_{\text{cr}} = 2\sqrt{mk} \approx 6.32$，系统振荡但振幅逐渐衰减

• 临界阻尼（绿色）：c=ccr​，系统最快地回到平衡位置且不振荡——这是实际工程中（如汽车减震器）追求的最优阻尼

中图和右图是相图（phase portrait）——将位移 x 作为横坐标、速度 v 作为纵坐标，展示系统状态在相空间中的轨迹：

• 无阻尼相图：闭合的椭圆轨道——系统能量守恒，状态在相空间中循环往复

• 欠阻尼相图：向内螺旋的曲线——每次振荡都损失能量，最终收敛到原点 (0,0)

相图是分析非线性系统的核心工具——后续章节中的捕食者-猎物模型、传染病模型都要用相图来理解系统的全局行为。

2.3 单摆方程

2.3.1 从线性到非线性

单摆的运动方程是经典的二阶非线性 ODE：

其中 θ 为摆角，g 为重力加速度，L 为摆长，b 为阻尼系数。

小角度近似下 sinθ≈θ，方程退化为线性谐振子。但大角度时必须保留非线性项——这正是数值方法的优势所在：解析解只对线性方程容易获得，而非线性方程可以直接数值求解。

降阶后的一阶方程组：

def pendulum(t, state, g, L, b):
    theta, omega = state
    return [omega, -(g / L) * np.sin(theta) - b * omega]

g, L, b = 9.8, 1.0, 0.2
sol = solve_ivp(pendulum, [0, 20], [np.pi / 3, 0], args=(g, L, b))

2.3.2 不同初始角度的行为

左图展示了从不同初始角度释放的单摆运动（有阻尼 b=0.2）：

• 小角度（30°、60°）：运动接近简谐振动，周期几乎不变

• 大角度（90°、162°）：周期明显变长，波形偏离正弦——这就是非线性效应

• 162° 的红色曲线在 t=0 附近有较大的初始速度变化——摆从接近倒立的位置释放，开始阶段运动很快

中图是有阻尼相图：所有轨迹都螺旋收敛到 (0,0)——无论初始条件如何，摆最终都会静止在最低点。

右图是无阻尼相图：闭合轨道，系统能量守恒。值得注意的是：

• 小角度轨道接近圆形（近似简谐振动）

• 大角度轨道呈蛋形——非线性导致相空间轨道变形

• 最外层的红色轨道接近分界线（separatrix）——这是摆从振荡运动过渡到旋转运动的能量边界。超过这条边界的轨迹将变成绕定点旋转而非摆动

2.4 受迫振动与共振

2.4.1 模型

在阻尼谐振子基础上增加外部驱动力 $F(t) = F_0\cos(\omega t)$：

降阶后：

​​关键点：右端函数显含时间 t——这在 solve_ivp 中是完全允许的，因为函数签名本身就是 f(t, y, *args)。

def forced_oscillator(t, state, m, c, k, F0, omega):
    x, v = state
    return [v, (-c * v - k * x + F0 * np.cos(omega * t)) / m]

sol = solve_ivp(forced_oscillator, [0, 50], [0, 0], args=(1.0, 0.2, 4.0, 1.0, 2.0))

2.4.2 共振现象

左图展示了四种驱动频率下的响应——注意初始阶段的暂态过程（前约 10 个时间单位）和之后的稳态振荡：

• 低频驱动（0.5ωn​，蓝色）：响应振幅较小，相位滞后小

• 接近共振（0.9ωn​，绿色）：振幅显著增大

• 共振（ωn​，红色）：振幅达到最大——驱动力频率与系统固有频率匹配，能量输入最大化

• 高频驱动（1.1ωn​，紫色）：振幅减小，相位反转

右图是频率响应曲线：横坐标为驱动频率 ω，纵坐标为稳态振幅。数值解（蓝色实线）与理论值（红色虚线）完美吻合：

​​峰值出现在 $\omega_{\text{peak}} = \sqrt{\frac{k}{m} - \frac{c^2}{2m^2}}$，略小于固有频率 $\omega_n = \sqrt{k/m}$​——阻尼越大，共振峰越低、越宽。

建模应用：在数学建模中，共振分析常用于：

• 机械结构设计（避免共振频率与外部激励匹配）

• 电路谐振分析

• 地震工程中建筑物固有频率与地震波频率的关系

2.5 降阶的一般步骤总结

代码模板：

def higher_order_system(t, Y, *params):
    """
    将 n 阶 ODE 转化为一阶方程组的通用模板
    
    Y = [y_1, y_2, ..., y_n] 对应 [y, y', ..., y^{(n-1)}]
    返回 [y_1', y_2', ..., y_n'] = [y_2, y_3, ..., f(t, Y, *params)]
    """
    y1, y2, y3 = Y  # 对三阶方程
    return [
        y2,                        # y_1' = y_2
        y3,                        # y_2' = y_3
        f(t, y1, y2, y3, *params)  # y_3' = f(t, Y)
    ]

2.6 本节小结

• solve_ivp 只能求解一阶方程——高阶方程必须降阶为一阶方程组

• 降阶公式：对 n 阶方程，令 yi​=y(i−1)，得到 n 个一阶方程

• 相图（y vs y′）是理解系统行为的直观工具：闭合曲线 = 周期解，螺旋线 = 衰减解

• 非线性方程（如单摆）无法解析求解时，数值方法依然可以直接处理

• 右端函数可以显含时间 t——受迫振动就是典型例子

• 共振分析在工程建模中极为重要：固有频率与驱动频率的匹配决定振幅大小

第3章 常微分方程组

3.1 方程组的基本形式

solve_ivp 的核心能力是求解一阶方程组：

​其中 y=(y1​,y2​,…,yn​) 是 n 维向量。右端函数 f 接收标量 t 和数组 y，返回同样长度的数组 dy/dt。

from scipy.integrate import solve_ivp
import numpy as np

# 最简示例：两个耦合方程
def coupled_system(t, y, a, b):
    y1, y2 = y
    dy1 = a * y1 - b * y1 * y2
    dy2 = b * y1 * y2 - a * y2
    return [dy1, dy2]

sol = solve_ivp(coupled_system, [0, 20], [1.0, 2.0], args=(1.0, 0.5))

# 解的索引：sol.y[0] 是 y1 的解，sol.y[1] 是 y2 的解

3.2 Lotka-Volterra 捕食者-猎物模型

3.2.1 模型推导

Lotka-Volterra 模型描述了捕食者-猎物两个种群的相互作用：

​其中：

四项的含义：

• αx：猎物在无捕食者时的指数增长

• −βxy：捕食导致的猎物减少（与两个种群数量的乘积成正比）

• δxy：捕食者因捕食而增长

• −γy：捕食者无猎物时的自然衰减

def lotka_volterra(t, state, alpha, beta, delta, gamma):
    x, y = state
    dx = alpha * x - beta * x * y
    dy = delta * x * y - gamma * y
    return [dx, dy]

alpha, beta, delta, gamma = 1.0, 0.1, 0.1, 1.5
sol = solve_ivp(lotka_volterra, [0, 30], [10, 5],
                args=(alpha, beta, delta, gamma), dense_output=True)

3.2.2 周期性振荡

左图（时间序列）揭示了两个种群的相位差：

• 猎物（蓝色）先增长，捕食者（红色）滞后约 1/4 周期才达到峰值——因果链：猎物增多 → 捕食者食物充足 → 捕食者增多 → 猎物被大量捕食 → 猎物减少 → 捕食者饿死 → 猎物恢复

• 每个周期中，猎物的峰值约为捕食者峰值的 2 倍——捕食者需要大量猎物才能维持种群

• 周期约为 5-6 个时间单位

右图（相图）展示了闭合轨道——这是保守系统的标志：

• 系统存在守恒量（首次积分）：$V(x, y) = \delta x - \gamma \ln x + \beta y - \alpha \ln y$，在每条轨道上为常数

• 不同初始条件会产生不同大小的闭合轨道——振幅取决于初始状态

• 中心点 $(\gamma/\delta, \alpha/\beta) = (15, 10)$ 是非零平衡点

建模技巧：在数学建模中，Lotka-Volterra 型模型可以扩展到：

• 三种或更多物种的食物链

• 加入环境承载力（logistic 增长）

• 加入季节性变化（参数随时间周期性变化）

• 加入随机扰动（SDE 模型）

3.3 SIR 传染病模型

3.3.1 模型推导

SIR 模型将人群分为三个舱室：​

其中 S+I+R=1（归一化人群比例）。

关键参数：基本再生数 R0​=β/γ

• R0​>1：每个感染者平均传染超过 1 人，疫情爆发

• R0​<1：疫情自行消退

• R0​ 是衡量传染病传播能力的核心指标——新冠原始毒株 R0​≈2.5−3，Delta 约 5-8

def sir_model(t, state, beta_sir, gamma_sir):
    S, I, R = state
    dS = -beta_sir * S * I
    dI = beta_sir * S * I - gamma_sir * I
    dR = gamma_sir * I
    return [dS, dI, dR]

beta_sir, gamma_sir = 0.5, 0.1
sol = solve_ivp(sir_model, [0, 160], [0.99, 0.01, 0.0],
                args=(beta_sir, gamma_sir))

3.3.2 SIR 模型的动态行为

左图（时间序列）展示了典型的疫情曲线——与新冠疫情期间大家熟悉的曲线一致：

• 感染者 I（红色）：先指数增长，约第 15 天达到峰值（约 48% 人口同时感染），然后快速下降，约第 80 天基本清零

• 易感者 S（蓝色）：持续下降——越来越多的人被感染后移出易感群体

• 康复者 R（绿色）：持续增长，最终接近 100%——疫情结束后几乎所有人都感染过

右图（I vs S 相图）揭示了 SIR 模型的重要特性：

• 轨迹从 $(S, I) = (0.99, 0.01)$开始，I 先上升后下降——这是因为 $dI/dt = \beta SI - \gamma I = (\beta S - \gamma)I$，当 $S > \gamma/\beta = 1/R_0$ 时 I 增长，当 $S < 1/R_0$​ 时 I 减少

• 峰值感染点（红星）恰好出现在 $S = 1/R_0 = 0.2$ 处——这是 SIR 模型的精确理论结果

• 轨迹不闭合——SIR 不是保守系统，最终收敛到 (S∞​,0)

建模应用：SIR 模型是数学建模竞赛中的高频题型，常见的扩展方向：

• SEIR：加入潜伏期（Exposed 舱室）

• 加入疫苗接种：$dS/dt = -\beta SI - vS$

• 加入隔离措施：将 β 设为时间的分段函数

• 多人群模型：不同年龄组的耦合 SIR

3.4 参数敏感性分析

3.4.1 Lotka-Volterra 的捕食率影响

左图展示了捕食率 β 对猎物动态的影响：

• β=0.05（蓝色）：捕食压力小，猎物峰值更高、周期更长

• β=0.1（绿色）：基准情形

• β=0.2（红色）：捕食压力大，猎物峰值更低、周期更短——两个种群更紧密地耦合

3.4.2 感染率对疫情峰值的影响

中图展示了感染率 β 对 SIR 模型中感染者峰值的影响：

• β=0.2（R0​=2，蓝色）：峰值约 12%，疫情发展缓慢，持续约 100 天

• β=0.5（R0​=5，绿色）：峰值约 48%，爆发迅猛

• β=0.8（R0​=8，红色）：峰值超过 60%，疫情在 20 天内达到高峰

建模启示：这解释了为什么减少接触（降低 β）是疫情防控的核心——即使 R0​ 仍然大于 1，降低 β 也能显著降低峰值感染比例，减轻医疗系统压力。这就是"压平曲线"（flatten the curve）的数学原理。

3.4.3 R₀ 与感染峰值的定量关系

右图是基本再生数 R0​ 与最大感染比例的定量关系：

• R0​≤1（红线左侧）：疫情不会爆发，I 单调下降

• R0​=2：峰值约 12%

• R0​=5：峰值约 48%

• R0​=10：峰值超过 60%

关系呈高度非线性——R0​ 从 2 增加到 5，峰值从 12% 飙升到 48%。这提醒我们：微小的参数变化可能导致灾难性的后果。

3.5 建模实战：自定义方程组模板

def custom_odes(t, Y, *params):
    """
    自定义 ODE 方程组的通用模板
    
    参数：
        t : float - 当前时间
        Y : array - 状态向量 [y1, y2, ..., yn]
        params : tuple - 模型参数
    
    返回：
        list - 各变量的导数 [dy1/dt, dy2/dt, ..., dyn/dt]
    """
    y1, y2, y3 = Y  # 根据变量个数调整
    p1, p2, p3 = params
    
    dy1 = ...  # 根据模型方程填写
    dy2 = ...
    dy3 = ...
    
    return [dy1, dy2, dy3]

# 求解
sol = solve_ivp(custom_odes, [0, 100], [y1_0, y2_0, y3_0],
                args=(p1_val, p2_val, p3_val), dense_output=True)

3.6 本节小结

• 方程组的写法与单个方程完全相同，只是 y 变为数组，返回值为数组

• Lotka-Volterra：捕食者-猎物的周期性振荡，相图呈闭合轨道

• SIR：疫情从爆发到消退的完整过程，R0​=β/γ 决定是否会爆发

• 相图是理解系统全局行为的利器：闭合 = 周期，螺旋 = 衰减，开放曲线 = 非保守系统

• 参数敏感性分析是数学建模的核心环节——微小的参数变化可能导致质的改变

• SIR 相图中峰值感染点满足 S=1/R0​——这是重要的理论结果

第4章 求解器选择与刚性方程

4.1 什么是"刚性"问题

在求解 ODE 时，solve_ivp 的默认方法是 RK45（显式 Runge-Kutta 4(5) 阶），它在大多数情况下都能很好地工作。但对于刚性问题（stiff problems），显式方法会遭遇严重困难。

刚性的直观理解

刚性系统的特征是系统中存在多个差异极大的时间尺度。快速衰减的瞬态分量强迫求解器取极小的步长以保证数值稳定，但慢速变化分量又需要积分很长时间才能看到有意义的结果。显式方法因为数值稳定性限制（CFL 条件），步长被快速分量"绑架"，导致计算成本爆炸。

形象地说：显式方法像一个只靠视觉判断路况的司机——前方有一点颠簸，就猛踩刹车减速；而隐式方法像一个有预知能力的司机——它"看穿"了颠簸的本质，知道减速幅度应该多大。

数学定义

线性系统 $\frac{dy}{dt} = Ay$ 的刚性由 Jacobian 矩阵 A 的特征值比衡量：

​比值越大，问题越刚性。刚性比达到 $10^3 \sim 10^6$ 或更高时，显式方法就难以胜任了。

4.2 Robertson 问题——经典的刚性测试

Robertson 问题是刚性 ODE 的标准测试案例，描述一个三组分化学反应：

​​参数 $k_1 = 0.04,\ k_2 = 3\times 10^7,\ k_3 = 10^4$ 跨越了 9 个数量级——这就是刚性的根源。

from scipy.integrate import solve_ivp
import numpy as np

def robertson(t, y, k1, k2, k3):
    y1, y2, y3 = y
    dy1 = -k1 * y1 + k3 * y2 * y3
    dy2 = k1 * y1 - k2 * y2**2 - k3 * y2 * y3
    dy3 = k2 * y2**2
    return [dy1, dy2, dy3]

k1, k2, k3 = 0.04, 3e7, 1e4
y0 = [1.0, 0.0, 0.0]

# 使用 BDF（隐式方法）求解刚性问题
sol = solve_ivp(robertson, [0, 1e11], y0, args=(k1, k2, k3),
                method='BDF', rtol=1e-6, atol=1e-10)

求解结果

左图展示了 Robertston 反应的浓度变化（横轴为对数刻度，覆盖 t∈[1,10¹¹]）：

• y₁（蓝色）：反应物 A，从 1.0 单调衰减，约在 t=10⁵ 时接近零

• y₂（红色）：中间物 B，浓度始终极低（约 10⁻⁴量级）——它生成极快也消耗极快，是典型的刚性瞬态分量

• y₃（绿色）：产物 C，逐渐累积至接近 1.0

右图对比了三种适用于刚性问题的求解器：

RK45 和 RK23 呢？ 它们在 Robertson 问题上会尝试取步长约 10⁻¹⁰ 才能维持稳定——在 [0,10¹¹] 的积分区间内需要约 10²¹ 步，实际上永远无法完成。这就是刚性的威力。

4.3 非刚性问题：所有方法都能胜任

对于非刚性问题，各方法都能给出准确解，但效率差异明显。以简谐振子 $\ddot{x} + x = 0$ 为例：

左图（误差对比）：

• DOP853（紫色）精度最高——8 阶方法在相同容差下产生的截断误差最小

• RK45（蓝色）精度适中，是精度与效率的最佳平衡

• BDF（红色）和 LSODA（橙色）作为隐式方法也能工作，但在非刚性问题上并不比显式方法更优

• RK23（绿色）误差最大——2/3 阶方法精度不足

右图（函数调用次数）：

• DOP853 仅需 797 次——高阶方法步长可以更大

• RK45 需要 1778 次——中等开销

• RK23 高达 6677 次——低阶导致步长过小，反而最慢

结论：对于非刚性问题，优先使用 RK45（默认）或 DOP853（高精度需求）。不要用 RK23——它并不快。

4.4 Van der Pol 振荡器：从非刚性到刚性的过渡

Van der Pol 方程是一个经典的非线性振荡器：

降阶后：

​参数 μ 控制非线性强度，也决定了问题的刚性程度。

μ 较小时：非刚性

左图展示了 μ=0.5,1,5 时的时间响应。当 μ 较小时，系统表现为平滑的准正弦振荡——这是非刚性的，RK45 可以高效求解。

极限环

中图展示了不同 μ 值下的相图（极限环）：

• μ 越小，极限环越接近圆形（接近线性谐振子）

• μ 越大，极限环越趋向矩形——系统在慢变和快变两种模式间切换，这就是刚性的来源

μ 很大时：刚性显现

右图对比了 μ=1000 时 RK45 与 BDF 的表现：

当 μ 很大时，Van der Pol 振荡器的解在"快速跳跃"和"缓慢漂移"之间交替（右图中几乎垂直的跳变部分）。显式方法必须在快速阶段取极小步长，即使缓慢阶段可以用大步长也做不到——因为步长一旦放大就不稳定。

这就是刚性的典型表现。此时切换到 BDF、Radau 或 LSODA 是必须的。

4.5 求解器选择速查表

切换方法很简单

# 默认方法（非刚性）
sol = solve_ivp(fun, [0, 10], y0)  # 默认 method='RK45'

# 发现警告或速度极慢 → 改用刚性求解器
sol = solve_ivp(fun, [0, 10], y0, method='BDF')

# 不确定 → 让系统自动判断
sol = solve_ivp(fun, [0, 10], y0, method='LSODA')

如何识别刚性问题

以下信号提示你遇到了刚性问题：

1. 求解极慢——函数调用次数异常多

2. 警告信息——"Required step size is less than spacing between numbers"

3. 解出现非物理振荡——数值不稳定导致的伪振荡

4. 模型包含快慢两个时间尺度——如化学反应中快反应和慢反应共存

4.6 实用技巧

from scipy.integrate import solve_ivp

# 技巧1：限制最大步长（防止错过快速瞬态）
sol = solve_ivp(fun, [0, 100], y0, max_step=0.1)

# 技巧2：收紧容差提高精度
sol = solve_ivp(fun, [0, 100], y0, rtol=1e-8, atol=1e-10)

# 技巧3：检查求解状态
if not sol.success:
    print(f"求解失败: {sol.message}")

# 技巧4：监控函数调用次数（性能指标）
print(f"函数调用: {sol.nfev} 次")
print(f"成功: {sol.success}, 状态码: {sol.status}")

4.7 本节小结

• 刚性源于系统中存在差异极大的时间尺度，显式方法被迫取极小步长

• RK45/RK23/DOP853 是显式方法，适用于非刚性问题

• BDF/Radau 是隐式方法，适用于刚性问题

• LSODA 能自动切换显式/隐式模式，是不确定时的安全选择

• Van der Pol 方程展示了从非刚性到刚性的连续过渡：μ 越大越刚性

• Robertson 问题是经典刚性测试：三个速率常数跨越 9 个数量级

• 遇到求解极慢或警告时，第一反应应该是换用 BDF 或 LSODA

第5章 事件检测与密集输出

5.1 事件检测（Events）

solve_ivp 最强大的功能之一是事件检测——可以在积分过程中自动检测特定条件（如变量穿过零点），并精确地定位事件发生的时间点。

事件函数的基本用法

事件函数是一个签名为 event(t, y, *args) → float 的函数。当返回值为 0 时表示事件发生：

from scipy.integrate import solve_ivp
import numpy as np

def my_event(t, y):
    return y[0] - 1.0  # 当 y[0] = 1 时触发

sol = solve_ivp(fun, [0, 100], y0, events=my_event)

每个事件函数有两个可选属性：

my_event.terminal = True   # 事件发生时终止积分（默认 False）
my_event.direction = -1    # 检测方向：-1=仅下降穿过, 1=仅上升穿过, 0=双向（默认）

示例1：弹跳球

一个经典的物理事件场景：球从高处落下，每次触地反弹，恢复系数为 0.8。

左图展示了弹跳球的高度-时间曲线。每次触地（y = 0）时，事件函数被触发：

• 红色虚线标记事件触发时刻——solve_ivp 精确地定位了每次撞击的时间，而不是简单地在某个时间步"发现"y 已变为负值

• 红色圆点是每次弹跳的撞击点——事件检测精度远高于求解器的步长精度

• 通过设置 terminal = True，积分在每次撞击时暂停，我们修改速度（乘以恢复系数）后重新启动积分

右图展示了每次弹跳的最大高度——按照物理规律，每次反弹后高度为前一次的 $e^2 = 0.64$ 倍（恢复系数 0.8 的平方）。

def ground_event(t, state, g, restitution):
    """事件函数：球触地"""
    y, v = state
    return y  # y = 0 时触发

ground_event.terminal = True   # 触地时停止积分
ground_event.direction = -1    # 只检测 y 下降穿过 0（排除反弹瞬间）

# 每次触地后，链式调用 solve_ivp
for bounce in range(6):
    sol = solve_ivp(bouncing_ball, [t_start, t_start + 50], y_current,
                    args=(g, restitution), events=ground_event)

    if sol.t_events[0].size > 0:
        t_bounce = sol.t_events[0][0]
        v_after = -restitution * v_before  # 反弹
        y_current = [0.0, v_after]
        t_start = t_bounce

示例2：抛体运动的射程计算

事件检测的另一个经典应用是自动停止积分——不需要预先知道精确的积分区间。

左图展示了以不同角度发射的抛体轨迹。每个轨迹都设置了地面事件（y = 0），积分在精确的落地时刻自动停止——不需要猜测积分时间。图例中直接显示了用事件检测精确计算的射程。

右图通过遍历所有角度，绘制了射程-角度曲线。红色圆点标记了最优发射角度：

• 最优角度 ≈ 45°，与理论结果 $v_0^2/g = 2500/9.8 \approx 255$ 米完全吻合

• 30° 和 60° 的射程相同——这是抛体运动的经典对称性（互补角射程相等）

def ground_hit(t, state, g):
    return state[1]  # y = 0 时触发

ground_hit.terminal = True
ground_hit.direction = -1  # 只检测落地（y 从正变零）

# 求解器在精确的落地时刻停止
sol = solve_ivp(projectile, [0, 20], y0, args=(g,), events=ground_hit)

# 精确的射程
range_val = sol.t_events[0][0] * v0 * np.cos(theta)

5.2 密集输出（Dense Output）

solve_ivp 返回的 sol.t 是求解器自适应选择的时间点——这些点可能不均匀，且数量不一定满足你的需求。dense_output=True 提供了一个连续的插值函数，可以在任意时间点查询解。

密集输出的工作原理

求解器在每个内部步长上构建一个插值多项式，使得：

• 在步的端点处精确匹配解值

• 插值多项式的阶数与求解方法匹配

sol = solve_ivp(fun, [0, 40], y0, dense_output=True)

# 在任意时间点查询（甚至不是求解器的内部步长点）
t_query = np.linspace(0, 40, 2000)
y_at_query = sol.sol(t_query)  # 连续插值，返回形状 (n_eq, len(t_query))

示例：Lorenz 混沌系统

左图对比了求解器的离散输出点（红色圆点）和密集输出插值（蓝色曲线）：

• Lorenz 系统是经典的混沌系统，变化剧烈且不规则

• 求解器只用了约 600 个离散点就覆盖了 [0, 40] 的积分区间

• 通过 dense_output=True，我们可以在任意时间点（这里是 2000 个均匀点）获得平滑的插值结果

• 插值精度由求解器的容差保证——对于 rtol=1e-8 的设置，插值误差也在相同量级

中图和右图展示了 Lorenz 吸引子的两个投影——混沌系统的标志性"蝴蝶"形状。密集输出确保了相图曲线的平滑性。

def lorenz(t, state, sigma, rho, beta):
    x, y, z = state
    return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]

# 高精度求解
sol = solve_ivp(lorenz, [0, 40], [1,1,1],
                args=(10, 28, 8.0/3.0),
                rtol=1e-8, atol=1e-10, dense_output=True)

# 密集输出用于平滑的相图
t_dense = np.linspace(0, 40, 2000)
x, y, z = sol.sol(t_dense)

# 相图投影
plt.plot(x, y, 'b-', lw=0.8)  # x-y 投影
plt.plot(y, z, 'r-', lw=0.8)  # y-z 投影

5.3 多事件检测

solve_ivp 可以同时监测多个事件函数：

def event1(t, y):
    return y[0] - 1.0  # 当 y[0] = 1 时触发

def event2(t, y):
    return y[1] - 0.5  # 当 y[1] = 0.5 时触发

# 将多个事件函数放在列表中
sol = solve_ivp(fun, [0, 100], y0, events=[event1, event2])

# 每个事件的结果分别存储在 t_events 列表中
print(f"事件1触发时间: {sol.t_events[0]}")
print(f"事件2触发时间: {sol.t_events[1]}")

方向控制的实际应用

direction 属性在物理场景中非常有用：

# 弹簧振子：检测每次经过平衡位置
def zero_crossing(t, y):
    return y[0]  # x = 0

zero_crossing.direction = 1  # 只检测从左向右穿过

# 弹跳球：只检测落地（y 从正变零），不检测反弹
def ground(t, y):
    return y[0]  # y = 0

ground.direction = -1  # 只检测下降穿过
ground.terminal = True  # 触地即停

5.4 实战技巧

技巧1：链式事件求解

对于需要反复处理同一事件的情况（如多次碰撞、周期性触发），链式调用 solve_ivp：

t_start = 0.0
y_current = y0
results = []

while t_start < t_final:
    sol = solve_ivp(fun, [t_start, t_final], y_current,
                    events=my_event, dense_output=True)

    results.append(sol)

    if sol.t_events[0].size == 0:
        break  # 没有更多事件了

    t_start = sol.t_events[0][0]
    y_current = handle_event(sol.y[:, -1])  # 自定义事件处理

技巧2：事件函数的精度

事件检测的精度独立于求解器的步长——solve_ivp 会用二分法精确找到事件发生的时刻，精度通常远优于 rtol 和 atol 设置的容差。

sol = solve_ivp(fun, [0, 100], y0, events=my_event, rtol=1e-3)

# 事件检测精度通常比 rtol 高几个数量级
# 即使 rtol=1e-3，事件触发时间的精度可能达到 1e-8 或更高

技巧3：密集输出的使用场景

5.5 本节小结

• 事件函数返回零时触发事件，通过 terminal 和 direction 控制行为

• terminal=True 在事件发生时停止积分——适用于碰撞、越界等场景

• direction=-1/1/0 控制检测穿越方向——避免重复触发

• 多个事件函数可以同时监测，结果分别存储在 sol.t_events 列表中

• 密集输出（dense_output=True）提供连续插值函数 sol.sol(t)

• 事件检测精度通常远高于求解器步长精度——二分法定位

• Lorenz 混沌系统的相图是密集输出的经典应用——确保曲线平滑

• 链式调用 solve_ivp 可以处理反复发生的事件序列

第6章 偏微分方程简介

6.1 从 ODE 到 PDE：方法线（Method of Lines）

scipy.integrate 没有内置的 PDE 求解器，但有限差分法配合 solve_ivp 可以高效求解许多常见的偏微分方程。

核心思想——方法线（Method of Lines）：

1. 将空间域离散化为网格点

2. 用有限差分近似空间导数

3. 得到一组耦合的 ODE（每个网格点一个方程）

4. 用 solve_ivp 求解这个 ODE 系统

import numpy as np
from scipy.integrate import solve_ivp

# 空间离散
N = 50
L = 1.0
dx = L / (N - 1)
x = np.linspace(0, L, N)

# 将 PDE 转化为 ODE 系统
def pde_to_ode(t, u, alpha, dx):
    dudt = np.zeros_like(u)
    # 内部点用中心差分近似 d²u/dx²
    dudt[1:-1] = alpha * (u[2:] - 2*u[1:-1] + u[:-2]) / dx**2
    # 边界条件
    dudt[0] = 0
    dudt[-1] = 0
    return dudt

# 用 solve_ivp 求解
sol = solve_ivp(pde_to_ode, [0, T_final], u0, args=(alpha, dx))

这种方法的优势是复用 solve_ivp 的全部功能——自适应步长、刚性求解器、事件检测——全部直接可用。

6.2 热传导方程

一维热传导方程是最经典的抛物型 PDE：

​其中 α 是热扩散系数，u(x,t) 是温度分布。

数值离散

用二阶中心差分近似空间二阶导数：

​​得到 ODE 系统：

def heat_eq_ode(t, u, alpha, dx):
    dudt = np.zeros_like(u)
    dudt[1:-1] = alpha * (u[2:] - 2*u[1:-1] + u[:-2]) / dx**2
    # 边界条件：u(0,t) = u(L,t) = 0
    dudt[0] = 0
    dudt[-1] = 0
    return dudt

N = 50
x = np.linspace(0, 1.0, N)
dx = 1.0 / (N - 1)
alpha = 0.01

# 初始条件：高斯脉冲
u0 = np.exp(-((x - 0.5) / 0.1)**2)

sol = solve_ivp(heat_eq_ode, [0, 0.5], u0, args=(alpha, dx),
                rtol=1e-6, atol=1e-8, dense_output=True)

数值结果

左图展示了高斯初始温度分布随时间的扩散过程：

• t = 0（深紫色）：初始高斯脉冲，中心在 x = 0.5，宽度 σ = 0.1

• t = 0.02 ~ 0.20：脉冲逐渐变宽、变低——热量从高温区向低温区扩散

• t = 0.50（浅黄绿色）：脉冲已显著展宽，峰值降低约一半——但曲线下的面积（总热量）保持不变

右图是时空热力图，更直观地展示了扩散过程：

• 横轴是空间坐标 x，纵轴是时间 t

• 颜色越亮温度越高——初始的高亮区域随时间向两侧扩展

• 扩散前沿呈抛物线型扩展——符合热传导方程的标度律 x∼αt​

建模应用：热传导方程的离散形式可以类比到：

• 污染物在水体中的扩散

• 种群的空间扩散

• 金融中的 Black-Scholes 方程（经变量变换后等价于热传导方程）

6.3 波动方程

一维波动方程是经典的双曲型 PDE：

​其中 c 是波速。这是二阶时间导数方程，需要先降阶为一阶系统：

令 $v = \frac{\partial u}{\partial t}$，则：

​​离散后得到 2N 个 ODE——这是高阶方程降阶技巧（第2章）和空间离散的结合。

def wave_eq_ode(t, state, c, dx):
    N = len(state) // 2
    u = state[:N]  # 位移
    v = state[N:]  # 速度

    dudt = v.copy()
    dvdt = np.zeros(N)
    dvdt[1:-1] = c**2 * (u[2:] - 2*u[1:-1] + u[:-2]) / dx**2
    dvdt[0] = 0  # 边界
    dvdt[-1] = 0

    return np.concatenate([dudt, dvdt])

# 初始条件：高斯脉冲，零初速度
u0 = np.exp(-((x - 0.3) / 0.05)**2)
v0 = np.zeros(N)
state0 = np.concatenate([u0, v0])

sol = solve_ivp(wave_eq_ode, [0, 2.0], state0, args=(1.0, dx))

波动现象

左图展示了脉冲在固定端点条件下的传播：

• t = 0（深紫色）：初始脉冲位于 x = 0.3 处

• t = 0.25（紫色）：脉冲分裂为两个子脉冲，分别向左和向右传播——这是波动方程的特征解结构

• t = 0.50（粉紫色）：左行波到达左端点并反射，反射后波形反转（固定端点的反射）

• t = 0.75（粉色）：右行波到达右端点并反射

• t = 1.0（橙色）：两个反射波在中心相遇并叠加——产生干涉

右图的时空热力图清晰地展示了波的传播轨迹：

• 两条交叉的对角线——脉冲以恒定速度 c 向两侧传播

• 到达边界后折返——斜率不变但符号反转（反射）

• 中心交叉点——两列波相遇叠加

与热传导的对比：

• 热传导是耗散性的——脉冲不断展宽、振幅衰减

• 波动是保守性的——脉冲形状基本保持（有色散），能量守恒

6.4 Fisher-KPP 方程——反应扩散

Fisher-KPP 方程结合了扩散和非线性反应：

这是最简单的反应-扩散方程，描述了一个物种在空间中扩散并 logistic 增长的过程。

def fisher_kpp(t, u, D, r):
    dudt = np.zeros_like(u)
    # 扩散项 + 反应项
    dudt[1:-1] = D * (u[2:] - 2*u[1:-1] + u[:-2]) / dx**2 + r * u[1:-1] * (1 - u[1:-1])
    dudt[0] = 0
    dudt[-1] = 0
    return dudt

D, r = 0.001, 1.0
u0 = np.zeros(N)
u0[N//2] = 1.0  # 点源

sol = solve_ivp(fisher_kpp, [0, 5.0], u0, args=(D, r), dense_output=True)

行波解

左图展示了点源初始条件下的演化——形成了向右传播的行波前（traveling wave front）：

• t = 0（深紫色）：只在 x = 0.5 处有值 1，其余为零

• t = 0.5 ~ 2.0：波形向右扩展，形成一个陡峭的前沿

• t = 3.0 ~ 5.0：前沿持续向右传播，波形形状趋于稳定——这就是行波解

中图（热力图）更直观地展示了行波前：红色区域从左向右推进，形成一条斜线。

右图定量分析了波前传播速度：

• 蓝色曲线：数值计算得到的波前位置（u = 0.5 的位置）随时间的变化

• 红色虚线：理论波速 $c = 2\sqrt{rD} = 2\sqrt{1 \times 0.001} \approx 0.063$

• 数值结果与理论预测高度吻合——波前以恒定速度传播

建模应用：Fisher-KPP 方程可以描述：

• 入侵物种的空间扩散

• 化学反应前沿的传播

• 神经脉冲沿轴突的传播

• 肿瘤生长前沿

6.5 方法线的通用模板

import numpy as np
from scipy.integrate import solve_ivp

def solve_pde_1d(rhs_spatial, u0, t_span, x=None, N=50, L=1.0,
                 bc_left=None, bc_right=None, **kwargs):
    """
    一维 PDE 求解通用模板（方法线）

    参数：
        rhs_spatial : callable - 空间离散后的右端函数
                      签名 f(t, u) -> du/dt
        u0 : array - 初始条件
        t_span : tuple - 时间区间
        x : array - 空间网格（可选）
        N : int - 网格点数
        L : float - 空间域长度
        **kwargs : 传递给 solve_ivp 的额外参数

    返回：
        sol : Bunch - solve_ivp 的返回值
    """
    if x is None:
        x = np.linspace(0, L, N)

    sol = solve_ivp(rhs_spatial, t_span, u0, **kwargs)
    return sol, x

# 使用示例
N = 50
x = np.linspace(0, 1, N)
dx = x[1] - x[0]

def my_pde(t, u):
    dudt = np.zeros_like(u)
    dudt[1:-1] = ...  # 填入空间离散公式
    dudt[0] = ...     # 左边界
    dudt[-1] = ...    # 右边界
    return dudt

sol, x = solve_pde_1d(my_pde, u0, [0, 10], x=x)

6.6 注意事项与局限性

稳定性限制

显式时间积分方法（RK45）对空间离散有CFL 稳定性条件：

如果网格太细，RK45 的自适应步长可能变得非常小。此时应改用隐式方法：

# 细网格时改用 BDF
sol = solve_ivp(heat_eq_ode, [0, T], u0, args=(alpha, dx),
                method='BDF', rtol=1e-6, atol=1e-8)

二维及更高维度

方法线可以推广到二维：

# 2D 热传导: u_t = alpha * (u_xx + u_yy)
u = U.reshape(Nx, Ny)  # 将一维状态向量重塑为二维
laplacian = (np.roll(u, 1, axis=0) + np.roll(u, -1, axis=0) +
             np.roll(u, 1, axis=1) + np.roll(u, -1, axis=1) - 4*u) / dx**2
dudt = alpha * laplacian.flatten()

但当维度增加时，ODE 系统的规模变为 N² 或 N³，计算量急剧增加——此时应考虑专门的 PDE 求解器（如 FEniCS、FiPy）。

刚性问题

PDE 离散后的 ODE 系统通常是刚性的——空间离散引入了快慢两种时间尺度。推荐始终使用 BDF 或 LSODA：

# 推荐做法
sol = solve_ivp(pde_ode, [0, T], u0, method='BDF',
                rtol=1e-6, atol=1e-8)

6.7 本节小结

• PDE 可以通过方法线（Method of Lines）转化为 ODE 系统，用 solve_ivp 求解

• 热传导方程：抛物型 PDE，扩散过程，脉冲展宽

• 波动方程：双曲型 PDE，需降阶为两个一阶方程，脉冲传播

• Fisher-KPP：反应-扩散方程，产生行波解，波速 $c = 2\sqrt{rD}$​

• PDE 离散后的 ODE 系统通常是刚性的——优先使用 BDF 或 LSODA

• 空间网格越细，CFL 条件越严格——必要时减小 max_step

• 方法线适用于一维和简单的二维问题——复杂 PDE 需要专门的求解器

附录：完整代码获取

本教程所有代码均可通过以下链接下载：

ode.zip