python - 使用solve_ivp和odeint的解决方案的曲线拟合差异

标签 python scipy curve-fitting differential-equations numerical-integration

这是一个基本示例方程,我试图将其拟合到示例数据。

\frac{dx}{dt} = -k x

目标是假设数据遵循上述方程,找到最适合我的数据的 k。一个明显的方法是对方程进行数值积分,然后使用曲线拟合方法最小化最小二乘法并得到 k

使用 odeint 进行数值积分和 ivp_solve并在 curve_fit 上使用它们产生了相当大的差异。与新版 solve_ivp 相比,旧版 odeint 的拟合效果更好。 k 的最佳拟合值也有很大不同。

Best fit k using ivp = [2.74421966]
Best fit k using odeint = [161.82220545]

这背后的原因是什么?

Discrepancy

## SciPy, NumPy etc imports here

N_SAMPLES = 20
rnd = np.random.default_rng(12345)
times = np.sort(rnd.choice(rnd.uniform(0, 1, 100), N_SAMPLES))
vals = np.logspace(10, 0.1, N_SAMPLES) + rnd.uniform(0, 1, N_SAMPLES)


def using_ivp(t, k):
    eqn = lambda t, x, k: -k * x
    y0 = vals[0]
    sol = solve_ivp(eqn, t, y0=[y0], args=(k, ),
                    dense_output=True)
    return sol.sol(t)[0]

def using_odeint(t, k):
    eqn = lambda x, t: -k * x
    y0 = vals[0]
    sol = odeint(eqn, y0, t)
    return sol[:,0]
    
tfit = np.linspace(min(times), max(times), 100)


#Fitting using ivp
k_fit1, kcov1 = curve_fit(using_ivp, times, vals, p0=1.3)
fit1 = using_ivp(tfit, k_fit1)


#Fitting using odeint
k_fit2, kcov2 = curve_fit(using_odeint, times, vals, p0=1.3)
fit2 = using_odeint(tfit, k_fit2)

plt.figure(figsize=(5, 5))
plt.plot(times, vals, 'ro', label='data')
plt.plot(tfit, fit1, 'r-', label='using ivp')
plt.plot(tfit, fit2, 'b-', label='using odeint')
plt.legend(loc='best');

print('Best fit k using ivp = {}\n'\
      'Best fit k using odeint = {}\n'.\
      format(k_fit1, k_fit2))

最佳答案

再次检查solve_ivp的输入参数是什么。积分间隔由 t_span 参数中的前两个数字给出,因此在您的应用程序中 sol.sol(t) 中的大多数值都是通过狂野外推获得的。

通过将间隔指定为 [min(t),max(t)] 来更正该问题。

为了获得更兼容的计算,请显式设置误差容限,因为默认值不必相等。例如 atol=1e-22, rtol=1e-9 这样只有相对容差才会产生影响。

很好奇你如何使用args机制。它最近才被引入到 solve_ivp 中,以便与 odeint 更加兼容。在这两种情况下我都不会使用它,因为参数的定义及其使用包含在一个 3 行 block 中。它有其用途,其中正确的封装和与其他代码的隔离是真正关心的问题。

关于python - 使用solve_ivp和odeint的解决方案的曲线拟合差异,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/69513659/

相关文章:

python - 巴特沃斯滤波器 - 输出 x (-1)?

python - 替换numpy数组值python

python - 使用 curve_fit 获取 r 平方值

matlab - 在Matlab中移除图像偏移(二维基线)

Python:子进程可以暂停/恢复父进程吗

python - Python 中的 __set__ 和 __setattr__ 有什么区别,应该在什么时候使用?

iphone - 使用 Python 解析 .strings 文件

python - 将 GET 参数转换为 Django REST Framework 中 Request 对象上的 POST 数据

python - 如何使用 "scipy.optimize.curve_fit"平滑地拟合我的数据点?

python - 如何使用数据集来拟合 3D 表面?