De differentiaalvergelijkingen voor de enkele en dubbele (gekoppelde) slinger kunnen alleen in de harmonische benadering exact opgelost worden. Gebruikmakend van numerieke methoden kunnen deze echter met willekeurig kleine fout opgelost worden voor elke waarde van $ \alpha$ en $ \beta$.

De methode die in de simulatie gebruikt wordt om numeriek te integreren heet het velocity Verlet algoritme. Dit algoritme is waarschijnlijk het meest gebruikte algoritme. Eerst wordt het gewone Verlet algoritme besproken.

De bewegingsvergelijking van Newton voor een deeltje met massa $ m$ onderhevig aan een kracht $ F$ is:

$\displaystyle \ddot r=\frac{F(r)}{m}, \textrm{ with }r(t=0)=r_0,\, \dot r(t=0)=v_0 \nonumber $

De Taylor benaderingen voor $ r(t+h)$ en $ r(t-h)$ zijn:

$\displaystyle r(t+h)$ $\displaystyle = r(t) + h \dot r(t) + \frac{h^2}{2}\ddot r(t) + \ldots$ (3)
$\displaystyle r(t-h)$ $\displaystyle = r(t) - h \dot r(t) + \frac{h^2}{2}\ddot r(t) + \ldots$ (4)

Optellen van (3) en (4) geeft:

$\displaystyle r(t+h) + r(t-h) = 2r(t) + h^2\ddot r(t)+ O(h^4) $

herschikken:

$\displaystyle r(t+h)$ $\displaystyle = 2r(t) - r(t-h) + h^2\ddot r(t) + O(h^4)$    
  $\displaystyle = 2r(t) - r(t-h) + h^2 \frac{F(r)}{m} + O(h^4)$ (5)

Dit is het verlet algoritme.

De snelheid ( $ \dot r(t) = v(t)$), kan uitgerekend worden met:

$\displaystyle v(t)=\frac{r(t+h)-r(t-h)}{2h}+O(h^2)$ (6)

Het gebruik van vergelijking (6) heeft een nadeel. Het heeft een fout van $ O(h^2)$ in plaats van de $ O(h^4)$ fout van (5). Als de snelheid erg belangrijk is (bijvoorbeeld om de kinetische energie te bepalen om energie behoud te kunnen controlleren) is het gebruik van het zogenaamde velocity Verlet algoritme aan te raden.

Eerst worden $ r(t+h)$ en $ v(t+h/2)$ bepaald met de krachten op tijd $ t$:

$\displaystyle r(t+h)$ $\displaystyle = r(t) +hv(t)+\frac{h^2}{2}\frac{F(r(t))}{m}$    
$\displaystyle v(t+h/2)$ $\displaystyle = v(t) +\frac{h}{2} \frac{F(r(t))}{m}$    

daarna wordt de kracht op $ t+h$ berekend met de nieuwe positie en wordt deze gebruikt om $ v(t+\delta t)$ te berekenen:

$\displaystyle v(t+h) = v(t+h/2) + \frac{h}{2} \frac{F(r(t+h))}{m} \nonumber $

Deze vergelijkingen laten zich eenvoudig generaliseren, zodat ze ook te gebruiken zijn in systemen van bewegingsvergelijkingen.

Velocity Verlet met aanpassende stapgrootte

De stapgrootte moet klein genoeg gekozen worden om energiebehoud te garanderen gedurende de duur van het experiment. Het velocity Verlet algoritme met aanpassende stapgrootte wordt gebruikt om een geschikte stapgrootte te schatten. Dit algoritme schat de fout per integratiestap en past zijn stapgrootte hieraan aan.

In deze simulatie zal het algoritme toegepast worden op 1 periode van een enkele slinger. De kleinste stap die het algoritme kiest tijdens deze periode, zal gebruikt worden als stapgrootte voor de gehele simulatie (met de gekoppelde slingers). Het algoritme kan als volgt worden afgeleid:

Een simulatie van een deeltje met massa $ m$ onderhevig aan een kracht $ F$ start op tijd $ t_0$ met een bekende positie $ r(t_0)$. Eerst wordt 1 integratie stap gedaan om de positie $ r(t_1=t_0+h)$ met een fout $ O(h^4)$ te schatten, deze benaderde positie wordt $ r_h(t_1)$ genoemd. De fout in deze stap is gedefinieerd als:

$\displaystyle \epsilon_h = \vert r(t_1)-r_h(t_1)\vert $

Aangezien de orde van de fout $ h^4$ is, zal het nemen van een half zo grote stap resulteren in een fout die ongeveer $ 2^4=16$ keer kleiner is. Twee stappen zijn dan nodig om $ r_h(t_1)$ te schatten, dus de fout is:

$\displaystyle \epsilon_h \approx 8 \,\epsilon_{h/2} $

en ook:

$\displaystyle \epsilon_h \approx \frac{8}{7}(\epsilon_h - \epsilon_{h/2})= \frac{8}{7} \vert r_h(t_1)-r_{h/2}(t_1)\vert $

Er is dus een benadering voor $ \epsilon_h$ gevonden, ten koste van 2 extra stappen van $ h/2$.

Kies een maximale fout $ \epsilon_{\textrm{max}}$ voor de gehele simulatie. Beginnend op $ t_0$ wordt er 1 integratiestap $ h$ gedaan om $ r(t_1)$ te benaderen. De positie $ r(t_1)$ wordt ook benaderd door 2 stappen $ h/2$ te doen. Beide resultaten worden vergeleken om een schatting voor $ \epsilon_h$ te krijgen. Met deze fout wordt de maximaal te nemen stapgrootte $ h_{\textrm{max}}$ behorende bij $ \epsilon_{\textrm{max}}$ geschat:

$\displaystyle \left\vert \frac{\epsilon_{\textrm{max}}}{\epsilon_h}\right\vert ... ...} = \left\vert \frac{\epsilon_{\textrm{max}}}{\epsilon_h}\right\vert^{1/4} h $


(vanwege de $ O(h^4)$ fout).

Als $ h_{\textrm{max}}<h/2$, dan zal $ h$ worden vervangen door $ 0.9\times 2h_{\textrm{max}}$, en zal de integratie worden herhaald op $ t_0$. De nieuwe stapgrootte wordt vermenigvuldigd met 0.9 om er zeker van te zijn dat de stapgrootte altijd klein genoeg is.

Als $ h_{\textrm{max}}\geq h/2$ dan zal de stapgrootte klein genoeg zijn en zal de simulatie verder gaan naar $ t_0+h$ (waar $ h$ vervangen kan zijn bij de vorige stap).

In de meeste toepassingen van dit algoritme wordt de stapgrootte ook vergroot als deze te klein is vergeleken met $ h_{\textrm{max}}$. Bij deze simulatie wordt dit niet gedaan omdat we juist de kleinste stap $ h_{\textrm{max}}$ willen hebben die gebruikt is tijdens de simulatie. Dit zorgt er namelijk voor dat de simulatie (in het geval van de gekoppelde slingers) altijd wordt gedaan met een fout kleiner of gelijk aan $ \epsilon_{\textrm{max}}$. Elke keer dat parameters van de simulatie veranderd zijn voor het starten ervan, wordt het algoritme gebruikt om de stapgrootte te schatten.

De stapgrootte die nodig is om zeker te zijn van redelijk energiebehoud is zo klein dat de simulatie niet per integratiestap naar het scherm wordt getekend. De simulatie wordt ongeveer 10 keer per seconde getekend.