常微分方程（十一）数值方法
Chen Kai BOSS
  2019-05-29 15:30:00 2019-05-29 15:30      5.4k 字  25 分钟  

很多实际问题中的微分方程无法用解析方法求解——这时候数值方法就成了我们最强大的武器。从 Euler 的简单想法到现代自适应算法，数值方法让我们能够"近似地"解决几乎任何微分方程。本章将深入探讨各种数值方法的原理、实现和误差分析。

为什么需要数值方法？

回顾前面的章节，我们学习了许多解析方法：分离变量、积分因子、常数变易法、 Laplace 变换等。这些方法很优美，但有一个致命的限制——只对特定类型的方程有效。

欧拉方法的几何解释

考虑这个看似简单的方程：

试着用之前学过的任何方法解它——你会发现根本行不通！这个方程没有初等函数形式的解。

更一般地，工程和科学中遇到的大多数微分方程都无法解析求解：

  • 非线性方程（如 Navier-Stokes 方程）
  • 变系数方程
  • 高维系统
  • 带有复杂边界条件的问题

数值方法的核心思想：既然无法得到精确解，不如在离散点 上计算近似值。

Euler 方法：最简单的数值方法

基本思想

Runge-Kutta方法对比

Euler 方法是所有数值方法的鼻祖，由 Leonhard Euler 在 1768 年提出。它的思想极其简单：用切线近似曲线。

考虑初值问题：

在点 处，解曲线的切线斜率是。沿着这条切线走一小步，我们得到：

这就是向前 Euler 公式（ Forward Euler）。

重复这个过程：

几何直觉

想象你在迷宫中行走，但只能看到脚下一小块地面。你的策略是：

  1. 看看当前位置的地形坡度（导数）
  2. 沿着这个方向走一小步
  3. 重复

这就是 Euler 方法的本质——局部线性化。

Python 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
import numpy as np
import matplotlib.pyplot as plt

def euler_method(f, x0, y0, x_end, h):
    """
    向前 Euler 方法
    
    参数:
        f: 函数 f(x, y)，表示 dy/dx = f(x, y)
        x0, y0: 初始条件
        x_end: 积分终点
        h: 步长
    
    返回:
        x_values, y_values: 数值解
    """
    x_values = [x0]
    y_values = [y0]
    
    x, y = x0, y0
    while x < x_end:
        y = y + h * f(x, y)
        x = x + h
        x_values.append(x)
        y_values.append(y)
    
    return np.array(x_values), np.array(y_values)

# 示例： dy/dx = -2y, y(0) = 1
# 解析解： y = exp(-2x)
f = lambda x, y: -2 * y
x0, y0 = 0, 1
x_end = 3

# 不同步长的比较
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左图：不同步长的 Euler 解
step_sizes = [0.5, 0.2, 0.1, 0.05]
x_exact = np.linspace(0, x_end, 100)
y_exact = np.exp(-2 * x_exact)

axes[0].plot(x_exact, y_exact, 'k-', linewidth=2, label='精确解')
for h in step_sizes:
    x_num, y_num = euler_method(f, x0, y0, x_end, h)
    axes[0].plot(x_num, y_num, 'o-', markersize=4, label=f'h = {h}')

axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('Euler 方法：不同步长的比较')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 右图：误差随步长的变化
errors = []
step_sizes_fine = np.logspace(-3, 0, 20)
for h in step_sizes_fine:
    x_num, y_num = euler_method(f, x0, y0, x_end, h)
    y_exact_at_end = np.exp(-2 * x_end)
    error = abs(y_num[-1] - y_exact_at_end)
    errors.append(error)

axes[1].loglog(step_sizes_fine, errors, 'bo-', label='实际误差')
axes[1].loglog(step_sizes_fine, step_sizes_fine, 'r--', label='O(h) 参考线')
axes[1].set_xlabel('步长 h')
axes[1].set_ylabel('终点误差')
axes[1].set_title('Euler 方法的误差阶数')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('euler_method.png', dpi=150, bbox_inches='tight')
plt.show()

