Posts on Physics

A proof of Earnshaw's theorem


Posted by Diego Assencio on 2014.09.19 under Physics (Electromagnetism)

Consider a set of $n$ fixed point charges $q_1, q_2, \ldots, q_n$ in vacuum and let ${\bf E}({\bf x})$ be the electric field produced by these charges (see figure 1). Can we find a point ${\bf z}$ in space such that if a charge $q$ is placed there, it will remain in stable equilibrium?

Fig. 1: A set of $n$ fixed charges and a very small Gaussian surface $S$ centered at a point ${\bf z}$. For ${\bf z}$ to be a point of stable equilibrium (assuming we would place a positive charge $q$ there), ${\bf E}({\bf z})$ must be zero and ${\bf E}({\bf x})$ must point towards ${\bf z}$ at all points ${\bf x}$ of $S$.

The answer is no, and this fact is referred to as the Earnshaw's theorem. We will prove this assuming $q \gt 0$, but the proof is similar for $q \lt 0$.

Let's first think about what do we mean by stable equilibrium for a charge $q$ placed at ${\bf z}$:

1.the electric force $q{\bf E}({\bf z})$ acting on $q$ must be zero, i.e, ${\bf E}({\bf z}) = {\bf 0}$
2.for a small displacement $\delta{\bf z}$ of $q$ around ${\bf z}$, $q{\bf E}({\bf z} + \delta{\bf z})$ must point towards ${\bf z}$.

The second property means the electric force on $q$ must be restoring, meaning if we move $q$ a little bit away from ${\bf z}$, this force will will tend to bring $q$ back to ${\bf z}$. Since we are assuming $q \gt 0$, ${\bf E}({\bf z} + \delta{\bf z})$ itself must then point towards ${\bf z}$ for the electric force on $q$ to be restoring.

Consider now a very small Gaussian surface $S$ shaped like a sphere of radius $r$ and centered at ${\bf z}$ (as shown in figure 1). Since ${\bf E}({\bf x})$ points towards $q$ (i.e., inwards) for any point ${\bf x}$ on the surface of this sphere, we have that: $$ \int_{S} {\bf E}\cdot d{\bf A} \lt 0 \label{post_bc04395b103021d338b4e30a061bfc74_surface_int} $$ where $d{\bf A}$ is an outward-oriented infinitesimal surface area element. But using the divergence theorem and Gauss's law on equation \eqref{post_bc04395b103021d338b4e30a061bfc74_surface_int} implies we must have a negative charge $Q_{\textrm{enc}}$ enclosed inside our sphere: $$ \int_{S} {\bf E}\cdot d{\bf A} = \int_{V}\nabla\cdot{\bf E}\,dV = \displaystyle\frac{Q_{\textrm{enc}}}{\epsilon_0} \lt 0 $$ This goes against our initial assumption that ${\bf z}$ is a point in vacuum where no charge is placed ($q$ itself does not count since it does not contribute to the electric field ${\bf E}$).

One of the consequences of Earnshaw's theorem is the fact that a classical model for the atom cannot be stable if it is based on a static (fixed) distribution of charges. For such an atom, the charge configuration would not be in stable equilibrium and therefore unable to resist even the slightest disturbance.

Comments (5) Direct link

The double pendulum: Hamiltonian formulation


Posted by Diego Assencio on 2014.05.17 under Physics (Mechanics)

Consider the double pendulum shown on figure 1. A double pendulum is formed by attaching a pendulum directly to another one. Each pendulum consists of a bob connected to a massless rigid rod which is only allowed to move along a vertical plane. The pivot of the first pendulum is fixed to a point $O$. All motion is frictionless.

Fig. 1: A double pendulum.

