Recap - Thermolab

import numpy as np
import matplotlib.pyplot as plt

Recap en uitbreiding¶

Je hebt al gezien hoe je een botsend deeltjesmodel maakt. Je startte toen met 1 of 2 deeltjes en hebt dit uitgebreid naar ~100 deeltjes. Daarbij was het logisch om gebruik te maken van classes.

Voor het maken van een simulatie van gasmoleculen in een gesloten volume zal je logischerwijs de volgende stappen doorlopen:

bepalen van initiële condities
definiëren van de deeltjes
definiëren van volume met randvoorwaarden
simulaties bestaande uit:
- bepaling van positie van deeltjes
- controleren op onderlinge botsingen
- controleren op botsingen met wanden
- aanpassen van eigenschappen zoals totaal volume of temperatuur
- opslaan van gegevens

De initiële condities¶

In onze simulatie gaan we uit van een aantal aannames. Bij deze aannames wil je dat de simulatie een voorspellende waarde heeft voor een reëel systeem, maar ook dat de benodigde rekenkracht voor het uitvoeren van de simulatie op te brengen is. We kiezen daarom voor de volgende condities:

Ons model heeft een beperkt en constant aantal deeltjes.
We maken gebruik van een 2D simulatie.
De deeltjes voelen geen onderlinge kracht en ondergaan alleen elastische botsingen.
(tijdens simuleren) De tijdstap is voldoende klein. (d.w.z. de afgelegde weg $v_0*dt$ is klein in vergelijking met de diameter van de deeltjes, zodat we geen botsingen missen)

Hieronder geven we een aantal constanten voor onze code. Om de simulatie straks te kunnen interpreteren, vergelijken we deze eerst met een realistische situatie:

Solution to Exercise 1

