De Maxwell-Boltzmann snelheidsverdeling

Doelen¶

We hebben al gezien in het vorige werkblad dat de snelheid van de deeltjes en de temperatuur aan elkaar gerelateerd zijn. In dit werkblad gaan we dat verder bekijken.

Eerst nemen we alle delen over van de code die we opnieuw moeten gebruiken:

class voor particles
functies voor een lijst aan deeltjes

Daarna voegen we de code toe voor het bekijken van de dynamiek van de deeltjes:

We kijken naar de verschillende richtingen
We bepalen de kansverdeling van de snelheden van de deeltjes

Zie voor de verdere inhoudelijke achtergrond de Feynman lectures on Physics.

Laden van eerdere code¶

Eerst nemen we de pakketten van Python en de constanten voor de simulatie over. Voer je eigen getallen in, die je ook in de vorige werkbladen hebt gebruikt.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.optimize import curve_fit

BOX_SIZE_0 = 10              # Hoogte en breedte startvolume in nm
N = 40                         # Aantal deeltjes
V_0 = 410                        # Startsnelheid van deeltjes in nm/ns
RADIUS = 0.3                    # Straal van moleculen in nm
DT = 7.5e-5                        # Tijdstap om geen botsing te missen in ns
M = 28 * 1.66053873e-27    # Massa 'lucht molecuul' in kg

Dan maken we weer gebruik van de klasse voor het deeltje:

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 * self.v**2
    
def collide_detection(p1: ParticleClass, p2: ParticleClass) -> bool:
    """ Geeft TRUE als de deeltjes overlappen """
    return np.linalg.norm(p1.r - p2.r) < (p1.R + p2.R)


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

En we laten de deeltjes met de wanden botsen, zodat er sprake is van een druk en een temperatuur in een gesloten volume.

def box_collision(particle: ParticleClass):
    ''' botsing met harde wanden '''
    if abs(particle.r[0]) + particle.R > BOX_SIZE_0/2: 
        particle.v[0] = -particle.v[0]                                        # Omdraaien van de snelheid
        particle.r[0] = np.sign(particle.r[0]) * (BOX_SIZE_0/2 - particle.R)  # Zet terug net binnen box                 
    if abs(particle.r[1]) + particle.R > BOX_SIZE_0/2: 
        particle.v[1] = -particle.v[1]     
        particle.r[1] = np.sign(particle.r[1]) * (BOX_SIZE_0/2 - particle.R)

Functies aan een lijst van deeltjes¶

Waarbij we al deze functies uitvoeren en samennemen over een lijst met deeltjes:

def create_particles(particles):
    """ Leegmaken en opnieuw aanmaken van deeltjes  in lijst """
    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=M, v=[vx, vy], r=pos, R=RADIUS))
        
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])

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

def take_time_step(particles):
    """ zet tijdstap voor een lijst deeltjes en verwerk alle botsingen onderling en met wanden """
    for p in particles:
        p.update_position()
    handle_collisions(particles)
    handle_walls(particles)

Eerste simulatie ter controle¶

Zoals we inmiddels gewend zijn draaien we eerst een korte simulatie om te verifiëren of alle code correct is overgenomen:

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=10)
    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()

Equipartitiebeginsel¶

Dit beginsel is heel belangrijk voor de thermodynamica en stelt dat de energie in thermodynamisch evenwicht gelijk wordt verdeeld over de toegankelijke vrijheidsgraden. Laten we dit eerst verifiëren in onze simulatie.

# jouw antwoord

def create_uniform_particles(particles):
    """ Leegmaken en opnieuw aanmaken van deeltjes met uniforme snelheid in lijst """
    particles.clear()
    for i in range(N):
        vy = np.random.choice([-1, 1]) * V_0  #|vy|=v_0, omdat er geen snelheid in de x richting is en vy krijgt een radom positive of negatieve richting    
        pos = np.random.uniform(-BOX_SIZE_0/2 + RADIUS, BOX_SIZE_0/2 - RADIUS, 2) #random positie in de doos met een afstand minimaal de radius tot de wand
        particles.append(ParticleClass(m=1.0, v=[0, vy], r=pos, R=RADIUS)) #particle aanmaken

particles = []
timesteps = int(((BOX_SIZE_0/V_0)*2)/DT) #aantal tijdstappen berekenen om de deeltjes 2x op en neer te laten gaan
E_k = np.zeros((2,timesteps))
t_array = np.arange(0,timesteps*DT,DT)

