This webpage uses the
Runge-Kutta-Fehlberg fourth-order method with fifth-order error checking (RKF45) to approximate the solution to the problem of the
double pendulum (the equations were taken from
Science World): \[ \begin{aligned} \dfrac{d \theta_1}{dt} &= \dfrac{l_2 p_{\theta_1} - l_1 p_{\theta_2} \cos{\left(\theta_1-\theta_2\right)}}{l_1^2 l_2 \left(m_1 + m_2 \sin^2{\left(\theta_1 - \theta_2\right)}\right)} \\ \dfrac{d\theta_2}{dt} &= \dfrac{l_1 (m_1 + m_2)p_{\theta_2} - l_2 m_2 p_{\theta_1} \cos{\left(\theta_1-\theta_2\right)}}{l_1 l_2^2 m_2 \left(m_1 + m_2 \sin^2{\left(\theta_1-\theta_2\right)}\right)} \\ \dfrac{d p_{\theta_1}}{dt} &= -(m_1+m_2) gl_1 \sin{\theta_1} - C_1 + C_2 \\ \dfrac{d p_{\theta_2}}{dt} &= -m_2 gl_2 \sin{\theta_2} + C_1 - C_2 \end{aligned} \] where: \[ \begin{aligned} C_1 &= \dfrac{p_{\theta_1} p_{\theta_2} \sin{\left(\theta_1-\theta_2\right)}}{l_1 l_2 \left(m_1 + m_2 \sin^2{\left(\theta_1-\theta_2\right)}\right)} \\ C_2 &= \dfrac{l_2^2 m_2 p_{\theta_1}^2 + l_1^2 (m_1+m_2)p_{\theta_2}^2 - l_1 l_2 m_2 p_{\theta_1} p_{\theta_2} \cos{\left(\theta_1-\theta_2\right)}}{2l_1^2 l_2^2 \left(m_1 + m_2 \sin^2{\left(\theta_1-\theta_2\right)}\right)^2} \sin{2\left(\theta_1-\theta_2\right)} \end{aligned} \] and \(\theta_1\) is the angle (in radians) the first pendulum rod makes with the negative \(y\)-axis and \(\theta_2\) is the angle (in radians) the second pendulum rod makes with the negative \(y\)-axis. \(l_1\) and \(m_1\) are the length (in metres) and mass (in kilograms) of the first pendulum, respectively, and \(l_2\) (in metres) and \(m_2\) (in kilograms) are the length and mass of the second pendulum, respectively. \(g\) is the acceleration due to gravity in metres per second squared.