Deeltjes model - Thermolab

Zoom je heel ver in, dan zie je deeltjes rond vliegen. Elk met een eigen massa, een eigen snelheid en richting. De deeltjes botsen onderling, wisselen energie uit. We zouden op basis van een botsingsmodel van deeltjes iets moeten kunnen leren over thermodynamica. In de thermodynamica gaat het dan om heel veel deeltjes. Maar laten we beginnen met twee botsende ‘deeltjes’.

Je raadt het misschien al... deze simulatie gaan we na bouwen! Daarbij maken we gebruik van Python classes uit het vorige hoofdstuk en de basics van simulaties geleerd in Q1. Zorg er dus voor dat je weet hoe dit werkt! We maken ook gebruik van plotten, en daarbovenop een animatie. Hoe de animatie precies werkt en hoe je die zelf maakt hoef je niet te weten. Je zou wel in staat moeten zijn de code te lezen.

In de onderstaande cell maken we de ParticleClass aan, en geven we enkele parameters van onze simulatie op.

# Importeren van libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Maken van de class
class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                         
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = np.array(R, dtype=float)  

    def update_position(self, dt):
        self.r += self.v * dt

# Simulation parameters
dt = 0.1                                                            # tijd stap
num_steps = 500                                                     # aantal te nemen stappen
particle = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0], R=1.0)    # het maken van ons deeltje

We hebben nu een deeltje met massa, een snelheid, een begin positie en een straal. We hebben ook al de stapgrootte bepaald!

We willen de beweging van dat deeltje straks bestuderen en moeten dus een plot maken:

# Creeer een figuur en de assen
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Animatie")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Toon het deeltje als een rode stip
dot, = ax.plot([], [], 'ro', markersize=10); # semicolon to suppress output

Bij het aanmaken van ons deeltje hebben we het deeltje een beginpositie en snelheid mee gegeven. Als we dan per tijdstap de positie bepalen en deze laten plotten en die plots achter elkaar plakken, dan krijgen we een animatie van het deeltje. Met FuncAnimation wordt die animatie voor ons gedaan.

# Initialisieren van de functie voor de animatie
def init():
    dot.set_data([], [])
    return dot,

# Updaten van de functie voor elk frame
def update(frame):
    particle.update_position(dt)
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,

# Creeer de animatie
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# Omdat we werken met Jup. Notebooks (en niet een .py file)
from IPython.display import HTML
HTML(ani.to_jshtml())

Loading...

Frames in FuncAnimation verwijst naar het totaal aantal frames dat gebruikt wordt. Interval naar de snelheid van de animatie, nl. 50 milliseconden ofwel 20 frames per seconde.

Merk op dat als we de laatste cel opnieuw runnen, het deeltje zich niet in de oorsprong bevindt. Dat is een ‘eigenaardigheid’ van Jupyter Notebooks. Het is nu beter om alle code in één cel te plaatsen (hieronder gedaan voor je), zodat we ervoor zorgen dat het deeltje altijd in de oorsprong begint.

# Simulation parameters
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 0], r=[1.0, 0.0],R=1.0)  

# Create the figure and axis
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Create the particle as a red dot
dot, = ax.plot([], [], 'ro', markersize=10)

# Initialization function for animation
def init():
    dot.set_data([], [])
    return dot,

# Update function for each frame
def update(frame):
    particle.update_position(dt)
    if particle.r[0] < -10:
        particle.r[0] = 10
    elif particle.r[0] > 10:
        particle.r[0] = -10
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,

# Create animation
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)
    
# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())

Loading...

Een van de dingen die je kunt opmerken, is dat het deeltje niet in zijn doos blijft.

# doorlopende doos

#your code/answer
# Simulation parameters
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0],R=1.0)  

# Create the figure and axis
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Create the particle as a red dot
dot, = ax.plot([], [], 'ro', markersize=10)

# Initialization function for animation
def init():
    dot.set_data([], [])
    return dot,

# Update function for each frame
def update(frame):
    particle.update_position(dt)
    if particle.r[0] < -10:
        particle.r[0] = 10
    elif particle.r[0] > 10:
        particle.r[0] = -10
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,

# Create animation
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())

Loading...

# doos met harde wanden

#your code/answer
# Simulation parameters
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0],R=1.0)  

# Create the figure and axis
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Create the particle as a red dot
dot, = ax.plot([], [], 'ro', markersize=10)

# Initialization function for animation
def init():
    dot.set_data([], [])
    return dot,

