如何避免Scipy solve_ivp大求解域误差增大及量化误差

一、如何避免长时间积分的误差累积 #

针对你遇到的t>55后数值解偏离解析解的问题，可从以下方向优化：

1. 收紧求解器精度参数 #

scipy.integrate.solve_ivp默认的相对误差rtol=1e-3、绝对误差atol=1e-6，长时间积分时局部截断误差的累积会逐渐放大。尝试将精度参数调至更严格级别，强制求解器使用更精细步长，减少每一步的截断误差：

sol = solve_ivp(fun, t_span, y0, method='RK45', rtol=1e-8, atol=1e-10)

2. 变量替换转化为稳定方程 #

原方程为：
$$h'(t) = -\frac{2}{h(t)-\lambda(t)}$$
令y(t) = h(t) - \lambda(t)，则y'(t) = h'(t) - \lambda'(t)，代入原方程可得：
$$y'(t) = -\frac{2}{y(t)} - \lambda'(t)$$
当\lambda(t)=t时，方程简化为y'(t) = -\frac{2}{y(t)} - 1，其解析解为常数y(t)=-2。求解常数解的数值稳定性远高于线性增长的h(t)=t-2，能大幅降低长时间积分的误差累积。

3. 换用高阶自适应步长方法 #

尝试DOP853方法（solve_ivp支持的高阶显式Runge-Kutta方法），它的截断误差阶数高于RK45，相同精度要求下步长可更大，长时间积分的误差累积更慢：

sol = solve_ivp(fun, t_span, y0, method='DOP853', rtol=1e-8, atol=1e-10)

二、无解析解时的误差量化方法 #

无解析解可对比时，通过以下方式量化求解误差：

1. 高精度参考解对比 #

用极小的rtol和atol（如rtol=1e-12，atol=1e-15）求解一个高精度结果作为“参考解”，将常规求解结果与参考解逐点对比，计算绝对误差或相对误差，评估解的准确性。

2. 残差验证 #

将数值解代入原微分方程，计算残差：
$$\text{res}(t) = h'{\text{num}}(t) + \frac{2}{h{\text{num}}(t)-\lambda(t)}$$
其中h'_{\text{num}}(t)可通过数值差分（如中心差分）或求解器返回的插值结果求导得到。残差绝对值越小，数值解越接近真实解。

3. 多方法交叉验证 #

用两种不同类型的求解器（如显式RK45和隐式Radau）求解同一问题，对比结果差异。若差异在预设误差范围内（如小于1e-6），则解的可靠性较高；若差异过大，需调整精度参数或检查方程建模。

4. 守恒量检验 #

若能从原方程推导得到守恒量（不随时间变化的量），可通过计算数值解的守恒量波动评估误差。比如变量替换后的方程，积分可得守恒量：
$$\frac{1}{2}y(t)^2 + \int_0^t y(\tau)\lambda'(\tau)d\tau + 2t = C$$
其中C由初始条件确定，数值解的守恒量与理论值的偏差，可作为误差量化指标。

内容的提问来源于stack exchange，提问作者nicholas-t-nguyen