误差分析

Euler 方法的误差分为两类：

局部截断误差（ LTE）：单步引入的误差

通过 Taylor 展开：

而 Euler 公式给出，所以：

全局截断误差（ GTE）：从 到 的累积误差

经过 步后：

这意味着 Euler 方法是一阶方法：步长减半，误差也减半。

稳定性分析

考虑测试方程，其中（衰减问题）。

Euler 公式给出：

要使数值解不发散，需要，即：

由于，稳定性条件是：

这称为条件稳定性——步长必须足够小才能保证稳定。

向后 Euler 方法

隐式格式

自适应步长控制

向后 Euler 方法（ Backward Euler）使用下一时刻的导数值：

注意 同时出现在等式两边——这是一个隐式方程！

为什么需要隐式方法？

对于刚性问题（ stiff problems），显式方法需要极小的步长才能稳定，而隐式方法可以使用大步长。

刚性问题的特点：系统中同时存在快变和慢变的成分。例如化学反应动力学中，有些反应在微秒内完成，有些需要几小时。

稳定性分析

对测试方程，向后 Euler 给出：

放大因子是。当 时，对任意都有放大因子。

这称为A-稳定（ A-stable）或无条件稳定。

Python 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def backward_euler(f, df_dy, x0, y0, x_end, h, tol=1e-10, max_iter=100):
    """
    向后 Euler 方法（使用 Newton 迭代求解隐式方程）
    
    参数:
        f: 函数 f(x, y)
        df_dy: f 对 y 的偏导数
        x0, y0: 初始条件
        x_end: 积分终点
        h: 步长
        tol: Newton 迭代容差
        max_iter: 最大迭代次数
    """
    x_values = [x0]
    y_values = [y0]
    
    x, y = x0, y0
    while x < x_end:
        x_new = x + h
        
        # Newton 迭代求解 y_new = y + h*f(x_new, y_new)
        y_new = y  # 初始猜测
        for _ in range(max_iter):
            g = y_new - y - h * f(x_new, y_new)
            dg = 1 - h * df_dy(x_new, y_new)
            y_new_next = y_new - g / dg
            
            if abs(y_new_next - y_new) < tol:
                y_new = y_new_next
                break
            y_new = y_new_next
        
        x, y = x_new, y_new
        x_values.append(x)
        y_values.append(y)
    
    return np.array(x_values), np.array(y_values)

# 刚性问题示例： y' = -1000(y - cos(x)) - sin(x), y(0) = 1
# 这个问题有快速衰减到 cos(x)的瞬态行为
def f_stiff(x, y):
    return -1000 * (y - np.cos(x)) - np.sin(x)

def df_dy_stiff(x, y):
    return -1000

# 比较显式和隐式方法
x0, y0 = 0, 1
x_end = 0.1

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 向前 Euler（可能不稳定）
try:
    h_forward = 0.001  # 小步长才稳定
    x_fwd, y_fwd = euler_method(f_stiff, x0, y0, x_end, h_forward)
    axes[0].plot(x_fwd, y_fwd, 'b-', label=f'向前 Euler (h={h_forward})')
except:
    pass

# 向后 Euler（大步长也稳定）
h_back = 0.01
x_bwd, y_bwd = backward_euler(f_stiff, df_dy_stiff, x0, y0, x_end, h_back)
axes[0].plot(x_bwd, y_bwd, 'ro-', label=f'向后 Euler (h={h_back})')

# 精确解
x_exact = np.linspace(x0, x_end, 200)
y_exact = np.cos(x_exact) + (y0 - 1) * np.exp(-1000 * x_exact)
axes[0].plot(x_exact, y_exact, 'k--', label='精确解')

axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('刚性问题：显式 vs 隐式方法')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