As shown on a previous post, if we define $\theta_1$ and $\theta_2$ as the angles which the first and second rods make with the vertical direction respectively, the following expression can be obtained for the Lagrangian of the system: $$ \begin{eqnarray} L =\frac{1}{2}(m_1 + m_2) l_1^2 \dot{\theta}_1^2 &+& \frac{1}{2}m_2 l_2^2 \dot{\theta}_2^2 + m_2l_1l_2\dot{\theta}_1\dot{\theta}_2 \cos(\theta_1 - \theta_2)\nonumber\\[3pt] &+&(m_1 + m_2) g l_1 \cos\theta_1 + m_2 g l_2\cos\theta_2 \end{eqnarray} $$ From the Lagrangian, one can obtain the canonical momenta of the system: $$ \begin{eqnarray} \displaystyle p_{\theta_1} &=& \frac{\partial L}{\partial \dot{\theta}_1} = (m_1 + m_2) l_1^2 \dot{\theta}_1 + m_2 l_1 l_2 \dot{\theta}_2\cos(\theta_1-\theta_2) \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_ptheta1} \\[5pt] \displaystyle p_{\theta_2} &=& \frac{\partial L}{\partial \dot{\theta}_2} = m_2 l_2^2 \dot{\theta}_2 + m_2 l_1 l_2 \dot{\theta}_1\cos(\theta_1-\theta_2) \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_ptheta2} \end{eqnarray} $$ The Hamiltonian of the system is given by: $$ H = \sum_{i=1}^2 \dot{\theta}_i p_{\theta_i} - L \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian} $$ From the Hamiltonian $H$, we can obtain a set of equations of motion for the system which are equivalent to the Euler-Lagrange equations: $$ \displaystyle \dot{\theta}_i = \frac{\partial H}{\partial p_{\theta_i}}, \quad\quad \dot{p}_{\theta_i} = -\frac{\partial H}{\partial \theta_i} \quad \textrm{ for } \quad i = 1,2 \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eqs} $$ Equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eqs} requires us to write $H$ as a function of the variables $\theta_1$, $\theta_2$, $p_{\theta_1}$ and $p_{\theta_2}$, so from equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian} we see that we must determine $\dot{\theta}_i$ and $L$ in terms of these variables. With that in mind, notice that equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_ptheta1} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_ptheta2} can be written in matrix form as shown below: $$ \left(\begin{matrix} p_{\theta_1} \\[4pt] p_{\theta_1} \end{matrix}\right) = B \left(\begin{matrix} \dot{\theta}_1 \\ \dot{\theta}_2 \end{matrix}\right) \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_p_func_dottheta} $$ where $B$ is a $2 \times 2$ matrix whose entries depend on $\theta_1$ and $\theta_2$: $$ B = \left(\begin{matrix} (m_1 + m_2) l_1^2 & m_2 l_1 l_2 \cos(\theta_1-\theta_2) \\ m_2 l_1 l_2 \cos(\theta_1-\theta_2) & m_2 l_2^2 \end{matrix}\right) $$ From equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_p_func_dottheta}, we can obtain the generalized velocities $\dot{\theta}_i$ in terms of of the canonical momenta $p_{\theta_i}$ and the angles $\theta_i$: $$ \left(\begin{matrix} \dot{\theta}_1 \\ \dot{\theta}_2 \end{matrix}\right) = B^{-1} \left(\begin{matrix} p_{\theta_1} \\[4pt] p_{\theta_1} \end{matrix}\right) \label{%INDEX_eq_dotheta_func_p} $$ The matrix $B$ is indeed invertible for all values of $\theta_1$ and $\theta_2$ since: $$ \begin{eqnarray} \det(B) &=& m_1m_2l_1^2 l_2^2 + m_2^2l_1^2 l_2^2\left[1 - \cos^2(\theta_1 - \theta_2)\right] \nonumber\\[5pt] &=& m_1m_2l_1^2 l_2^2 + m_2^2l_1^2 l_2^2\sin^2(\theta_1 - \theta_2) \nonumber\\[5pt] &\geq& m_1m_2l_1^2 l_2^2 \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_det_B} \end{eqnarray} $$ Being a $2 \times 2$ matrix, $B$ can be inverted directly: $$ B^{-1} = \displaystyle\frac{1}{\det(B)} \left( \begin{matrix} m_2 l_2^2 & -m_2 l_1 l_2 \cos(\theta_1-\theta_2) \\ -m_2 l_1 l_2 \cos(\theta_1-\theta_2) & (m_1 + m_2) l_1^2 \end{matrix} \right) \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_inv_B} $$ Using equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_inv_B} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_det_B} on equation \eqref{%INDEX_eq_dotheta_func_p}, we get (after canceling out common factors and rearranging some terms): $$ \begin{eqnarray} \displaystyle \dot{\theta}_1 &=& \frac{l_2 p_{\theta_1} - l_1 p_{\theta_2}\cos(\theta_1-\theta_2)} {l_1^2l_2\left[m_1 + m_2 \sin^2(\theta_1 - \theta_2)\right]} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta1_func_p}\\[5pt] \displaystyle \dot{\theta}_2 &=& \frac{-m_2 l_2 p_{\theta_1}\cos(\theta_1-\theta_2) + (m_1 + m_2)l_1 p_{\theta_2}} {m_2l_1 l_2^2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta2_func_p} \end{eqnarray} $$ Using equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta1_func_p} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta2_func_p} on equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian} yields the Hamiltonian $H$ in terms of $\theta_1$, $\theta_2$, $p_{\theta_1}$ and $p_{\theta_2}$ (after quite a bit of algebraic work): $$ \begin{eqnarray} H &=& \frac{m_2 l_2^2 p^2_{\theta_1} + (m_1 + m_2)l_1^2 p^2_{\theta_2} - 2m_2 l_1 l_2 p_{\theta_1}p_{\theta_2}\cos(\theta_1-\theta_2)} {2m_2 l_1^2l_2^2 \left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}\nonumber\\[3pt] & & \quad\quad- (m_1 + m_2) g l_1 \cos\theta_1 - m_2 g l_2\cos\theta_2 \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_right_vars} \end{eqnarray} $$ Equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_right_vars} can finally be used on equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eqs} to give us the Hamiltonian equations of motion for the double pendulum: $$ \begin{eqnarray} \displaystyle \dot{\theta}_1 &=& \;\;\frac{\partial H}{\partial p_{\theta_1}} &=& \frac{l_2 p_{\theta_1} - l_1 p_{\theta_2}\cos(\theta_1-\theta_2)} {l_1^2l_2\left[m_1 + m_2 \sin^2(\theta_1 - \theta_2)\right]} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_1} \\[5pt] \displaystyle \dot{\theta}_2 &=& \;\;\frac{\partial H}{\partial p_{\theta_2}} &=& \frac{-m_2 l_2 p_{\theta_1}\cos(\theta_1-\theta_2) + (m_1 + m_2)l_1 p_{\theta_2}} {m_2l_1 l_2^2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_2} \\[5pt] \dot{p}_{\theta_1} &=& -\frac{\partial H}{\partial \theta_1} &=& -(m_1 + m_2)gl_1 \sin\theta_1 - h_1 + h_2\sin\left[2(\theta_1 - \theta_2)\right] \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_3} \\[5pt] \dot{p}_{\theta_2} &=& -\frac{\partial H}{\partial \theta_2} &=& -m_2gl_2 \sin\theta_2 + h_1 - h_2\sin\left[2(\theta_1 - \theta_2)\right] \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_4} \end{eqnarray} $$ where: $$ h_1 = \frac{p_{\theta_1}p_{\theta_2}\sin(\theta_1 - \theta_2)} {l_1l_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]} \\[5pt] \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_h1} $$ $$ h_2 = \frac{m_2 l_2^2 p^2_{\theta_1} + (m_1 + m_2)l_1^2 p^2_{\theta_2} - 2 m_2 l_1l_2 p_{\theta_1}p_{\theta_2}\cos(\theta_1 - \theta_2)} {2l_1^2l_2^2 \left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]^2} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_h2} $$ Notice that equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_1} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_2} are the same as equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta1_func_p} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta2_func_p}, so in the process of determining the Hamiltonian in terms of the canonical momenta $p_{\theta_i}$ and the angles $\theta_i$, we ended up already obtaining two of the Hamiltonian equations.

Equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_1}-\eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_4} form a set of coupled first order differential equations on the variables $\theta_1$, $\theta_2$, $p_{\theta_1}$ and $p_{\theta_2}$. This set of equations can be solved numerically using a Runge-Kutta (RK) method. A simulator based on the fourth-order RK method can be found here.

References

[1]math24.net
[2]physics.buffalo.edu
Comments (1) Direct link

The double pendulum: Lagrangian formulation


Posted by Diego Assencio on 2014.02.28 under Physics (Mechanics)

Consider the double pendulum shown on figure 1. A double pendulum is formed by attaching a pendulum directly to another one. Each pendulum consists of a bob connected to a massless rigid rod which is only allowed to move along a vertical plane. The pivot of the first pendulum is fixed to a point $O$. All motion is frictionless.

Fig. 1: A double pendulum.

In our discussion, the fixed point $O$ will be taken as the origin of the Cartesian coordinate system with the $x$ axis pointing along the horizontal direction and the $y$ axis pointing vertically upwards. Let $\theta_1$ and $\theta_2$ be the angles which the first and second rods make with the vertical direction respectively. As can be seen on figure 1, the positions of the bobs are given by: $$ \begin{eqnarray} x_1 &=& l_1 \sin\theta_1 \quad\quad & y_1 &=& -l_1 \cos\theta_1\\[5pt] x_2 &=& l_1 \sin\theta_1 + l_2 \sin\theta_2 \quad\quad & y_2 &=& -l_1\cos\theta_1 -l_2\cos\theta_2 \end{eqnarray} $$ Differentiating the quantities above with respect to time, we obtain the velocities of the bobs: $$ \begin{eqnarray} \dot{x}_1 &=& l_1 \dot{\theta}_1\cos\theta_1 \quad\quad & \dot{y}_1 &=& l_1 \dot{\theta}_1\sin\theta_1\\[5pt] \dot{x}_2 &=& l_1 \dot{\theta}_1\cos\theta_1 + l_2 \dot{\theta}_2\cos\theta_2 \quad\quad & \dot{y}_2 &=& l_1 \dot{\theta}_1\sin\theta_1 + l_2\dot{\theta}_2\sin\theta_2 \end{eqnarray} $$

