Single pendulum

This webpage uses the Runge-Kutta-Fehlberg fourth-order method with fifth-order error checking (RKF45) to approximate the solution to the problem of a single pendulum. To calculate the equations of motion, we must first determine the Lagrangian which is defined as the kinetic energy minus the potential energy. The kinetic energy of the bob is simply its mass times its velocity squared over two. The kinetic energy of the rod has a term like this, too, (although the velocity we will focus on will be the velocity of its centre of mass) but we must also consider its rotational velocity. As for the potential energy of the bob and rod, it will simply be the mass of them times by the yy-coordinate of its centre of mass.

L=mb2vb2+12Icmω2+mr2vr,cm2mbglsinθmrglsinθ2.\begin{aligned} \mathcal{L} &= \dfrac{m_b}{2} |\vec{v}_b|^2 + \dfrac{1}{2} I_{\mathrm{cm}}\omega^2 + \dfrac{m_r}{2} |\vec{v}_{r, \mathrm{cm}}|^2 - m_bg l\sin{\theta} -\dfrac{m_r gl\sin{\theta}}{2}. \end{aligned}

Where gg is the acceleration due to gravity in metres per second squared, ll is the length of the pendulum in metres, θ\theta is the angle from the positive xx axis (in radians). mbm_b is the mass of the bob in kilograms. vb\vec{v}_b is the velocity vector for the bob in metres per second. vr,cm\vec{v}_{r,\mathrm{cm}} is the velocity of the rod's centre of mass in metres per second. mrm_r is the mass of the rod in kilograms. tt is the time in seconds. IcmI_{\mathrm{cm}} is the moment of inertia around the centre of mass of the rod. It would be mrl212\dfrac{m_r l^2}{12}.[1] Velocity squared, for the bob, in polar coordinates is given by l˙2+l2θ˙2\dot{l}^2 + l^2\dot{\theta}^2 but here ll is a constant and hence its time derivative is zero. As for the velocity of the rod's centre of mass, we will assume the rod is uniform in mass and hence its centre of mass is in its centre, which we will assume is midway between the bob and the origin. Hence its velocity will be half that of the bob.

L=mb2l2θ˙2+mr24l2θ˙2+mr8l2θ˙2(mb+mr2)glsinθ=[mb2+mr6]l2θ˙2(mb+mr2)glsinθ.\begin{aligned} \mathcal{L} &= \dfrac{m_b}{2} l^2\dot{\theta}^2 + \dfrac{m_r}{24}l^2\dot{\theta}^2 + \dfrac{m_r}{8}l^2\dot{\theta}^2 - \left(m_b+\dfrac{m_r}{2}\right)gl\sin{\theta} \\ &= \left[\dfrac{m_b}{2} + \dfrac{m_r}{6}\right]l^2\dot{\theta}^2 - \left(m_b+\dfrac{m_r}{2}\right)gl\sin{\theta}. \end{aligned}

Hence its generalized momentum, momentum's time derivative and force are:

Lθ˙=[mb+mr3]l2θ˙ddtLθ˙=[mb+mr3]l2θ¨Lθ=(mb+mr2)glcosθ.\begin{aligned} \dfrac{\partial \mathcal{L}}{\partial \dot{\theta}} &= \left[m_b + \dfrac{m_r}{3}\right]l^2\dot{\theta}\\ \dfrac{d}{dt}\dfrac{\partial \mathcal{L}}{\partial \dot{\theta}} &= \left[m_b + \dfrac{m_r}{3}\right]l^2\ddot{\theta} \\ \dfrac{\partial \mathcal{L}}{\partial \theta} &= - \left(m_b+\dfrac{m_r}{2}\right)gl\cos{\theta}. \end{aligned}

Next we will define δfδqi\dfrac{\delta' f}{\delta' q_i} as the negative functional derivative of f(qi,q˙i,t)f(q_i, \dot{q}_i, t) with respect to qiq_i.

δLδθ=ddtLθ˙Lθ=[mb+mr3]l2θ¨+(mb+mr2)glcosθ.\begin{aligned} \dfrac{\delta' \mathcal{L}}{\delta' \theta} &= \dfrac{d}{dt}\dfrac{\partial \mathcal{L}}{\partial \dot{\theta}} - \dfrac{\partial \mathcal{L}}{\partial \theta} \\ &= \left[m_b + \dfrac{m_r}{3}\right]l^2\ddot{\theta} + \left(m_b+\dfrac{m_r}{2}\right)gl\cos{\theta}. \end{aligned}

Next we will calculate the generalized dissipation force QiQ_i which is defined as being the dot product of the dissipation force vector with generalized basis vector for each component of a physical system. Here there are two components of the system: the rod and the bob. The dissipation force vector will be assumed to be FD,j=(bj+cjvj)vj\vec{F}_{D,j} = -(b_j + c_j|\vec{v}_j|)\vec{v}_j.

