如何利用Adams来强迫枕叶进行集成？

Question

我正在使用scipy.integrate.odeint，当它切换到bdf方法时遇到了麻烦(因为系统看起来很僵硬)。不幸的是，bdf需要计算雅可比。我的系统有成千上万的方程，所以这会杀死Python (分割错误！！)

我有充分的理由相信亚当斯在这个系统上会做得很好。所以我想强迫枕木使用它。我认为这在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。

Answer 1

scipy.ode在表单f(t, y, *f_args)中需要一个右手边，奥丁在表单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实现的是亚当斯-穆尔顿，而不是亚当斯-巴什福思。也就是说，您仍然有一个隐式方法。如果集成继续失败，请尝试在integrate.ode中实现Dormand/Prince方法。