# Two-dimensional collisions

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. Particles 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.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
@x.setter
def x(self, value):
self.r = value
@property
def y(self):
return self.r
@y.setter
def y(self, value):
self.r = value
@property
def vx(self):
return self.v
@vx.setter
def vx(self, value):
self.v = value
@property
def vy(self):
return self.v
@vy.setter
def vy(self, value):
self.v = value

def overlaps(self, other):
"""Does the circle of this Particle overlap that of other?"""

def draw(self, ax):
"""Add this Particle's Circle patch to the Matplotlib Axes ax."""

return circle

"""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.vx = -self.vx
if self.x + self.radius > 1:
self.vx = -self.vx
if self.y - self.radius < 0:
self.vy = -self.vy
if self.y + self.radius > 1:
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.

"""

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.

"""

"""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:
except TypeError:
# r isn't iterable: turn it into a generator that returns the
# same value n times.
for i in range(n):

self.n = n
self.particles = []
# 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.
# 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.

"""

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])

"""Advance the animation by dt, returning the updated Circles list."""

for i, p in enumerate(self.particles):
self.circles[i].center = p.r
self.handle_collisions()
return self.circles

for i, p in enumerate(self.particles):
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."""

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
styles = {'edgecolor': 'C0', 'linewidth': 2, 'fill': None}
sim.do_animation(save=False)

Current rating: 4.5 #### Lazaro 3 years ago

Excellent!!!!

Currently unrated #### christian 3 years ago

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

Current rating: 5 #### jodenda 3 years ago

Thanks a lot ! Wonderful

Currently unrated #### TAllen 2 years, 11 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!

Currently unrated #### christian 2 years, 11 months ago

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

Currently unrated #### fhubbell 1 year, 10 months ago

Great example, just what I was looking for.

Currently unrated #### Felipe 1 year ago

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

Currently unrated #### christian 1 year 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}

Currently unrated #### phy 10 months, 2 weeks ago

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

Currently unrated #### christian 10 months, 1 week ago

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

%matplotlib notebook

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

Currently unrated #### orztrickster 8 months, 3 weeks ago

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