I wrote the following code for a simple chemical network to be solved by 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)
As I'm expanding the size of the network, the number of equations and the length of each equation increases, as some species have many possible reactions. I would like to automate the process of creating these equations, i.e. in pseudocode:
#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])
But I'm not quite sure how to do this, or if it's even possible. I'm still a bit of a python newbie, but it seems that this should be more straightforward than it is.
You could automate the construction of the equation but it is much easier to automate the evaluation without explicitly constructing the equation.
Modifying your example:
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)
With linear relationships like this, you can also construct sparse interaction matrices and use matrix products to do your updates.