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 three-dimensional 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 so-called 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 three-dimensional space, this matrix would result in:

Position matrix

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 multi-dimensional arrays, which can be understood as matrices of matrices:

In previous figure, each axis is devoted to a particular task:

  • The axis-0 of the ndarray is devoted to store the different position matrices at each one of the time steps.
  • The axis-1 refers to the i-th particle.
  • And finally, the axis-2 is used to store the position coordinates.

But how do you create these multi-dimensional 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 low-level method (ndarray(…)) for instantiating an array.


In this first post of the series, the basic elements of particle systems and an effective way of storing simulation data via multi-dimensional 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