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):
        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()
    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())

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[5], line 17
     15 # Omdat we werken met Jup. Notebooks (en niet een .py file)
     16 from IPython.display import HTML
---> 17 HTML(ani.to_jshtml())

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\animation.py:1376, in Animation.to_jshtml(self, fps, embed_frames, default_mode)
   1372         path = Path(tmpdir, "temp.html")
   1373         writer = HTMLWriter(fps=fps,
   1374                             embed_frames=embed_frames,
   1375                             default_mode=default_mode)
-> 1376         self.save(str(path), writer=writer)
   1377         self._html_representation = path.read_text()
   1379 return self._html_representation

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\animation.py:1122, in Animation.save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
   1119 for data in zip(*[a.new_saved_frame_seq() for a in all_anim]):
   1120     for anim, d in zip(all_anim, data):
   1121         # TODO: See if turning off blit is really necessary
-> 1122         anim._draw_next_frame(d, blit=False)
   1123         if progress_callback is not None:
   1124             progress_callback(frame_number, total_frames)

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\animation.py:1158, in Animation._draw_next_frame(self, framedata, blit)
   1156 self._pre_draw(framedata, blit)
   1157 self._draw_frame(framedata)
-> 1158 self._post_draw(framedata, blit)

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\animation.py:1183, in Animation._post_draw(self, framedata, blit)
   1181     self._blit_draw(self._drawn_artists)
   1182 else:
-> 1183     self._fig.canvas.draw_idle()

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\backend_bases.py:1893, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
   1891 if not self._is_idle_drawing:
   1892     with self._idle_draw_cntx():
-> 1893         self.draw(*args, **kwargs)

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\backends\backend_agg.py:382, in FigureCanvasAgg.draw(self)
    379 # Acquire a lock on the shared font cache.
    380 with (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
    381       else nullcontext()):
--> 382     self.figure.draw(self.renderer)
    383     # A GUI class may be need to update a window using this draw, so
    384     # don't forget to call the superclass.
    385     super().draw()

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     92 @wraps(draw)
     93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94     result = draw(artist, renderer, *args, **kwargs)
     95     if renderer._rasterizing:
     96         renderer.stop_rasterizing()

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\figure.py:3257, in Figure.draw(self, renderer)
   3254             # ValueError can occur when resizing a window.
   3256     self.patch.draw(renderer)
-> 3257     mimage._draw_list_compositing_images(
   3258         renderer, self, artists, self.suppressComposite)
   3260     renderer.close_group('figure')
   3261 finally:

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\axes\_base.py:3226, in _AxesBase.draw(self, renderer)
   3223 if artists_rasterized:
   3224     _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
   3227     renderer, self, artists, self.get_figure(root=True).suppressComposite)
   3229 renderer.close_group('axes')
   3230 self.stale = False

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\axis.py:1404, in Axis.draw(self, renderer)
   1401     return
   1402 renderer.open_group(__name__, gid=self.get_gid())
-> 1404 ticks_to_draw = self._update_ticks()
   1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
   1407 for tick in ticks_to_draw:

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\axis.py:1288, in Axis._update_ticks(self)
   1286     tick.label1.set_text(label)
   1287     tick.label2.set_text(label)
-> 1288 minor_locs = self.get_minorticklocs()
   1289 minor_labels = self.minor.formatter.format_ticks(minor_locs)
   1290 minor_ticks = self.get_minor_ticks(len(minor_locs))

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\axis.py:1532, in Axis.get_minorticklocs(self)
   1530 minor_locs = np.asarray(self.minor.locator())
   1531 if self.remove_overlapping_locs:
-> 1532     major_locs = self.major.locator()
   1533     transform = self._scale.get_transform()
   1534     tr_minor_locs = transform.transform(minor_locs)

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\ticker.py:2230, in MaxNLocator.__call__(self)
   2228 def __call__(self):
   2229     vmin, vmax = self.axis.get_view_interval()
-> 2230     return self.tick_values(vmin, vmax)

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\ticker.py:2238, in MaxNLocator.tick_values(self, vmin, vmax)
   2235     vmin = -vmax
   2236 vmin, vmax = mtransforms.nonsingular(
   2237     vmin, vmax, expander=1e-13, tiny=1e-14)
