使用numpy數組與scipy odeint

Question

我使用scipy來求解一個常微分方程組。爲簡單起見，把我的代碼是：

import scipy as sp 
import numpy as np 
from scipy.integrate import odeint 
from numpy import array 

def deriv(y,t): # return derivatives of the array y 
a = -2.0 
b = -0.1 
return array([ y[1], a*y[0]+b*y[1] ]) 
time = np.linspace(0.0,10.0,100) 
yinit = array([0.0005,0.2]) # initial values 
y = odeint(deriv,yinit,time) 

但現在我想解決這個系統的常數「a」的幾個值。因此，例如，我不想只有一個= -2.0，我想有：

a = array([[-2.0,-1.5,-1.,-0.5]]) 

並解決系統的每個值的一個。有沒有辦法做到這一點，而不必循環數組的每個元素？我可以一次完成嗎？

Answer 1

您可以通過爲每個a的值編寫兩個方程來重新表達方程組。做到這一點的方法之一是

from scipy.integrate import odeint 
from numpy import array,linspace,tile,empty_like 
a = array([-2.0,-1.5,-1.,-0.5]) 
b = array([-.1,-.1,-.1,-.1]) 

yinit = tile(array([0.0005,0.2]),a.size) 

def deriv(y,_): 
    dy = empty_like(y) 
    dy[0::2]=y[1::2] 
    dy[1::2]=y[::2]*a+b*y[1::2] 
    return dy 
time = linspace(0.0,10.0,100) 
y = odeint(deriv,yinit,time) 

你會y[:,0]找到y的解決方案a=-2，在y[:,2]解決方案爲a=-1.5，等等等等與y[:,-1]是的y爲a=-.5的衍生物。

無論如何，如果你想矢量化這個以獲得速度，你可能會對將你的python代碼轉換成c的庫感興趣，然後編譯它。我個人使用pydelay，因爲我需要模擬時間延遲，並會推薦它。你甚至不需要處理C代碼，翻譯是完全自動化的。