Trouble getting odeINT to accept equations generated by a loop

Question

I want to use Python's odeINT to integrate multiple sets of equations that are generated by a loop. The equations are all coupled and so must be integrated simultaneously, via one call to odeINT. The problem is that the initial conditions ("y0") must be a list of lists or a matrix (rather than simply a list). odeINT gives this error: "Initial condition y0 must be one-dimensional". I am wondering how to work around this. Here is an example of the code; thanks so much for any ideas.

class network:
    def __init__(self):
        self.i_range = 3
        ## INITIAL CONDITIONS WILL BE A LIST OF LISTS. 
        ## THIS IS THE SOURCE OF odeINT's ERROR.
        self.init = [[] for i in range(self.i_range)]  
        for i in range(0,self.i_range):
            self.init[i].append(-50.+0.1*(random.random()))
            self.init[i].append(1.+1.*(random.random()))

        self.Tfinal = 10  # final time   
        self.dt = 1.    # time step

    def eqns(self, x, t):
        a, b = x
        dadt = zeros_like(a)
        dbdt = zeros_like(b)
        for i in range (0,i_range):
            dadt[i] = np.cos(b[i])
            dbdt[i] = np.sin(a[i])
        return dadt, dbdt

    def run(self):
        self.times = sp.arange(0,self.Tfinal,self.dt)
        self.sim = odeint(self.eqns,self.init,self.times)

Answer 1

Use

a,b = reshape(x,(2,-1))

to split the flat vector and

return reshape((dadt,dbdt), -1)

to recombine them into a flat array. See https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html

---

More concretely, and playing with some of the list functionality of python, in the initialization do

a = [ -50. + 0.1*random.random() for _ in range(self.i_range) ]
b = [   1. + 1.0*random.random() for _ in range(self.i_range) ]
self.init = np.reshape((a,b), -1)

and the ode function becomes

def eqns(self, x, t):
    a, b = np.reshape(x,(2,-1))
    dadt = [ np.cos(bb) for aa,bb in zip(a,b) ]
    dbdt = [ np.sin(aa) for aa,bb in zip(a,b) ]
    return np.reshape((dadt,dbdt), -1)

This should be sufficient to get some result. You might want to transform self.sim into a list of pairs of lists to regain the structure of the problem, however there you have to work, as there are now 2 "large" dimensions to provide.

np.reshape(self.sim, (self.times.size,2,-1))

could work as a first step..