# How best to use scipy.integrate's ODE solvers when state is not naturally stored as a vector?

I have a large system of ODEs. For various reasons, it is natural to store the values of the dependent variables in a multidimensional array. For example, these values might represent the solution of a PDE at a point in a 2D or 3D grid, after semi-discretization.

The function scipy.integrate.odeint takes a vector (i.e., a one-dimensional array) of initial values, and returns back vectors of the same size. If also requires a right-hand-side function $f$ that takes and returns vectors. So one must flatten the state vector in order to use scipy.integrate.odeint. Meanwhile, inside $f$ it might be convenient to transform the state back to a multidimensional array, compute the derivative, and flatten the derivative array to a vector.

Is this flattening(for passing to f)->unflattening(for working within f)->flattening(for returning something from f) cycle the correct way to do things? Is there some other, better way to set things up (that might also happen to be efficient by avoiding the flattening/unflattening cycle?

The scipy.integrate.odeint function does not take a vector as input. It takes a function f, an initial state y0 and a vector of points in time t (and possibly a bunch of other optional arguments to control the integration strategy). So you will need to provide it with a function that calculates the right-hand-side of your ODE system, which will boil down to flattening (i.e. looping over all your 'bodies' and building a vector ft which is the RHS evaluated at time t. There is no way around this.
• y0 is a vector though -- it's only scalar when you have a single ODE (i.e. not a system of ODEs)? – user89 Oct 1 '14 at 4:14