I am attempting to write a Python code to simulate many particles in a confined box. These particles behave in such a way that they move in the box in straight lines with a slight angular noise (small changes in the direction of the particle path). They should interact by acknowledging the other particle and 'shuffle/squeeze' past each other and continue on their intended path, much like humans on a busy street. Eventually, the particles should cluster together when the density of particles (or packing fraction) reaches a certain value.

However, I have a feeling there are parts of the code that are inefficient or which could be either sped up or written more conveniently.

If anyone has any improvements for the code speed or ideas which may help with the interactions and/or angular noise that would be much appreciated. I will also leave an example of an animation which is my aim: https://warwick.ac.uk/fac/sci/physics/staff/research/cwhitfield/abpsimulations

The above link shows the animation I am looking for, although I don't need the sliders, just the box, and moving particles. The whole code is shown below:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def set_initial_coordinates():
    x_co = [np.random.uniform(0, 2) for i in range(n_particles)]
    y_co = [np.random.uniform(0, 2) for i in range(n_particles)]
    return x_co, y_co

def set_initial_velocities():
    x_vel = np.array([np.random.uniform(-1, 1) for i in range(n_particles)])
    y_vel = np.array([np.random.uniform(-1, 1) for i in range(n_particles)])
    return x_vel, y_vel

def init():
    ax.set_xlim(-0.05, 2.05)
    ax.set_ylim(-0.07, 2.07)
    return ln,

def update(dt):
    xdata = initialx + vx * dt
    ydata = initialy + vy * dt
    fx = np.abs((xdata + 2) % 4 - 2)
    fy = np.abs((ydata + 2) % 4 - 2)

    for i in range(n_particles):
        for j in range(n_particles):
            if i == j:

            dx = fx[j] - fx[i]  # distance in x direction
            dy = fy[j] - fy[i]  # distance in y direction
            dr = np.sqrt((dx ** 2) + (dy ** 2))  # distance between x

            if dr <= r:
                force = k * ((2 * r) - dr)  # size of the force if distance is less than or equal to radius

                # Imagine a unit vector going from i to j
                x_comp = dx / dr  # x component of force
                y_comp = dy / dr  # y component of force

                fx[i] += -x_comp * force  # x force
                fy[i] += -y_comp * force  # y force
    ln.set_data(fx, fy)
    return ln,

# theta = np.random.uniform(0, 2) for i in range(n_particles)
n_particles = 10
initialx, initialy = set_initial_coordinates()
vx, vy = set_initial_velocities()
fig, ax = plt.subplots()
x_co, y_co = [], []
ln, = plt.plot([], [], 'bo', markersize=15)  # radius 0.05
plt.xlim(0, 2)
plt.ylim(0, 2)

k = 1
r = 0.1
t = np.linspace(0, 10, 1000)

ani = FuncAnimation(fig, update, t, init_func=init, blit=True, repeat=False)
  • 2
    \$\begingroup\$ @SᴀᴍOnᴇᴌᴀ Hi, the code supplied does work for the functions I have mentioned above and I would like to improve the lines above. In the other post, I am asking questions about lines that I have attempted but not completed. I would like to both improve this code and add new features but my question here is only asking for review/improvement, hopefully, this can keep it on topic? \$\endgroup\$
    – CMTT
    Dec 29, 2020 at 18:29
  • 1
    \$\begingroup\$ Okay - thanks for explaining that. I have retracted my close vote. I hope you receive valuable feedback. \$\endgroup\$ Dec 29, 2020 at 18:50
  • 1
    \$\begingroup\$ @SᴀᴍOnᴇᴌᴀ No worries, apologies for the confusion \$\endgroup\$
    – CMTT
    Dec 29, 2020 at 19:38
  • \$\begingroup\$ Has it to be matplotlib? The OO structure there has a huge penalty in time complexity, these lags are noticeable in animations. I did once a similar simulation in an answer to a question using pygame with a much smoother feeling: stackoverflow.com/questions/29374247/… \$\endgroup\$ Jan 16, 2021 at 20:20

2 Answers 2


Coordinates as a matrix


x_co = [np.random.uniform(0, 2) for i in range(n_particles)]
y_co = [np.random.uniform(0, 2) for i in range(n_particles)]
return x_co, y_co

is a little baffling. You call into Numpy to generate uniformly distributed numbers (good) but then put it in a list comprehension (not good) and separate x and y into separate lists (not good).

Consider instead

return np.random.uniform(low=0, high=2, size=(2, n_particles))

and similarly for set_initial_velocities:

return np.random.uniform(low=-1, high=1, size=(2, n_particles))