-> 2238 locs = self._raw_ticks(vmin, vmax)
   2240 prune = self._prune
   2241 if prune == 'lower':

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\ticker.py:2168, in MaxNLocator._raw_ticks(self, vmin, vmax)
   2166 if self._nbins == 'auto':
   2167     if self.axis is not None:
-> 2168         nbins = np.clip(self.axis.get_tick_space(),
   2169                         max(1, self._min_n_ticks - 1), 9)
   2170     else:
   2171         nbins = 9

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\axis.py:2581, in XAxis.get_tick_space(self)
   2580 def get_tick_space(self):
-> 2581     ends = mtransforms.Bbox.unit().transformed(
   2582         self.axes.transAxes - self.get_figure(root=False).dpi_scale_trans)
   2583     length = ends.width * 72
   2584     # There is a heuristic here that the aspect ratio of tick text
   2585     # is no more than 3:1

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\transforms.py:466, in BboxBase.transformed(self, transform)
    462 """
    463 Construct a `Bbox` by statically transforming this one by *transform*.
    464 """
    465 pts = self.get_points()
--> 466 ll, ul, lr = transform.transform(np.array(
    467     [pts[0], [pts[0, 0], pts[1, 1]], [pts[1, 0], pts[0, 1]]]))
    468 return Bbox([ll, [lr[0], ul[1]]])

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\transforms.py:1496, in Transform.transform(self, values)
   1493 values = values.reshape((-1, self.input_dims))
   1495 # Transform the values
-> 1496 res = self.transform_affine(self.transform_non_affine(values))
   1498 # Convert the result back to the shape of the input values.
   1499 if ndim == 0:

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\transforms.py:2411, in CompositeGenericTransform.transform_affine(self, values)
   2409 def transform_affine(self, values):
   2410     # docstring inherited
-> 2411     return self.get_affine().transform(values)

File c:\Users\omloo\AppData\Local\Programs\Python\Python312\Lib\site-packages\matplotlib\transforms.py:2437, in CompositeGenericTransform.get_affine(self)
   2435     return self._b.get_affine()
   2436 else:
-> 2437     return Affine2D(np.dot(self._b.get_affine().get_matrix(),
   2438                            self._a.get_affine().get_matrix()))

KeyboardInterrupt:

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=[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()
    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
# 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):
    if particle.r[0] == 10.0:
        particle.r[0] = -10.0
    particle.update_position()
    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

# doorlopende doos
# 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):
    if particle.r[0] == 10.0 or particle.r[0] == -10.0:
        particle.v[0] = -particle.v[0]
    particle.update_position()
    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

# doorlopende doos
# Simulation parameters
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 3.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):
    if particle.r[0] > 10.0 or particle.r[0] < -10.0:
        particle.v[0] = -particle.v[0]
    if particle.r[1] > 10.0 or particle.r[1] < -10.0:
        particle.v[1] = -particle.v[1]
    particle.update_position()
    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[1] += a*dt

# Simulation parameters
dt = 0.1         # time step
num_steps = 500  # number of time steps
particle = ParticleClass(m=1.0, v=[5.0, 3.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):
    if particle.r[0] > 10.0 or particle.r[0] < -10.0:
        particle.v[0] = -particle.v[0]
    if particle.r[1] > 10.0 or particle.r[1] < -10.0:
        particle.v[1] = -particle.v[1]
    particle.update_position()
    particle.update_velocity(-9.81)
    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

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=[0, 0], 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.

# 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  


particle2 = ParticleClass(m=1.0, v=[3, 0], r=[0.0, 0.0],R=1.0)
particle = ParticleClass(m=1.0, v=[0, 0], r=[0.0, 0.0],R=1.0)
a = np.array([0.0, -5.0])  

track_x = []
track_y = []
track_x2 = []
track_y2 = []

# 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)
track_line2, = ax.plot([], [], 'b--', linewidth=1)

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

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

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

    particle2.update_position()
    particle2.update_velocity(a)

    track_x.append(particle.r[0])
    track_y.append(particle.r[1])
    track_line.set_data(track_x, track_y)

    track_x2.append(particle2.r[0])
    track_y2.append(particle2.r[1])
    track_line2.set_data(track_x2, track_y2)


    dot.set_data([particle.r[0]], [particle.r[1]])
    if (particle.r[0]+particle.R)**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]+particle.R)**2>100: # Check if particle is outside the bounds, np.abs could be used but is slower
        particle.v[1] = -particle.v[1]
    
    dot2.set_data([particle2.r[0]], [particle2.r[1]])
    if particle2.r[0]**2>100: # Check if particle is outside the bounds, np.abs could be used but is slower
        particle2.v[0] = -particle2.v[0]
    if particle2.r[1]**2>100: # Check if particle is outside the bounds, np.abs could be used but is slower
        particle2.v[1] = -particle2.v[1]
    return dot, track_line, dot2, track_line2

