优化 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])

但我不太确定该怎么做，或者是否可能。我仍然是一个 python 新手，但看起来这应该比现在更直接。

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)

有了这样的线性关系，您还可以构建稀疏交互矩阵并使用矩阵乘积来进行更新。

优化 odeint 的书写方程

Optimizing writing equations for odeint

python

scipy

odeint