This small Python project is a physical simulation of two-dimensional physics. The animation is carried out using Matplotlib's `FuncAnimation`

method and is implemented by the class `Simulation`

. Each "particle" of the simulation is represented by an instance of the `Particle`

class and depicted as a circle with a fixed radius which undergoes elastic collisions with other particles.

A `Particle`

's position and velocity vectors are represented by NumPy arrays `r`

and `v`

, but their components can be retrieved and set by using the shortcuts `Particle.x`

, `Particle.y`

etc., implemented as class properties. `Particle`

s move within the domain $0 \le x < 1$, $0 \le y < 1$ and are reflected (elastically) off the edges (walls) of this domain. The elastic collision between two particles occurs when their centres are separated by the sum of their radii and the collision event changes their velocities to:

$$ \boldsymbol{u}_1 = \boldsymbol{v}_1 - \frac{2m_2}{M}\frac{(\boldsymbol{v}_1 - \boldsymbol{v}_2)\cdot(\boldsymbol{r}_1 - \boldsymbol{r}_2)}{|\boldsymbol{r}_1 - \boldsymbol{r}_2|^2}(\boldsymbol{r}_1 - \boldsymbol{r}_2) $$

and

$$ \boldsymbol{u}_2 = \boldsymbol{v}_2 - \frac{2m_1}{M}\frac{(\boldsymbol{v}_2 - \boldsymbol{v}_1)\cdot(\boldsymbol{r}_2 - \boldsymbol{r}_1)}{|\boldsymbol{r}_1 - \boldsymbol{r}_2|^2}(\boldsymbol{r}_2 - \boldsymbol{r}_1), $$

as dictated by the requirements of conservation of momentum and kinetic energy. The Wikipedia article on elastic collisions is a good place to learn more about this. Here, $M=m_1+m_2$ is the total mass of the colliding pair, and the particle's mass is considered to be just the square of its radius (*i.e.* its density is $\pi^{-1}$).

The most up-to-date code is on GitHub but as initially committed is given below.