改进的 Euler 方法（ Heun 方法）

预测-校正思想

刚性方程的数值求解

Euler 方法用起点的斜率，但曲线的斜率是在变化的。能不能用平均斜率？

Heun 方法（也称改进的 Euler 方法）：

  1. 预测：用 Euler 公式得到

  2. 校正：用起点和预测点的平均斜率

几何解释

想象你要从 A 点走到 B 点：

  • Euler 方法：按 A 点的方向一直走
  • Heun 方法：先按 A 点方向试探到 A'，然后取 A 和 A'方向的平均

误差阶数

Heun 方法是二阶方法：

步长减半，误差减为四分之一！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def heun_method(f, x0, y0, x_end, h):
    """
    Heun 方法（改进的 Euler 方法）
    """
    x_values = [x0]
    y_values = [y0]
    
    x, y = x0, y0
    while x < x_end:
        k1 = f(x, y)
        y_pred = y + h * k1
        k2 = f(x + h, y_pred)
        
        y = y + h * (k1 + k2) / 2
        x = x + h
        
        x_values.append(x)
        y_values.append(y)
    
    return np.array(x_values), np.array(y_values)

Runge-Kutta 方法

RK4：工程师的瑞士军刀

数值稳定性分析

Runge-Kutta 方法家族是数值 ODE 求解的主力。其中最著名的是四阶 Runge-Kutta 方法（ RK4）：

直觉解释

RK4 在区间 上采样 4 个斜率：

-：左端点斜率 -：中点斜率（用预测） -：中点斜率（用预测，更准确） -：右端点斜率

然后加权平均： 这类似于Simpson 积分公式的权重！

Python 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
def rk4_method(f, x0, y0, x_end, h):
    """
    经典四阶 Runge-Kutta 方法
    """
    x_values = [x0]
    y_values = [y0]
    
    x, y = x0, y0
    while x < x_end:
        k1 = f(x, y)
        k2 = f(x + h/2, y + h*k1/2)
        k3 = f(x + h/2, y + h*k2/2)
        k4 = f(x + h, y + h*k3)
        
        y = y + h * (k1 + 2*k2 + 2*k3 + k4) / 6
        x = x + h
        
        x_values.append(x)
        y_values.append(y)
    
    return np.array(x_values), np.array(y_values)

# 比较 Euler 、 Heun 和 RK4
f = lambda x, y: np.sin(x * y)  # 非线性方程
x0, y0 = 0, 1
x_end = 5
h = 0.2

# 参考解（用很小的步长）
from scipy.integrate import solve_ivp
sol = solve_ivp(lambda x, y: np.sin(x * y[0]), [x0, x_end], [y0], 
                dense_output=True, method='RK45', rtol=1e-10)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

x_ref = np.linspace(x0, x_end, 200)
y_ref = sol.sol(x_ref)[0]

# 左图：三种方法的解
x_euler, y_euler = euler_method(f, x0, y0, x_end, h)
x_heun, y_heun = heun_method(f, x0, y0, x_end, h)
x_rk4, y_rk4 = rk4_method(f, x0, y0, x_end, h)

axes[0].plot(x_ref, y_ref, 'k-', linewidth=2, label='参考解')
axes[0].plot(x_euler, y_euler, 'b.-', label=f'Euler (h={h})')
axes[0].plot(x_heun, y_heun, 'g.-', label=f'Heun (h={h})')
axes[0].plot(x_rk4, y_rk4, 'r.-', label=f'RK4 (h={h})')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('三种方法的比较')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 右图：误差阶数比较
step_sizes = np.logspace(-2.5, -0.5, 15)
euler_errors, heun_errors, rk4_errors = [], [], []

