04_数值模拟：微分方程求解

数值模拟：微分方程求解 (Numerical Simulation)

---

章节概述

科学计算的核心任务之一是求解微分方程——从牛顿力学到量子力学，从电路分析到流行病模型，微分方程无处不在。本章从常微分方程（ODE）的 scipy.integrate.solve_ivp 入手，覆盖 Runge-Kutta 方法家族；随后使用 NumPy 有限差分法求解一维热方程（PDE）实例；最后介绍蒙特卡洛模拟（π 估计、期权定价）作为随机模拟的入门。我们将对比 C 语言中手写 RK4 的代码量与 SciPy 一行调用的差距，并简要提及可视化输出的去向。

核心理念：数值模拟是”理论”与”实验”之外的第三条科研道路。Python 科学栈让你在几分钟内搭建起在 C 中需要几天才能写好的仿真原型——先理解你的模型行为，再决定是否需要 C 级别的性能优化。

---

第一节：常微分方程（ODE）数值解

1.1 solve_ivp 基础

SciPy 的 solve_ivp 求解初值问题 dy/dt = f(t, y)，提供多种 Runge-Kutta 方法。

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
 
# 经典 ODE：指数衰减 dy/dt = -k*y, y(0) = y0
def decay(t, y, k):
 return -k * y
 
k, y0 = 0.5, 10.0
sol = solve_ivp(decay, [0, 10], [y0], args=(k,),
 method='RK45', dense_output=True)
 
print(f"求解状态: {sol.message}")
print(f"时间点数: {len(sol.t)}")
print(f"y(5) ≈ {sol.sol(5)[0]:.4f}") # 解析解: 10*exp(-0.5*5) ≈ 0.8208
 
# 设置固定的评估点
t_eval = np.linspace(0, 10, 100)
sol = solve_ivp(decay, [0, 10], [y0], args=(k,),
 t_eval=t_eval, method='RK45')

1.2 Runge-Kutta 方法家族

方法	类型	精度	适用场景
RK45（默认）	显式，自适应步长	5(4) 阶	非刚性 ODE
RK23	显式，自适应步长	3(2) 阶	低精度需求
DOP853	显式，8 阶	高精度	高精度需求
Radau	隐式	5 阶	刚性 ODE
BDF	隐式，多步	可变阶	严重刚性 ODE
LSODA	自动切换	可变	未知刚性时

刚性（stiffness） 是数值 ODE 求解的核心挑战：当系统中存在差异巨大的时间尺度时，显式方法必须使用极小的步长才能稳定。隐式方法（Radau, BDF）通过求解非线性方程组来保持稳定性，虽然每步更贵但总体更高效。

1.3 刚性 ODE 示例

import numpy as np
from scipy.integrate import solve_ivp
import time
 
# Van der Pol 振荡器（在 μ 很大时刚性）
def vanderpol(t, y, mu):
 return [y[1], mu * (1 - y[0]**2) * y[1] - y[0]]
 
mu = 1000 # 大 μ → 刚性
y0 = [2, 0]
 
# 显式 RK45 — 极慢（步长被限制）
t0 = time.perf_counter()
sol_rk45 = solve_ivp(vanderpol, [0, 3000], y0, args=(mu,),
 method='RK45')
t1 = time.perf_counter()
print(f"RK45: {t1-t0:.2f}s, 步数: {len(sol_rk45.t)}")
 
# 隐式 Radau — 快很多
t0 = time.perf_counter()
sol_radau = solve_ivp(vanderpol, [0, 3000], y0, args=(mu,),
 method='Radau')
t1 = time.perf_counter()
print(f"Radau: {t1-t0:.2f}s, 步数: {len(sol_radau.t)}")
# Radau 步数远少于 RK45，整体时间大幅缩短

1.4 事件（Events）——提前终止

import numpy as np
from scipy.integrate import solve_ivp
 
# 抛体运动：y'' = -g，求解何时落地
def projectile(t, state):
 x, y, vx, vy = state
 return [vx, vy, 0, -9.81]
 
# 定义事件：y = 0（触地），检测下降过程中穿过零点
def hit_ground(t, state):
 return state[1] # y 坐标
 
hit_ground.terminal = True # 触地时停止积分
hit_ground.direction = -1 # 仅在下降时检测
 
