In Q2 werken we aan een deeltjes model. Dat model bestaat uiteindelijk uit een stuk of 100 deeltjes die in een afgesloten volume bewegen. Om voor elk deeltje apart de eigenschappen te bepalen, wordt enorm veel code. Omdat de deeltjes dezelfde eigenschappen hebben (massa, radius, snelheid, positie) kunnen we gebruik maken van classes.
Een class is een blauwdruk voor het creëren van objecten (zoals deeltjes, atomen, planeten, auto’s, studenten... noem maar op!). Het stelt je in staat om data (zoals positie, massa) en gedrag (zoals bewegen, botsen) te bundelen tot één overzichtelijk geheel. Precies wat we nodig hebben. Laten we de anatomie van een deeltje creëren.
import numpy as np
import matplotlib.pyplot as plt# Definieer een class voor een deeltje.
# Een class beschrijft alleen het type eigenschappen dat een ding heeft, bijvoorbeeld: in de beschrijving van de class zeg je:
# - een deeltje heeft een positie.
# - een deeltje heeft een functie die, wanneer aangeroepen, de positie van dat deeltje bijwerkt.
# Merk op dat de class zelf geen deeltje is!
# Een object dat tot een bepaalde class behoort, wordt een instant van die class genoemd.
# Je kunt meerdere instants van een class hebben.
class ParticleClass:
def __init__(self, m, v, r, R):
self.m = m # massa van het deeltje
self.v = np.array(v, dtype=float) # snelheids vector
self.r = np.array(r, dtype=float) # positie vector
self.R = np.array(R, dtype=float) # radius van het deeltje
def update_position(self,dt):
"""Werk de positie van het deeltje bij op basis van zijn snelheid en tijdstap dt."""
self.r += self.v * dt
print(ParticleClass)<class '__main__.ParticleClass'>
We hebben alleen de class aangemaakt, maar nog niet een object zelf. We hebben ook de mogelijkheid gemaakt om een nieuwe positie uit te rekenen: de volgende positie is de huidige positie + de snelheid maal de tijdstap, ofwel:
Laten we nu een ook echt een deeltje creëren en kijken wat we ermee kunnen doen.
# een enkel deeltje:
this_particle = ParticleClass(m=1.0, v=[5.0, 0], r=[0.0, 0.0], R=1.0)
print(this_particle)
%whos<__main__.ParticleClass object at 0x00000243D1E6ECF0>
Variable Type Data/Info
--------------------------------------------
N int 100000
ParticleClass type <class '__main__.ParticleClass'>
array ndarray 100000: 100000 elems, type `float64`, 800000 bytes (781.25 kb)
eind float 1763125506.4847884
eind_new float 1763125506.4861414
i int 99999
lengte float 0.765779972076416
lengte2 float 0.0015728473663330078
lengte_new float 0.0011203289031982422
np module <module 'numpy' from 'C:\<...>ges\\numpy\\__init__.py'>
particle int 9
particle_array list n=10
particle_object ParticleClass <__main__.ParticleClass o<...>ct at 0x00000243D3FF6580>
plt module <module 'matplotlib.pyplo<...>\\matplotlib\\pyplot.py'>
random module <module 'random' from 'C:<...>kfra8p0\\Lib\\random.py'>
random_array ndarray 100000: 100000 elems, type `float64`, 800000 bytes (781.25 kb)
random_float float 2.6615755401860985
random_float_2 float 0.7664111351594481
shorter float 683.5313896573739
start float 1763125505.7190084
start_new float 1763125506.485021
start_x int 9
that_particle ParticleClass <__main__.ParticleClass o<...>ct at 0x00000243D1D8D310>
this_particle ParticleClass <__main__.ParticleClass o<...>ct at 0x00000243D1E6ECF0>
time module <module 'time' (built-in)>
Nu we een object hebben, laten we eens kijken of we de eigenschappen van het object kunnen opvragen:
print("Massa van het deeltje " + str(this_particle.m) + " kg")
print("Huidige positie is " + str(this_particle.r))
# oproepen van de functie 'update_position' met dt = 1.0 seconde
this_particle.update_position(1.0)
print("Nieuwe positie " + str(this_particle.r))Massa van het deeltje 1.0 kg
Huidige positie is [0. 0.]
Nieuwe positie [5. 0.]
Let op! Run je bovenstaande cell opnieuw, dan zal de positie ook steeds veranderen!
Laten we nog een deeltje maken om te laten zien dat het aanroepen van update op het ene deeltje geen invloed heeft op het andere deeltje: het zijn onafhankelijke objecten!
that_particle = ParticleClass(m=1.0, v=[2.0, 0], r=[5.0, 0.0],R=1.0)print("Huidige positie van this_particle " + str(this_particle.r))
print("Huidige position van that_particle " + str(that_particle.r))
# oproepen van de functie 'update_position' met dt = 2.0 seconde
that_particle.update_position(2.0)
print("Nieuwe positie van this_particle " + str(this_particle.r))
print("Nieuwe positie van that_particle " + str(that_particle.r))Huidige positie van this_particle [5. 0.]
Huidige position van that_particle [5. 0.]
Nieuwe positie van this_particle [5. 0.]
Nieuwe positie van that_particle [9. 0.]
We kunnen ook een array van objecten maken. Dat zou eenvoudig moeten kunnen als het gaat om dezelfde deeltjes (zelfde massa, snelheid, maar andere startpositie).
particle_array = []
for start_x in range(4):
particle_array.append(ParticleClass(m=1.0, v=[2.0, 0], r=[start_x, 0.0],R=1.0))print("Hoe het array eruit ziet volgens python: " + str(particle_array))
print("Hoe een element van het array eruit ziet volgens python: " + str(particle_array[0]))
print("Hoe een eigenschap (in dit geval: 'r') van het element van het array eruit ziet volgens python: " + str(particle_array[0].r))Hoe het array eruit ziet volgens python: [<__main__.ParticleClass object at 0x00000243D1D8D310>, <__main__.ParticleClass object at 0x00000243D3F87A80>, <__main__.ParticleClass object at 0x00000243D3F87820>, <__main__.ParticleClass object at 0x00000243D1DA9FD0>]
Hoe een element van het array eruit ziet volgens python: <__main__.ParticleClass object at 0x00000243D1D8D310>
Hoe een eigenschap (in dit geval: 'r') van het element van het array eruit ziet volgens python: [0. 0.]
plt.figure()
for particle in particle_array:
plt.plot(particle.r[0],particle.r[1],'k.')
plt.show()
Wat eerder op moet zijn gevallen, en wat hierboven ook weer duidelijk wordt is dat de positie een 2D vector is. We kunnen dus onafhankelijk in de x en y richting bewegen!
Door gebruik te maken van een array waarin alle deeltjes zijn opgeslagen, kun je voor elk deeltje dezelfde bewerking uitvoeren, bijvoorbeeld allemaal een stukje updaten in de tijd! Hierbij kan je mooi gebruik maken van hoe Python loops maakt: als je één voor één de verschillende elementen van array af wilt gaan en daar dezelfde bewerking op doen, gaat dat zo:
for particle, particle_object in enumerate(particle_array):
print("Deeltje " + str(particle) + ", huidige positie: " + str(particle_object.r))
print("oproepen van de functie 'update_position' met dt = 1.0 second ")
particle_object.update_position(1.0)
print("Volgende positie " + str(particle_object.r))
if particle < len(particle_array) - 1:
print("Naar het volgende deeltje \n")
Deeltje 0, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [2. 0.]
Naar het volgende deeltje
Deeltje 1, huidige positie: [1. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [3. 0.]
Naar het volgende deeltje
Deeltje 2, huidige positie: [2. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [4. 0.]
Naar het volgende deeltje
Deeltje 3, huidige positie: [3. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [5. 0.]
Let op! In bovenstaande code maken we handig gebruik van twee programmeerconcepten.
enumerate hangt een nummer (counter) aan elk item in de array. Zo kun je bijhouden met welk deeltje je bezig bent.
Omdat onze code steeds aan het eind aangeeft dat het naar het volgende deeltje gaat, moet die code alleen stoppen bij het laatste deeltje.
We hebben hier gebruik gemaakt van het if statement.
Deze manier van werken is misschien niet heel efficient - het kost rekentijd - maar kan enorm helpen bij debuggen van code: je ziet snel welke tekst wel en welke tekst niet geprint wordt.
import random
# Clearing the array
particle_array = []
N = 10 # Amount of particles
for start_x in range(N):
random_float = random.uniform(-5.0, 5.0)
random_float_2 = random.uniform(0,1)
particle_array.append(ParticleClass(m=1.0, v=[random_float*np.cos(2*np.pi*random_float_2), random_float*np.sin(2*np.pi*random_float_2)], r=[0.0, 0.0],R=1.0))
plt.plot(particle_array[-1].r[0],particle_array[-1].r[1],'k.')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.title('Particles in the beginning')
plt.show()
# Since the speed is supposed to be between -5 and 5, and the velocity of the Particles is a 2D vector (so with a x and y component)
# The randomised velocity is put in Euler Coordinates
for particle, particle_object in enumerate(particle_array):
print("Deeltje " + str(particle) + ", huidige positie: " + str(particle_object.r))
print("oproepen van de functie 'update_position' met dt = 1.0 second ")
particle_object.update_position(1.0)
print("Volgende positie " + str(particle_object.r))
plt.plot(particle_object.r[0],particle_object.r[1],'k.')
if particle < len(particle_array) - 1:
print("Naar het volgende deeltje \n")
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.title('Particles Updated Plot')
plt.show()
Deeltje 0, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [-0.94316909 -1.25381056]
Naar het volgende deeltje
Deeltje 1, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [2.21245953 3.54931546]
Naar het volgende deeltje
Deeltje 2, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [0.82117006 0.3017578 ]
Naar het volgende deeltje
Deeltje 3, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [-1.62999136 1.24398312]
Naar het volgende deeltje
Deeltje 4, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [-0.77345394 2.51606872]
Naar het volgende deeltje
Deeltje 5, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [3.57405709 1.41096206]
Naar het volgende deeltje
Deeltje 6, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [-0.80787269 -0.00879554]
Naar het volgende deeltje
Deeltje 7, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [-2.00971428 3.68232068]
Naar het volgende deeltje
Deeltje 8, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [ 0.1835452 -0.0549715]
Naar het volgende deeltje
Deeltje 9, huidige positie: [0. 0.]
oproepen van de functie 'update_position' met dt = 1.0 second
Volgende positie [ 4.22992162 -0.37955853]
Onze simulaties worden snel “zwaar”. Er moeten veel berekeningen gedaan worden en dat kost nu eenmaal tijd. Nu is het wel zo dat sommige berekeningen meer tijd kosten dan anderen. Om te kijken of een stukje code geoptimaliseerd kan worden, moeten we een timer maken. Dat kan met onderstaande code:
import time
import numpy as np
start = time.time()
array = np.array([])
N = int(1e5)
for i in range(N):
array = np.append(array,np.random.rand(1))
eind = time.time()
lengte = eind - start
print("Dit kostte ", lengte, "seconde!")
#------------Nieuwe methode-----------------------
start_new = time.time()
%time random_array = np.random.rand(N)
print(random_array)
eind_new = time.time()
lengte_new = eind_new - start_new
print("Dit kostte ", lengte_new, "seconde!")
shorter= lengte/lengte_new
print('De nieuwe methode is ', shorter, ' keer korter.')
Dit kostte 1.733433485031128 seconde!
CPU times: total: 0 ns
Wall time: 1 ms
[0.38823666 0.44730676 0.80081593 ... 0.67386637 0.14540919 0.82981584]
Dit kostte 0.0010001659393310547 seconde!
De nieuwe methode is 1733.1458879618594 keer korter.
Idealiter draaien we de bovenstaande code een aantal keer zodat we een goede schatting hebben van de rekentijd.