Double elastic pendulum problem solver

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 elastic pendulum.

The equations that will be integrated here are derived in this article.

Figure 1: Diagram of the double elastic pendulum.

The ordinary differential equation system being solved is

\begin{aligned} \begin{bmatrix} \ddot{r}_1 \\ \ddot{r}_2 \\ \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} &= \begin{bmatrix} 1 & \dfrac{m_2\cos{\Delta}}{m_1+m_2} & 0 & -\dfrac{m_2r_2\sin{\Delta}}{m_1+m_2} \\ \cos{\Delta} & 1 & r_1\sin{\Delta} & 0 \\ 0 & \dfrac{m_2\sin{\Delta}}{(m_1+m_2)r_1} & 1 & \dfrac{m_2r_2\cos{\Delta}}{(m_1+m_2)r_1} \\ -\dfrac{\sin{\Delta}}{r_2} & 0 & \dfrac{r_1\cos{\Delta}}{r_2} & 1 \end{bmatrix}^{-1} \begin{bmatrix} r_1\dot{\theta}_1^2-g\sin{\theta_1} + \dfrac{m_2}{m_1+m_2}\left(r_2\dot{\theta}_2^2\cos{\Delta} + 2\dot{r}_2\dot{\theta}_2\sin{\Delta}\right) + \dfrac{Q_{r_1}-k_1(r_1-l_1)}{m_1+m_2}\\ r_2\dot{\theta}_2^2-g\sin{\theta_2} + r_1\dot{\theta}_1^2\cos{\Delta} - 2\dot{r}_1\dot{\theta}_1\sin{\Delta} + \dfrac{Q_{r_2}-k_2(r_2-l_2)}{m_2} \\ -\dfrac{2\dot{r}_1\dot{\theta}_1}{r_1} - \dfrac{g\cos{\theta_1}}{r_1} - \dfrac{m_2}{(m_1+m_2)r_1}\left[2\dot{r}_2\dot{\theta}_2\cos{\Delta} -r_2\dot{\theta}_2^2\sin{\Delta}\right] + \dfrac{Q_{\theta_1}}{(m_1+m_2)r_1^2} \\ -\dfrac{2\dot{r}_2\dot{\theta}_2}{r_2}- \dfrac{g\cos{\theta_2}}{r_2} - \dfrac{2\dot{r}_1\dot{\theta}_1\cos{\Delta}}{r_2} - \dfrac{r_1\dot{\theta}_1^2\sin{\Delta}}{r_2} + \dfrac{Q_{\theta_2}}{m_2r_2^2} \end{bmatrix}. \end{aligned}

Simulation parameter form.
Parameter	Value	Explanation
$g$ :		Gravitational acceleration in metres (m) per second (s) squared ( $\mathrm{m}\cdot \mathrm{s}^{-2}$ ).
$l_1$ :		Rest length (m) of pendulum 1.
$l_2$ :		Rest length (m) of pendulum 2.
$m_{1b}$ :		Mass (kilograms or kg) of pendulum bob 1.
$m_{1r}$ :		Mass (kilograms or kg) of pendulum rod 1.
$m_{2b}$ :		Mass of pendulum bob 2 in kilograms.
$m_{2r}$ :		Mass of pendulum rod 2 in kilograms.
$b_{1b}$ :		Linear dissipation coefficient of pendulum bob 1.
$b_{1r}$ :		Linear dissipation coefficient of pendulum rod 1.
$c_{1b}$ :		Quadratic dissipation coefficient of pendulum bob 1.
$c_{1r}$ :		Quadratic dissipation coefficient of pendulum rod 1.
$b_{2b}$ :		Linear dissipation coefficient for pendulum bob 2.
$b_{2r}$ :		Linear dissipation coefficient for pendulum rod 2.
$c_{2b}$ :		Quadratic dissipation coefficient for pendulum bob 2.
$c_{2r}$ :		Quadratic dissipation coefficient for pendulum rod 2.
$t_f$ :		End time (s) for the simulation.
$r_{1,0}$ :		Initial value of $r_1$ (m).
$\dot{r}_{1,0}$ :		Initial value of $\dot{r}_1$ in $\mathrm{m}\cdot \mathrm{s}^{-1}$ .
$r_{2,0}$ :		Initial value of $r_2$ (m).
$\dot{r}_{2,0}$ :		Initial value of $\dot{r}_2$ in $\mathrm{m}\cdot \mathrm{s}^{-1}$ .
$\theta_{1,0}$ :		Initial value of $\theta_1$ in radians (r).
$\dot{\theta}_{1,0}$ :		Initial value of $\dot{\theta}_1$ in $\mathrm{r}\cdot \mathrm{s}^{-1}$ .
$\theta_{2,0}$ :		Initial value of $\theta_2$ in r.
$\dot{\theta}_{2,0}$ :		Initial value of $\dot{\theta}_2$ in $\mathrm{r}\cdot \mathrm{s}^{-1}$ .
$\epsilon$ :		Error tolerance.
$\mathrm{Tol}_{\mathrm{type}}$ :		Tolerance type, can be either absolute (0) or relative (1).
$h_{\mathrm{Initial}}$ :		Initial step size.
$h_{\mathrm{Min}}$ :		Minimum allowed step size.
$\Delta t$ :		Time increment for skipping ahead in animation.
$t_1$ :		Time you want to skip ahead to in animation when you press the skip button.
$\mathrm{Width}$ :		Width (in px) of Plotly windows used for plotting and animation below.
$\mathrm{Height}$ :		Height (in px) of Plotly windows used for plotting and animation below.
$t_\mathrm{Scale}$ :		Proportion of animation time passed per real time. $t_{\mathrm{Scale}}=1.0$ means animation and real time match. $t_{\mathrm{Scale}}<1.0$ means the animation is going more slowly than real time. $t_{\mathrm{Scale}}>1.0$ means it is going more rapidly.

© Brenton Horne. I do not claim to be infallible, so feel free to open issues or pull requests on GitHub if you notice flaws or things that can be improved.
I favour Oxford spelling, so if you notice I generally use different spelling to Australian spelling, that is probably why.
Homepage. Source code: https://github.com/fusion809/fusion809.github.io.
Website built with Franklin.jl.

Double elastic pendulum problem solver

Mathematical solver and related pages

Partial differential equation solver (non-interactive)

Ordinary differential equation solvers (interactive)

Complex analysis pages

Numerical mathematics

Lagrangian mechanics pages

Other mathematics

Education, my personal opinion largely

Linux