All of the rest of your code, such as the calculations for xdata/ydata should similarly vectorize over the new length-2 dimension eliminating distinct code for x and y.

  • \$\begingroup\$ Hi @Reinderien, thanks for your advice. Quick question about vectorising; should I keep my x_co and y_co lines below the 'def...' line and replace my current return line with your return line like this "def set_initial_coordinates(): x_co = [np.random.uniform(0, 2) for i in range(n_particles)] y_co = [np.random.uniform(0, 2) for i in range(n_particles)] return np.random.uniform(low=0, high=2, size=(2, n_particles))" or do this "def set_initial_coordinates(): return np.random.uniform(low=0, high=2, size=(2, n_particles))" or are both of those not what you meant? \$\endgroup\$
    – CMTT
    Jan 1, 2021 at 10:49
  • \$\begingroup\$ The second one is what I meant. \$\endgroup\$
    – Reinderien
    Jan 1, 2021 at 15:25
  • \$\begingroup\$ Awesome, thanks! \$\endgroup\$
    – CMTT
    Jan 1, 2021 at 20:29
  • \$\begingroup\$ @LutzLehmann You're entirely right. The premise of that entire section was wrong so I just deleted it. \$\endgroup\$
    – Reinderien
    Jan 16, 2021 at 17:33

On the implementation of the dynamic

The (first) value that animate passes to the update function as argument is the frame number, if you use that as dt you can expect some strangely accelerating particles. As your dynamic does not depend on time, the frame number should not enter the computation. The time step is a global constant related to the frame rate.

On second glance,

xdata = initialx + vx * dt

is completely wrong. That is to say, yes, you can use it, but the resulting dynamic is rather unnatural. What you intend to do is to implement the (symplectic) Euler method for

dx/dt = vx,  dy/dt = vy
dvx/dt = force_x(x), dvy/dt = force_y(x)

so you should increase the particle positions incrementally as you do with the velocities.

As dynamic on the torus, dx and dy should also be reduced by the torus period 4 using the same centering formula, dx = (dx+2)%4-2 etc. Else you miss interactions that cross the wrap-around lines.

Just as you have the time step in the position update, you should also have it in the velocity update

fx[i] += -x_comp * force * dt  # x force

Otherwise you get a substantially different dynamic if you change dt.

Overall, I got following code implementing the corrected dynamic and some numpy tricks

from matplotlib.animation import FuncAnimation

def dynamic_gen(n_particles, r_particles, elasticity, dt):
    x,y = np.random.uniform(low=0, high=2, size=(2, n_particles))
    vx,vy = np.random.uniform(low=-1, high=1, size=(2, n_particles))
    yield x,y
    r = r_particles
    k = elasticity
    while True:
        x += vx*dt; x = (x+2)%4-2
        y += vy*dt; y = (y+2)%4-2
        # dx[i,j] = x[i]-x[j]
        dx = x[:,None]-x[None,:]; dx = (dx+2)%4-2
        dy = y[:,None]-y[None,:]; dy = (dy+2)%4-2

        dr = np.hypot(dx,dy)
        dr = np.ma.masked_where( ~((0<dr)&(dr<r)), dr )
        force = k*(2*r-dr)
        vx += np.ma.sum(force*dx/dr,axis=1)*dt
        vy += np.ma.sum(force*dy/dr,axis=1)*dt
        yield x,y
def update(t,dynamic,scatter):
    x,y = next(dynamic)
    return scatter,

n_particles = 400
k = 20
r = 0.2
t = np.linspace(0, 4, 1000)
dynamic = dynamic_gen(n_particles,r,k,t[1]-t[0])
x,y = next(dynamic)

fig, ax = plt.subplots(1,1,figsize=(8,8))
scatter, = ax.plot(x,y,'ob',ms=5)
ax.set_xlim(-2,2); ax.set_ylim(-2,2)

ani = FuncAnimation(fig, update, frames=t, fargs=(dynamic,scatter), blit=True, repeat=False)

Guiding ideas:

  • separate computation and plotting code by providing a stepper for the dynamic progress in the form of a generator function, so that in the update function the only steps are "get new model data" and "modify plot scene data"

  • In the physics simulation, avoid explicit python loops, shift everything to implicit numpy loops behind the scenes. This is at the cost of creating some large objects in memory, but should be significantly faster.

  • Use the broadcasting trick dx = x[:,None]-x[None,:], so that `dx[i,j]=x[i]-x[j]

  • Use a masked arrays from numpy.ma to select the pairs of particles with non-zero forces and restrict the computation of the forces to them.

  • Use the appropriate sum function, sum over axis=1, that is, over the j index.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.