create_uniform_particles(particles) #particles aanmaken

#simulatie uitvoeren en E_k voor iedere timestep berekenen
for i in range(timesteps):
    take_time_step(particles)
    E_k[0,i] = sum(p.kin_energy[0] for p in particles)
    E_k[1,i] = sum(p.kin_energy[1] for p in particles)
    
#E_k plotten tegen de tijd
plt.figure(figsize=(9,5))
plt.plot(t_array,E_k[0,:],label='x direction')
plt.plot(t_array,E_k[1,:],label='y direction')
plt.xlabel('$t \mathrm{(s)}$')
plt.ylabel('$E_k \mathrm{(J)}$')
plt.legend()
plt.show()

Aan het resultaat van je simulatie kan je zien dat de deeltjes zich inderdaad gedragen volgens equipartitiebeginsel. Afwijkingen van het gemiddelde zijn het gevolg van de statistiek en worden voor grotere aantallen deeltjes relatief kleiner.

De snelheidsverdeling¶

In het vorige werkblad hebben we al gezien dat er een relatie is tussen de snelheid van de deeltjes en de temperatuur. Je kan hierdoor al vermoeden dat de snelheid van de deeltjes wordt gegeven door een verdeling die van de temperatuur afhangt. De vraag is dus of we deze functie kunnen vinden.

Om de snelheidsverdeling van de deeltjes te bepalen kan je gebruik maken van de functie histogram van python. Laten we daarom een simulatie draaien waarin we vanuit een willekeurige beginsituatie 100 tijdstappen zetten en daarna de snelheidsverdeling plotten:

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

speeds = [np.linalg.norm(p.v) for p in particles]
counts, bins = np.histogram(speeds, bins=10, density='True')
plt.figure()
plt.xlabel('Speed')
plt.ylabel('Frequency')
plt.stairs(counts, bins, fill='True')
plt.show()

Voldoende statistiek¶

Als je de simulatie hierboven een aantal keer uitvoert, dan zal je zien dat de statistiek onvoldoende is om een reproduceerbaar antwoord te krijgen. We kunnen het aantal deeltjes toe laten nemen om meer statistiek te krijgen, maar dat kost heel veel rekenkracht. Het is een goedkopere oplossing om de statistiek te bepalen op verschillende momenten in de tijd en deze statistische resultaten te middelen. De onderstaand code laat de deeltjes 5000 tijdstappen zetten en noteert de snelheid van alle deeltjes op elke 250e tijdstap. Het histogram wordt dan bepaald door de snelheden bij elke 250e tijdstap samen te nemen.

Door de simulatie een aantal keer te draaien zie je dat de verdeling al een stuk stabieler wordt.

particles = []
num_bins = 10
tot_counts = np.zeros(num_bins, dtype=float)

plt.figure()
plt.xlabel('Speed')
plt.ylabel('Frequency')

create_uniform_particles(particles)

for i in range(5000):
    take_time_step(particles)
    if i % 250 == 0:
        speeds = [np.linalg.norm(p.v) for p in particles]
        partial_counts, bins = np.histogram(speeds, bins=num_bins, density='True', range=[0,3*V_0])
        tot_counts += partial_counts

norm_counts = tot_counts / 20
plt.stairs(norm_counts, bins, fill='True')
plt.show()

De mathematische vorm van de snelheidsverdeling¶

De meest algemene vorm van de snelheidsverdeling heeft de vorm $f_{2D}(\vec{v})$ . Dat wil zeggen dat er bij elke unieke combinatie van $x$ en $y$ -component van de snelheid een specifieke kans hoort. We weten echter al dat dit niet het geval kan zijn. De natuur heeft helemaal geen voorkeur voor richting en wij kunnen zelf kiezen hoe ons assenstelsel georiënteerd is. De kans is dus alleen afhankelijk van de modulus van de snelheid en onafhankelijk van de richting. We kunnen de verdeling daarom weergeven als $f_{2D}(v)$ (zonder pijl voor de vector).

We kunnen de verdeling nog scherper definiëren door te stellen dat de componenten van de snelheid onderling onafhankelijk zijn. De functie $f_{2D}$ is daarom te splitsen in aparte functies voor de $x$ en $y$ -richting. Combineren we dit met onze conclusie van de vorige paragraaf dan moet dus gelden dat we de functie kunnen splitsen in twee functies voor $x$ en $y$ die onderling precies hetzelfde zijn en ook hetzelfde als $f(v)$ :

