欧拉方法的截断误差分析与龙格-库塔法的阶数提升

数值求解常微分方程初值问题 y' = f(t, y)，初始条件 y(t0) = y0，是科学计算的基础。欧拉方法（Euler's Method）是最直观的数值解法，但其精度往往不够高。本文将手把手分析欧拉方法的误差来源，并展示如何通过龙格-库塔法（Runge-Kutta Methods）系统地提升解的精度阶数。

1. 欧拉方法：理解与误差来源

理解欧拉方法的核心思想

欧拉方法的基本原理是用离散的线性近似来替代连续的曲线。在点 (tn, yn) 处，我们利用微分方程给出的信息——即该点的斜率 f(tn, yn)——来预测下一个点 (tn+1, yn+1) 的值。

1. 写出微分方程的标准形式：y' = f(t, y)，并明确初始条件 y(t0) = y0。
2. 确定步长 h，即时间 t 每次前进的固定长度。
3. 初始化：令 t0 和 y0 为初始值。
4. 迭代：对于从 n=0 开始的每一步，计算斜率 k1 = f(tn, yn)，然后更新：yn+1 = yn + h * k1，tn+1 = tn + h。

分析局部截断误差

局部截断误差是指，假设前一个点 yn 是精确的，在单步计算 yn+1 时所产生的误差。我们可以通过泰勒展开来精确分析。

假设精确解为 y(t)，在 t = tn 处进行泰勒展开：

y(tn+1) = y(tn + h) = y(tn) + h y'(tn) + (h^2 / 2) y''(ξn), 其中 ξn 在 tn 和 tn+1 之间。

由于 y'(tn) = f(tn, y(tn))，且欧拉方法的计算值为 yn+1 = yn + h f(tn, yn)。如果我们假设 yn = y(tn)（即上一步是精确的），那么将精确解与数值解相减，得到局部截断误差 τn：

τn = y(tn+1) - [yn + h f(tn, yn)]
τn = [y(tn) + h y'(tn) + (h^2 / 2) y''(ξn)] - [y(tn) + h y'(tn)]
τn = (h^2 / 2) y''(ξn)

核心结论：局部截断误差与步长 h 的平方成正比。我们称欧拉方法为一阶方法，因为其局部误差是 O(h^2)。这意味着，全局误差（所有步累积的误差）将是 O(h)，即误差随 h 的减小而线性减小。要显著提高精度，需要减小步长，但计算量会大幅增加。

2. 龙格-库塔法：提升精度的系统方法

既然误差来源于对解曲线的低阶（线性）近似，那么一个自然的改进思路就是使用更高阶的多项式来近似。龙格-库塔法正是这样一类方法，它们通过在单个步长 h 内计算多个点的斜率，然后进行加权平均，来获得一个更精确的“平均斜率”，从而实现更高阶的精度。

理解经典四阶龙格-库塔法 (RK4)

四阶龙格-库塔法是最常用的方法之一，其局部截断误差为 O(h^5)，全局误差为 O(h^4)，是一个四阶方法。

1. 计算第一个点处的斜率：k1 = f(tn, yn)
2. 计算中点处的斜率（使用 k1 预测）：k2 = f(tn + h/2, yn + (h/2) * k1)
3. 计算中点处的另一个斜率（使用 k2 预测）：k3 = f(tn + h/2, yn + (h/2) * k2)
4. 计算终点处的斜率（使用 k3 预测）：k4 = f(tn + h, yn + h * k3)
5. 加权平均四个斜率，以获得一个更优的估计斜率，然后更新：yn+1 = yn + (h/6) * (k1 + 2k2 + 2k3 + k4)，tn+1 = tn + h。

为什么它是四阶？

对比欧拉法只使用一个斜率 k1，RK4 通过精心设计的四个中间步骤和加权方案，实际上对精确解 y(tn+h) 进行了更高阶的泰勒展开匹配。其更新公式可以理解为对 y(tn+h) 的四阶泰勒多项式的一个优秀近似。

通用 s 级龙格-库塔格式

所有龙格-库塔法都可以用一个通用格式表示，其中 s 代表在一步内计算斜率的次数（称为级）。计算 s 个斜率 ki：

