Double pendulum 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 pendulum. See this article for the derivation of the equations of motion.

Figure 1: Diagram of the double pendulum.

The equations being solved in this webpage are

θ¨1=1(m1r12+m1b+m2b)l1[m2bl2(θ¨2cos(θ1θ2)+θ˙22sin(θ1θ2))+gcosθ1(m1r2+m2r+m1b+m2b)(b1b+c1bl1θ˙1)l1θ˙1+(b2b+c2bl12θ˙12+l22θ˙22+2l1l2θ˙1θ˙2cos(θ1θ2))(l1θ˙1+l2θ˙2cos(θ1θ2))+(b1r+c1rl1θ˙12)l1θ˙14+(b2r+c2rl12θ˙12+l22θ˙224+l2θ˙1θ˙2cos(θ1θ2))(l1θ˙1+l2θ˙2cos(θ1θ2)2)]θ¨2=1(m2r12+m2b)l2m2b2l2cos2(θ1θ2)(m1r12+m1b+m2b)[m2bcos(θ1θ2)(m1r12+m1b+m2b)[m2bl2θ˙22sin(θ1θ2)+gcosθ1(m1r2+m2r+m1b+m2b)(b1b+c1bl1θ˙1)l1θ˙1+(b2b+c2bl12θ˙12+l22θ˙22+2l1l2θ˙1θ˙2cos(θ1θ2))(l1θ˙1+l2θ˙2cos(θ1θ2))+(b1r+c1rl1θ˙12)l1θ˙14+(b2r+c2rl12θ˙12+l22θ˙224+l2θ˙1θ˙2cos(θ1θ2))(l1θ˙1+l2θ˙2cos(θ1θ2)2)]+m2b(l1θ˙12sin(θ1θ2)gcosθ2)14(b2r+c2rl12θ˙12+l22θ˙224+l1l2θ˙1θ˙2cos(θ1θ2))(2l1θ˙1cos(θ1θ2)+l2θ˙2)(b2b+c2bl12θ˙12+l22θ˙22+2l1l2θ˙1θ˙2cos(θ1θ2))(l1θ˙1cos(θ1θ2)+l2θ˙2)].\begin{aligned} &\ddot{\theta}_1 = \dfrac{-1}{\left(\dfrac{m_{1r}}{12} + m_{1b}+m_{2b}\right)l_1} \left[m_{2b}l_2 ( \ddot{\theta}_2\cos{(\theta_1-\theta_2)} +\dot{\theta}_2^2\sin{(\theta_1-\theta_2)}) + g \cos{\theta_1}\left(\dfrac{m_{1r}}{2} +m_{2r} +m_{1b} + m_{2b}\right) -(b_{1b} + c_{1b} l_1 \dot{\theta}_1)l_1 \dot{\theta}_1 \right.\\ &\left.+\left(b_{2b}+c_{2b}\sqrt{l_1^2 \dot{\theta}_1^2 + l_2^2 \dot{\theta}_2^2 +2l_1 l_2\dot{\theta}_1 \dot{\theta}_2 \cos{(\theta_1-\theta_2)}}\right)(l_1 \dot{\theta}_1 + l_2 \dot{\theta}_2 \cos{(\theta_1-\theta_2)}) +\left(b_{1r} + \dfrac{c_{1r}l_1 \dot{\theta}_1}{2}\right) \dfrac{l_1 \dot{\theta}_1}{4} \right.\\ &\left.+\left(b_{2r} + c_{2r}\sqrt{l_1^2 \dot{\theta}_1^2 + \dfrac{l_2^2 \dot{\theta}_2^2}{4} + l_2 \dot{\theta}_1 \dot{\theta}_2 \cos{(\theta_1 -\theta_2)}}\right)\left(l_1 \dot{\theta}_1 + \dfrac{l_2\dot{\theta}_2 \cos{\left(\theta_1 - \theta_2\right)}}{2}\right)\right]\\ &\ddot{\theta}_2 = \dfrac{1}{\left(\dfrac{m_{2r}}{12} + m_{2b}\right)l_2 - \dfrac{m_{2b}^2l_2\cos^2{(\theta_1-\theta_2)}}{\left(\dfrac{m_{1r}}{12} + m_{1b}+m_{2b}\right)}}\left[\dfrac{m_{2b}\cos{(\theta_1-\theta_2)}}{\left(\dfrac{m_{1r}}{12} + m_{1b}+m_{2b}\right)}\left[m_{2b}l_2\dot{\theta}_2^2\sin{(\theta_1-\theta_2)} + g \cos{\theta_1}\left(\dfrac{m_{1r}}{2} +m_{2r} +m_{1b} + m_{2b}\right) -(b_{1b} + c_{1b} l_1 \dot{\theta}_1)l_1 \dot{\theta}_1 \right.\right.\\ &\quad\left.\left.+\left(b_{2b}+c_{2b}\sqrt{l_1^2 \dot{\theta}_1^2 + l_2^2 \dot{\theta}_2^2 +2l_1 l_2\dot{\theta}_1 \dot{\theta}_2 \cos{(\theta_1-\theta_2)}}\right)(l_1 \dot{\theta}_1 + l_2 \dot{\theta}_2 \cos{(\theta_1-\theta_2)})+\left(b_{1r} + \dfrac{c_{1r}l_1 \dot{\theta}_1}{2}\right) \dfrac{l_1 \dot{\theta}_1}{4} +\left(b_{2r} + c_{2r}\sqrt{l_1^2 \dot{\theta}_1^2 + \dfrac{l_2^2 \dot{\theta}_2^2}{4} + l_2 \dot{\theta}_1 \dot{\theta}_2 \cos{(\theta_1 -\theta_2)}}\right)\left(l_1 \dot{\theta}_1 \right.\right.\right. \\ &\quad\left.\left.\left.+ \dfrac{l_2\dot{\theta}_2 \cos{\left(\theta_1 - \theta_2\right)}}{2}\right)\right] + m_{2b}(l_1\dot{\theta}_1^2\sin{(\theta_1-\theta_2)}-g\cos{\theta_2}) -\dfrac{1}{4}\left(b_{2r} + c_{2r}\sqrt{l_1^2 \dot{\theta}_1^2 + \dfrac{l_2^2 \dot{\theta}_2^2}{4} + l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos{(\theta_1 -\theta_2)}}\right)(2l_1 \dot{\theta}_1 \cos{(\theta_1-\theta_2)}+ l_2 \dot{\theta}_2)\right.\\ &\quad\left.-\left(b_{2b}+c_{2b}\sqrt{l_1^2 \dot{\theta}_1^2 + l_2^2 \dot{\theta}_2^2 +2l_1l_2\dot{\theta}_1 \dot{\theta}_2 \cos{(\theta_1-\theta_2)}}\right)(l_1 \dot{\theta}_1 \cos{(\theta_1-\theta_2)} + l_2 \dot{\theta}_2)\right]. \end{aligned}
Simulation parameter form.
Parameter Value Explanation
Gravitational acceleration (in ms2\mathrm{m}\cdot \mathrm{s}^{-2}).
Length of pendulum 1 in metres.
Length of pendulum 2 in metres.
Mass of pendulum bob 1 in kilograms.
Mass of pendulum rod 1 in kilograms.
Mass of pendulum bob 2 in kilograms.
Mass of pendulum rod 2 in kilograms.
Linear dissipation coefficient of pendulum bob 1.
Linear dissipation coefficient of pendulum rod 1.
Quadratic dissipation coefficient of pendulum bob 1.
Quadratic dissipation coefficient of pendulum rod 1.
Linear dissipation coefficient for pendulum bob 2.
Linear dissipation coefficient for pendulum rod 2.
Quadratic dissipation coefficient for pendulum bob 2.
Quadratic dissipation coefficient for pendulum rod 2.
End time for the simulation in seconds.
Initial value of θ1\theta_1 in radians.
Initial value of θ˙1\dot{\theta}_1 in radians per second.
Initial value of θ2\theta_2 in radians.
Initial value of θ˙2\dot{\theta}_2 in radians per second.
Error tolerance.
Tolerance type, can be either absolute (0) or relative (1).
Initial step size.
Minimum allowed step size.
Time increment for skipping ahead in animation.
Time you want to skip ahead to in animation when you press the skip button.
Width (in px) of Plotly windows used for plotting and animation below.
Height (in px) of Plotly windows used for plotting and animation below.
Proportion of animation time passed per real time. tScale=1.0t_{\mathrm{Scale}}=1.0 means animation and real time match. tScale<1.0t_{\mathrm{Scale}}<1.0 means the animation is going more slowly than real time. tScale>1.0t_{\mathrm{Scale}}>1.0 means it is going more rapidly.