Pendulum motion
Considering a single pendulum:
A few assumptions:
- the rod has no mass with length L;
- the mass (represented by the ball of finite radius) is {m};
- all pivots are assumed to be without friction.
The an xy-coordinate system is set, with the origin of which coincides with suspension point of the rod. The coordinates of the pendulums are then in this coordinates system:
\begin{aligned} {x} & = L\sin\left(\theta\right) \\ {y} & = - L\cos\left(\theta\right) \end{aligned}
The kinetic and potential energy of the pendulums (respectively T and V) are expressed by the formulas:
\begin{aligned} T & = \frac{mv^2}{2} \\ & = \frac{m\left( \dot x^2 + \dot y^2 \right)}{2}\\ V & = m\,g\,y \end{aligned}
Then the Lagrangian can be written as:
\mathcal{L} = T - V = \frac{m}{2}\left( \dot x^2 + \dot y^2 \right) - m\,g\,y
Take into account that:
\begin{aligned} \dot x & = L\cos\left(\theta\right) \cdot \dot \theta \\ \dot y & = L\sin\left(\theta\right) \cdot \dot \theta \end{aligned}
Consequently:
\begin{aligned} {T} & = \frac{m}{2}\left( {\dot x^2 + \dot y^2} \right) = \frac{m}{2}\left( {L^2\dot \theta^2\cos^2\left(\theta\right) + L^2\dot \theta^2\sin^2\left(\theta\right)} \right) = \frac{m}{2}L^2\dot \theta^2 \\ {V} & = {m}g{y} = -{m}g{L}\cos\left(\theta\right) \end{aligned}
As a result, the Lagrangian of the system takes the following form:
\mathcal{L} = T - V = \frac{m}{2} L^2\dot \theta^2 + m\,g\,L\cos\left(\theta\right)
Now we can write the Euler-Lagrange equations:
\frac{\mathrm d}{\mathrm dt}\frac{\partial \mathcal{L}}{\partial \dot \theta } - \frac{\partial \mathcal{L}}{\partial \theta} = 0
The partial derivatives in these equations are expressed by the following formulas:
\begin{aligned} \frac{\mathrm d}{\mathrm dt}\frac{\partial L\mathcal{L}}{ \partial \dot \theta} & = ml^2 \ddot \theta \\[10px] \frac{\partial \mathcal{L}}{\partial \theta} & = - m g l \sin\left(\theta\right) \end{aligned}
The Euler-Lagrange becomes:
\begin{aligned} & mL^2 \ddot \theta = -m \,g \,L \sin\left(\theta\right) \\ & \ddot \theta = -\omega_0^2 \sin\left(\theta\right), \quad \omega = \sqrt{\frac{g}{L}} \end{aligned}
Defining \omega = \dot \theta, it is possible to transform into a system of two first order Ordinal Differential Equations (ODE):
\begin{aligned} & \dot \theta = \omega \\ & \dot \omega = \omega_0^2 \sin\left(\theta\right) \end{aligned}
Integrating this system with Runge-Kutta, it is possible to display the angular position \theta, the angular velocity \dot \theta and show that the energy that should be constant (as the system is ideal) is not constant (so the integrator is introducing a numerical error). Depending from the initial conditions, this numerical error it could be diverging to \pm \infty as time passes.