sol = solve_ivp(projectile, [0, 10], [0, 0, 10, 20],
 events=hit_ground, max_step=0.01)
 
print(f"飞行时间: {sol.t[-1]:.3f} s")
print(f"落地点 x: {sol.y[0, -1]:.3f} m")

---

第二节：偏微分方程（PDE）— 热方程

2.1 一维热方程的有限差分解

热方程 ∂u/∂t = α ∂²u/∂x² 使用显式 Euler 时间推进 + 中心差分空间离散：

import numpy as np
 
# 参数
alpha = 0.01 # 热扩散系数
L = 1.0 # 杆的长度
nx = 50 # 空间网格点
dx = L / (nx - 1)
dt = 0.001 # 时间步长
nt = 5000 # 时间步数
 
# 稳定性条件（显式格式的 CFO 条件）
cfl = alpha * dt / dx**2
print(f"CFL = {cfl:.4f} （必须 < 0.5 以保证稳定）")
 
# 初始条件：高斯热脉冲
x = np.linspace(0, L, nx)
u = np.exp(-((x - 0.5)**2) / 0.01)
 
# 时间推进 — 向量化（无 Python 循环）
for n in range(nt):
 # 二阶中心差分：∂²u/∂x² ≈ (u[i+1] - 2*u[i] + u[i-1]) / dx²
 u[1:-1] += alpha * dt / dx**2 * (u[2:] - 2*u[1:-1] + u[:-2])
 
 # Dirichlet 边界条件：两端固定为 0
 u[0] = u[-1] = 0.0
 
 if n % 1000 == 0:
 print(f"步骤 {n:5d}, 总热量: {u.sum() * dx:.6f}")
 
print(f"最终总热量: {u.sum() * dx:.6f}")
# 热量应几乎完全耗散（因两端边界为 0）

CFL 稳定性条件：显式格式要求 α·dt/dx² < 0.5，否则数值解振动发散。隐式格式（如 Crank-Nicolson）无条件稳定，但每步需解线性方程组。

2.2 与 C 实现的对比

对比同样用 C 语言实现显式热方程求解器（伪代码）：

// C 语言手写热方程（约 40 行，需管理内存、边界、CFL 检查）
double *u = malloc(nx * sizeof(double));
double *u_new = malloc(nx * sizeof(double));
 
// 初始化...
for (int n = 0; n < nt; n++) {
 for (int i = 1; i < nx-1; i++) {
 u_new[i] = u[i] + alpha * dt / (dx*dx)
 * (u[i+1] - 2*u[i] + u[i-1]);
 }
 // 交换指针、处理边界...
 double *tmp = u; u = u_new; u_new = tmp;
}
// free / 错误检查 ...

NumPy 版本的代码行数约为 C 版本的 1/5，且无需管理内存、边界检查。但 C 版本在编译器 -O3 优化下可能比纯 NumPy 向量化还快（因为无需 Python 对象）。权衡：开发速度 vs 运行时性能。

关于性能对比的深入讨论，参见 性能对比章节。

---

第三节：蒙特卡洛模拟

3.1 π 估计

蒙特卡洛方法的核心思想是用随机采样逼近确定性量：

import numpy as np
 
N = 1_000_000
 
# 在 [-1,1] 正方形中均匀采样
x = np.random.uniform(-1, 1, N)
y = np.random.uniform(-1, 1, N)
 
# 点在单位圆内的比例 → π/4
inside = (x**2 + y**2) <= 1
pi_estimate = 4 * inside.sum() / N
 
print(f"π 估计值: {pi_estimate:.6f}")
print(f"真实值: {np.pi:.6f}")
print(f"误差: {abs(pi_estimate - np.pi):.6f}")
# 蒙特卡洛误差 ∝ 1/√N，约 0.001 量级

3.2 欧式期权定价

Black-Scholes 模型中，股票价格遵循几何布朗运动：

ST = S0 * exp((r - σ²/2)*T + σ*√T*ε) 其中 ε ~ N(0,1)

import numpy as np
 
S0, K, r, sigma, T = 100, 105, 0.05, 0.2, 1.0
N = 500_000
 
np.random.seed(42)
Z = np.random.randn(N)
ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
 
# 看涨期权 payoff = max(ST - K, 0)
payoff = np.maximum(ST - K, 0)
 
