Оптимизация написания уравнений для odeint

Я написал следующий код для простой химической сети, которую должен решить 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, но кажется, что это должно быть проще, чем есть.


python scipy odeint
person SteelAngel    schedule 18.05.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Вы можете автоматизировать построение уравнения, но гораздо проще автоматизировать оценку без явного построения уравнения.

Изменение вашего примера:

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)

С такими линейными отношениями вы также можете создавать разреженные матрицы взаимодействия и использовать матричные произведения для своих обновлений.

person Neapolitan    schedule 18.05.2016