4. 初值敏感与病态初值问题的求解诊断

问题方程与初值扰动设置

本节选用一个经典的对初值高度敏感的线性常微分方程：

y' = λ·y,  y(0) = y₀,  λ = 8

其解析解为 y(t) = y₀·e^(8t)。取两组仅相差 10⁻³ 的初值：

• 基准轨线：y₀ = 1.000，解析解 y(t) = e^(8t)
• 扰动轨线：y₀ = 1.001，解析解 y(t) = 1.001·e^(8t)

两条解析轨线的相对分离量为 |Δy|/|y| = 0.001/1.000 ≈ 10⁻³，且不随时间增长——这说明方程本身对初值扰动的相对误差是稳定的（λ>0 时绝对误差放大，但相对误差恒定）。真正需要警惕的是：当 λ 很大或为负且绝对值很大时，数值求解器的局部截断误差会被同样放大，导致数值解与解析解严重偏离。下面用前向欧拉法在固定步长 h=0.1 下求解，观察数值轨线的分离。

数值实验：不同初值下的轨线分离

用前向欧拉法 y_{n+1} = y_n + h·λ·y_n = (1 + hλ)·y_n 求解，放大因子 R = 1 + 0.1×8 = 1.8。解析放大因子为 e^(0.8) ≈ 2.2255。两者每步相差约 0.4255 倍，误差随步数指数累积。

import numpy as np

lam, h, T = 8.0, 0.1, 1.0
n_steps = int(T / h)
t = np.linspace(0, T, n_steps + 1)

def euler(y0):
    y = np.zeros(n_steps + 1)
    y[0] = y0
    for n in range(n_steps):
        y[n+1] = y[n] + h * lam * y[n]
    return y

y_base = euler(1.000)
y_pert = euler(1.001)
y_exact_base = np.exp(lam * t)
y_exact_pert = 1.001 * np.exp(lam * t)

# 数值解相对误差
err_base = np.abs(y_base - y_exact_base) / y_exact_base
err_pert = np.abs(y_pert - y_exact_pert) / y_exact_pert
print("t=1.0 数值解:", y_base[-1], y_pert[-1])
print("t=1.0 解析解:", y_exact_base[-1], y_exact_pert[-1])
print("相对误差:", err_base[-1], err_pert[-1])

预期结果（示例计算，非实测）：t=1.0 时数值解约为 1.8¹⁰ ≈ 357.05 与 357.41，解析解约为 e⁸ ≈ 2980.96 与 2983.94。数值解严重偏低，相对误差约 88%，且两组初值的数值解几乎重合——求解器误差远大于初值扰动带来的差异，这正是需要区分的核心现象。

区分方程不稳定与求解器误差放大

判断依据可以归纳为三条可操作的检查：

1. 解析相对分离量：对线性方程 y'=λy，解析相对差恒为 |Δy₀/y₀|，与 t 无关。若数值解的相对差随时间显著增长，说明误差来自离散化而非方程。
2. 步长减半检验：将 h 从 0.1 减到 0.05，前向欧拉放大因子变为 1+0.05×8=1.4，与 e^0.4≈1.4918 更接近。若数值解随 h 减小而收敛到解析解，则误差是求解器截断误差。
3. 方法切换对照：改用四阶龙格-库塔（RK4）或隐式方法，若结果显著改善，进一步确认是显式低阶方法的稳定性/精度不足。

def rk4(y0):
    y = np.zeros(n_steps + 1)
    y[0] = y0
    for n in range(n_steps):
        k1 = lam * y[n]
        k2 = lam * (y[n] + 0.5*h*k1)
        k3 = lam * (y[n] + 0.5*h*k2)
        k4 = lam * (y[n] + h*k3)
        y[n+1] = y[n] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
    return y

y_rk4 = rk4(1.000)
print("RK4 t=1.0:", y_rk4[-1], "解析:", y_exact_base[-1])
print("RK4 相对误差:", abs(y_rk4[-1]-y_exact_base[-1])/y_exact_base[-1])

预期 RK4 在 h=0.1 时相对误差约 10⁻⁴ 量级，远优于前向欧拉的 88%。这确认了原问题中数值解偏离的主因是显式欧拉法在 hλ=0.8 时的截断误差累积，而非方程对初值敏感。

可信求解的判定标准