# 现值贴现
option_price = np.exp(-r * T) * payoff.mean()
std_error = np.exp(-r * T) * payoff.std() / np.sqrt(N)
 
print(f"期权价格: {option_price:.4f}")
print(f"95% 置信区间: [{option_price - 1.96*std_error:.4f}, {option_price + 1.96*std_error:.4f}]")
 
# Black-Scholes 解析解（对比用）
from scipy.stats import norm
d1 = (np.log(S0/K) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
d2 = d1 - sigma*np.sqrt(T)
bs_price = S0*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2)
print(f"BS 解析价格: {bs_price:.4f}")

蒙特卡洛的收敛速度是 O(1/√N)——要想提高一位小数精度，需要 100 倍的计算量。对于低维问题，确定性积分方法（如 Simpson 规则）更高效；蒙卡的优势在于高维积分（> 4 维），其中确定性方法遭遇”维数诅咒”。

3.3 方差缩减技巧

import numpy as np
 
# 对偶变量法：用一对相反的随机数，降低方差
N_half = N // 2
Z = np.random.randn(N_half)
ST1 = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
ST2 = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*(-Z))
 
payoff_avg = 0.5 * (np.maximum(ST1 - K, 0) + np.maximum(ST2 - K, 0))
price_antithetic = np.exp(-r*T) * payoff_avg.mean()
print(f"对偶变量法价格: {price_antithetic:.4f}")

---

第四节：C 实现对比

4.1 手写 RK4 对比

以下是在 C 中手写一个单步的经典四阶 Runge-Kutta（RK4）所需的核心函数：

// C: 手写 RK4 的 dydt（导数函数）
void dydt(double t, double *y, double *dydt_out, void *params) {
 double k = *(double*)params;
 dydt_out[0] = -k * y[0]; // dy/dt = -k*y
}
 
// C: 通用 RK4 一步
void rk4_step(void (*f)(double, double*, double*, void*),
 double t, double *y, double h, int n_eq,
 double *y_new, void *params) {
 double *k1 = malloc(n_eq * sizeof(double));
 double *k2 = malloc(n_eq * sizeof(double));
 double *k3 = malloc(n_eq * sizeof(double));
 double *k4 = malloc(n_eq * sizeof(double));
 double *ytmp = malloc(n_eq * sizeof(double));
 
 f(t, y, k1, params);
 
 for (int i = 0; i < n_eq; i++)
 ytmp[i] = y[i] + 0.5 * h * k1[i];
 f(t + 0.5*h, ytmp, k2, params);
 
 for (int i = 0; i < n_eq; i++)
 ytmp[i] = y[i] + 0.5 * h * k2[i];
 f(t + 0.5*h, ytmp, k3, params);
 
 for (int i = 0; i < n_eq; i++)
 ytmp[i] = y[i] + h * k3[i];
 f(t + h, ytmp, k4, params);
 
 for (int i = 0; i < n_eq; i++)
 y_new[i] = y[i] + (h/6.0) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);
 
 free(k1); free(k2); free(k3); free(k4); free(ytmp);
}

对比 Python 等价代码：

from scipy.integrate import solve_ivp
 
def dydt(t, y, k):
 return -k * y
 
sol = solve_ivp(dydt, [0, 10], [10.0], args=(0.5,), method='RK45')

Python 版本不仅更短，而且具有：

• 自适应步长（C 版本需要额外数百行代码实现）
• 事件检测
• 刚性自动检测（LSODA）
• 稠密插值输出

4.2 何时该用 C

场景	推荐
原型开发、实验	Python（SciPy）
单次运行、低维	Python
大规模 Monte Carlo	Python + Numba 或 C
实时嵌入式系统	C（无 Python 解释器可用）
大规模 PDE（如 CFD）	C/Fortran/C++（或 Python 调用 C 后端）
参数扫描（同一模型多次运行）	C（或 Python + 多进程/MPI）

关于 Python → C 性能桥接的完整讨论，见 性能对比和 pybind11 与 Cython。

---

练习

以下题目用于验证本章所学内容：

题号	题目	链接	涉及知识点
50	Pow(x, n)	https://leetcode.cn/problems/powx-n/	数值计算、迭代
69	x 的平方根	https://leetcode.cn/problems/sqrtx/	牛顿迭代、数值方法
509	斐波那契数	https://leetcode.cn/problems/fibonacci-number/	递推、数值序列