for h in step_sizes:
    _, y_e = euler_method(f, x0, y0, x_end, h)
    _, y_h = heun_method(f, x0, y0, x_end, h)
    _, y_r = rk4_method(f, x0, y0, x_end, h)
    
    y_true = sol.sol(x_end)[0]
    euler_errors.append(abs(y_e[-1] - y_true))
    heun_errors.append(abs(y_h[-1] - y_true))
    rk4_errors.append(abs(y_r[-1] - y_true))

axes[1].loglog(step_sizes, euler_errors, 'b.-', label='Euler O(h)')
axes[1].loglog(step_sizes, heun_errors, 'g.-', label='Heun O(h ²)')
axes[1].loglog(step_sizes, rk4_errors, 'r.-', label='RK4 O(h ⁴)')
axes[1].loglog(step_sizes, step_sizes, 'b--', alpha=0.5)
axes[1].loglog(step_sizes, step_sizes**2, 'g--', alpha=0.5)
axes[1].loglog(step_sizes, step_sizes**4, 'r--', alpha=0.5)
axes[1].set_xlabel('步长 h')
axes[1].set_ylabel('误差')
axes[1].set_title('误差阶数比较')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rk_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

RK4 的推导

RK 方法的一般形式是：

其中：

参数、、通过阶条件确定——要求数值解的 Taylor 展开与真解匹配到某一阶。

对于 RK4，需要满足的条件有 11 个，而自由参数只有 13 个，所以解不唯一。经典 RK4 只是其中一种选择。

多步法

思想：利用历史信息

Runge-Kutta 方法是单步法——每一步只用到。但我们已经计算了那么多历史点，为什么不利用呢？

多步法使用多个历史点来提高精度或效率。

Adams-Bashforth 方法（显式）

二步 Adams-Bashforth（ AB2）：

四步 Adams-Bashforth（ AB4）：

Adams-Moulton 方法（隐式）

二步 Adams-Moulton（ AM2）：

这是隐式的，因为 包含未知的。

预测-校正法

实践中常用预测-校正组合：

  1. 用 Adams-Bashforth预测${n+1}计算{n+1} = f(x_{n+1}, {n+1})用校正得到y{n+1}$ 
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    def adams_bashforth_moulton(f, x0, y0, x_end, h):
        """
        Adams-Bashforth-Moulton 预测-校正法（ 4 阶）
        使用 RK4 启动
        """
        # 用 RK4 计算前 4 个点
        x_values = [x0]
        y_values = [y0]
        f_values = [f(x0, y0)]
        
        x, y = x0, y0
        for _ in range(3):
            k1 = f(x, y)
            k2 = f(x + h/2, y + h*k1/2)
            k3 = f(x + h/2, y + h*k2/2)
            k4 = f(x + h, y + h*k3)
            
            y = y + h * (k1 + 2*k2 + 2*k3 + k4) / 6
            x = x + h
            
            x_values.append(x)
            y_values.append(y)
            f_values.append(f(x, y))
        
        # Adams-Bashforth-Moulton
        while x < x_end:
            # AB4 预测
            y_pred = y_values[-1] + h/24 * (
                55*f_values[-1] - 59*f_values[-2] + 
                37*f_values[-3] - 9*f_values[-4]
            )
            
            x_new = x_values[-1] + h
            f_pred = f(x_new, y_pred)
            
            # AM4 校正
            y_corr = y_values[-1] + h/24 * (
                9*f_pred + 19*f_values[-1] - 
                5*f_values[-2] + f_values[-3]
            )
            
            x_values.append(x_new)
            y_values.append(y_corr)
            f_values.append(f(x_new, y_corr))
        
        return np.array(x_values), np.array(y_values)

启动问题

多步法需要多个初始点，但初值问题只给出一个。解决方案：

  1. 用单步法（如 RK4）计算前几个点
  2. 然后切换到多步法

自适应步长控制

问题：如何选择步长？

  • 步长太大：误差不可接受
  • 步长太小：计算太慢

自适应步长的思想：让算法自己根据局部误差调整步长。

嵌入式 Runge-Kutta 方法

