使用solve_ivp与odeint的解进行曲线拟合出现明显差异的原因分析

问题原因 #

两个核心原因共同导致了结果差异：

• 求解器适配性问题：你拟合的方程当k值很大时属于刚性常微分方程，odeint默认使用LSODA求解器，支持自动在刚性/非刚性求解算法之间切换，对刚性问题适配性好；而solve_ivp默认使用RK45求解器，属于显式非刚性求解器，处理大k的刚性方程时步长受稳定性限制，默认精度设置下积分误差极大，最终拟合出来的参数完全偏离真实值。
• 精度设置差异：odeint默认的相对误差阈值rtol=1.4e-8、绝对误差阈值atol=1.4e-8，而solve_ivp默认rtol=1e-3、atol=1e-6，精度设置差了几个数量级，进一步放大了结果偏差。

修正方案 #

修改using_ivp函数，指定求解器为LSODA，同时对齐精度阈值即可，修正后的代码如下：

def using_ivp(t, k):
    eqn = lambda t, x, k: -k * x
    y0 = vals[0]
    # 显式指定求解器、精度，以及积分起止区间
    sol = solve_ivp(eqn, t_span=(t.min(), t.max()), y0=[y0], args=(k, ),
                    method='LSODA', rtol=1.4e-8, atol=1.4e-8,
                    dense_output=True)
    return sol.sol(t)[0]

修正后两者拟合得到的k值会基本一致，拟合曲线也会重合。

内容的提问来源于stack exchange，提问作者cmbfast