I’m working with discrete-time simulations of ODEs with time varying parameters. I have time series of various data (e.g. time series of state vectors generated by solve_ivp, time series of system matrices generated by my control algorithm, time series of system matrices in modal form, and so on).

My question: in what order should I place the indices? My intuition is that since numpy arrays are (by default) stored in row-major order, and I want per-item locality, each row should contain the “item” (i.e. a vector or matrix), and so the number of rows is the number of time points, and the number of columns is the dimension of my vector, e.g.:

x_k = np.array((5000, 4)) # a time series of 5000, 4-vectors
display(x_k[25]) # the 26th timepoint

Or for matrices I might use:

A_k = np.array((5000, 4, 4)) # a time series of 5000, 4x4-matrices

However, solve_ivp appears to do the opposite and returns a row-major array with the time series in columns (sol.y shape is (4, 5000)). Furthermore, transposing the result with .T just flips a flag to column-major so it is not really clear what the developers of solve_ivp and numpy intend me to do to write cache efficient code.

What are the conventions? Should I use the first index for the time index, as in my examples above, or last index as solve_ivp does?