Runge-Kutta-Fehlberg（ RKF45）使用同一组 值同时计算 4 阶和 5 阶近似：

误差估计：

步长调整策略

给定容差，如果，减小步长并重算；如果，可以增大步长。

新步长的估计：

实际中会加一个安全系数（如 0.9）和上下限。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def rkf45_adaptive(f, x0, y0, x_end, tol=1e-6, h_init=0.1):
    """
    Runge-Kutta-Fehlberg 4(5)自适应步长方法
    """
    # Butcher 表系数
    a2, a3, a4, a5, a6 = 1/4, 3/8, 12/13, 1, 1/2
    
    b21 = 1/4
    b31, b32 = 3/32, 9/32
    b41, b42, b43 = 1932/2197, -7200/2197, 7296/2197
    b51, b52, b53, b54 = 439/216, -8, 3680/513, -845/4104
    b61, b62, b63, b64, b65 = -8/27, 2, -3544/2565, 1859/4104, -11/40
    
    # 4 阶和 5 阶的权重
    c4 = [25/216, 0, 1408/2565, 2197/4104, -1/5, 0]
    c5 = [16/135, 0, 6656/12825, 28561/56430, -9/50, 2/55]
    
    x_values = [x0]
    y_values = [y0]
    
    x, y = x0, y0
    h = h_init
    
    while x < x_end:
        if x + h > x_end:
            h = x_end - x
        
        # 计算 k 值
        k1 = f(x, y)
        k2 = f(x + a2*h, y + h*b21*k1)
        k3 = f(x + a3*h, y + h*(b31*k1 + b32*k2))
        k4 = f(x + a4*h, y + h*(b41*k1 + b42*k2 + b43*k3))
        k5 = f(x + a5*h, y + h*(b51*k1 + b52*k2 + b53*k3 + b54*k4))
        k6 = f(x + a6*h, y + h*(b61*k1 + b62*k2 + b63*k3 + b64*k4 + b65*k5))
        
        ks = [k1, k2, k3, k4, k5, k6]
        
        # 4 阶和 5 阶近似
        y4 = y + h * sum(c4[i]*ks[i] for i in range(6))
        y5 = y + h * sum(c5[i]*ks[i] for i in range(6))
        
        # 误差估计
        err = abs(y5 - y4)
        
        if err < 1e-15:
            err = 1e-15  # 避免除零
        
        # 步长调整
        if err <= tol:
            # 接受这一步
            x = x + h
            y = y5  # 用更高阶的结果
            x_values.append(x)
            y_values.append(y)
        
        # 计算新步长
        h_new = 0.9 * h * (tol / err) ** 0.2
        h = min(max(h_new, h/10), h*2)  # 限制变化幅度
    
    return np.array(x_values), np.array(y_values)

微分方程组的数值解

向量化

高阶方程和方程组都可以写成一阶向量方程：

所有前面的方法都可以直接推广。

示例：洛伦兹系统

洛伦兹系统是混沌理论的经典例子：

参数,, 时产生著名的"蝴蝶效应"。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
from mpl_toolkits.mplot3d import Axes3D

def lorenz_system(state, t, sigma=10, rho=28, beta=8/3):
    x, y, z = state
    return np.array([
        sigma * (y - x),
        x * (rho - z) - y,
        x * y - beta * z
    ])

def rk4_system(f, t0, y0, t_end, h, args=()):
    """
    向量版 RK4
    """
    t_values = [t0]
    y_values = [np.array(y0)]
    
    t = t0
    y = np.array(y0)
    
    while t < t_end:
        k1 = f(y, t, *args)
        k2 = f(y + h*k1/2, t + h/2, *args)
        k3 = f(y + h*k2/2, t + h/2, *args)
        k4 = f(y + h*k3, t + h, *args)
        
        y = y + h * (k1 + 2*k2 + 2*k3 + k4) / 6
        t = t + h
        
        t_values.append(t)
        y_values.append(y.copy())
    
    return np.array(t_values), np.array(y_values)