Een stikstofmolecuul heeft een lengte van ongeveer 109 picometer. Een zuurstofmolecuul heeft een lengte van ongeveer 120 picometer. Voor de volume van een gasmolecuul wordt daarom genomen: (0.8109 + 0.2120)^3 = 1375036.928 pm^3 = 1.375x10^(-30) m^3 met de aanname dat lucht in een fietsband ongeveer voor 80% uit stikstofmoleculen en 20 % uit zuurstofmoleculen bestaat.```

‘’’ De simulatie die we gaan maken bevat maar twee dimensies. We kunnen het oppervlak per deeltje in deze simulatie daarom kiezen als het zojuist berekende $V_{molecuul}^{2/3}$ .

Voor de gemiddelde snelheid, $v$ , van een molecuul geldt dat

\frac{1}{2}m\left< v^2\right> = \frac{f}{2} k_B T.

(1)

Hierbij is $m$ de massa van het deeltje, $f$ het aantal dimensies, $k_B$ de constante van Boltzmann en $T$ de temperatuur.‘’’

Onze simulatie modelleert dus een klein systeem zonder rekening te houden met de kracht tussen de luchtmoleculen. Je zult zien dat we toch al behoorlijk wat van de processor vragen om deze berekeningen te doen. Om een meer realistische modelsimulatie te maken is daarom een serieuze uitdaging, die veel wiskunde en programmeerervaring vraagt. Op dit ogenblik wordt dit bijvoorbeeld gevraagd bij een van de eerste opdrachten van het vak ‘Computational Physics’ van de master Applied Physics.

# Ruimte voor uitwerking
BOX_SIZE_0 = 10                                         # Hoogte en breedte startvolume (in nanometer) 
N = 40                                                  # Aantal deeltjes 
V_0 = 1 #400 m/s                                         # Startsnelheid van deeltjes (in nanometer/2.5 picoseconden)
RADIUS = 0.3                                            # Straal van moleculen (in nanometer)
DT = 0.1*RADIUS/V_0                                  # Tijdstap om geen botsing te missen (in picoseconden)

#your code/answer

Particle class¶

Onze ParticleClass definieert van een deeltje de massa $m$ , snelheid $v$ , positie $r$ , en straal $R$ . De positie moet per tijdstap bepaald worden.

Nieuw aan ParticleClass zijn twee properties, nl. momentum en energy. Een property is eigenlijk een variabele die samenhangt met de waarde van andere variabelen binnen de klasse. Om de interne consistentie te bewaren moet je de juiste waarde van deze variabele dus telkens afleiden en dat doe je met een functie. Zo geeft kin_energy de kinetische energie terug van het deeltje door die af te leiden van de variabelen m en v.

Buiten ParticleClass voegen we twee functies toe, die je ook al hebt gezien. De eerste heet collide_detection en detecteert of twee deeltjes onderling overlappen en dus botsen. De tweede heet particle_collision en berekent de nieuwe snelheden van de twee deeltjes ten gevolge van hun botsing.

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

    def update_position(self):
        """ verandert positie voor één tijdstap """
        self.r += self.v * DT 
            
    @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)
    
def collide_detection(p1: ParticleClass, p2: ParticleClass) -> bool:
    """ Geeft TRUE als de deeltjes overlappen """
    """ Geeft TRUE als de deeltjes overlappen """
    dx = p1.r[0]-p2.r[0]
    dy = p1.r[1]-p2.r[1]
    rr = p1.R + p2.R
    return  dx**2+dy**2 < rr**2 
    #Note: np.linalg.norm(p1.r - p2.r) < (p1.R + p2.R) had ook gekund, maar is veel langzamer. Controleer zelf!

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

Volume en randvoorwaarden¶

Voor het volume maken we gebruik van het oppervlak van een pyplot, net zoals we dat in Q1 gedaan hebben. Daar is op dit ogenblik niets anders voor nodig dan een functie die de deeltjes laat botsen met de wanden.

Je hebt hiervoor al een voorbeeld van gezien, maar hieronder gebruiken we een andere vorm. Als je code ontwikkelt is het zaak om regelmatig en structureel te controleren of je code klopt. Het maakt daarbij niet uit of de code door jezelf is gemaakt, door een ander (zoals de docent) is aangeboden of van een generatieve AI komt. Hoe langer je code toevoegt zonder een controle uit te voeren, hoe groter de kans is dat je fout op fout stapelt en niet meer kan traceren wat er precies mis is.

We maken gebruik van een truc die in veel simulaties wordt gebruikt: we maken gebruik van periodieke randvoorwaarden. In dit geval wil dit zeggen dat een deeltje dat botst met de wand aan de ene kant van het volume, aan de andere kant verschijnt met exact dezelfde snelheid. Op deze manier verandert de snelheid van een deeltje niet bij het botsen met de wand en blijven alle behoudswetten gelden. Die behoudswetten kunnen we straks goed controleren.

def box_collision(particle: ParticleClass):
    ''' botsing met wanden geeft periodieke randvoorwaarden '''
    for i in range(2):
        if particle.r[i] > BOX_SIZE_0 / 2:
            particle.r[i] -= BOX_SIZE_0 
        elif particle.r[i] < -BOX_SIZE_0 / 2:
            particle.r[i] += BOX_SIZE_0

Maken van de simulatie¶

Nu dat we alle benodigde functies hebben, kunnen we ze in de juiste structuur zetten om de simulatie uit te voeren.

Bepaling positie deeltjes¶

Eerst moeten we alle deeltjes aanmaken en de startwaardes van hun variabelen invoeren:

def create_particles(particles):
    ''' Leegmaken en opnieuw aanmaken van deeltjes '''
    particles.clear
    for i in range(N):
        vx = np.random.uniform(-V_0, V_0)
        vy = np.random.choice([-1, 1]) * np.sqrt(V_0**2 - vx**2)        
        pos = np.random.uniform(-BOX_SIZE_0/2 + RADIUS, BOX_SIZE_0/2 - RADIUS, 2)
        particles.append(ParticleClass(m=1.0, v=[vx, vy], r=pos, R=RADIUS))

Controleren op onderlinge botsingen¶

Om de botsingen van alle deeltjes goed te verwerken moeten we controleren welke paren deeltjes met elkaar moeten botsen en deze botsingen één keer uitvoeren.

def handle_collisions(particles):
    global count
    """ alle onderlinge botsingen afhandelen voor deeltjes in lijst """
    num_particles = len(particles)
    count = 0
    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])
                count+=1
                print(count)

Controleren op botsingen met wanden¶

Voor het bepalen of deeltjes botsen met de wanden maken we gebruik van de functie die we net hebben aangepast.

def handle_walls(particles):
    ''' botsing met wanden controleren voor alle deeltjes in volume '''
    for p in particles:
        box_collision(p)

Integreren in tijdstap¶

Nu moeten we alle functies op de juiste manier combineren om de simulatie een tijdstap te laten maken:

def take_time_step(particles):
    for p in particles:
        p.update_position()
        
    handle_collisions(particles)
    handle_walls(particles)

Uitvoeren simulatie en tonen resultaat¶

Zoals aangegeven, kunnen we een animatie maken van de positie en snelheid als functie van de tijd, maar we kunnen ook het eindresultaat tonen en interpreteren:

particles = []
create_particles(particles)
for i in range(100):
    take_time_step(particles)

plt.figure()
plt.xlabel('x')
plt.ylabel('y')
plt.gca().set_aspect('equal')
plt.xlim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
plt.ylim(-BOX_SIZE_0/2, BOX_SIZE_0/2)

for p in particles:
    plt.plot(p.r[0], p.r[1], 'k.', ms=25)
    plt.arrow(p.r[0], p.r[1], p.v[0], p.v[1], 
              head_width=0.05, head_length=0.1, color='red')
plt.show()

Toevoegen animatie¶

Zoals gezegd: het kost veel meer rekenkracht om een animatie te maken over de tijd. Soms kan het wel helpen om te zien of de simulatie aan alle verwachtingen voldoet.

from matplotlib.animation import FuncAnimation
from IPython.display import HTML

fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_ylim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_aspect('equal')
dot, = ax.plot([], [], 'ro', ms=13);


def init():
    dot.set_data([], [])
    return dot,

def update(frame):
    take_time_step(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())

Loading...

Controle code¶

Om te controleren of onze code werkt passen we nu de code aan.

#your code/answer
kinenergy = []
for i in range(100):
    take_time_step(particles)
    kinenergy.append(sum([p.kin_energy for p in particles]))

print(kinenergy)
plt.plot(kinenergy)
plt.xlabel('tijd')
plt.ylabel('kinetische energie')
print(kinenergy)

1
1
1
1
1
1
1
1
1
2
1
1
1
1
1
1
1
1
1
1
2
1
1
1
1
1
1
1
1
1
2
1
1
1
2
1
1
2
1
2
1
2
1
1
2
1
2
1
1
1
1
1
1
1
1
1
[np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.000000000000004), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000007), np.float64(20.000000000000007), np.float64(20.000000000000007), np.float64(20.000000000000007), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(20.0), np.float64(20.0), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996)]
[np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.000000000000004), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000007), np.float64(20.000000000000007), np.float64(20.000000000000007), np.float64(20.000000000000007), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.000000000000004), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(20.0), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(20.0), np.float64(20.0), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996), np.float64(19.999999999999996)]

#your code/answer

momentum = []
for i in range(100):
    take_time_step(particles)
    momentum.append(sum([np.linalg.norm(p.momentum) for p in particles]))
   # momentum.append(sum(p.momentum for p in particles))

print(momentum)

plt.plot(momentum)
# run meerdere keren bij niet behouden

1
1
1
1
1
1
1
2
1
1
1
2
1
1
1
1
1
1
2
1
2
1
1
1
1
1
1
1
1
2
1
1
1
1
1
1
1
1
2
1
1
1
1
1
1
[np.float64(36.97577563918075), np.float64(36.9601337367073), np.float64(36.9601337367073), np.float64(36.9601337367073), np.float64(36.9601337367073), np.float64(37.231824847010834), np.float64(37.231824847010834), np.float64(37.231824847010834), np.float64(37.231824847010834), np.float64(37.231824847010834), np.float64(37.231824847010834), np.float64(37.231824847010834), np.float64(37.23500246498563), np.float64(37.3193533259608), np.float64(37.3193533259608), np.float64(37.33198975562693), np.float64(37.33198975562693), np.float64(37.2830607293605), np.float64(37.2830607293605), np.float64(36.860310075412826), np.float64(36.860310075412826), np.float64(36.860310075412826), np.float64(36.87953516535438), np.float64(36.87953516535438), np.float64(36.87953516535438), np.float64(36.45806124870188), np.float64(36.30823647810613), np.float64(36.281225958777426), np.float64(36.12297914111312), np.float64(36.10992845591167), np.float64(36.10992845591167), np.float64(35.51683645536223), np.float64(35.51683645536223), np.float64(35.77980522433056), np.float64(35.77980522433056), np.float64(35.76519668363307), np.float64(35.76519668363307), np.float64(35.76519668363307), np.float64(35.76519668363307), np.float64(35.93170412421624), np.float64(35.93170412421624), np.float64(35.93170412421624), np.float64(35.885559263477525), np.float64(35.885559263477525), np.float64(35.885559263477525), np.float64(35.885559263477525), np.float64(35.885559263477525), np.float64(35.885559263477525), np.float64(35.885559263477525), np.float64(36.06057638950535), np.float64(36.06057638950535), np.float64(36.06057638950535), np.float64(36.06057638950535), np.float64(35.98661076139177), np.float64(35.98661076139177), np.float64(35.98661076139177), np.float64(35.98661076139177), np.float64(35.98661076139177), np.float64(35.98661076139177), np.float64(36.201101427388764), np.float64(36.201101427388764), np.float64(36.159160917966005), np.float64(36.16820595064458), np.float64(36.16820595064458), np.float64(36.16820595064458), np.float64(36.2316786513024), np.float64(36.2316786513024), np.float64(36.2316786513024), np.float64(36.2316786513024), np.float64(36.2316786513024), np.float64(35.7088908736522), np.float64(35.698202302774945), np.float64(35.698202302774945), np.float64(35.70177253021195), np.float64(35.66635111348457), np.float64(35.66635111348457), np.float64(35.66635111348457), np.float64(36.04200093083378), np.float64(36.039564230319016), np.float64(36.039564230319016), np.float64(36.039564230319016), np.float64(36.039564230319016), np.float64(36.23525752052969), np.float64(36.23525752052969), np.float64(36.131472681771676), np.float64(36.131472681771676), np.float64(36.131472681771676), np.float64(36.131472681771676), np.float64(36.131472681771676), np.float64(36.0764282418637), np.float64(36.0764282418637), np.float64(36.0764282418637), np.float64(36.0764282418637), np.float64(36.07056225733888), np.float64(36.10833309713487), np.float64(36.10833309713487), np.float64(36.10833309713487), np.float64(36.10833309713487), np.float64(36.10833309713487), np.float64(36.10833309713487)]

Voor simulaties waarbij de kans op botsingen heel groot is (richting vloeistof), is er een handige methode om de snelheid van je simulatie op te voeren. In onderstaande code wordt daar een voorbeeld van gegeven.

def handle_collisions(particles):
    """ 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])

Exercise 11 (Controleren van beperkingen 🌶)

Leg uit wat deze code precies doet.

Deze code definieert de actie die genomen wordt bij een botsing. Er wordt gekeken naar een botsing tussen twee deeltjes. Hierbij wordt gebruik gemaakt van een for loop binnen een for loop: eerst wordt gekeken naar het eerste deeltje in de lijst in een botsing met elk ander deeltje in de lijst, daarna wordt gekeken naar het tweede deeltje in de lijst in een botsing met elk ander deeltje, en zo voort. Hierbij wordt elke keer dat er een botsing plaatsvindt, een functie uitgevoerd om de deeltjes te laten botsen door de richtingen om te draaien zodat de deeltjes van elkaar wegbewegen.

Maak een script waarbij twee deeltjes met dezelfde (absolute) snelheid op een derde deeltje afbewegen en tegelijk zouden botsen.

Laat zien dat ‘slechts’ twee deeltjes een interactie met elkaar aan gaan.
Leg uit waarom dat een acceptabele benadering is wanneer je veel deeltjes botsingen per tijdseenheid hebt.

Dit is een acceptabele benadering, omdat er zoveel botsingen zijn per tijdseenheid dat in zo een situatie, de verandering in snelheid van het derde deeltje zo klein is ten opzichte van de andere veranderingen in snelheid en richting die het de hele tijd ondervindt (van andere botsingen) dat het verwaarloosbaar is.

Voor wie een uitdaging wil: maak een grafiek met het aantal botsingen per tijdseenheid. Dit om bovenstaande benadering verder te verantwoorden.

BOX_SIZE_0 = 10
N = 40
V_0 = 1
RADIUS = 0.3
dt = 0.1*RADIUS/V_0

particleA = ParticleClass(m=1.0, v = [V_0, 0], r=[-2, 0], R=RADIUS)
particleB = ParticleClass(m=1.0, v = [-V_0, 0], r = [2, 0], R=RADIUS)
particleC = ParticleClass(m=1.0, v = [0, 0], r=[0, 0], R=RADIUS)
particle_array = [particleA, particleB, particleC]

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(p1: ParticleClass, p2: ParticleClass) -> bool:
    dx = p1.r[0]-p2.r[0]
    dy = p1.r[1]-p2.r[1]
    rr = p1.R + p2.R
    return dx**2+dy**2 < rr**2

def particle_collision(p1: ParticleClass, p2: ParticleClass):
    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)

    if dot_product >= 0:
        return
    
    distance_squared = np.dot(delta_r, delta_r)

    p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
    p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_r       

def box_collision(particle: ParticleClass):
    ''' botsing met wanden geeft periodieke randvoorwaarden '''
    for i in range(2):
        if particle.r[i] > BOX_SIZE_0 / 2:
            particle.r[i] -= BOX_SIZE_0 
        elif particle.r[i] < -BOX_SIZE_0 / 2:
            particle.r[i] += BOX_SIZE_0




def handle_collisions(particles):
    global count
    """ alle onderlinge botsingen afhandelen voor deeltjes in lijst """
    num_particles = len(particle_array)
    count = 0
    for i in range(num_particles):
        for j in range(i+1, num_particles):
            if collide_detection(particles[i], particles[j]):
                count+=1
                particle_collision(particles[i], particles[j])



fig, ax = plt.subplots()

ax.set_xlim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_ylim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_aspect('equal')

dot, = ax.plot([], [], 'ro', ms=13);
dotA, = ax.plot([], [], 'ro', ms=13)
dotB, = ax.plot([], [], 'bo', ms=13)
dotC, = ax.plot([], [], 'go', ms=13)

def init():
    dot.set_data([], [])
    return dot,

def update(frame):
    particleA.update_position()
    particleB.update_position()
    particleC.update_position()

    dotA.set_data([particleA.r[0]], [particleA.r[1]])
    dotB.set_data([particleB.r[0]], [particleB.r[1]])
    dotC.set_data([particleC.r[0]], [particleC.r[1]])

    handle_collisions(particle_array)
    return dot, 

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

# Deze animatie moet ook twee keer gerund worden soms voordat het werkt. A.U.B run dubbel.

Loading...

counter = []

for i in range(100):
    take_time_step(particles)
    counter.append(count)

print(counter)

plt.figure()
plt.plot(counter)

1
2
1
2
1
1
1
1
1
1
1
2
1
1
1
1
1
1
1
1
1
1
1
1
2
1
1
1
1
1
1
2
3
1
1
2
1
1
1
1
1
1
2
1
1
1
1
1
1
1
1
2
1
2
[2, 0, 0, 0, 0, 2, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 2, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 3, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 2, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 2, 2]

Nog sneller

Er is een manier om je script nog sneller te maken. Die methode gaat wel voorbij aan wat we verwachten van eerstejaars, maar het idee is de moeite waard om te bekijken. Bij $N$ deeltjes zou je $\frac{N(N-1)}{2}$ mogelijke botsingen moeten bekijken. Door gebruik te maken van een ignore list maak je die mogelijk aantal botsingen al een stuk minder. Maar je zou ook er van uit kunnen gaan dat botsingen slechts plaats vinden met deeltjes die in de buurt zijn. Als je een grid maakt en elk deeltje in een grid zet, hoeft elk deeltje slechts in de naastgelegen velden te “kijken” of daar een deeltje aanwezig is waarmee een botsing kan plaats vinden.

In onderstaande figuur is dat weergegeven. Het lichtrode deeltje heeft potentieel maar negen mogelijke botsingen. Dat script kan verder geoptimaliseerd worden door van links naar rechts per grid te lopen, dan hoeven zelfs niet alle velden bekeken te worden.

In de volgende notebooks maken we veelal gebruik van numpy functies. In sommige gevallen is dat sneller, maar lang niet altijd. Het is good-practice om na te gaan wat sneller is. Aan de andere kant moet de tekst ook leesbaar blijven...

Als onderdeel om excellent te scoren verwachten we dat je een aantal checks tussendoor doet om code sneller te maken. In onderstaande testcode, gemaakt door Josh Sleijfer wordt een van die functies getest. Wat blijkt? Voor “kleine” vectoren is het sneller om geen gebruik te maken van Numpy maar voor grotere vectoren wel!

import numpy as np
import time

L = 2                # lengte van de vector
repeats = 1_000_000  # aantal herhalingen

rnd_vec1 = np.random.rand(L, repeats)
 
# test numpy dot product
t0 = time.perf_counter()
for vec1 in rnd_vec1.T:
    dot = np.dot(vec1, vec1)
t1 = time.perf_counter()

print(f"Numpy: Time taken for {repeats} dot products of length {L}: {t1 - t0:.4f} seconds")

# test no numpy dot product 
t0 = time.perf_counter()
for vec1 in rnd_vec1.T:
    dot = vec1[0] * vec1[0] + vec1[1] * vec1[1] # + vec1[2] * vec1[2] + vec1[3] * vec1[3] # don't forget to adjust!
t1 = time.perf_counter()

print(f"No Numpy: Time taken for {repeats} dot products of length {L}: {t1 - t0:.4f} seconds")

# test fully vectorized 
t0 = time.perf_counter()
dot = np.sum(rnd_vec1[0] * rnd_vec1[0] + rnd_vec1[1] * rnd_vec1[1])
t1 = time.perf_counter()
print(f"Fully vectorized: Time taken for {repeats} dot products of length {L}: {t1 - t0:.4f} seconds")