爲odeint

Question

優化寫作公式我寫了下面的代碼簡單的化學網絡由odeint解決：

def chemnet(y,t): 
    assoc=0.1,0.001,1 
    oxy=0.001,0.1,0.001 
    f=zeros(6,float) 
    f[0]= -2*assoc[0]*y[0]**2-assoc[1]*y[0]*y[1]-assoc[2]*y[0]*y[2]-oxy[0]*y[0]*y[4]+oxy[1]*y[1]*y[4] 
    f[1]= +assoc[0]*y[0]**2-assoc[1]*y[0]*y[1]-oxy[1]*y[1]*y[4]+oxy[2]*y[2]*y[4] 
    f[2]= +assoc[1]*y[0]*y[1]-assoc[2]*y[0]*y[2]-oxy[2]*y[2]*y[4] 
    f[3]= +assoc[2]*y[0]*y[2] 
    f[4]= -oxy[0]*y[0]*y[4]-oxy[1]*y[1]*y[4]-oxy[2]*y[2]*y[4] 
    f[5]= +oxy[0]*y[0]*y[4]+oxy[1]*y[1]*y[4]+oxy[2]*y[2]*y[4] 
    return f 
time = linspace(0.0,10.0,1000) 
yinit = array([1,0,0,0,0.25,0]) 
y = odeint(chemnet,yinit,time) 

因爲我擴大網絡規模，方程的數量和長度每個方程都增加，因爲一些物種有許多可能的反應。我想在僞代碼的過程中自動創建這些方程，即：

#reaction 0+n->n+1 
for n in range(len(assoc)): 
    f[0].append(-assoc[n]*y[0]*y[n]) 
    f[n].append(-assoc[n]*y[0]*y[n]) 
    f[n+1].append(+assoc[n]*y[0]*y[n]) 

但我不太清楚如何做到這一點，或者如果它甚至有可能。我仍然是一個蟒蛇新手，但似乎這應該是比它更直接。

Answer 1

您可以自動構建等式，但在未明確構造等式的情況下自動執行評估會容易得多。

修改你的例子：

def chemnet(y,t): 
    assoc=0.1,0.001,1 
    oxy=0.001,0.1,0.001 
    f=zeros(6,float) 
    f[0]= -oxy[0]*y[0]*y[4]+oxy[1]*y[1]*y[4] 
    f[1]= -oxy[1]*y[1]*y[4]+oxy[2]*y[2]*y[4] 
    f[2]= -oxy[2]*y[2]*y[4] 
    f[3]= 0. 
    f[4]= -oxy[0]*y[0]*y[4]-oxy[1]*y[1]*y[4]-oxy[2]*y[2]*y[4] 
    f[5]= +oxy[0]*y[0]*y[4]+oxy[1]*y[1]*y[4]+oxy[2]*y[2]*y[4] 

    #reaction 0+n->n+1 
    for n,assoc_n in enumerate(assoc): 
     f[0] -= assoc_n*y[0]*y[n] 
     f[n] -= assoc_n*y[0]*y[n] 
     f[n+1] += assoc_n*y[0]*y[n] 
    return f 

time = linspace(0.0,10.0,1000) 
yinit = array([1,0,0,0,0.25,0]) 
y = odeint(chemnet,yinit,time) 

有了這樣的線性關係，還可以構建稀疏矩陣的交互和使用矩陣產品做你的更新。