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.

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 inputs

Parameter	Value	Explanation
\(g\):		Acceleration due to gravity in \(\mathrm{m}\cdot \mathrm{s}^{-2}\).
\(l_1\):		Rest length of pendulum 1 in metres.
\(l_2\):		Rest length of pendulum 2 in metres.
\(m_1\):		Mass of pendulum bob 1 in kilograms.
\(m_2\):		Mass of pendulum bob 2 in kilograms.
\(k_1\):		Coefficient for first pendulum spring.
\(k_2\):		Coefficient for second pendulum spring.
\(b_1\):		Linear dissipation coefficient for pendulum 1.
\(b_2\):		Linear dissipation coefficient for pendulum 2.
\(c_1\):		Quadratic dissipation coefficient for pendulum 1.
\(c_2\):		Quadratic dissipation coefficient for pendulum 2.
\(t_0\):		Starting time for the simulation in seconds (s).
\(t_f\):		End time for the simulation in seconds.
\(r_{1,0}\):		Value of \(r_1\) at \(t_0\) in metres.
\(\dot{r}_{1,0}\):		Value of \(\dot{r}_1\) at \(t_0\) in \(\mathrm{m}\cdot \mathrm{s}^{-1}\).
\(r_{2,0}\):		Value of \(r_2\) at \(t_0\) in metres.
\(\dot{r}_{2,0}\):		Value of \(\dot{r}_2\) at \(t_0\) in \(\mathrm{m}\cdot \mathrm{s}^{-1}\).
\(\theta_{1,0}\):		Value of \(\theta_1\) at \(t_0\) in radians.
\(\dot{\theta}_{1,0}\):		Value of \(\dot{\theta}_1\) at \(t_0\) in radians per second.
\(\theta_{2,0}\):		Value of \(\theta_2\) at \(t_0\) in radians.
\(\dot{\theta}_{2,0}\):		Value of \(\dot{\theta}_2\) at \(t_0\) in radians per second.
\(\epsilon\):		Absolute error tolerance. Equation solver can fail if too small.
dtInitial:		Initial guess for step size.
dtMin:		Minimum limit for step size.
\(\Delta t\)		Time you want to add to animation.
\(t_1\)		Time you want to skip ahead to in the animation.

Tabulate the solution	Remove the solution table
Generate pendulum position plots	Remove pendulum position plots
Generate a time plot of all variables (+time derivatives)	Remove time plot of all variables (+time derivatives)
Generate a \(r_1\) against time plot	Remove \(r_1\) against time plot
Generate a \(\dot{r}_1\) against time plot	Remove \(\dot{r}_1\) against time plot
Generate a phase plot of \(\dot{r}_1\) vs \(r_1\)	Remove phase plot of \(\dot{r}_1\) vs \(r_1\)
Generate a \(r_2\) against time plot	Remove \(r_2\) against time plot
Generate a \(\dot{r}_2\) against time plot	Remove \(\dot{r}_2\) against time plot
Generate a phase plot of \(\dot{r}_2\) vs \(r_2\)	Remove phase plot of \(\dot{r}_2\) vs \(r_2\)
Generate a \(r_2\) against \(r_1\) phase plot	Remove \(r_2\) against \(r_1\) phase plot
Generate a \(\theta_1\) against time plot	Remove \(\theta_1\) against time plot
Generate a \(\dot{\theta}_1\) against time plot	Remove \(\dot{\theta}_1\) against time plot
Generate a phase plot of \(\dot{\theta}_1\) vs \(\theta_1\)	Remove phase plot of \(\dot{\theta}_1\) vs \(\theta_1\)
Generate a \(\theta_2\) against time plot	Remove \(\theta_2\) against time plot
Generate a \(\dot{\theta}_2\) against time plot	Remove \(\dot{\theta}_2\) against time plot
Generate a phase plot of \(\dot{\theta}_2\) vs \(\theta_2\)	Remove phase plot of \(\dot{\theta}_2\) vs \(\theta_2\)
Generate a \(\theta_2\) against \(\theta_1\) phase plot	Remove \(\theta_2\) against \(\theta_1\) phase plot
Generate all solution plots	Remove all plots
Generate an animation of the system	Remove system animation
Generate a \(r_1\) phase plot animation	Remove \(r_1\) phase plot animation
Generate a \(r_2\) phase plot animation	Remove \(r_2\) phase plot animation
Generate a \(\theta_1\) phase plot animation	Remove \(\theta_1\) phase plot animation
Generate a \(\theta_2\) phase plot animation	Remove \(\theta_2\) phase plot animation
Generate all animations	Remove all animations

The buttons below pertain to the animation immediately below them. Pause Restart Add time \(\Delta t\) to animation Skip to \(t_1\) Pause Restart Add time \(\Delta t\) to animation Skip to \(t_1\) Pause Restart Add time \(\Delta t\) to animation Skip to \(t_1\) Pause Restart Add time \(\Delta t\) to animation Skip to \(t_1\) Pause Restart Add time \(\Delta t\) to animation Skip to \(t_1\)