f_{2D}(v)=f(v_x)f(v_y)

(1)

Laten we nu de snelheidsverdeling beschouwen langs een willekeurige richting $r$ , die een lineaire combinatie van de $x$ - en $y$ -richting is. Omdat dit een deelverzameling is van de twee-dimensionale snelheidsverdeling moet bovenstaande relatie hiervoor nog steeds gelden. Om de vorm van deze functie te vinden stellen we nu een differentiaalvergelijking op door te differentiëren naar $v_x$ :

\frac{d}{dv_x}f(v_r)=\frac{d}{dv_x}f(v_x)f(v_y)

(2)

Om de linkerkant uit te werken, passen we de kettingregel toe:

\left(\frac{\partial f(v)}{\partial v_r}\right)\left(\frac{\partial v_r}{\partial v_x}\right)=f(v_y)\frac{df(v_x)}{dv_x}

(3)

Met de Stelling van Pythagoras wordt dit:

\frac{v_x}{v_r}\left(\frac{\partial f(v_r)}{\partial v_r}\right)=f(v_y) \frac{df(v_x)}{dv_x}

(4)

Om deze differentiaalvergelijking op te lossen willen we een scheiding van variabelen toepassen. Daarvoor willen we eerst af van de variabele $v_y$ . Dat kan door de vergelijking te delen door $f(v_r)=f(v_x)f(v_y)$ en de $v_x$ naar de andere kant te brengen:

\frac{1}{v_r f(v_r)}\left(\frac{\partial f(v_r)}{\partial v_r}\right)=\frac{1}{v_x f(v_x)} \frac{df(v_x)}{dv_x}

(5)

Je kan precies dezelfde redenering ook opzetten vanuit de $y$ -coördinaat in plaats van de $x$ -coördinaat. De termen die je hier gevonden hebt zijn daarmee functies van verschillende en onderling onafhankelijke variabelen die toch hetzelfde antwoord geven. Ze moeten daarom constant zijn. Die constante noemen we $-2\alpha$ , omdat dit de formules verderop vereenvoudigt:

\frac{1}{v_r f(v_r)}\left(\frac{\partial f(v_r)}{\partial v_r}\right)=\frac{1}{v_x f(v_x)} \frac{df(v_x)}{dv_x}=\frac{1}{v_y f(v_y)} \frac{df(v_y)}{dv_y} = -2\alpha

(6)

Misschien herken je hier al de vorm die $f(v_x)$ moet hebben om hieraan te voldoen, maar om dit makkelijker te maken, kunnen we de vergelijking een beetje herschrijven:

\frac{df(v_x)}{f(v_x)}=-2\alpha v_x dv_x

(7)

We kunnen beide zijden nu rustig integreren, zodat de oplossing wordt gegeven door:

f(v_x)=A \exp\left(-\alpha v_x^2\right)

(8)

Bij het vak Multivariabele Analyse zal je deze integraal tegenkomen aan het einde van dit blok. Daar wordt bewezen dat de integraal onder deze functie precies “1” is als (we zeggen ook wel: de functie is genormaliseerd):

A = \sqrt{\frac{\alpha}{\pi}}

(9)

Om ook de waarde van $\alpha$ te bepalen, kunnen we nu eisen dat deze formule de juiste waarde moet geven voor het gemiddelde van het kwadraat van de snelheid in de $x$ -richting:

\left< v_x^2 \right> = \frac{kT}{m} = \sqrt{\frac{\alpha}{\pi}} \int_{-\infty}^{\infty} v_x^2 \exp\left(-\alpha v_x^2\right)dv_x

(10)

Partieel integreren levert op:

\alpha = \frac{m}{2kT},

(11)

zodat de snelheidsverdeling voor twee dimensies uiteindelijk de vorm heeft:

f_{2D}(v_x,v_y)=f(v_x)f(v_y)=\frac{m}{2\pi kT} \exp\left( -\frac{m(v_x^2+v_y^2)}{2kT} \right)

(12)

Om hieruit de snelheidsverdeling $f_{2D}(v)$ te bepalen, moeten we nog een extra stap nemen. Het aantal combinaties van $v_x$ en $v_y$ dat overeenkomt met de snelheid $v$ is namelijk afhankelijk van $v$ . Dit wordt gegeven door de cirkelomtrek met middelpunt $(v_x=0,v_y=0)$ en straal $v$ . Zodoende is de snelheidsverdeling voor de modus van de snelheid gegeven door:

f_{2D}(v)=\frac{mv}{kT} \exp\left( -\frac{mv^2}{2kT} \right)

(13)

#your code/answer
k = 1.380649e-23 #

#fitfunctie aanmaken op basis van f_2d(v) die hierboven gegeven is
def f_2d(v,T):
    return (M*v)/(k*T)*np.exp(-(M*v**2)/(2*k*T))

bins_center = np.convolve(bins, [0.5, 0.5], mode='valid') #het centrum van de bins uitrekenen

var, cov = curve_fit(f_2d, bins_center, norm_counts, p0=300, maxfev=10000) #curve fit uitvoeren

#arrays maken om de curve fit te plotten
v_test = np.linspace(0,max(bins),1000)
f_2d_fit = f_2d(v_test,*var)

#simulatie en curvefit plotten
plt.figure()
plt.xlabel('Speed')
plt.ylabel('Frequency')
plt.stairs(norm_counts, bins, fill='True', label='sim')
plt.plot(v_test,f_2d_fit,'r--',label=f'T={var[0]:.2e}')
plt.legend()
plt.show()

def create_particles(particles, m_array):
    """ Leegmaken en opnieuw aanmaken van deeltjes  in lijst """
    particles.clear()
    for m in m_array: #2x de onderstaande loop uitvoeren, 1x met een zware massa en 1x met een lichte massa
        for i in range(int(N/len(m_array))): #de eerste of de tweede helft van de particles selecteren
            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=m, v=[vx, vy], r=pos, R=RADIUS))

particles = []
num_bins = 10
tot_counts = np.zeros(num_bins, dtype=float)
light_mass = M
heavy_mass = 4*light_mass
tot_counts_light = 0
tot_counts_heavy = 0

create_particles(particles,[light_mass,heavy_mass]) #zware en lichte particles aanmaken

#simulatie uitvoeren
for i in range(5000):
    take_time_step(particles)
    if i % 250 == 0: #iedere 250 stappen uitvoeren
        #particles splitsen in die met een lichte of een zware massa
        particles_light = list(filter(lambda p: p.m == light_mass, particles))
        particles_heavy = list(filter(lambda p: p.m == heavy_mass, particles))
        #snelheid van de lichte en zwarre massas gerekenen
        speeds_light = [np.linalg.norm(p.v)*np.sqrt(light_mass) for p in particles_light]
        speeds_heavy = [np.linalg.norm(p.v)*np.sqrt(heavy_mass) for p in particles_heavy]
        #histogram maken van de snelheden van de lichte en zwarre massas
        partial_counts_light, bins_light = np.histogram(speeds_light, bins=num_bins, density='True', range=[0,3*V_0*np.sqrt(M)])
        partial_counts_heavy, bins_heavy = np.histogram(speeds_heavy, bins=num_bins, density='True', range=[0,3*V_0*np.sqrt(M)])
        #som van de histogrammen
        tot_counts_light += partial_counts_light
        tot_counts_heavy += partial_counts_heavy

#delen door het aantal keer dat de histogrammen berekend zijn om een gemidelde te berekenen
norm_counts_light = tot_counts_light / 20
norm_counts_heavy = tot_counts_heavy / 20

#histogrammen voor de snelheden van de lichte en de zware deeltjes over elkaar heen plotten
plt.figure()
plt.xlabel('Speed*sqrt(mass)')
plt.ylabel('Frequency')
plt.stairs(norm_counts_light, bins_light, fill='True', color='g', label='light')
plt.stairs(norm_counts_heavy, bins_heavy, fill='True', color='b', label='heavy')
plt.legend()
plt.show()

in plaats van botsingen met elkaar zijn het gelijke ladingen die elkaar afstoten

import itertools

#SI eenheden
BOX_SIZE_0 = 10e-9               # Hoogte en breedte startvolume in m
N = 40                         # Aantal deeltjes
V_0 = 410                        # Startsnelheid van deeltjes in m/s
RADIUS = 0.3e-9                     # Straal van moleculen in m
DT = 1e-15                         # Tijdstap om geen botsing te missen in s
M = 28 * 1.66053873e-27           # massa van een 'lucht atoom' in kg
Q = 1.602e-19                           #lading van de deeltjes