# 求解洛伦兹系统
t0 = 0
y0 = [1, 1, 1]
t_end = 50
h = 0.01

t, y = rk4_system(lorenz_system, t0, y0, t_end, h)

# 3D 绘图
fig = plt.figure(figsize=(12, 5))

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot(y[:, 0], y[:, 1], y[:, 2], lw=0.5)
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('Z')
ax1.set_title('洛伦兹吸引子')

# 敏感性测试
y0_perturbed = [1.001, 1, 1]  # 微小扰动
_, y_pert = rk4_system(lorenz_system, t0, y0_perturbed, t_end, h)

ax2 = fig.add_subplot(122)
ax2.plot(t, y[:, 0], 'b-', label='原始')
ax2.plot(t, y_pert[:, 0], 'r--', label='扰动 (Δ x ₀=0.001)')
ax2.set_xlabel('时间')
ax2.set_ylabel('x(t)')
ax2.set_title('蝴蝶效应：初值敏感性')
ax2.legend()
ax2.set_xlim([0, 30])

plt.tight_layout()
plt.savefig('lorenz.png', dpi=150, bbox_inches='tight')
plt.show()

辛积分与几何积分

问题：长时间模拟

对于 Hamiltonian 系统（如天体力学），传统方法会导致能量漂移：

辛 Euler 方法

显式辛 Euler：

注意这里用 更新——这个"半隐式"结构保持了辛结构。

St ö rmer-Verlet 方法

更常用的是蛙跳法（ Leapfrog）或Verlet 方法：

这是二阶方法，且保持辛结构。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def verlet_harmonic_oscillator(omega, q0, p0, t_end, h):
    """
    Verlet 方法求解简谐振子
    H = p ²/2 + ω² q ²/2
    """
    t_values = [0]
    q_values = [q0]
    p_values = [p0]
    energy_values = [p0**2/2 + omega**2*q0**2/2]
    
    q, p = q0, p0
    t = 0
    
    while t < t_end:
        # 半步更新动量
        p_half = p - h/2 * omega**2 * q
        # 全步更新位置
        q = q + h * p_half
        # 半步更新动量
        p = p_half - h/2 * omega**2 * q
        
        t = t + h
        t_values.append(t)
        q_values.append(q)
        p_values.append(p)
        energy_values.append(p**2/2 + omega**2*q**2/2)
    
    return (np.array(t_values), np.array(q_values), 
            np.array(p_values), np.array(energy_values))

# 比较 RK4 和 Verlet 的能量守恒
omega = 1
q0, p0 = 1, 0
t_end = 100
h = 0.1

# Verlet
t_v, q_v, p_v, E_v = verlet_harmonic_oscillator(omega, q0, p0, t_end, h)

# RK4
def harmonic_ode(state, t, omega):
    q, p = state
    return np.array([p, -omega**2 * q])

t_rk, y_rk = rk4_system(harmonic_ode, 0, [q0, p0], t_end, h, args=(omega,))
E_rk = y_rk[:, 1]**2/2 + omega**2 * y_rk[:, 0]**2/2

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 相空间轨迹
axes[0].plot(q_v, p_v, 'b-', alpha=0.7, label='Verlet')
axes[0].plot(y_rk[:, 0], y_rk[:, 1], 'r--', alpha=0.7, label='RK4')
circle = plt.Circle((0, 0), 1, fill=False, color='k', linestyle=':')
axes[0].add_patch(circle)
axes[0].set_xlabel('q')
axes[0].set_ylabel('p')
axes[0].set_title('相空间轨迹')
axes[0].legend()
axes[0].set_aspect('equal')

