我使用的是scipy.integrate.odeint,当切换到bdf方法时遇到了麻烦(因为系统看起来很僵硬)。不幸的是,bdf需要计算雅可比行列式。我的系统有成千上万个方程,因此这会杀死Python(分段错误!)

我有充分的理由相信亚当斯在这个系统上会做的很好。所以我想强迫scipy使用它。我认为odeint无法实现,因此我正在尝试使用ode。我已经编写了自己的函数,该函数应该与odeint相同。看起来我的代码不正确:

def my_odeint(dfunc, V0, times, args=()):
    r = integrate.ode(dfunc)
    r.set_integrator('vode', method='adams')
    r.set_initial_value(V0,times[0]).set_f_params(*args)
    V=[V0]
    for time in times[1:]:
        V.append(r.integrate(time))
    V = scipy.array(V)
    return V


当我测试它时,我得到一个错误:

def dx(X, t):
    x=X[0]
    y=X[1]
    return (-x, x+y)
X0 = [1,1]
times = scipy.linspace(0,2,10)
my_odeint(dx, X0, times)
> capi_return is NULL
> Call-back cb_f_in_dvode__user__routines failed.


然后

TypeError                                 Traceback (most recent call last)
<ipython-input-6-cc286a9956a2> in <module>()
----> 1 my_odeint(dx, X0, times)

in my_odeint(dfunc, V0, times, args)
   1352     V=[V0]
   1353     for time in times[1:]:
-> 1354         V.append(r.integrate(time))

in integrate(self, t, step, relax)
    406             self._y, self.t = mth(self.f, self.jac or (lambda: None),
    407                                 self._y, self.t, t,
--> 408                                 self.f_params, self.jac_params)

 in run(self, f, jac, y0, t0, t1, f_params, jac_params)
    863         args = ((f, jac, y0, t0, t1) + tuple(self.call_args) +
    864                 (f_params, jac_params))
--> 865         y1, t, istate = self.runner(*args)

<ipython-input-4-1770d3974ec9> in dx(X, t)
      1 def dx(X, t):
----> 2     x=X[0]
      3     y=X[1]

> TypeError: 'float' object has no attribute '__getitem__'


所以我的问题是-如何修复我的代码,以使我的函数接受与odeint相同的输入,但使用adams并且从不切换到bdf

最佳答案

scipy.ode要求右侧为f(t, y, *f_args)形式,odeint要求其右侧为f(y, t, *f_args)形式。

所以这将工作:

def my_odeint(dfunc, V0, times, args=()):
    r = integrate.ode(lambda t, X: scipy.array(dfunc(X, t, *args)))
    r.set_integrator('vode', method='adams')
    r.set_initial_value(V0,times[0])
    V=[V0]
    for time in times[1:]:
        V.append(r.integrate(time))
    V = scipy.array(V)
    return V


请注意,根据文档,vode实现的是Adams-Moulton,而不是Adams-Bashforth。也就是说,您仍然有一个隐式方法。如果集成继续失败,请尝试在integrate.ode中也实现的Dormand / Prince方法。

关于python - 如何使用Adams强制scipy进行整合?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/40317096/

10-14 17:44
查看更多