首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >使用solve_ivp和odeint的解决方案进行曲线拟合时存在差异

使用solve_ivp和odeint的解决方案进行曲线拟合时存在差异
EN

Stack Overflow用户
提问于 2021-10-10 08:44:42
回答 1查看 101关注 0票数 2

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

我们的目标是找到最适合我的数据的k,假设数据遵循上面的等式。一种显而易见的方法是对方程进行数值积分,然后使用曲线拟合方法来最小化最小二乘,从而获得k

使用odeintivp_solve进行数值集成,并在curve_fit上使用它们会产生相当大的差异。与较新的solve_ivp相比,较旧的odeint具有更好的适应性。k的最佳拟合值也有很大不同。

代码语言:javascript
运行
复制
Best fit k using ivp = [2.74421966]
Best fit k using odeint = [161.82220545]

这背后的原因可能是什么?

代码语言:javascript
运行
复制
## 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))
EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2021-10-10 09:00:18

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

通过将间隔设置为[min(t),max(t)]来纠正此错误。

要获得更兼容的计算,请显式设置误差容差,因为默认值不必相等。例如atol=1e-22, rtol=1e-9,因此只有相对公差才有效。

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

票数 3
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/69513659

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档