# Update function for each frame
def update(frame):
    particle.update_position(dt)
    if particle.r[0] < -10+particle.R:
        particle.v[0] *= -1
    elif particle.r[0] > 10-particle.R  :
        particle.v[0] *= -1
    if particle.r[1] < -10+particle.R:
        particle.v[1] *= 1
    elif particle.r[1] > 10-particle.R  :
        particle.v[1] *= -1
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,

# Create animation
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())

Loading...


#your code/answer
# Simulation parameters
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 3], r=[0.0, 0.0],R=1.0)  

# Create the figure and axis
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Create the particle as a red dot
dot, = ax.plot([], [], 'ro', markersize=10)

# Initialization function for animation
def init():
    dot.set_data([], [])
    return dot,

# Update function for each frame
def update(frame):
    particle.update_position(dt)
    if particle.r[0] < -10+particle.R:
        particle.v[0] *= -1
    elif particle.r[0] > 10-particle.R:
        particle.v[0] *= -1
    if particle.r[1] < -10+particle.R:
        particle.v[1] *= -1
    elif particle.r[1] > 10-particle.R:
        particle.v[1] *= -1
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,

# Create animation
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())

Loading...

Laten we teruggaan naar ons deeltje. Er is een functie om de positie bij te werken, hoewel de snelheid hetzelfde lijkt te blijven... kunnen we de snelheid veranderen door (bijvoorbeeld) de versnelling door de zwaartekracht?

# Maken van de class met versnelling
class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                  # mass of the particle
        self.v = np.array(v, dtype=float)  # velocity vector
        self.r = np.array(r, dtype=float)  # position vector
        self.R = np.array(R, dtype=float)  # radius of the particle

    def update_position(self):
        self.r += self.v * dt
    
    def update_velocity(self, a):
        self.v += a*dt
        
a = np.array([0, -9.81])

#your code/answe

# Simulation parameters
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 3], r=[0.0, 0.0],R=1.0)  

# Create the figure and axis
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Create the particle as a red dot
dot, = ax.plot([], [], 'ro', markersize=10)

# Initialization function for animation
def init():
    dot.set_data([], [])
    return dot,

# Update function for each frame
def update(frame):
    particle.update_position()
    if particle.r[0] < -10 + particle.R:
        particle.v[0] *= -1
    elif particle.r[0] > 10 - particle.R:
        particle.v[0] *= -1
    if particle.r[1] < -10 + particle.R:
        particle.v[1] *= -1
    elif particle.r[1] > 10 - particle.R:
        particle.v[1] *= -1
    particle.update_velocity(a)
    dot.set_data([particle.r[0]], [particle.r[1]])
    return dot,

# Create animation
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())

Loading...

Solution to Exercise 6

Jouw antwoord hier;

De “Sinking” Error (Overshoot): Omdat we werken met vaste tijdstappen ( $dt$ ), “schiet” het deeltje in de stap dat hij de vloer raakt vaak een klein stukje door de vloer heen voordat de code de botsing opmerkt. In de volgende stap wordt de snelheid wel omgeklapt, maar het deeltje begint zijn weg omhoog vanaf een lagere positie ( $r$ ) dan waar de vloer eigenlijk is. Hierdoor is de potentiële energie ( $m \cdot g \cdot h$ ) aan het begin van de stijging lager dan hij zou moeten zijn.

Bij een vrije val (zonder botsing) zorgt Forward Euler er paradoxaal genoeg vaak voor dat de energie juist een klein beetje toeneemt. Maar zodra er botsingen bij komen kijken, zorgt de “overshoot” (het deeltje dat in de muur zakt) er meestal voor dat er netto energie verloren gaat.

Je merkt op dat de energie niet behouden blijft. Het model is een benadering en door de “discretisatie” (het opdelen in stapjes $dt$ ) ontstaan er fouten in de energiebalans. In een ideaal model zou de kinetische energie voor en na de botsing exact gelijk moeten zijn, maar hier “lekt” er energie weg door de manier waarop de tijdstappen werken.

Een optie om de simulatie te verbeteren is de tijdstap $\Delta t$ kleiner te maken, maar dan hebben we ook meer geduld meer nodig - het aantal berekeningen schaalt met $\frac{1}{\Delta t}$ . Een tweede optie is een meer directe oplossing: We weten dat $a = const.$ en daarom weten we ook de bewegingsvergelijking van het deeltje!

Daarnaast willen we graag weten waar het deeltje is geweest, dat is in onderstaande code toegevoegd.