# 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())

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 21
     17 dt = 0.1         
     18 num_steps = 500  
---> 21 particle2 = ParticleClass(m=1.0, v=[3, 0], r=[0.0, 0.0],R=1.0)
     22 particle = ParticleClass(m=1.0, v=[0, 0], r=[0.0, 0.0],R=1.0)
     23 a = np.array([0.0, -5.0])  

Cell In[2], line 5, in ParticleClass.__init__(self, m, v, r, R)
      3 def __init__(self, m, v, r, R):
      4     self.m = m                  
----> 5     self.v = np.array(v, dtype=float)  
      6     self.r = np.array(r, dtype=float)  
      7     self.R = np.array(R, dtype=float)

NameError: name 'np' is not defined

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.

# Define a class for a particle

# Maken van de class met botsing
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):
        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 

# Simulation parameters
dt = 0.1         # time step
num_steps = 200  # number of time steps

particleA = ParticleClass(m=1.0, v=[2.5, 0], r=[-2.0, 0.0],R=0.45)  
particleB = ParticleClass(m=1.0, v=[-1, 0], r=[0.0, 0.0],R=0.45)  

track_x = []
track_y = []


# Creer de plot

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)  

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

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

# Updaten van de functie per frame
def update(frame):
    particleA.update_position()
    particleB.update_position()

    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]])

    # botsing tussen de deeltjes onderling
    if particleA.collide_detection(particleB):
        particleA.v[0] = -particleA.v[0]
        particleB.v[0] = -particleB.v[0]


    # botsing met de wand
    if particleA.r[0]**2>100:
        particleA.v[0] = -particleA.v[0]
    if particleA.r[1]**2>100: 
        particleA.v[1] = -particleA.v[1]

    dot.set_data([particleB.r[0]], [particleB.r[1]])
    if particleB.r[0]**2>100: 
        particleB.v[0] = -particleB.v[0]
    if particleB.r[1]**2>100: 
        particleB.v[1] = -particleB.v[1]

    return dot, track_line

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

# Voor Jupyter notebook:
from IPython.display import HTML
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                  # 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 collide_detection(self, other):
        return np.linalg.norm(self.r - other.r) <= (self.R + other.R)

# Simulation parameters
dt = 0.1         # time step
num_steps = 530  # number of time steps
m1 = 1.0
m2 = 1000.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)  


# 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")
dot, = ax.plot([], [], 'ro', markersize=10); # semicolon to suppress output

counter = 0

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

counter_text = ax.text(-9.5, 9, "")

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

# Update function for each frame
def update(frame):
    global counter
    particleA.update_position()
    particleB.update_position()

    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}")

    #collision detection and response
    if particleA.collide_detection(particleB):
        vA, vB, mA, mB, rA, rB = particleA.v, particleB.v, particleA.m, particleB.m, particleA.r, particleB.r
        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
#your code/answer


    # wall collision detection and response
    if particleA.r[0]**2>100: # Check if particle is outside the bounds, np.abs could be used but is slower
        particleA.v[0] = -particleA.v[0]
#your code/answer
    if particleA.r[1]**2>100: 
        particleA.v[1] = -particleA.v[1]
#your code/answer

    dot.set_data([particleB.r[0]], [particleB.r[1]])
    if particleB.r[0]**2>100: 
        particleB.v[0] = -particleB.v[0]
#your code/answer
    if particleB.r[1]**2>100: 
        particleB.v[1] = -particleB.v[1]
#your code/answer

    return dot, counter_text

# 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...