# 能量随时间
E0 = p0**2/2 + omega**2*q0**2/2
axes[1].plot(t_v, (E_v - E0)/E0, 'b-', label='Verlet')
axes[1].plot(t_rk, (E_rk - E0)/E0, 'r-', label='RK4')
axes[1].set_xlabel('时间')
axes[1].set_ylabel('相对能量误差')
axes[1].set_title('能量守恒性比较')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('symplectic.png', dpi=150, bbox_inches='tight')
plt.show()

收敛性与稳定性理论

收敛性

定义：数值方法收敛，如果当 时，。

Lax 等价定理：对于线性差分格式，一致性 + 稳定性 = 收敛性。

稳定性区域

对测试方程，定义。

方法的稳定性区域是使得 的集合，其中 是放大因子。

  • Euler 方法：，稳定区域是圆盘
  • RK4：

A-稳定性

A-稳定：整个左半平面都在稳定区域内。

遗憾的是，不存在 A-稳定的显式 Runge-Kutta 方法（ Dahlquist 障碍）。

使用 SciPy 求解 ODE

实际工作中，我们通常使用成熟的库：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
from scipy.integrate import solve_ivp, odeint

# 方法 1： odeint（旧 API，但简单）
def f(y, t):
    return -2 * y

t = np.linspace(0, 5, 100)
y0 = 1
sol_odeint = odeint(f, y0, t)

# 方法 2： solve_ivp（新 API，功能更强）
def f_ivp(t, y):  # 注意参数顺序不同！
    return -2 * y

sol_ivp = solve_ivp(f_ivp, [0, 5], [y0], t_eval=t, method='RK45')

# 可用方法：
# 'RK45': 自适应 Runge-Kutta 4(5)
# 'RK23': 自适应 Runge-Kutta 2(3)
# 'DOP853': 高阶 Dormand-Prince
# 'Radau': 隐式 Runge-Kutta（刚性问题）
# 'BDF': 后向差分公式（刚性问题）
# 'LSODA': 自动检测刚性

# 刚性问题示例
def van_der_pol(t, y, mu=1000):
    x, v = y
    return [v, mu * (1 - x**2) * v - x]

sol_stiff = solve_ivp(van_der_pol, [0, 3000], [2, 0], 
                       method='Radau', dense_output=True)

t_plot = np.linspace(0, 3000, 1000)
y_plot = sol_stiff.sol(t_plot)

plt.figure(figsize=(12, 4))
plt.subplot(121)
plt.plot(t_plot, y_plot[0])
plt.xlabel('t')
plt.ylabel('x')
plt.title('Van der Pol 振子 (μ=1000)')

plt.subplot(122)
plt.plot(y_plot[0], y_plot[1])
plt.xlabel('x')
plt.ylabel('dx/dt')
plt.title('相图')

plt.tight_layout()
plt.show()

实际应用：疫情模型的数值模拟

SIR 模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def sir_model(t, y, beta, gamma, N):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return [dSdt, dIdt, dRdt]

# 参数设置（以 COVID-19 早期数据为参考）
N = 1e6  # 总人口
I0 = 1   # 初始感染者
R0_ratio = 2.5  # 基本再生数
gamma = 1/14  # 康复率（平均 14 天康复）
beta = R0_ratio * gamma  # 传染率

# 求解
sol = solve_ivp(sir_model, [0, 200], [N-I0, I0, 0],
                args=(beta, gamma, N),
                t_eval=np.linspace(0, 200, 500),
                method='RK45')

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].plot(sol.t, sol.y[0]/1e3, 'b-', label='易感者 S')
axes[0].plot(sol.t, sol.y[1]/1e3, 'r-', label='感染者 I')
axes[0].plot(sol.t, sol.y[2]/1e3, 'g-', label='康复者 R')
axes[0].set_xlabel('天数')
axes[0].set_ylabel('人数 (千人)')
axes[0].set_title(f'SIR 模型 (R ₀ = {R0_ratio})')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 不同 R0 的比较
R0_values = [1.5, 2.0, 2.5, 3.0]
for R0 in R0_values:
    beta = R0 * gamma
    sol = solve_ivp(sir_model, [0, 200], [N-I0, I0, 0],
                    args=(beta, gamma, N),
                    t_eval=np.linspace(0, 200, 500))
    axes[1].plot(sol.t, sol.y[1]/1e3, label=f'R ₀ = {R0}')

