How to simulate particle systems with NumPy and SciPy  Part I
Modeling particle systems usually involves the usage of some kind of list or array objects in order to store the different vector quantities such us positions, velocities or accelerations. Most of these simulations consider a threedimensional space and the physical system to evolve over time. Therefore, how can we efficiently track all these vector quantities? In the first post of this series, a closer look at the general particle simulation elements is introduced together with how to effectively store simulation data.
Fundamentals of particle simulations
No matter the purpose of your particle simulation, there are some common elements between them, as the figure below depicts:
Each one of the previous elements except the animation is important from the point of view of the simulation computations:

The number of particles is directly related to the workload of the simulation. These particles are usually assumed to be points. Particles live within a physical space, usually \(\Re^2\) or \(\Re^3\). The computation workload is thus proportional to the number of particles, \(N_p\) and the number of dimensions, named \(N_d\).

Any particles simulator needs to solve the socalled equations of motion of the system. To do so, a span of time called here “time of integration” is declared. A numerical integrator is required to solve the equations of motion. Small equal frames of time, named integration steps, are used to feed the integrator so it can compute the problem. The computation workload increases with the time of integration and the number of integration steps. A balance must be found between the time it takes for the computer to solve the problem and the accuracy of the results.

Finally, although the animation is not mandatory, it provides a better understanding on how the system evolves. Complex animations might be difficult to produce, requiring computer graphics knowledge.
Effectively storing simulation data
A simulation will consider a system composed of \(N_p\) particles living in a vector space of \(\Re^{N_d}\) dimensions. At a given instant of time \(t\), the position of the particles can be described using a matrix who’s shape is \((N_p, N_d)\). For the particular case of a threedimensional space, this matrix would result in:
However, previous option is only useful for a given instant of time. What about
storing every position for each particle over a set of time steps? In that case,
the best way to store our simulation data is using NumPy
ndarray
objects, which are the final objects Numpy manages in the background, even when
declaring np.array
ones. The np.ndarray
class models multidimensional
arrays, which can be understood as matrices of matrices:
In previous figure, each axis is devoted to a particular task:
 The
axis0
of thendarray
is devoted to store the different position matrices at each one of the time steps.  The
axis1
refers to the ith particle.  And finally, the
axis2
is used to store the position coordinates.
But how do you create these multidimensional arrays? Well, it is discouraged to
instantiate ndarray
objects directly, as from the official NumPy
documentation:
Arrays should be constructed using array, zeros or empty (refer to the See Also section below). The parameters given here refer to a lowlevel method (ndarray(…)) for instantiating an array.
Conclusion
In this first post of the series, the basic elements of particle systems and an effective way of storing simulation data via multidimensional arrays was introduced. In the next post of the series, we will see how to declare these elements and index them for retrieving all necessary data.
Link to the second post of the series: Effectively simulating particle systems  Part II