class ParticleClass_new:
    def __init__(self, m, v, r, R, q, F):
        """ maakt een deeltje (constructor) """
        self.m = m                         
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = R
        self.q = q #lading
        self.F = np.array(F, dtype='float') #netto kracht die op het deeltje staat

    def update_position(self):
        """ verandert positie voor één tijdstap """
        #positie en snelheid berekenen volgens de verlet intergratie methode om de uitdemping/versterking van de euler intergratie methode te verminderen
        a = self.F / self.m
        self.v += 0.5*a*DT
        self.r += self.v*DT + 0.5*a*DT**2
        self.v += 0.5*a*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 box_collision(particle: ParticleClass_new):
    ''' botsing met harde wanden '''
    if abs(particle.r[0]) + particle.R > BOX_SIZE_0/2: 
        particle.v[0] = -particle.v[0]                                        # Omdraaien van de snelheid
        particle.r[0] = np.sign(particle.r[0]) * (BOX_SIZE_0/2 - particle.R)  # Zet terug net binnen box                 
    if abs(particle.r[1]) + particle.R > BOX_SIZE_0/2: 
        particle.v[1] = -particle.v[1]     
        particle.r[1] = np.sign(particle.r[1]) * (BOX_SIZE_0/2 - particle.R)

def create_particles(particles):
    """ Leegmaken en opnieuw aanmaken van deeltjes  in lijst """
    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_new(m=M, v=[vx, vy], r=pos, R=RADIUS, q=Q, F=[0,0])) #particle aanmaken met lading Q en beginkracht (0,0)

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

#functie om de kracht op iedere particle te berekenen
def handle_F(particles):
    #voor alle particles F reseten naar (0,0)
    for p in particles:
        p.F = [0,0]
    #de kracht op 2 particles berekenen (een positief en een negatief) voor iedere mogenlijke combinatie van 2 particles
    for p in list(itertools.combinations(particles,2)):
        F = 8.988e9 * (p[0].r-p[1].r) * p[0].q*p[1].q/((np.linalg.norm(p[0].r - p[1].r))**3) #F berekenen volgens die schaalt met delta_r^-2
        p[0].F += F #kracht in de positieve richting voor het ene deeltje
        p[1].F -= F #kracht in de negatieve richting voor het andere deeltje
        

def take_time_step(particles):
    """ zet tijdstap voor een lijst deeltjes en verwerk alle botsingen onderling en met wanden """
    handle_F(particles) #als eerst voor alle deeltjes de netto kracht berekenen
    for p in particles:
        p.update_position()
    handle_walls(particles)  
    
#functie om de potentiele energie van ieder deeltje te berekenen
def potential_energy(particles):
    U = 0.0
    for p0, p1 in itertools.combinations(particles, 2):
        dr = p0.r - p1.r
        dist = np.linalg.norm(dr)
        U += 8.988e9 * p0.q * p1.q / dist #de potentiele energie berekenen door middel van de eerste afgeleide van de kracht tegen de afstand en zonder richting
    return U

particles = []
timesteps = int(((BOX_SIZE_0/V_0)*2)/DT)

E_k = np.zeros(timesteps)
U = np.zeros(timesteps)
E_tot = np.zeros(timesteps)
t_array = np.arange(0,timesteps*DT,DT)

create_particles(particles)

for i in range(timesteps):
    take_time_step(particles)
    #totale potentiele en kinetische energie berekenen, en de totale hoeveelheid mechanische energie
    E_k[i] = sum(p.kin_energy for p in particles)
    U[i] = potential_energy(particles)
    E_tot[i] = U[i] + E_k[i]

#totale, potentiele en kinetische energie plotten
plt.figure(figsize=(9,5))
plt.plot(t_array,E_k,'g',label='kin')
plt.plot(t_array,U,'b',label='pot')
plt.plot(t_array,E_tot,'r',label='tot')
plt.xlabel('$t \mathrm{(s)}$')
plt.ylabel('$E \mathrm{(J)}$')
plt.legend()
plt.show()

Het afnemen van de totale energie is waarschijnlijk door het gebruik van numirieke intergratie intergratie, andere intergratie methodes, zoals RK4 zullen beter zijn, maar die vereisen veel meer rekenwerk. Dit bewijst ook warom de potentiele energie redenlijk constant blijft, maar de kinetische energie wel afneemt, omdat die afhankelijk is van de snelheid.