axes[1].set_xlabel('天数')
axes[1].set_ylabel('感染者数 (千人)')
axes[1].set_title('基本再生数的影响')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('sir_model.png', dpi=150, bbox_inches='tight')
plt.show()

练习题

基础题

练习 1：用 Euler 方法求解, 在 上的近似解。取，与精确解比较。

练习 2：实现 Heun 方法，求解,。验证它是二阶方法。

练习 3：证明 RK4 方法对（ 为常数）是精确的。

练习 4：考虑刚性方程,。 (a) 用向前 Euler 方法，找到保证稳定的最大步长 (b) 用向后 Euler 方法，验证对任意步长都稳定

练习 5：实现 Adams-Bashforth 二步法，并与 RK4 比较计算效率。

中级题

练习 6：实现自适应步长 RK4 方法，使局部误差保持在指定容差内。测试在 上的表现。

练习 7：用 Verlet 方法模拟双摆系统。观察混沌行为。

练习 8：考虑 Lotka-Volterra 捕食者-猎物模型：(a) 用 RK4 方法模拟种群动态 (b) 绘制相图，观察周期轨道

练习 9：用数值方法验证：对于简谐振子，辛积分器（如 Verlet）长时间模拟的能量守恒性优于 RK4 。

练习 10：实现 Dormand-Prince 方法（ RK5(4)），并与 SciPy 的solve_ivp比较。

高级题

练习 11：研究数值方法在刚性问题上的表现。考虑：

精确解是。 (a) 用显式方法，探索稳定性限制 (b) 实现隐式梯形法，观察无条件稳定性

练习 12：考虑二体问题（开普勒问题）：

用不同的方法（ RK4 、辛方法）长时间模拟轨道，比较： (a) 轨道形状的保持 (b) 能量的守恒 (c) 角动量的守恒

练习 13：实现指数积分器（ exponential integrator）用于求解，其中 是刚性线性部分。

练习 14：研究误差的积累。对于混沌系统（如洛伦兹系统），误差会指数增长。设计实验量化这种增长率。

练习 15：实现一个简单的自动刚性检测算法：根据 Jacobian 的特征值判断问题是否刚性，并自动选择合适的求解器。

总结

本章我们学习了常微分方程数值方法的完整体系：

方法 阶数 类型 特点
向前 Euler 1 显式 简单，条件稳定
向后 Euler 1 隐式 A-稳定，适合刚性问题
Heun 2 显式 预测-校正
RK4 4 显式 工程首选
RKF45 4(5) 显式 自适应步长
Adams-Bashforth 可变 显式多步 高效
Adams-Moulton 可变 隐式多步 更稳定
Verlet 2 辛 保结构

选择方法的原则：

  1. 一般问题： RK45（自适应）
  2. 刚性问题： BDF 或 Radau
  3. 长时间 Hamilton 系统：辛积分器
  4. 高精度需求：高阶嵌入式 RK

实际建议：

  • 先用scipy.integrate.solve_ivp试试
  • 如果太慢或不稳定，考虑换方法
  • 总是检查收敛性：减小步长看结果是否变化

下一章，我们将学习边值问题——当条件不仅在初始点，还在终点给出时如何求解。

  • 本文标题：常微分方程（十一）数值方法
  • 本文作者：Chen Kai
  • 创建时间：2019-05-29 15:30:00
  • 本文链接：https://www.chenk.top/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%EF%BC%88%E5%8D%81%E4%B8%80%EF%BC%89%E6%95%B0%E5%80%BC%E6%96%B9%E6%B3%95/
  • 版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
 评论