2014-01-21 38 views
0

我的问题是关于当前的scipy颂歌求解器。从scipy doc page,它们的用法是:矢量化SciPy颂歌求解器

# A problem to integrate and the corresponding jacobian: 

from scipy.integrate import ode 
y0, t0 = [1.0j, 2.0], 0 
def f(t, y, arg1): 
    return [1j*arg1*y[0] + y[1], -arg1*y[1]**2] 
def jac(t, y, arg1): 
    return [[1j*arg1, 1], [0, -arg1*2*y[1]]] 

# The integration: 
r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True) 
r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0) 
t1 = 10 
dt = 1 
while r.successful() and r.t < t1: 
    r.integrate(r.t+dt) 
    print("%g %g" % (r.t, r.y)) 

我的问题是:它使用了很多蟒蛇循环(while循环)基本上减慢程序下来。我可以尝试编写C代码并使用ctypes使其更快,但我无法在scipy中访问像dopri5这样的漂亮算法(除非我自己实现它)。

是否有任何矢量化的编码方式,使其更快?

谢谢!

回答

2

'矢量化'意味着一次完成并行计算。是的,详细的实现将涉及迭代,但它在C和顺序并不重要你,程序员Python

但是这样的解决方案本质上是一个串行操作。您必须在时间t解决问题,然后在时间t+dt解决它。您无法通过时间向量化解决方案。你所能做的最好的是选择一个ode求解器,它可以为时间步骤(dt)提供智能选择,在可能的情况下采取大步骤,在需要捕获快速变化时采取小步骤。

一个很好的颂歌求解器可以让您矢量化空间维度 - 即平行求解10个颂。您可能还可以矢量化计算雅可比矩阵,一次返回y+dyy-dy。基本上你想尽可能快地计算fjac

+0

非常感谢!这对我来说是完全意义上的。 –

3

这有点晚,但我认为你需要的是scipy.integrate.odeint。它需要一系列时间点并计算每个时间点的解决方案。这使得循环在Python之外进行计算(在FORTRAN中)。

+0

永不嫌晚!这是我最终使用的,而其他人在未来看这篇文章可能会发现它有帮助。非常感谢! –