The double pendulum is a very interesting system as it is very simple but can show chaotic behavior for certain initial conditions. In this regime, slightly changing the initial values of the angles ($\theta_1,\theta_2$) and angular velocities ($\dot{\theta}_1,\dot{\theta}_2$) makes the trajectories of the bobs become very different from the original ones.

The Lagrangian for the double pendulum is given by $L = T - V$, where $T$ and $V$ are the kinetic and potential energies of the system respectively. The kinetic energy $T$ is given by: $$ \begin{eqnarray} T &=& \displaystyle\frac{1}{2}m_1v_1^2 + \frac{1}{2}m_2v_2^2 \nonumber \\[5pt] &=& \frac{1}{2}m_1(\dot{x}_1^2 + \dot{y}_1^2) + \frac{1}{2}m_2(\dot{x}_2^2 + \dot{y}_2^2) \nonumber \\[5pt] &=& \frac{1}{2}m_1 l_1^2 \dot{\theta}_1^2 + \frac{1}{2}m_2\left[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] \end{eqnarray} $$ where above we used the fact that $\cos\theta_1\cos\theta_2 + \sin\theta_1\sin\theta_2 =$ $\cos(\theta_1-\theta_2)$. The potential energy $V$ is given by: $$ \begin{eqnarray} V &=& m_1 g y_1 + m_2gy_2 \nonumber \\[5pt] &=& -m_1 g l_1 \cos\theta_1 - m_2 g (l_1 \cos\theta_1 + l_2 \cos\theta_2) \nonumber\\[5pt] &=& -(m_1 + m_2) g l_1 \cos\theta_1 - m_2 g l_2\cos\theta_2 \end{eqnarray} $$ The Lagrangian of the system is then: $$ \begin{eqnarray} L =\frac{1}{2}(m_1 + m_2) l_1^2 \dot{\theta}_1^2 &+& \frac{1}{2}m_2 l_2^2 \dot{\theta}_2^2 + m_2l_1l_2\dot{\theta}_1\dot{\theta}_2 \cos(\theta_1 - \theta_2)\nonumber\\[3pt] &+&(m_1 + m_2) g l_1 \cos\theta_1 + m_2 g l_2\cos\theta_2 \end{eqnarray} $$ The canonical momenta associated with the coordinates $\theta_1$ and $\theta_2$ can be obtained directly from $L$: $$ \begin{eqnarray} \displaystyle p_{\theta_1} &=& \frac{\partial L}{\partial \dot{\theta}_1} = (m_1 + m_2) l_1^2 \dot{\theta}_1 + m_2 l_1 l_2 \dot{\theta}_2\cos(\theta_1-\theta_2) \\[5pt] \displaystyle p_{\theta_2} &=& \frac{\partial L}{\partial \dot{\theta}_2} = m_2 l_2^2 \dot{\theta}_2 + m_2 l_1 l_2 \dot{\theta}_1\cos(\theta_1-\theta_2) \end{eqnarray} $$ The equations of motion of the system are the Euler-Lagrange equations: $$ \displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}_i}\right) - \frac{\partial L}{\partial \theta_i} = 0 \Longrightarrow \frac{d p_{\theta_i}}{dt} - \frac{\partial L}{\partial \theta_i} = 0 \quad \textrm{ for } \quad i = 1,2 \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq} $$ Since: $$ \begin{eqnarray} \displaystyle \frac{d p_{\theta_1}}{dt} &=& (m_1 + m_2) l_1^2 \ddot{\theta}_1 + m_2 l_1 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2) \nonumber \\ &\;&\quad -\, m_2 l_1 l_2 \dot{\theta}_2\dot{\theta}_1\sin(\theta_1 - \theta_2) + m_2 l_1 l_2 \dot{\theta}_2^2\sin(\theta_1 - \theta_2) \\[5pt] \displaystyle \frac{d p_{\theta_2}}{dt} &=& m_2 l_2^2 \ddot{\theta}_2 + m_2 l_1 l_2 \ddot{\theta}_1 \cos(\theta_1 - \theta_2) \nonumber \\ &\;&\quad -\, m_2 l_1 l_2 \dot{\theta}_1^2 \sin(\theta_1 - \theta_2) + m_2 l_1 l_2 \dot{\theta}_1\dot{\theta}_2 \sin(\theta_1 - \theta_2)\\[5pt] \displaystyle \frac{\partial L}{\partial \theta_1} &=& - m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1-\theta_2) - (m_1 + m_2) g l_1 \sin\theta_1 \\[5pt] \displaystyle \frac{\partial L}{\partial \theta_2} &=& m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1-\theta_2) - m_2 g l_2 \sin\theta_2 \end{eqnarray} $$ then equation \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq} yields (after dividing by $l_1$ when $i = 1$ and by $m_2 l_2$ when $i = 2$): $$ \begin{eqnarray} (m_1 + m_2) l_1 \ddot{\theta}_1 &+& m_2 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2) \nonumber\\ &+& m_2 l_2 \dot{\theta}_2^2\sin(\theta_1 - \theta_2) + (m_1 + m_2) g \sin\theta_1 = 0 \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1} \end{eqnarray} $$ $$ l_2 \ddot{\theta}_2 + l_1 \ddot{\theta}_1 \cos(\theta_1 - \theta_2) - l_1 \dot{\theta}_1^2 \sin(\theta_1 - \theta_2) + g \sin\theta_2 = 0 \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2} $$ Equations \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1} and \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2} form a system of coupled second-order nonlinear differential equations. Dividing equation \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1} by $(m_1 + m_2)l_1$ and equation \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2} by $l_2$ and also moving all terms which do not involve $\ddot{\theta}_1$ and $\ddot{\theta}_2$ to the right-hand side, we obtain: $$ \begin{eqnarray} \ddot{\theta}_1 + \alpha_1(\theta_1,\theta_2) \ddot{\theta}_2 &=& f_1(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1_simple}\\[5pt] \ddot{\theta}_2 + \alpha_2(\theta_1,\theta_2) \ddot{\theta}_1 &=& f_2(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2_simple} \end{eqnarray} $$ where: $$ \begin{eqnarray} \alpha_1(\theta_1,\theta_2) &:=& \displaystyle\frac{l_2}{l_1}\left(\frac{m_2}{m_1 + m_2}\right)\cos(\theta_1 - \theta_2) \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_alpha1}\\[5pt] \alpha_2(\theta_1,\theta_2) &:=& \frac{l_1}{l_2}\cos(\theta_1-\theta_2) \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_alpha2} \end{eqnarray} $$ and: $$ \begin{eqnarray} \displaystyle f_1(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) &:=& -\frac{l_2}{l_1}\left(\frac{m_2}{m_1+m_2}\right) \dot{\theta}_2^2\sin(\theta_1 - \theta_2) - \frac{g}{l_1} \sin\theta_1 \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_f_1}\\[5pt] \displaystyle f_2(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) &:=& \frac{l_1}{l_2}\dot{\theta}_1^2\sin(\theta_1-\theta_2) - \frac{g}{l_2} \sin\theta_2 \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_f_2} \end{eqnarray} $$ Interestingly, $f_1$ does not depent on $\dot{\theta}_1$ and $f_2$ does not depend on $\dot{\theta}_2$. Equations \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1_simple} and \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2_simple} can be combined into a single equation: $$ A \left( \begin{matrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{matrix} \right) = \left( \begin{matrix} 1 & \alpha_1 \\[3pt] \alpha_2 & 1 \end{matrix} \right) \left( \begin{matrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{matrix} \right) = \left( \begin{matrix} f_1 \\ f_2 \end{matrix} \right) \label{%INDEX_eq_lin_sys_ddtheta} $$ where the matrix $A$ depends on $\theta_1$ and $\theta_2$ since $\alpha_1$ and $\alpha_2$ depend on these variables. Being a $2 \times 2$ matrix, $A$ can be inverted directly: $$ A^{-1} = \displaystyle\frac{1}{\det(A)} \left( \begin{matrix} 1 & -\alpha_1 \\ -\alpha_2 & 1 \end{matrix} \right) = \displaystyle\frac{1}{1 - \alpha_1\alpha_2} \left( \begin{matrix} 1 & -\alpha_1 \\ -\alpha_2 & 1 \end{matrix} \right) \label{post_1500c66ae7ab27bb0106467c68feebc6_inv_A} $$ Before we continue, notice that $A$ is always invertible since: $$ \det(A) = 1 - \alpha_1\alpha_2 = 1 - \left(\frac{m_2}{m_1 + m_2}\right)\cos^2(\theta_1 - \theta_2) \gt 0 $$ because $m_2 / (m_1 + m_2) \lt 1$ and $\cos^2(x) \leq 1$ for all real values of $x$. From equations \eqref{%INDEX_eq_lin_sys_ddtheta} and \eqref{post_1500c66ae7ab27bb0106467c68feebc6_inv_A} we obtain: $$ \left( \begin{matrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{matrix} \right) = A^{-1} \left( \begin{matrix} f_1 \\[3pt] f_2 \end{matrix} \right) = \displaystyle\frac{1}{1 - \alpha_1\alpha_2} \left( \begin{matrix} f_1 - \alpha_1 f_2\\[3pt] -\alpha_2 f_1 + f_2 \end{matrix} \right) $$ Finally, letting $\omega_1 := \dot{\theta}_1$ and $\omega_2 := \dot{\theta}_2$, we can write the equations of motion of the double pendulum as a system of coupled first order differential equations on the variables $\theta_1$, $\theta_2$, $\omega_1$, $\omega_2$: $$ \displaystyle\frac{d}{dt} \left( \begin{matrix} \theta_1 \\[1pt] \theta_2 \\[1pt] \omega_1 \\[1pt] \omega_1 \end{matrix} \right) = \left( \begin{matrix} \omega_1 \\ \omega_2 \\ g_1(\theta_1,\theta_2,\omega_1,\omega_2) \\ g_2(\theta_1,\theta_2,\omega_1,\omega_2) \end{matrix} \right) \label{post_1500c66ae7ab27bb0106467c68feebc6_first_order_eq_theta_omega} $$ $$ g_1 := \displaystyle\frac{f_1 - \alpha_1 f_2}{1 - \alpha_1\alpha_2} \quad\quad g_2 := \displaystyle\frac{-\alpha_2 f_1 + f_2}{1 - \alpha_1\alpha_2} \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_g1_g2} $$ where $\alpha_i = \alpha_i(\theta_1,\theta_2)$ and $f_i = f_i(\theta_1,\theta_2,\omega_1,\omega_2)$ for $i = 1,2$ are given on equations \eqref{post_1500c66ae7ab27bb0106467c68feebc6_eq_alpha1}-\eqref{post_1500c66ae7ab27bb0106467c68feebc6_eq_f_2}.

Equation \eqref{post_1500c66ae7ab27bb0106467c68feebc6_first_order_eq_theta_omega} can be solved numerically using a Runge-Kutta (RK) method. A simulator based on the fourth-order RK method can be found here.

References

[1]myphysicslab.com
[2]math24.net
[3]scienceworld.wolfram.com
[4]sophia.dtp.fmph.uniba.sk
Comments (5) Direct link