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.
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?