1. 常微分方程初值问题求解失败：从报错定位到方法选择

问题复现：一个“看起来很简单”的初值问题

考虑一阶常微分方程初值问题：

y' = -1000·(y - cos t) - sin t,   y(0) = 1,   t ∈ [0, 2]

该方程的解析解为 y(t) = cos t（可代入验证：右端 = -1000·(cos t - cos t) - sin t = -sin t = y'）。因此任何数值解与 cos t 的偏差都可以直接量化。

若直接用显式欧拉法或经典 RK4 以固定步长 h = 0.01 求解，常见现象是：解在 t 很小时就迅速偏离 cos t，甚至出现数值爆炸或求解器报“步长过小、无法满足容差”的错误。这不是代码写错，而是方法选择与方程特性不匹配。

三步定位：连续性、刚性、区间

第一步，检查右端函数连续性。 右端 f(t,y) = -1000(y - cos t) - sin t 对 t 和 y 都是多项式/三角函数组合，在全平面连续且对 y 满足 Lipschitz 条件（Lipschitz 常数 L = 1000）。所以报错不是间断点或不满足存在唯一性条件导致的，可以排除“右端函数病态”这一类原因。

第二步，判断刚性特征。 对 y 求偏导得 ∂f/∂y = -1000，即雅可比矩阵（一维情形就是标量）特征值为 λ = -1000。求解区间长度 T = 2，则 |λ|·T = 2000 ≫ 1，且 λ 实部为大幅负值——这是典型的刚性信号：存在一个衰减极快的瞬态分量 e^{-1000t}，它在 t ≈ 0.005 内就衰减殆尽，但显式方法的稳定性却要求步长满足 h ≤ 2/|λ| = 0.002（欧拉法稳定域边界），否则该瞬态分量会被数值放大。

第三步，检查求解区间。 区间 [0, 2] 相对瞬态时间尺度 1/1000 极长，意味着绝大部分区间上解已进入光滑的慢变阶段（cos t），但显式方法仍被稳定性约束“锁死”在极小步长上，计算代价与精度收益严重失衡。

可复现的修正方案

下面用 Python 的 scipy.integrate.solve_ivp 演示：先用显式方法复现失败，再切换到隐式方法得到正确解。代码自包含，可直接运行。

import numpy as np
from scipy.integrate import solve_ivp

def f(t, y):
    return -1000.0 * (y - np.cos(t)) - np.sin(t)

y0 = [1.0]
t_span = (0.0, 2.0)
t_eval = np.linspace(0.0, 2.0, 201)
exact = np.cos(t_eval)

# 1) 显式方法：固定步长 RK45 在刚性下会被迫压缩步长
sol_rk = solve_ivp(f, t_span, y0, method='RK45',
                   t_eval=t_eval, rtol=1e-6, atol=1e-9)
print("RK45 成功:", sol_rk.success, "| 消息:", sol_rk.message)
print("RK45 步数:", sol_rk.t.size, "| 最大误差:",
      np.max(np.abs(sol_rk.y[0] - exact)) if sol_rk.success else "N/A")

# 2) 隐式方法：BDF（A-稳定，适合刚性）
sol_bdf = solve_ivp(f, t_span, y0, method='BDF',
                    t_eval=t_eval, rtol=1e-6, atol=1e-9)
print("BDF 成功:", sol_bdf.success, "| 步数:", sol_bdf.t.size)
print("BDF 最大误差:", np.max(np.abs(sol_bdf.y[0] - exact)))

# 3) 隐式 Radau（L-稳定，对极刚性更稳）
sol_rad = solve_ivp(f, t_span, y0, method='Radau',
                    t_eval=t_eval, rtol=1e-6, atol=1e-9)
print("Radau 成功:", sol_rad.success, "| 步数:", sol_rad.t.size)
print("Radau 最大误差:", np.max(np.abs(sol_rad.y[0] - exact)))

预期现象（示例运行结果，具体数值随 SciPy 版本略有差异）：RK45 会报告步长被压缩到极小、步数远超 BDF/Radau，或在严格容差下直接失败；BDF 与 Radau 则以远少的步数完成求解，最大误差在 1e-6 量级，与 cos t 吻合。这验证了“刚性 → 隐式方法”的判断。

方法选择决策表

验证与通过条件

• 成功标志：求解器返回 success=True，且 max|y_num - cos t| 在设定容差（如 1e-5）以内。
• 步数对比：隐式方法步数应显著少于显式方法在同等精度下的步数，说明稳定性约束已解除。
• 区间端点检查：在 t=0 附近解应快速贴合 cos t，不应出现数值振荡或发散。
• 反例确认：若把 h 固定为 0.01 用显式欧拉，应观察到解在数步内偏离——这反向印证刚性判断正确。