Qθ=(br+crvr)vrrrθ(bb+cbvb)vbrbθ=(br+crlθ˙2)lθ˙2[sinθcosθ]l2[sinθcosθ](bb+cblθ˙)lθ˙[sinθcosθ]l[sinθcosθ]=(br+crlθ˙2)l2θ˙4(bb+cblθ˙)l2θ˙.\begin{aligned} Q_{\theta} &= -(b_r+c_r|\vec{v}_r|)\vec{v}_r \cdot \dfrac{\partial \vec{r}_r}{\partial \theta} -(b_b+c_b|\vec{v}_b|)\vec{v}_b \cdot \dfrac{\partial \vec{r}_b}{\partial \theta}\\ &= -\left(b_r+c_r\dfrac{l|\dot{\theta}|}{2}\right)\dfrac{l\dot{\theta}}{2}\begin{bmatrix} -\sin{\theta} \\ \cos{\theta} \end{bmatrix} \cdot \dfrac{l}{2} \begin{bmatrix} -\sin{\theta} \\ \cos{\theta} \end{bmatrix} -(b_{b} +c_{b} l|\dot{\theta}|)l\dot{\theta}\begin{bmatrix} -\sin{\theta} \\ \cos{\theta} \end{bmatrix} \cdot l \begin{bmatrix} -\sin{\theta} \\ \cos{\theta} \end{bmatrix} \\ &= -\left(b_r+c_r \dfrac{l|\dot{\theta}|}{2}\right)\dfrac{l^2\dot{\theta}}{4} -(b_b+c_b l|\dot{\theta}|)l^2\dot{\theta}. \end{aligned}

Hence the Euler-Lagrange equation with dissipative forces is

δLδθ=Qθ[mb+mr3]l2θ¨+(mb+mr2)glcosθ=(br+crlθ˙2)l2θ˙4(bb+cblθ˙)l2θ˙.\begin{aligned} \dfrac{\delta' \mathcal{L}}{\delta' \theta} &= Q_{\theta} \\ \left[m_b + \dfrac{m_r}{3}\right]l^2\ddot{\theta} + \left(m_b+\dfrac{m_r}{2}\right)gl\cos{\theta} &= -\left(b_r+c_r \dfrac{l|\dot{\theta}|}{2}\right)\dfrac{l^2\dot{\theta}}{4} -(b_b+c_b l|\dot{\theta}|)l^2\dot{\theta}. \end{aligned}

Dividing by the coefficient of θ¨\ddot{\theta}, namely [mb+mr3]l2\left[m_b + \dfrac{m_r}{3}\right]l^2 yields

θ¨+mb+mr2[mb+mr3]lgcosθ=(br+crlθ˙2)θ˙4+(bb+cblθ˙)θ˙mb+mr3.\begin{aligned} \ddot{\theta} + \dfrac{m_b+\dfrac{m_r}{2}}{\left[m_b + \dfrac{m_r}{3}\right]l}g\cos{\theta} &= -\dfrac{\left(b_r+c_r \dfrac{l|\dot{\theta}|}{2}\right)\dfrac{\dot{\theta}}{4} +(b_b+c_b l|\dot{\theta}|)\dot{\theta}}{m_b + \dfrac{m_r}{3}}. \end{aligned}

Isolating θ¨\ddot{\theta} yields

θ¨=mb+mr2[mb+mr3]lgcosθ(br+crlθ˙2)θ˙4+(bb+cblθ˙)θ˙mb+mr3=6mb+3mr[6mb+2mr]lgcosθ(6br+24bb)θ˙+(3cr+24cb)lθ˙θ˙24mb+8mr.\begin{aligned} \ddot{\theta} &= -\dfrac{m_b+\dfrac{m_r}{2}}{\left[m_b + \dfrac{m_r}{3}\right]l}g\cos{\theta} -\dfrac{\left(b_r+c_r \dfrac{l|\dot{\theta}|}{2}\right)\dfrac{\dot{\theta}}{4} +(b_b+c_b l|\dot{\theta}|)\dot{\theta}}{m_b + \dfrac{m_r}{3}} \\ &= -\dfrac{6m_b+3m_r}{\left[6m_b + 2m_r\right]l}g\cos{\theta} -\dfrac{(6b_r+24b_b)\dot{\theta}+(3c_r+24c_b)l|\dot{\theta}|\dot{\theta}}{24m_b + 8m_r}. \end{aligned}

Below you can specify the various parameters for the problem we will solve.

References

[1] Dourmashkin, P (2022). 16.3 Rotational Kinetic Energy and Moment of Inertia in Classical Mechanics (Dourmashkin). Massachusetts Institute of Technology.

Simulation parameter form.
Parameter Value Explanation
Acceleration due to gravity in metres (m) per second (s) squared (ms2\mathrm{m} \cdot \mathrm{s}^{-2}).
Length (m) of the pendulum rod.
Mass (kilograms or kg) of rod.
Mass (kg) of bob.
Linear dissipation coefficient of rod.
Linear dissipation coefficient of bob.
Quadratic dissipation coefficient of rod.
Quadratic dissipation coefficient of bob.
End time (s) for the simulation. The default tft_f value is equal to four times the period TT of the problem.
Initial angle (radians or r) from the positive xx-axis.
Initial first derivative of θ\theta with respect to tt (rs1\mathrm{r}\cdot \mathrm{s}^{-1}).
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.