Een van de bekendste voorbeelden van botsende deeltjes in de natuur is Brownian motion. Fijn gemalen pollen in water lijken te dansen in willekeurige richting. Dit komt doordat de pollen worden geraakt door watermoleculen die in alle richtingen bewegen. Omdat de pollen veel zwaarder zijn dan watermoleculen, dus de beweging van de pollen is veel langzamer en minder “intens” dan die van de watermoleculen. Dit proces van willekeurige beweging door botsingen met kleinere deeltjes wordt Brownian motion genoemd en kunnen we simuleren op basis van ons (premature) botsingsmodel. Daarbij kunnen we ook gebruik maken van de zojuist geleerde manier van tracking van deeltjes, waarbij we een zowel het zware bolletjes als een enkel deeltje kunnen volgen.
Let op! We bestuderen hier nog geen thermische effecten, deze opdrachten zijn met name bedoeld om beter te begrijpen hoe het botsingsmodel in elkaar zit.
import numpy as np
import matplotlib.pyplot as plt# Maken van de ParticleClass
class ParticleClass:
# Het maken van het deeltje
def __init__(self, m, v, r, R, c):
self.m = m
self.v = np.array(v, dtype=float)
self.r = np.array(r, dtype=float)
self.R = np.array(R, dtype=float)
self.c = c
# Het updaten van de positie, eventueel met zwaartekracht
def update_position(self):
self.r += self.v * dt #+ 1/2 * a * dt**2
# Harde wand
def boxcollision(self):
if abs(self.r[0]) + self.R > Box_length:
self.v[0] = -self.v[0] # Omdraaien van de snelheid
self.r[0] = np.sign(self.r[0]) * (Box_length - self.R) # Zet terug net binnen box
return True
if abs(self.r[1]) + self.R > Box_length:
self.v[1] = -self.v[1]
self.r[1] = np.sign(self.r[1]) * (Box_length - self.R)
return True
@property
def momentum(self):
return self.m * self.v
@property
def kin_energy(self):
return 1/2 * self.m * np.dot(self.v, self.v)# Aanmaken van de randvoorwaarden en initiele condities
Box_size_0 = 10
Box_length_0 = Box_size_0/2
Box_length = Box_length_0 # De grootte van de box kan wijzigen!
# Particles
dt = 0.1
particles = []
N = 40
v_0 = 1
dt = 0.04
botsingen = []
counter = 0
# Aanmaken van deeltjes
for i in range(N-1):
vx = np.random.uniform(-v_0,v_0)
vy = np.random.choice([-1, 1])*np.sqrt(v_0**2-vx**2)
pos = Box_length_0*np.random.uniform(-1,1,2)
particles.append(ParticleClass(m=1.0, v=[vx, vy], r = pos, R=.5,c='blue'))
#extra particle toevoegen met speciale waardes SWA
particles.append(ParticleClass(m=20.0, v=[0, 0], r = [0, 0], R=.5,c='red'))
Er is een doos vol met deeltjes op willekeurige positie aangemaakt. We willen kijken waar de deeltjes zijn terechtgekomen. Hieronder staat dit weergegeven.
# Inspecteren van beginsituatie
plt.figure()
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-Box_length_0,Box_length_0)
plt.ylim(-Box_length_0,Box_length_0)
for particle, particle_object in enumerate(particles):
plt.plot(particle_object.r[0],particle_object.r[1],color=particle_object.c,marker='.')
plt.arrow(particle_object.r[0],particle_object.r[1],
particle_object.v[0],particle_object.v[1],
head_width=0.05, head_length=0.1, color='green')
plt.show()
{solution} ex-brownian-2 Laat de snelheid en de richting van elk deeltje zien.
dit wordt gedaan met de code:
De geeft een random float tussen en .
geeft een absolute waarden van en maakt het random of de snelheid negatief of positief is.
wordt gegeven door dus, is gelijk bij elk deeltje alleen met verscillende en componenten.
We gaan nu de functies van de simulatie weer aanroepen:
# Het bepalen of er een botsing plaats vindt
def collide_detection(self, other):
dx = self.r[0] - other.r[0]
dy = self.r[1] - other.r[1]
rr = self.R + other.R
return dx**2+dy**2 < rr**2
def particle_collision(p1: ParticleClass, p2: ParticleClass):
""" past snelheden aan uitgaande van overlap """
m1, m2 = p1.m, p2.m
delta_r = p1.r - p2.r
delta_v = p1.v - p2.v
dot_product = np.dot(delta_r, delta_v)
# Als deeltjes van elkaar weg bewegen dan geen botsing
if dot_product >= 0: # '='-teken voorkomt ook problemen als delta_r == \vec{0}
return
distance_squared = np.dot(delta_r, delta_r)
# Botsing oplossen volgens elastische botsing in 2D
p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_r
return True
def handle_collisions(particles):
#your code/answer
""" alle onderlinge botsingen afhandelen voor deeltjes in lijst """
num_particles = len(particles)
for i in range(num_particles):
for j in range(i+1, num_particles):
if collide_detection(particles[i], particles[j]):
particle_collision(particles[i], particles[j])
global counter
counter += 1
In onderstaande code geven we de code voor de simulatie en volgen we de positie van het zware deeltje.
# tracken van het zware deeltje
track_x = []
track_y = []
track_light_x = []
track_light_y = []
timedcounter = []
#your code/answer
#your code/answer
for i in range(400):
#your code/answer
for p in particles:
p.update_position() # Update positie
p.boxcollision() # Wandbotsing werkt per deeltje
#for n in range(len(particles)):
#if p.boxcollision or particle_collision == True:
# botsingen[n] += 1
handle_collisions(particles)
timedcounter.append(counter)
counter = 0
#your code/answer
track_x.append(particles[N-1].r[0])
track_y.append(particles[N-1].r[1])
track_light_x.append(particles[N-2].r[0])
track_light_y.append(particles[N-2].r[1])
#your code/answer
plt.figure()
plt.plot(track_x,track_y,'r')
plt.plot(track_light_x,track_light_y, 'green')
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.show(){solution} ex-brownian-3 Het zware deeltje, beweegt veel inder snel en heeft ook dus minder botsingen, het is minder chaotisch. Het lichte deeltje bevindt zich altijd op een willekeurige plek, en de beweging van beide deeltjes zijn random. Als niet alles wordt gerunt maar alleen de bovenstaande Code dan is het eindpunt van de vorige itteratie het startpunt van de volgende.
We zouden gevoel willen krijgen voor het aantal botsingen dat per tijdseenheid plaatsvindt. Elke keer dat er een botsing plaatsvindt, zou de counter met 1 omhoog moeten gaan. Idealiter wordt het aantal botsingen opgeslagen in een array zodat je het aantal botsingen als functie van de tijd kunt weergeven.
time = np.arange(0,400,1)
for i in range(400):
#your code/answer
for p in particles:
p.update_position() # Update positie
p.boxcollision() # Wandbotsing werkt per deeltje
handle_collisions(particles)
plt.figure()
plt.plot(time,timedcounter)
plt.xlim(0, 400)
plt.xlabel("time")
plt.ylabel("aantal botsingen")
plt.show()In zulke fysica modellen is de afgelegde weg (afstand tussen begin en eindpunt) van belang. Deze afgelegde weg zegt iets over de snelheid van difussie. Idealiter bekijken we een histogram. Maar voor een histogram hebben we veel deeltjes nodig.
# Maken van de ParticleClass
class ParticleClass:
# Het maken van het deeltje
def __init__(self, m, v, r, R, c):
self.m = m
self.v = np.array(v, dtype=float)
self.r = np.array(r, dtype=float)
self.R = np.array(R, dtype=float)
self.c = c
# Het updaten van de positie, eventueel met zwaartekracht
def update_position(self):
self.r += self.v * dt #+ 1/2 * a * dt**2
# Harde wand
def boxcollision(self):
if abs(self.r[0]) + self.R > Box_length:
self.v[0] = -self.v[0] # Omdraaien van de snelheid
self.r[0] = np.sign(self.r[0]) * (Box_length - self.R) # Zet terug net binnen box
return True
if abs(self.r[1]) + self.R > Box_length:
self.v[1] = -self.v[1]
self.r[1] = np.sign(self.r[1]) * (Box_length - self.R)
return True
@property
def momentum(self):
return self.m * self.v
@property
def kin_energy(self):
return 1/2 * self.m * np.dot(self.v, self.v)
# Aanmaken van de randvoorwaarden en initiele condities
Box_size_0 = 50
Box_length_0 = Box_size_0/2
Box_length = Box_length_0 # De grootte van de box kan wijzigen!
# Particles
dt = 0.1
particles = []
N = 362
v_0 = 1
dt = 0.04
botsingen = []
counter = 0
startpos = []
endpos = []
# Aanmaken van deeltjes
for i in range(N-1):
vx = np.random.uniform(-v_0,v_0)
vy = np.random.choice([-1, 1])*np.sqrt(v_0**2-vx**2)
pos = Box_length_0*np.random.uniform(-1,1,2)
startpos.append(pos)
particles.append(ParticleClass(m=1.0, v=[vx, vy], r = pos, R=.5,c='blue'))
#extra particle toevoegen met speciale waardes SWA
particles.append(ParticleClass(m=20.0, v=[0, 0], r = [0, 0], R=.5,c='red'))
# Inspecteren van beginsituatie
plt.figure()
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-Box_length_0,Box_length_0)
plt.ylim(-Box_length_0,Box_length_0)
for particle, particle_object in enumerate(particles):
plt.plot(particle_object.r[0],particle_object.r[1],color=particle_object.c,marker='.')
plt.arrow(particle_object.r[0],particle_object.r[1],
particle_object.v[0],particle_object.v[1],
head_width=0.25, head_length=1, color='green')
plt.show()
# Het bepalen of er een botsing plaats vindt
def collide_detection(self, other):
dx = self.r[0] - other.r[0]
dy = self.r[1] - other.r[1]
rr = self.R + other.R
return dx**2+dy**2 < rr**2
def particle_collision(p1: ParticleClass, p2: ParticleClass):
""" past snelheden aan uitgaande van overlap """
m1, m2 = p1.m, p2.m
delta_r = p1.r - p2.r
delta_v = p1.v - p2.v
dot_product = np.dot(delta_r, delta_v)
# Als deeltjes van elkaar weg bewegen dan geen botsing
if dot_product >= 0: # '='-teken voorkomt ook problemen als delta_r == \vec{0}
return
distance_squared = np.dot(delta_r, delta_r)
# Botsing oplossen volgens elastische botsing in 2D
p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_r
return True
def handle_collisions(particles):
#your code/answer
""" alle onderlinge botsingen afhandelen voor deeltjes in lijst """
num_particles = len(particles)
for i in range(num_particles):
for j in range(i+1, num_particles):
if collide_detection(particles[i], particles[j]):
particle_collision(particles[i], particles[j])
global counter
counter += 1
# tracken van het zware deeltje
track_x = []
track_y = []
track_light_x = []
track_light_y = []
timedcounter = []
#your code/answer
#your code/answer
for i in range(400):
#your code/answer
for p in particles:
p.update_position() # Update positie
p.boxcollision() # Wandbotsing werkt per deeltje
handle_collisions(particles)
timedcounter.append(counter)
counter = 0
#your code/answer
track_x.append(particles[N-1].r[0])
track_y.append(particles[N-1].r[1])
track_light_x.append(particles[N-2].r[0])
track_light_y.append(particles[N-2].r[1])
#your code/answer
plt.figure()
plt.plot(track_x,track_y, color=particles[N-1].c, marker=',')
plt.plot(track_light_x,track_light_y, color=particles[N-2].c, marker=',')
plt.xlim(-Box_length_0,Box_length_0)
plt.ylim(-Box_length_0,Box_length_0)
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.show()
time = np.arange(0,400,1)
plt.figure()
plt.plot(time,timedcounter)
plt.xlim(0, 400)
plt.xlabel("time")
plt.ylabel("aantal botsingen")
plt.show()
for i in range(N-1):
endpos.append(particles[i].r)
# print(particles[i].r)
distance = []
for i in range(N-1):
distance.append(np.sqrt((endpos[i][0] - startpos[i][0])**2 + (endpos[i][1] - startpos[i][1])**2))
distance_heavy = np.sqrt((particles[N-1].r[0] - 0)**2 + (particles[N-1].r[1] - 0)**2)
plt.figure()
plt.hist(distance, bins=40, density=True, alpha=0.7, label="Light particles")
plt.axvline(
distance_heavy,
color='red',
linestyle='--',
linewidth=2,
label='Heavy particle'
)
plt.xlabel("Distance")
plt.ylabel("Probability density")
plt.legend()
plt.show()
En nu we toch bezig zijn met twee verschillende deeltjes....
We kunnen twee “groepen” van deeltjes aanmaken, elk met een andere massa. Als we dan de zwaartekracht aan zetten, dan zouden we verwachten dat de lichtere deeltjes boven komen “drijven”.
valt even snel??
# Maken van de ParticleClass
class ParticleClass:
# Het maken van het deeltje
def __init__(self, m, v, r, R, c):
self.m = m
self.v = np.array(v, dtype=float)
self.r = np.array(r, dtype=float)
self.R = np.array(R, dtype=float)
self.c = c
# Het updaten van de positie, eventueel met zwaartekracht
def update_position(self):
self.r += self.v * dt + 0.5 * a * dt**2
self.v += a * dt
# Harde wand
def boxcollision(self):
if abs(self.r[0]) + self.R > Box_length_0:
self.v[0] = -self.v[0] # Omdraaien van de snelheid
self.r[0] = np.sign(self.r[0]) * (Box_length_0 - self.R) # Zet terug net binnen box
return True
if abs(self.r[1]) + self.R > Box_length:
self.v[1] = -self.v[1]
self.r[1] = np.sign(self.r[1]) * (Box_length - self.R)
return True
@property
def momentum(self):
return self.m * self.v
@property
def kin_energy(self):
return 1/2 * self.m * np.dot(self.v, self.v)
# Aanmaken van de randvoorwaarden en initiele condities
Box_size_0 = 50
Box_length_0 = Box_size_0/2
Box_length = Box_length_0 * 2 # De grootte van de box kan wijzigen!
# Particles
particles = []
N = 41
v_0 = 0.001
dt = 0.0004
botsingen = []
a = np.array([0, -50])
# Aanmaken van deeltjes
for i in range(N-1):
vx = np.random.uniform(-v_0,v_0)
vy = np.random.choice([-1, 1])*np.sqrt(v_0**2-vx**2)
pos = Box_length_0*np.random.uniform(-1,1,2)
startpos.append(pos)
particles.append(ParticleClass(m=1.0, v=[vx, vy], r = pos, R=.5,c='blue'))
#extra particles toevoegen met speciale waardes SWA
for i in range(N-1):
vx = np.random.uniform(-v_0,v_0)
vy = np.random.choice([-1, 1])*np.sqrt(v_0**2-vx**2)
pos = Box_length_0*np.random.uniform(-1,1,2)
startpos.append(pos)
particles.append(ParticleClass(m=1000.0, v=[vx, vy], r = pos, R=.5,c='red'))
# Inspecteren van beginsituatie
plt.figure(figsize=(5, 10))
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-Box_length_0,Box_length_0)
plt.ylim(-Box_length,Box_length)
for particle, particle_object in enumerate(particles):
plt.plot(particle_object.r[0],particle_object.r[1],color=particle_object.c,marker='.')
plt.arrow(particle_object.r[0],particle_object.r[1],
particle_object.v[0],particle_object.v[1],
head_width=0.5, head_length=1, color='green')
plt.show()
# Het bepalen of er een botsing plaats vindt
def collide_detection(self, other):
dx = self.r[0] - other.r[0]
dy = self.r[1] - other.r[1]
rr = self.R + other.R
return dx**2+dy**2 < rr**2
def particle_collision(p1: ParticleClass, p2: ParticleClass):
""" past snelheden aan uitgaande van overlap """
m1, m2 = p1.m, p2.m
delta_r = p1.r - p2.r
delta_v = p1.v - p2.v
dot_product = np.dot(delta_r, delta_v)
# Als deeltjes van elkaar weg bewegen dan geen botsing
if dot_product >= 0: # '='-teken voorkomt ook problemen als delta_r == \vec{0}
return
distance_squared = np.dot(delta_r, delta_r)
# Botsing oplossen volgens elastische botsing in 2D
p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_r
return True
def handle_collisions(particles):
#your code/answer
""" alle onderlinge botsingen afhandelen voor deeltjes in lijst """
num_particles = len(particles)
for i in range(num_particles):
for j in range(i+1, num_particles):
if collide_detection(particles[i], particles[j]):
particle_collision(particles[i], particles[j])
for i in range(10000):
#your code/answer
for p in particles:
p.update_position() # Update positie
p.boxcollision() # Wandbotsing werkt per deeltje
handle_collisions(particles)
plt.figure(figsize=(5, 10))
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-Box_length_0,Box_length_0)
plt.ylim(-Box_length,Box_length)
for particle, particle_object in enumerate(particles):
plt.plot(particle_object.r[0],particle_object.r[1],color=particle_object.c,marker='.')
# plt.arrow(particle_object.r[0],particle_object.r[1],
# particle_object.v[0],particle_object.v[1],
# head_width=0.5, head_length=1, color='green')
plt.show()
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
fig, ax = plt.subplots(figsize=(5, 10))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-Box_length_0,Box_length_0)
ax.set_ylim(-Box_length,Box_length)
dot, = ax.plot([], [], 'ro', ms=1);
def init():
dot.set_data([], [])
return dot,
def update(frame):
for i in range(100):
for p in particles:
p.update_position() # Update positie
p.boxcollision() # Wandbotsing werkt per deeltje
handle_collisions(particles)
dot.set_data([p.r[0] for p in particles], [p.r[1] for p in particles])
return dot,
ani = FuncAnimation(fig, update, frames=100, init_func=init, blit=True, interval=50)
HTML(ani.to_jshtml())