【问题标题】:Discrepancy in curve fit using solutions from solve_ivp and odeint使用 solve_ivp 和 odeint 解决方案的曲线拟合差异
【发布时间】:2021-11-29 11:07:49
【问题描述】:

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

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

使用odeintivp_solve 进行数值积分并在curve_fit 上使用它们产生了相当大的差异。与较新的 solve_ivp 相比,较旧的 odeint 更适合。 k 的最佳拟合值也非常不同。

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

这背后的原因是什么?

## 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))

【问题讨论】:

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


    【解决方案1】:

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

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

    要获得更兼容的计算,请显式设置容错,因为默认值不必相等。比如atol=1e-22, rtol=1e-9,那么只有相对容差有影响。

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

    【讨论】:

    • 我在using_ivp函数中使用了这一行的区间:sol = solve_ivp(eqn, [min(times), max(times)], y0=[y0], args=(k, ), dense_output=True)。我得到了与odeint 非常接近的结果。如果没有给出积分范围,确实solve_ivp 做了一些疯狂的推断!
    • 这里并不重要,因为它们在构造上是相同的,但是为了代码的一致性,最好使用本地定义的t 数组而不是全局times 数组。这也是一个封装/隔离问题。
    • 确实,在这种情况下使用t的min, max而不是全局的times是对的。
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 2021-01-24
    • 1970-01-01
    • 2017-07-19
    • 2023-04-09
    • 2014-12-05
    • 2021-11-20
    • 2021-12-24
    相关资源
    最近更新 更多