# Maken van de class met versnelling
class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                  
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = np.array(R, dtype=float)  

    def update_position(self):
        self.r += self.v * dt + 1/2 * a * dt**2  
    
    def update_velocity(self, a):
        """Update the particle's velocity."""
        self.v += a*dt

# Simulation parameters
dt = 0.1         
num_steps = 500  
particle = ParticleClass(m=1.0, v=[5, 3], r=[0.0, 0.0],R=1.0)  
a = np.array([0.0, -5.0])  

track_x = []
track_y = []

# creeeren van de plot en de assen
fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")
track_line, = ax.plot([], [], 'r--', linewidth=1)  

# creeeren van ons rode deeltje
dot, = ax.plot([], [], 'ro', markersize=10);

# initializeren van onze functie voor de animatie
def init():
    dot.set_data([], [])
    return dot,

# Update function for each frame
def update(frame):
    particle.update_position()
    particle.update_velocity(a)

    track_x.append(particle.r[0])
    track_y.append(particle.r[1])
    track_line.set_data(track_x, track_y)
    
    dot.set_data([particle.r[0]], [particle.r[1]])
    if particle.r[0]**2>100: # Check if particle is outside the bounds, np.abs could be used but is slower
        particle.v[0] = -particle.v[0]
    if particle.r[1]**2>100: # Check if particle is outside the bounds, np.abs could be used but is slower
        particle.v[1] = -particle.v[1]
    return dot, track_line

# Create animation
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())

Loading...

We hebben steeds slechts gewerkt met een enkel deeltje. Maar om de simulatie uit het filmpje te maken, hebben we twee deeltjes nodig.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m
        self.v = np.array(v, dtype=float)
        self.r = np.array(r, dtype=float)
        self.R = R

    def update_position(self, dt):
        self.r += self.v * dt

    def update_velocity(self, a, dt):
        self.v += a * dt

# Parameters
dt = 0.05
a = np.array([0, -9.81])
p1 = ParticleClass(m=1.0, v=[10.0, 5.0], r=[-5.0, 0.0], R=1.0)
p2 = ParticleClass(m=1.0, v=[-5.0, 15.0], r=[5.0, 0.0], R=1.0)

# Track data opslag
p1_x, p1_y = [], []
p2_x, p2_y = [], []

fig, ax = plt.subplots()
ax.set_xlim(-10, 10); ax.set_ylim(-10, 10)
ax.set_aspect('equal')

# Plot objecten
dot1, = ax.plot([], [], 'ro', markersize=8)
track1, = ax.plot([], [], 'r--', alpha=0.3)
dot2, = ax.plot([], [], 'ko', markersize=8)
track2, = ax.plot([], [], 'k--', alpha=0.3)

def init():
    dot1.set_data([], []); track1.set_data([], [])
    dot2.set_data([], []); track2.set_data([], [])
    return dot1, track1, dot2, track2

def update(frame):
    # Update Deeltje 1
    p1.update_position(dt)
    p1.update_velocity(a, dt)
    if abs(p1.r[0]) > 10 - p1.R: p1.v[0] *= -1
    if abs(p1.r[1]) > 10 - p1.R: p1.v[1] *= -1
    p1_x.append(p1.r[0]); p1_y.append(p1.r[1])
    dot1.set_data([p1.r[0]], [p1.r[1]])
    track1.set_data(p1_x, p1_y)

    # Update Deeltje 2
    p2.update_position(dt)
    p2.update_velocity(a, dt)
    if abs(p2.r[0]) > 10 - p2.R: p2.v[0] *= -1
    if abs(p2.r[1]) > 10 - p2.R: p2.v[1] *= -1
    p2_x.append(p2.r[0]); p2_y.append(p2.r[1])
    dot2.set_data([p2.r[0]], [p2.r[1]])
    track2.set_data(p2_x, p2_y)

    return dot1, track1, dot2, track2

ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=30)
HTML(ani.to_jshtml())

Loading...

We waren gebleven bij het maken van een botsingsmodel, waarbij we nu twee deeltjes hebben die onderhevig zijn aan zwaartekracht.

Laten we de zwaartekracht even vergeten en alleen 1D kijken.

Uitwerking Exercise 10: Voorwaarden voor botsen¶

1. Wiskundige voorwaarde Twee deeltjes botsen wanneer de afstand tussen hun middelpunten kleiner is dan of gelijk is aan de som van hun radii (). Om de rekenintensieve worteltrekking te vermijden, gebruiken we de gekwadrateerde afstand:

**2. Functie collide_detection** In de ParticleClass berekenen we het verschil tussen de positievectoren self.r en other.r. De functie retourneert True als de som van de gekwadrateerde verschillen in x en y kleiner is dan de gekwadrateerde som van de radii.

3. Verandering van richting Bij een botsing tussen twee deeltjes van gelijke massa wisselen we hun snelheidsvectoren (v) om. Dit simuleert een perfecte elastische botsing waarbij de deeltjes van richting veranderen en hun energie behouden.

4. Automatisering (ParticleArray) In plaats van losse variabelen gebruiken we een lijst (particle_array) om alle deeltjes-objecten in op te slaan. Met een geneste for-loop kunnen we efficiënt elk deeltje met elk ander deeltje in de lijst vergelijken, wat het model schaalbaar maakt voor grotere simulaties.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                                     # massa
        self.v = np.array(v, dtype=float)              # snelheidsvector
        self.r = np.array(r, dtype=float)              # positievector
        self.R = R                                     # straal (float)

    def update_position(self, dt):
        """Update de positie op basis van snelheid en tijdstap dt."""
        self.r += self.v * dt

    def collide_detection(self, other):
        """Check of dit deeltje botst met een ander deeltje."""
        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

# Simulatie parameters
dt = 0.1
num_steps = 200

# Aanmaken van deeltjes
particleA = ParticleClass(m=1.0, v=[2.5, 0.5], r=[-5.0, 0.0], R=0.45)
particleB = ParticleClass(m=1.0, v=[-1.0, 0.0], r=[0.0, 0.0], R=0.45)

# Voor de trackline
track_x, track_y = [], []

# Setup Figuur
fig, ax = plt.subplots()
ax.set_xlim(-10, 10); ax.set_ylim(-10, 10); ax.set_aspect('equal')
ax.set_title("Deeltjes Botsing Simulatie")

# Plot objecten
track_line, = ax.plot([], [], 'r--', linewidth=1, alpha=0.5)
dotA, = ax.plot([], [], 'ro', markersize=10)
dotB, = ax.plot([], [], 'bo', markersize=10)

def init():
    dotA.set_data([], []); dotB.set_data([], [])
    track_line.set_data([], [])
    return dotA, dotB, track_line

def update(frame):
    # Posities updaten
    particleA.update_position(dt)
    particleB.update_position(dt)

    # Botsing tussen deeltjes checken en afhandelen
    if particleA.collide_detection(particleB):
        v_temp = np.copy(particleA.v)
        particleA.v = particleB.v
        particleB.v = v_temp

    # Botsing met de wanden (reflectie)
    for p in [particleA, particleB]:
        if abs(p.r[0]) > 10 - p.R: p.v[0] *= -1
        if abs(p.r[1]) > 10 - p.R: p.v[1] *= -1

    # Visualisatie bijwerken
    track_x.append(particleA.r[0])
    track_y.append(particleA.r[1])
    track_line.set_data(track_x, track_y)
    
    dotA.set_data([particleA.r[0]], [particleA.r[1]])
    dotB.set_data([particleB.r[0]], [particleB.r[1]])

    return dotA, dotB, track_line

ani = FuncAnimation(fig, update, frames=range(num_steps), init_func=init, blit=True, interval=50)
HTML(ani.to_jshtml())

Loading...

Terug naar ons vraagstuk... we willen een simulatie waarbij een deeltje met een massa $m_1$ en snelheid $v_1$ op een andere stilstaand deeltje met massa $m_2$ botst. Deeltje twee beweegt naar een muur, botst tegen de muur en beweegt richting deeltje 1 en botst tegen dit deeltje. Hoe vaak vindt deze botsing plaats als functie van de massa verhouding $\frac{m_1}{m_2}$ ?

Daarvoor moeten we even terug naar het botsingsmodel zoals geleerd in Klassieke Mechanica. Bij elastische botsingen is zowel het impulsmoment als de kinetische energie behouden.

\sum_i \vec{p}_i^{before} = \sum_i \vec{p}_i^{after}

(1)

E_{kin, before} = E_{kin, after}

(2)

Voor een botsing met twee deeltjes levert dit een analytische oplossing:

\begin{align} \mathbf{v}'_1 &= \mathbf{v}_1-\frac{2 m_2}{m_1+m_2} \ \frac{\langle \mathbf{v}_1-\mathbf{v}_2,\,\mathbf{x}_1-\mathbf{x}_2\rangle}{\|\mathbf{x}_1-\mathbf{x}_2\|^2} \ (\mathbf{x}_1-\mathbf{x}_2), \\ \mathbf{v}'_2 &= \mathbf{v}_2-\frac{2 m_1}{m_1+m_2} \ \frac{\langle \mathbf{v}_2-\mathbf{v}_1,\,\mathbf{x}_2-\mathbf{x}_1\rangle}{\|\mathbf{x}_2-\mathbf{x}_1\|^2} \ (\mathbf{x}_2-\mathbf{x}_1) \end{align}

