The differential equations for the single pendulum and double pendulum can be solved exactly only in the harmonic approximation. Numerical methods can be used to solve the equations for any value of $ \alpha$ and $ \beta$.

The integration method used in this simulation is called the velocity Verlet algorithm. This algorithm is perhaps the most widely used algorithm for numerical integration of ordinary differential equation. First the (common) Verlet algorithm is discussed.

The Newton equation of motion for a particle of mass $ m$ subjected to a force $ F$ is:

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

The Taylor approximations for $ r(t+h)$ and $ r(t-h)$ are :

$\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)

Adding (3) and (4):

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

and rearranging terms:

$\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)

gives the Verlet algorithm.

The velocity ( $ \dot r(t) = v(t)$), if needed, can be calculated as:

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

Equation (6) has one drawback, it has an error $ O(h^2)$ instead of the $ O(h^4)$ error of (5). If the velocity is important (for example for determining the kinetic energy to study energy conservation) the so called velocity Verlet algorithm is preferable.

First $ r(t+h)$ and $ v(t+h/2)$ are calculated from the forces at time $ 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}$    

then the force at time $ t+h$ is calculated with this new position and used to calculate $ v(t+\delta t)$ as:

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

These equations are easily generalised to a system of interacting particles.

Velocity Verlet with adaptive stepsize

The step size should be chosen small enough to ensure energy conservation. To estimate a suitable stepsize, velocity Verlet with adaptive stepsize is used. It estimates the error of integration at every integration step and adjusts its stepsize accordingly.

In this simulation the algorithm is only used for one period of oscillation for a single pendulum. The smallest stepsize encountered by the adaptive stepsize algorithm will be used to perform the coupled pendulum simulation. The adaptive stepsize algorithm is derived as follows.

A simulation of a particle with mass $ m$ subjected to a force $ F$ starts at $ t_0$ with a known position $ r(t_0)$. One integration step is done to approximate the position $ r(t_1=t_0+h)$ with an error $ O(h^4)$, call this approximated position $ r_h(t_1)$. The error of this integration step is defined:

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

Because the order of the error is $ h^4$, taking half a stepsize results in an error approximately $ 2^4=16$ times smaller. Two steps are now needed to approximate $ r_h(t_1)$, so the error is:

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

and also:

$\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 $

So an approximation of $ \epsilon_h$ is found, at the cost of doing 2 extra steps of size $ h/2$.

Define a maximum error $ \epsilon_{\textrm{max}}$ for the simulation. Starting at $ t_0$ one integration step of size $ h$ is done to approximate $ r(t_1)$. The position $ r(t_1)$ is also approximated by doing 2 steps of size $ h/2$. Both results are compared to get an estimate of $ \epsilon_h$. With this error the maximum stepsize $ h_{\textrm{max}}$ belonging to $ \epsilon_{\textrm{max}}$ is estimated:

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

because of the $ O(h^4)$ of the error.

If $ h_{\textrm{max}}<h/2$, then $ h$ is replaced by $ 0.9\times 2h_{\textrm{max}}$, and the integration is redone starting at $ t_0$. The new stepsize is multiplied by 0.9, to assure that the stepsize is always small enough.

If $ h_{\textrm{max}}\geq h/2$ then the stepsize is small enough and the simulation proceeds to $ t_0+h$ (where $ h$ could have been changed by the previous step).

In most simulations the stepsize is also adusted when it is too small compared to $ h_{\textrm{max}}$. It is not needed in this simulation because the algorithm is only used to find the smallest $ h_{\textrm{max}}$. This ensures that the simulation with the two coupled pendula is performed with an error smaller than, or equal to $ \epsilon_{\textrm{max}}$. Everytime parameters of the simulation are changed, the adaptive stepsize algorithm is used to estimate the stepsize.

The step size needed to meet the energy conservation condition is so small, that the simulation cannot be drawn at every step. About ten frames per second are drawn.