ki = f(tn + ci*h, yn + h * Σ_{j=1}^{i-1} a_{ij} * kj), for i = 1, ..., s

最终更新为：

yn+1 = yn + h * Σ_{i=1}^{s} b_i * k_i

其中，系数 a_{ij}， b_i， c_i 构成了一个“布彻表”（Butcher tableau）。通过求解一组复杂的多项式方程来确定这些系数，可以使方法达到期望的阶数。RK4 对应的是 s=4 的一组特定系数。

3. 从欧拉法到高阶龙格-库塔法的实践步骤

假设你需要求解 y' = -y + t + 1，y(0) = 1，在 t = 1 处的值，并比较不同方法的精度。精确解为 y(t) = e^{-t} + t。

步骤一：实现并测试欧拉方法

1. 设置步长 h = 0.2。
2. 初始化 t = 0, y = 1, tfinal = 1。
3. 执行循环直到 t >= tfinal：
   • 计算 slope = -y + t + 1。
   • 更新 y = y + h * slope。
   • 更新 t = t + h。
4. 记录结果，并与精确解 y(1) = e^{-1} + 1 ≈ 1.367879 比较。你会发现误差较大。

步骤二：实现并测试四阶龙格-库塔法 (RK4)

使用相同的 h = 0.2 和初始条件，替换核心更新步骤为 RK4 的 k1 到 k4 计算和加权平均。再次记录结果并与精确解比较，误差将显著减小。

步骤三：对比误差收敛速度

1. 选择一系列递减的步长，如 h = 0.1, 0.05, 0.025。
2. 分别使用欧拉法和 RK4 对每个 h 进行计算，得到数值解 y_h(1)。
3. 计算误差 error = |y_h(1) - y_exact|。
4. 观察：当 h 减半时，欧拉法的误差大约减半（符合一阶 O(h)），而 RK4 的误差大约减少为原来的十六分之一（符合四阶 O(h^4)）。这直观证明了阶数越高，精度随 h 减小而提升的速度越快。

步骤四：理解阶数与计算量的权衡

虽然高阶方法每一步计算量更大（RK4 是欧拉法的4倍），但为达到相同精度，它所需的时间步数 N = (tfinal - t0) / h 要少得多。综合来看，高阶方法通常更高效。选择方法时，需要根据问题的光滑性、稳定性要求和精度需求进行权衡。

4. 代码示例：直观感受差异

以下是用 Python 实现欧拉法和 RK4 的核心函数片段。

def euler(f, t0, y0, h, steps):
    t, y = t0, y0
    trajectory = [(t, y)]
    for _ in range(steps):
        slope = f(t, y)
        y += h * slope
        t += h
        trajectory.append((t, y))
    return trajectory

def rk4(f, t0, y0, h, steps):
    t, y = t0, y0
    trajectory = [(t, y)]
    for _ in range(steps):
        k1 = f(t, y)
        k2 = f(t + h/2, y + (h/2) * k1)
        k3 = f(t + h/2, y + (h/2) * k2)
        k4 = f(t + h, y + h * k3)
        y += (h / 6) * (k1 + 2*k2 + 2*k3 + k4)
        t += h
        trajectory.append((t, y))
    return trajectory

# 使用示例
def my_func(t, y):
    return -y + t + 1

t0, y0, h = 0, 1, 0.2
steps = int(1.0 / h)

euler_traj = euler(my_func, t0, y0, h, steps)
rk4_traj = rk4(my_func, t0, y0, h, steps)

# 比较最终点 t=1 时的 y 值
euler_final = euler_traj[-1][1]
rk4_final = rk4_traj[-1][1]
exact_final = np.exp(-1) + 1

print(f"欧拉法结果: {euler_final:.6f}, 误差: {abs(euler_final - exact_final):.6f}")
print(f"RK4 结果:   {rk4_final:.6f}, 误差: {abs(rk4_final - exact_final):.6f}")

运行此代码，你将直接观察到 RK4 结果比欧拉法结果更接近精确解 1.367879。通过调整 h 重复运行，可以进一步验证误差收敛速度与理论阶数相符。