(3)

we vragen je natuurlijk niet om deze vergelijking zelf te schrijven in Python, maar die vergelijking operationaliseren we hieronder wel.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m                         
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = np.array(R, dtype=float)  

    def update_position(self):
        self.r += self.v * dt 

    def collide_detection(self, other):
        # Controleert of de afstand kleiner is dan de som van de radii
        return np.linalg.norm(self.r - other.r) <= (self.R + other.R)

# Simulatie parameters
dt = 0.1         
num_steps = 530  
m1 = 100.0       # Verander dit naar 100, 10000 etc. voor de Pi-proef
m2 = 1.0

particleA = ParticleClass(m=m1, v=[0, 0], r=[-4.0, 0.0], R=0.45)  
particleB = ParticleClass(m=m2, v=[-1, 0], r=[-2.0, 0.0], R=0.45)  

fig, ax = plt.subplots()
ax.set_xlim(-10, 10); ax.set_ylim(-10, 10); ax.set_aspect('equal')
ax.set_title("Particle Collision")

dotA, = ax.plot([], [], 'ro', markersize=10)
dotB, = ax.plot([], [], 'bo', markersize=10)
counter_text = ax.text(-9.5, 8.5, "Collisions: 0", fontsize=12, fontweight='bold')

counter = 0

def init():
    dotA.set_data([], []); dotB.set_data([], [])
    counter_text.set_text("Collisions: 0")
    return dotA, dotB, counter_text

def update(frame):
    global counter
    particleA.update_position()
    particleB.update_position()

    # Deeltje-deeltje botsing detectie en afhandeling
    if particleA.collide_detection(particleB):
        vA, vB, mA, mB, rA, rB = particleA.v, particleB.v, particleA.m, particleB.m, particleA.r, particleB.r
        # Analytische oplossing voor elastische botsing
        vA_new = vA - 2 * mB / (mA + mB) * np.dot(vA - vB, rA - rB) / (1e-12 + np.linalg.norm(rA - rB))**2 * (rA - rB)
        vB_new = vB - 2 * mA / (mA + mB) * np.dot(vB - vA, rB - rA) / (1e-12 + np.linalg.norm(rB - rA))**2 * (rB - rA)
        particleA.v = vA_new
        particleB.v = vB_new
        
        counter += 1 # OPHOGEN BIJ DEELTJESBOTSING

    # Botsing detectie met wanden (x en y)
    for p in [particleA, particleB]:
        if p.r[0]**2 > (10 - p.R)**2: 
            p.v[0] *= -1
            counter += 1 # OPHOGEN BIJ WANDBOTSING X
        if p.r[1]**2 > (10 - p.R)**2: 
            p.v[1] *= -1
            counter += 1 # OPHOGEN BIJ WANDBOTSING Y

    dotA.set_data([particleA.r[0]], [particleA.r[1]])
    dotB.set_data([particleB.r[0]], [particleB.r[1]])
    counter_text.set_text(f"Collisions: {counter}")

    return dotA, dotB, counter_text

ani = FuncAnimation(fig, update, frames=range(num_steps), init_func=init, blit=True, interval=20)

from IPython.display import HTML
HTML(ani.to_jshtml())

Loading...

Exercise 12¶

De uitspraak in de video klopt: wanneer de massaverhouding m1/m2 een macht van 100 is, dan is het totaal aantal botsingen gelijk aan de eerste N+1 cijfers van $\pi$ .

De grootste beperkende factor in ons huidige numerieke model is de tijdstap ( $\Delta t$ ).

Wanneer de massa erg groot is, wordt de snelheid van het kleine deeltje na een botsing extreem hoog. Als de tijdstap te groot is, verplaatst het deeltje zich in één stap zover dat het “door” de muur of “door” het andere deeltje heen springt zonder dat de collide_detection wordt geactiveerd.

Bij een zeer groot aantal botsingen stapelen kleine afrondingsfouten (floating point errors) zich op, waardoor de simulatie minder nauwkeurig wordt.

Omdat we een vaste gebruiken, moet de computer miljoenen stappen berekenen terwijl er soms lange tijd niets gebeurt, of juist te weinig stappen wanneer de snelheid hoog is.