```
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from matplotlib import animation
from itertools import combinations
class Particle:
"""A class representing a two-dimensional particle."""
def __init__(self, x, y, vx, vy, radius=0.01, styles=None):
"""Initialize the particle's position, velocity, and radius.
Any key-value pairs passed in the styles dictionary will be passed
as arguments to Matplotlib's Circle patch constructor.
"""
self.r = np.array((x, y))
self.v = np.array((vx, vy))
self.radius = radius
self.styles = styles
if not self.styles:
# Default circle styles
self.styles = {'edgecolor': 'b', 'fill': False}
# For convenience, map the components of the particle's position and
# velocity vector onto the attributes x, y, vx and vy.
@property
def x(self):
return self.r[0]
@x.setter
def x(self, value):
self.r[0] = value
@property
def y(self):
return self.r[1]
@y.setter
def y(self, value):
self.r[1] = value
@property
def vx(self):
return self.v[0]
@vx.setter
def vx(self, value):
self.v[0] = value
@property
def vy(self):
return self.v[1]
@vy.setter
def vy(self, value):
self.v[1] = value
def overlaps(self, other):
"""Does the circle of this Particle overlap that of other?"""
return np.hypot(*(self.r - other.r)) < self.radius + other.radius
def draw(self, ax):
"""Add this Particle's Circle patch to the Matplotlib Axes ax."""
circle = Circle(xy=self.r, radius=self.radius, **self.styles)
ax.add_patch(circle)
return circle
def advance(self, dt):
"""Advance the Particle's position forward in time by dt."""
self.r += self.v * dt
# Make the Particles bounce off the walls
if self.x - self.radius < 0:
self.x = self.radius
self.vx = -self.vx
if self.x + self.radius > 1:
self.x = 1-self.radius
self.vx = -self.vx
if self.y - self.radius < 0:
self.y = self.radius
self.vy = -self.vy
if self.y + self.radius > 1:
self.y = 1-self.radius
self.vy = -self.vy
class Simulation:
"""A class for a simple hard-circle molecular dynamics simulation.
The simulation is carried out on a square domain: 0 <= x < 1, 0 <= y < 1.
"""
def __init__(self, n, radius=0.01, styles=None):
"""Initialize the simulation with n Particles with radii radius.
radius can be a single value or a sequence with n values.
Any key-value pairs passed in the styles dictionary will be passed
as arguments to Matplotlib's Circle patch constructor when drawing
the Particles.
"""
self.init_particles(n, radius, styles)
def init_particles(self, n, radius, styles=None):
"""Initialize the n Particles of the simulation.
Positions and velocities are chosen randomly; radius can be a single
value or a sequence with n values.
"""
try:
iterator = iter(radius)
assert n == len(radius)
except TypeError:
# r isn't iterable: turn it into a generator that returns the
# same value n times.
def r_gen(n, radius):
for i in range(n):
yield radius
radius = r_gen(n, radius)
self.n = n
self.particles = []
for i, rad in enumerate(radius):
# Try to find a random initial position for this particle.
while True:
# Choose x, y so that the Particle is entirely inside the
# domain of the simulation.
x, y = rad + (1 - 2*rad) * np.random.random(2)
# Choose a random velocity (within some reasonable range of
# values) for the Particle.
vr = 0.1 * np.random.random() + 0.05
vphi = 2*np.pi * np.random.random()
vx, vy = vr * np.cos(vphi), vr * np.sin(vphi)
particle = Particle(x, y, vx, vy, rad, styles)
# Check that the Particle doesn't overlap one that's already
# been placed.
for p2 in self.particles:
if p2.overlaps(particle):
break
else:
self.particles.append(particle)
break
def handle_collisions(self):
"""Detect and handle any collisions between the Particles.
When two Particles collide, they do so elastically: their velocities
change such that both energy and momentum are conserved.
"""
def change_velocities(p1, p2):
"""
Particles p1 and p2 have collided elastically: update their
velocities.
"""
m1, m2 = p1.radius**2, p2.radius**2
M = m1 + m2
r1, r2 = p1.r, p2.r
d = np.linalg.norm(r1 - r2)**2
v1, v2 = p1.v, p2.v
u1 = v1 - 2*m2 / M * np.dot(v1-v2, r1-r2) / d * (r1 - r2)
u2 = v2 - 2*m1 / M * np.dot(v2-v1, r2-r1) / d * (r2 - r1)
p1.v = u1
p2.v = u2
# We're going to need a sequence of all of the pairs of particles when
# we are detecting collisions. combinations generates pairs of indexes
# into the self.particles list of Particles on the fly.
pairs = combinations(range(self.n), 2)
for i,j in pairs:
if self.particles[i].overlaps(self.particles[j]):
change_velocities(self.particles[i], self.particles[j])
def advance_animation(self, dt):
"""Advance the animation by dt, returning the updated Circles list."""
for i, p in enumerate(self.particles):
p.advance(dt)
self.circles[i].center = p.r
self.handle_collisions()
return self.circles
def advance(self, dt):
"""Advance the animation by dt."""
for i, p in enumerate(self.particles):
p.advance(dt)
self.handle_collisions()
def init(self):
"""Initialize the Matplotlib animation."""
self.circles = []
for particle in self.particles:
self.circles.append(particle.draw(self.ax))
return self.circles
def animate(self, i):
"""The function passed to Matplotlib's FuncAnimation routine."""
self.advance_animation(0.01)
return self.circles
def do_animation(self, save=False):
"""Set up and carry out the animation of the molecular dynamics.
To save the animation as a MP4 movie, set save=True.
"""
fig, self.ax = plt.subplots()
for s in ['top','bottom','left','right']:
self.ax.spines[s].set_linewidth(2)
self.ax.set_aspect('equal', 'box')
self.ax.set_xlim(0, 1)
self.ax.set_ylim(0, 1)
self.ax.xaxis.set_ticks([])
self.ax.yaxis.set_ticks([])
anim = animation.FuncAnimation(fig, self.animate, init_func=self.init,
frames=800, interval=2, blit=True)
if save:
Writer = animation.writers['ffmpeg']
writer = Writer(fps=100, bitrate=1800)
anim.save('collision.mp4', writer=writer)
else:
plt.show()
if __name__ == '__main__':
nparticles = 50
radii = np.random.random(nparticles)*0.03+0.02
styles = {'edgecolor': 'C0', 'linewidth': 2, 'fill': None}
sim = Simulation(nparticles, radii, styles)
sim.do_animation(save=False)
```

## Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

## Lazaro 4 years ago

Excellent!!!!

Link | Reply## christian 4 years ago

I'm glad you like it! You might also be interested in https://scipython.com/blog/the-maxwellboltzmann-distribution-in-two-dimensions/

Link | Reply## jodenda 4 years ago

Thanks a lot ! Wonderful

Link | Reply## TAllen 3 years, 10 months ago

Wow, I've been looking over your whole blog and I'm amazed at the plethora of topics that you've programmed about. It's amazing that you're able to learn about all of these things and to then program about them. I'm just starting out on my programming development and so I'm extremely grateful to have all of these resources of science&physics crossed with programming to help me learn more about pythons wide range of applicability. Thank you for these very helpful resources!

Link | Reply## christian 3 years, 10 months ago

That's very nice of you – thank you. I'm glad you find the site useful!

Link | Reply## fhubbell 2 years, 9 months ago

Great example, just what I was looking for.

Link | Reply## Felipe 2 years ago

How do I add different particle objects? With differentiated colors and sizes.

Link | Reply## christian 2 years ago

You pass the radius to the Particle class when you instantiate it, along with any custom styles if you want to change the defaults, e.g. styles = {'edgecolor': 'r', 'fill': True}

Link | Reply## phy 1 year, 9 months ago

It didn't work. when I start this, it appears the circles but not moving. why is that??

Link | Reply## christian 1 year, 9 months ago

If you are using Jupyter Notebook, you might need to add

Link | Reply%matplotlib notebook

to make animations work. To save the animation in mp4 format, you might also need to have FFmpeg installed.

## orztrickster 1 year, 8 months ago

Changing the bottom nparticles = 50 to nparticles = 100 seems to find that the circles overlap after a while

Link | Reply## New Comment