In this post, we will simulate a two-dimensional particle system in Python. This sort of simulation can be used to describe the movement of planets, gas particles or people.
This post is structured as follows:
Along the way, we apply the system we have built to some sample problems, such as earth-sun simulation and similar things.
Modelling a particle as a point in a two-dimensional space is part of basic mechanics. What one needs is:
The vector notation describes two-dimensional vectors. The equation of motion for a particle is:
$$ \vec{f} = m \, \vec{a} $$The quantity $\vec{a}$ is the acceleration of the particle. It can be integrated two times to yield the variation in velocity of the particle and in position. This is described by the following equations:
$$ \vec{v}(t+dt) = \vec{v}(t) + \frac{1}{m} \vec{f} dt $$$$ \vec{x}(t+dt) = \vec{x}(t) + \vec{v}(t+dt) dt $$The assumptions behind these equations are essentially that the timestep $dt$ is small.
In the next section, we implement this particle system.
As can be gathered from the previous equations, we need to store a couple of quantities for a given particle. We therefore decide to construct a class keeping all the properties of the particle together in one place.
#necessary imports
from pylab import *
%matplotlib inline
class Particle():
def __init__(self, mass, initial_position, initial_velocity):
self.mass = mass
self.position = initial_position
self.velocity = initial_velocity
def move(self, dt, external_force):
self.velocity += 1/self.mass * external_force * dt
self.position += self.velocity * dt
As one can see, this is pretty simple. What can be quite complicated is the computation of the external force that moves the particle. External forces can take several forms:
I would like to simulate the action of planets, so let's take a look at gravity for instance. If we lookup wikipedia, we can find that the gravitational pull felt by mass $i$ due to mass $j$ can be written as:
$$ \vec{F}_{ij} = \frac{G m_i m_j (\vec{x}_j - \vec{x}_i)}{\left\| \vec{x}_j - \vec{x}_i\right\|^3}, $$This can be implemented in Python with the following function:
def gravitation_force(particle_i, particle_j, gravitation_constant):
return gravitation_constant * particle_i.mass * particle_j.mass * (particle_j.position - particle_i.position) / norm((particle_j.position - particle_i.position))**3
We can now already implement the simulation of a sun-earth system:
sun = Particle(1.989e30, # mass of sun
zeros((2)), # sun is at origin of grid
zeros((2))) # intially not moving
earth = Particle(5.972e24, # mass of earth
array([149600000e3, 0]), #sun-earth distance
array([0, 29.8e3])) # speed of earth
def move_earth(dt):
pull = gravitation_force(earth, sun, 6.67e-11)
earth.move(dt, pull)
Finally before starting we still need a plotting routine to see what happens to our planets.
def plot_sun_earth():
plot(sun.position[0], sun.position[1], 'bo')
plot(earth.position[0], earth.position[1], 'ro')
plot_sun_earth()
We can also write a little tool to plot the pull exerted by the sun on the earth graphically:
def plot_pull(factor=1):
pull = gravitation_force(earth, sun, 6.67e-11)
pull /= norm(pull) * factor
plot(vstack((pull, zeros((2))))[:, 0], vstack((pull, zeros((2))))[:, 1], '->')
plot_sun_earth()
plot_pull(1000000000000000000)
move_earth(1)
plot_sun_earth()
As we can see, the earth has moved up by approximately 30000 meters. This seems correct, due to the fact that the earth is orbiting the sun at approximately 30 km per second.
We can now plot several positions of the earth as a function of time.
for i in range(10):
move_earth(120)
plot_sun_earth()
We can see that the earth radially escapes the sun. What happens if we iterate our simple model over a day, instead of 20 minutes?
for ind, t in enumerate(arange(0, 24 * 3600, 120)):
move_earth(120)
if ind % 100 == 0:
plot_sun_earth()
This is still not a very large movement, let's do our simulation over a half year and see what happens.
%%time
dt = 240.
for ind, t in enumerate(arange(0, 182.5 * 24 * 3600, dt)):
move_earth(dt)
if ind % (24 * 3600 / dt) == 0:
plot_sun_earth()
Wall time: 6.49 s
The earth circles around the sun. This is pretty nice!