3. 步长过大导致数值解振荡的排查与收敛性验证

问题复现：衰减方程为何出现锯齿振荡

取一个解析解已知的衰减型初值问题作为排障对象：

y' = -20·y,  y(0) = 1,  t ∈ [0, 1]
解析解：y(t) = e^(-20t)，单调衰减，无振荡

用经典四阶显式 Runge–Kutta（RK4）以固定步长 h = 0.25 求解，得到如下数值序列（示例计算值）：

数值解符号交替、幅值按约 2/3 比例衰减，与真实解的单调衰减完全不符。这不是方程本身的问题，而是步长越过了显式方法的稳定边界。

稳定性区域分析：定位步长越界

对测试方程 y' = λy（此处 λ = -20），RK4 的稳定函数为：

R(z) = 1 + z + z²/2 + z³/6 + z⁴/24,  z = λh

实轴上的稳定区间为 z ∈ [-2.785, 0]。代入 λ = -20：

• h = 0.25 → z = -5.0，落在稳定区间之外，|R(z)| ≈ 1.44 > 1，误差被逐点放大并交替变号，形成锯齿。
• 稳定上限 h_max = 2.785 / 20 ≈ 0.139。取 h = 0.125 → z = -2.5，|R(z)| ≈ 0.70 < 1，衰减正常。

缩小步长并做步长减半收敛性验证

把步长降到稳定区间内，并用“步长减半、比较相邻两级解”的方式验证收敛，而不是只看单次结果是否好看。以下脚本可直接运行（仅依赖标准库）：

import math

lam = -20.0
def f(t, y): return lam * y
def exact(t): return math.exp(lam * t)

def rk4(h, T=1.0):
    n = int(round(T / h))
    t, y = 0.0, 1.0
    ys = [y]
    for _ in range(n):
        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
        ys.append(y)
    return ys

def max_err(h):
    ys = rk4(h)
    n = len(ys) - 1
    return max(abs(ys[i] - exact(i*h)) for i in range(n+1))

for h in [0.25, 0.125, 0.0625, 0.03125]:
    e = max_err(h)
    print(f"h={h:<8} max|err|={e:.3e}")

# 步长减半收敛比
e1, e2 = max_err(0.125), max_err(0.0625)
print("收敛比 e(h)/e(h/2) =", round(e1/e2, 2))

预期输出（示例计算值）：h = 0.25 时误差量级约 10⁻¹，明显异常；h = 0.125 起误差迅速下降，相邻两级误差比接近 2⁴ = 16，符合 RK4 的四阶收敛特征。若收敛比远小于 16 或误差不降，说明步长仍在稳定边界附近或方程存在刚性，需要继续缩小步长或改用隐式方法。