2014-01-21 40 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

永不嫌晚!這是我最終使用的,而其他人在未來看這篇文章可能會發現它有幫助。非常感謝! –