04_Particle Dynamics.ipynb

 
 Particle Dynamics with Python

07 June 2021 | Jan H. Meinke
 import math
import random
import matplotlib.pyplot as plt
 
 %matplotlib widget
 
 Particle dynamics simulations are common in various scientific fields. They are used to simulate 
the formation of galaxies and the movements of molecules in a cell. Particles can have different
properties such as mass and charge and interact in different ways.
A classical particle dynamics code solves Newton's equation of motion:
 $F = m a,$ 
 $where F is the force, m the mass, and a the acceleration. F and a are vectors.$ In general, this problem is only solvable analytically for two particles. If there are more 
particles, we have to look for a numerical solution.
 $You may remember that you can calculate the velocity v of a particle as$  $v (t + d t) = v (t) + a (t) d t$ 
 $and the position r as$  $r (t + d t) = r (t) + v (t) d t + \frac{1}{2} a (t) d t^{2} .$ 
 $If we know all the positions, velocities and masses at time t and can calculate the forces, we can follow the motion of the particles over time.$ 
 Gravitational force
Let's assume our particles only interact via gravity. Then the force between two particles is given 
by
 $F_{ij} (t) = G \frac{m _{i} m _{j}}{r _{ij}^{2} ( t )} \hat{r}_{ij} (t),$ 
 $where F_{ij} (t) is the force on particle i due to particle j . r_{ij} (t) is the distance between particles i and j, and \hat{r}_{ij} (t) is the unit vector pointing from i to j .$  $To get the total force on particle i, we need to sum over all j \neq = i :$  $F_{i} (t) = \sum_{j \neq = i} F_{ij} (t) .$ 
 The algorithm
 
Calculate the force on each particle by summing up all the forces acting on it.
Integrate the equation of motion
 a) Calculate the position of each particle after a time step dt
 b) Calculate the velocity of each particle after a time step dt
 (Parallel) Patterns
In Think Vector, we got to know some patterns. Let's see how we can apply them here:
 $(i, j) -> F_{ij} : This is a map$  $$\mathbf F_{ij}$ -> : This is a reduction$ Calulate the new velocity:
    This is map
Calculate the new position:
    This is a map, too.
Now, let's try to express this in code.
    
 # Initialize positions and velocities
N = 50
L2 = 5  # half the length of the box size
epsilon = 1e-9  # softening factor
dt = 0.1  # time step
G = 1  # For simplicity we set the universal graviational constant to 1
m = 1  # This corresponds to 150 x 10^9 kg
x = [random.uniform(-L2, L2) for i in range(N)]
y = [random.uniform(-L2, L2) for i in range(N)]
z = [random.uniform(-L2, L2) for i in range(N)]
vx = [0 for i in range(N)]
vy = [0 for i in range(N)]
vz = [0 for i in range(N)]
 
 Calculating forces
 To calculate the force, we need the distance vector first. These are actually 3 maps (one for each component). The result is a distance matrix. As mentioned before maps are expressed as list generators:
 Dxx = [(i - j) for j in x for i in x]
Dyy = [(i - j) for j in y for i in y]
Dzz = [(i - j) for j in z for i in z]
D = [math.sqrt(i * i + j * j + k * k) for i, j, k in zip(Dxx, Dyy, Dzz)]
 
 Now that we have the vector components and the magnitude of the vector, we can calculate the forces.
 Fxx = [G * m * m * i / (d * d * d + epsilon) for i, d in zip(Dxx, D)]  # epsilon prevents a zero in the dominator.
Fyy = [G * m * m * i / (d * d * d + epsilon) for i, d in zip(Dyy, D)]
Fzz = [G * m * m * i / (d * d * d + epsilon) for i, d in zip(Dzz, D)]
Fx = [sum(Fxx[i * N: (i + 1) * N]) for i in range(N)]
Fy = [sum(Fyy[i * N: (i + 1) * N]) for i in range(N)]
Fz = [sum(Fzz[i * N: (i + 1) * N]) for i in range(N)]
 
 Let's visualize the forces on the particles:
 ax = plt.figure(figsize=(6, 6)).add_subplot(projection='3d')
ax.scatter3D(x, y, z)
ax.quiver(x, y, z, Fx, Fy, Fz)
 
 Integrating the equation of motion
 We are ready to update the positions and velocities of our particles:
 x = [i + v * dt + 0.5 * f / m * dt * dt for i, v, f in zip(x, vx, Fx)]
y = [i + v * dt + 0.5 * f / m * dt * dt for i, v, f in zip(y, vy, Fy)]
z = [i + v * dt + 0.5 * f / m * dt * dt for i, v, f in zip(z, vz, Fz)]
 
 vx = [v + f / m * dt for v, f in zip(vx, Fx)]
vy = [v + f / m * dt for v, f in zip(vy, Fy)]
vz = [v + f / m * dt for v, f in zip(vz, Fz)]
 
 Let's take a look at the particle positions and velocities:
 ax = plt.figure(figsize=(6, 6)).add_subplot(projection='3d')
ax.scatter3D(x, y, z)
ax.quiver(x, y, z, vx, vy, vz)
 
 That's it. By going back to the calculation of the forces, we can follow the motion of the particles over time.
 Exercise
Rewrite the program in a vectorized manner using ndarrays.
 
 
 Solution: