Double Pendulum

A double pendulum is undoubtedly an actual miracle of nature. The jump in complexity, which is observed at the transition from a simple pendulum to a double pendulum is amazing. The oscillations of a simple pendulum are regular. For small deviations from equilibrium, these oscillations are harmonic and can be described by sine or cosine function. In the case of nonlinear oscillations, the period depends on the amplitude, but the regularity of the motion holds. In other words, in the case of a simple pendulum, the approximation of small oscillations fully reflects the essential properties of the system.

Double pendulum "behaves" quite differently. In the regime of small oscillations, the double pendulum demonstrates the phenomenon of beats. The character of oscillations of the pendulums changes radically with increasing energy − the oscillations become chaotic. Despite the fact that the double pendulum can be described by a system of several ordinary differential equations, that is by a completely deterministic model, the appearance of chaos looks very unusual. This situation is reminiscent of the Lorenz system where a deterministic model of three equations also shows chaotic behavior. Try to experiment with the application below and watch the movement of the double pendulum at different mass ratios and initial angles.

Animation is available on desktop screens.

Next, we will build a mathematical model of the double pendulum in the form of a system of nonlinear differential equations. Let's start with the derivation of the Lagrange equations.

Lagrange Equations

In Lagrangian mechanics, evolution of a system is described in terms of the generalized coordinates and generalized velocities. In our case, the deflection angles of the pendulums \({\alpha _1},{\alpha _2}\) and the angular velocities \({\dot \alpha _1},{\dot \alpha _2}\) can be taken as the generalized variables. Using these variables, we construct the Lagrangian for the double pendulum and write the Lagrange differential equations.

A simplified model of the double pendulum is shown in Figure \(1.\)

We assume that the rods are massless. Their lengths are \({l_1}\) and \({l_2}.\) The point masses (they are represented by the balls of finite radius) are \({m_1}\) and \({m_2}.\) All pivots are assumed to be frictionless.

We introduce the \(xy\)-coordinate system, the origin \(O\) of which coincides with suspension point of the upper pendulum. The coordinates of the pendulums are defined by the following relationships:

\[{x_1} = {l_1}\sin {\alpha _1},\;\; {x_2} = {l_1}\sin {\alpha _1} + {l_2}\sin {\alpha _2},\;\; {y_1} = - {l_1}\cos{\alpha _1},\;\; {y_2} = - {l_1}\cos{\alpha _1} - {l_2}\cos{\alpha _2}.\]

The kinetic and potential energy of the pendulums (respectively \(T\) and \(V\)) are expressed by the formulas

\[T = \frac{{{m_1}v_1^2}}{2} + \frac{{{m_2}v_2^2}}{2} = \frac{{{m_1}\left( {\dot x_1^2 + \dot y_1^2} \right)}}{2} + \frac{{{m_2}\left( {\dot x_2^2 + \dot y_2^2} \right)}}{2},\;\; V = {m_1}g{y_1} + {m_2}g{y_2}.\]

Then the Lagrangian can be written as

\[L = T - V = {T_1} + {T_2} - \left( {{V_1} + {V_2}} \right) = \frac{{{m_1}}}{2}\left( {\dot x_1^2 + \dot y_1^2} \right) + \frac{{{m_2}}}{2}\left( {\dot x_2^2 + \dot y_2^2} \right) - {m_1}g{y_1} - {m_2}g{y_2}.\]

Take into account that

\[{{\dot x}_1} = {l_1}\cos{\alpha _1} \cdot {{\dot \alpha }_1},\;\; {{\dot x}_2} = {l_1}\cos{\alpha _1} \cdot {{\dot \alpha }_1} + {l_2}\cos{\alpha _2} \cdot {{\dot \alpha }_2},\]

\[{{\dot y}_1} = {l_1}\sin{\alpha _1} \cdot {{\dot \alpha }_1},\;\; {{\dot y}_2} = {l_1}\sin{\alpha _1} \cdot {{\dot \alpha }_1} + {l_2}\sin{\alpha _2} \cdot {{\dot \alpha }_2}.\]

Consequently,

\[{T_1} = \frac{{{m_1}}}{2}\left( {\dot x_1^2 + \dot y_1^2} \right) = \frac{{{m_1}}}{2}\left( {l_1^2\dot \alpha _1^2{{\cos }^2}{\alpha _1} + l_1^2\dot \alpha _1^2{\sin^2}{\alpha _1}} \right) = \frac{{{m_1}}}{2}l_1^2\dot \alpha _1^2,\]

\[ {T_2} = \frac{{{m_2}}}{2}\left( {\dot x_2^2 + \dot y_2^2} \right) = \frac{{{m_2}}}{2}\left[ {{{\left( {{l_1}{{\dot \alpha }_1}\cos {\alpha _1} + {l_2}{{\dot \alpha }_2}\cos {\alpha _2}} \right)}^2} + {{\left( {{l_1}{{\dot \alpha }_1}\sin{\alpha _1} + {l_2}{{\dot \alpha }_2}\sin{\alpha _2}} \right)}^2}} \right] = \frac{{{m_2}}}{2}\left[ {l_1^2\dot \alpha _1^2{{\cos }^2}{\alpha _1} + l_2^2\dot \alpha _2^2{{\cos }^2}{\alpha _2}} + 2{l_1}{l_2}{\dot \alpha _1}{\dot \alpha _2}\cos {\alpha _1}\cos {\alpha _2} + l_1^2\dot \alpha _1^2{\sin^2}{\alpha _1} + {l_2^2\dot \alpha _2^2{\sin^2}{\alpha _2} + 2{l_1}{l_2}{{\dot \alpha }_1}{{\dot \alpha }_2}\sin{\alpha _1}\sin{\alpha _2}} \right] = \frac{{{m_2}}}{2}\left[ {l_1^2\dot \alpha _1^2 + l_2^2\dot \alpha _2^2 + 2{l_1}{l_2}{{\dot \alpha }_1}{{\dot \alpha }_2}\cos \left( {{\alpha _1} - {\alpha _2}} \right)} \right],\]

\[{V_1} = {m_1}g{y_1} = - {m_1}g{l_1}\cos {\alpha _1},\]

\[{V_2} = {m_2}g{y_2} = - {m_2}g\left( {{l_1}\cos {\alpha _1} + {l_2}\cos {\alpha _2}} \right).\]

As a result, the Lagrangian of the system takes the following form:

\[ L = T - V = {T_1} + {T_2} - \left( {{V_1} + {V_2}} \right) = \left( {\frac{{{m_1}}}{2} + \frac{{{m_2}}}{2}} \right)l_1^2\dot \alpha _1^2 + \frac{{{m_2}}}{2}l_2^2\dot \alpha _2^2 + {m_2}{l_1}{l_2}{\dot \alpha _1}{\dot \alpha _2}\cos \left( {{\alpha _1} - {\alpha _2}} \right) + \left( {{m_1} + {m_2}} \right)g{l_1}\cos {\alpha _1} + {m_2}g{l_2}\cos {\alpha _2}.\]

Now we can write the Lagrange equations (sometimes they are called the Euler-Lagrange equations):

\[\frac{d}{{dt}}\frac{{\partial L}}{{\partial {{\dot \alpha }_i}}} - \frac{{\partial L}}{{\partial {\alpha _i}}} = 0,\;\; i = 1,2.\]

The partial derivatives in these equations are expressed by the following formulas:

\[\frac{{\partial L}}{{\partial {{\dot \alpha }_1}}} = \left( {{m_1} + {m_2}} \right)l_1^2{{\dot \alpha }_1} + {m_2}{l_1}{l_2}{{\dot \alpha }_2}\cos \left( {{\alpha _1} - {\alpha _2}} \right),\]

\[\frac{{\partial L}}{{\partial {\alpha _1}}} = - {m_2}{l_1}{l_2}{{\dot \alpha }_1}{{\dot \alpha }_2}\sin\left( {{\alpha _1} - {\alpha _2}} \right) - \left( {{m_1} + {m_2}} \right)g{l_1}\sin{\alpha _1},\]

\[\frac{{\partial L}}{{\partial {{\dot \alpha }_2}}} = {m_2}l_2^2{{\dot \alpha }_2} + {m_2}{l_1}{l_2}{{\dot \alpha }_1}\cos \left( {{\alpha _1} - {\alpha _2}} \right),\]

\[\frac{{\partial L}}{{\partial {\alpha _2}}} = {m_2}{l_1}{l_2}{{\dot \alpha }_1}{{\dot \alpha }_2}\sin\left( {{\alpha _1} - {\alpha _2}} \right) - {m_2}g{l_2}\sin{\alpha _2}.\]

Hence, the first Lagrange equation can be written as

\[\frac{d}{{dt}}\left[ {\left( {{m_1} + {m_2}} \right)l_1^2{{\dot \alpha }_1} + {m_2}{l_1}{l_2}{{\dot \alpha }_2}\cos \left( {{\alpha _1} - {\alpha _2}} \right)} \right] + {m_2}{l_1}{l_2}{\dot \alpha _1}{\dot \alpha _2}\sin \left( {{\alpha _1} - {\alpha _2}} \right) + \left( {{m_1} + {m_2}} \right)g{l_1}\sin {\alpha _1} = 0,\]

\[\Rightarrow \left( {{m_1} + {m_2}} \right)l_1^2{\ddot \alpha _1} + {m_2}{l_1}{l_2}{\ddot \alpha _2}\cos \left( {{\alpha _1} - {\alpha _2}} \right) + {m_2}{l_1}{l_2}\dot \alpha _2^2\sin\left( {{\alpha _1} - {\alpha _2}} \right) + \left( {{m_1} + {m_2}} \right)g{l_1}\sin {\alpha _1} = 0.\]

Cancelling \({l_1} \ne 0,\) we obtain:

\[\left( {{m_1} + {m_2}} \right){l_1}{\ddot \alpha _1} + {m_2}{l_2}{\ddot \alpha _2}\cos \left( {{\alpha _1} - {\alpha _2}} \right) + {m_2}{l_2}\dot \alpha _2^2\sin\left( {{\alpha _1} - {\alpha _2}} \right) + \left( {{m_1} + {m_2}} \right)g\sin {\alpha _1} = 0.\]

Similarly, we derive the second differential equation:

\[\frac{d}{{dt}}\left[ {{m_2}l_2^2{{\dot \alpha }_2} + {m_2}{l_1}{l_2}{{\dot \alpha }_1}\cos \left( {{\alpha _1} - {\alpha _2}} \right)} \right] - {m_2}{l_1}{l_2}{\dot \alpha _1}{\dot \alpha _2}\sin \left( {{\alpha _1} - {\alpha _2}} \right) + {m_2}g{l_2}\sin {\alpha _2} = 0,\]

\[\Rightarrow {m_2}l_2^2{\ddot \alpha _2} + {m_2}{l_1}{l_2}{\ddot \alpha _1}\cos \left( {{\alpha _1} - {\alpha _2}} \right) - {m_2}{l_1}{l_2}\dot \alpha _1^2\sin \left( {{\alpha _1} - {\alpha _2}} \right) + {m_2}g{l_2}\sin {\alpha _2} = 0.\]

After canceling \({m_2}{l_2} \ne 0\) the equation takes the form

\[{l_2}{\ddot \alpha _2} + {l_1}{\ddot \alpha _1}\cos \left( {{\alpha _1} - {\alpha _2}} \right) - {l_1}\dot \alpha _1^2\sin \left( {{\alpha _1} - {\alpha _2}} \right) + g\sin {\alpha _2} = 0.\]

Thus, the nonlinear system of two Lagrange differential equations can be written as

\[\left\{ {\begin{array}{*{20}{l}} {\left( {{m_1} + {m_2}} \right){l_1}{{\ddot \alpha }_1} + {m_2}{l_2}{{\ddot \alpha }_2}\cos \left( {{\alpha _1} - {\alpha _2}} \right) + {m_2}{l_2}\dot \alpha _2^2\sin \left( {{\alpha _1} - {\alpha _2}} \right) + \left( {{m_1} + {m_2}} \right)g\sin {\alpha _1} = 0}\\[0.4em] {{l_2}{{\ddot \alpha }_2} + {l_1}{{\ddot \alpha }_1}\cos \left( {{\alpha _1} - {\alpha _2}} \right) - {l_1}\dot \alpha _1^2\sin \left( {{\alpha _1} - {\alpha _2}} \right) + g\sin {\alpha _2} = 0} \end{array}.} \right.\]

Small Oscillations of the Double Pendulum

Assuming that the angles \({\alpha _1}\left( t \right),{\alpha _2}\left( t \right)\) are small, the oscillations of the pendulums near the zero equilibrium point can be described by a linear system of equations. To get such a system, let's return back to the original Lagrangian of the system:

\[L = T - V = \left( {\frac{{{m_1}}}{2} + \frac{{{m_2}}}{2}} \right)l_1^2\dot \alpha _1^2 + \frac{{{m_2}}}{2}l_2^2\dot \alpha _2^2 + {m_2}{l_1}{l_2}{\dot \alpha _1}{\dot \alpha _2}\cos \left( {{\alpha _1} - {\alpha _2}} \right) + \left( {{m_1} + {m_2}} \right)g{l_1}\cos {\alpha _1} + {m_2}g{l_2}\cos {\alpha _2}.\]

We write this Lagrangian in a simpler form, expanding it in a Maclaurin series and retaining the linear and quadratic terms. The trigonometric functions can be replaced by the following approximate expressions:

\[\cos {\alpha _1} \approx 1 - \frac{{\alpha _1^2}}{2},\;\; \cos {\alpha _2} \approx 1 - \frac{{\alpha _2^2}}{2},\;\; \cos \left( {{\alpha _1} - {\alpha _2}} \right) \approx 1 - \frac{{{{\left( {{\alpha _1} - {\alpha _2}} \right)}^2}}}{2} \approx 1.\]

Here we have taken into account that the term with \(\cos \left( {{\alpha _1} - {\alpha _2}} \right)\) contains the product of small quantities \({{\dot \alpha }_1}{{\dot \alpha }_2}\) and has the second order of smallness. Therefore, we can leave only the linear term in the cosine expansion.

Substituting this in the original Lagrangian and considering that the potential energy is defined up to a constant, we obtain the quadratic Lagrangian for the double pendulum in the form:

\[L = T - V = \left( {\frac{{{m_1}}}{2} + \frac{{{m_2}}}{2}} \right)l_1^2\dot \alpha _1^2 + \frac{{{m_2}}}{2}l_2^2\dot \alpha _2^2 + {m_2}{l_1}{l_2}{\dot \alpha _1}{\dot \alpha _2} + \left({m_1 + m_2}\right)g{l_1} - \left( {\frac{{{m_1}}}{2} + \frac{{{m_2}}}{2}} \right)g{l_1}\alpha _1^2 + {m_2}g{l_2} - \frac{{{m_2}}}{2}g{l_2}\alpha _2^2.\]

We derive the Lagrange differential equations for the given Lagrangian. They are written as

\[\frac{d}{{dt}}\frac{{\partial L}}{{\partial {{\dot \alpha }_1}}} - \frac{{\partial L}}{{\partial {\alpha _1}}} = 0,\;\; \frac{d}{{dt}}\frac{{\partial L}}{{\partial {{\dot \alpha }_2}}} - \frac{{\partial L}}{{\partial {\alpha _2}}} = 0.\]

Find the partial derivatives:

\[\frac{{\partial L}}{{\partial {{\dot \alpha }_1}}} = \left( {{m_1} + {m_2}} \right)l_1^2{{\dot \alpha }_1} + {m_2}{l_1}{l_2}{{\dot \alpha }_2},\]

\[\frac{{\partial L}}{{\partial {\alpha _1}}} = - \left( {{m_1} + {m_2}} \right)g{l_1}{\alpha _1},\]

\[\frac{{\partial L}}{{\partial {{\dot \alpha }_2}}} = {m_2}l_2^2{{\dot \alpha }_2} + {m_2}{l_1}{l_2}{{\dot \alpha }_1},\]

\[\frac{{\partial L}}{{\partial {\alpha _2}}} = - {m_2}g{l_2}{\alpha _2}.\]

We get the system of two differential equations

\[\left\{ {\begin{array}{*{20}{l}} {\frac{d}{{dt}}\left[ {\left( {{m_1} + {m_2}} \right)l_1^2{{\dot \alpha }_1} + {m_2}{l_1}{l_2}{{\dot \alpha }_2}} \right] + \left( {{m_1} + {m_2}} \right)g{l_1}{\alpha _1} = 0}\\[0.4em] {\frac{d}{{dt}}\left[ {{m_2}l_2^2{{\dot \alpha }_2} + {m_2}{l_1}{l_2}{{\dot \alpha }_1}} \right] + {m_2}g{l_2}{\alpha _2} = 0} \end{array}.} \right.\]

\[\left\{ {\begin{array}{*{20}{l}} {\left( {{m_1} + {m_2}} \right)l_1^2{{\ddot \alpha }_1} + {m_2}{l_1}{l_2}{{\ddot \alpha }_2} + \left( {{m_1} + {m_2}} \right)g{l_1}{\alpha _1} = 0}\\[0.4em] {{m_2}{l_1}{l_2}{{\ddot \alpha }_1} + {m_2}l_2^2{{\ddot \alpha }_2} + {m_2}g{l_2}{\alpha _2} = 0} \end{array}.} \right.\]

This system of equations can be written in a compact matrix form. We introduce the matrices

\[ \boldsymbol{\alpha} \left( t \right) = \left[ {\begin{array}{*{20}{c}} {{\alpha _1}\left( t \right)}\\ {{\alpha _2}\left( t \right)} \end{array}} \right],\;\; M = \left[ {\begin{array}{*{20}{c}} {\left( {{m_1} + {m_2}} \right)l_1^2}&{{m_2}{l_1}{l_2}}\\ {{m_2}{l_1}{l_2}}&{{m_2}l_2^2} \end{array}} \right],\;\; K = \left[ {\begin{array}{*{20}{c}} {\left( {{m_1} + {m_2}} \right)g{l_1}}&0\\ 0&{{m_2}g{l_2}} \end{array}} \right],\;\; \mathbf{0} = \left[ {\begin{array}{*{20}{c}} 0\\ 0 \end{array}} \right].\]

Then the system of differential equations can be represented as

\[M\boldsymbol{\ddot \alpha} + K\boldsymbol{\alpha} = \mathbf{0}.\]

In the case of one body, this equation describes the free undamped oscillations with a certain frequency. In the case of the double pendulum, the solution (as you will see below) will contain oscillations with two characteristic frequencies, which are called normal modes. The normal modes are the real part of the complex-valued vector function

where \({\mathbf{H}_1},\) \({\mathbf{H}_2}\) are the eigenvectors, \(\omega\) is the real frequency. The values of the normal frequencies \({\omega _{1,2}}\) are determined by solving the auxiliary equation

\[\det \left( {K - {\omega ^2}M} \right) = 0.\]

In the case of arbitrary masses \({m_1},\) \({m_2}\) and lengths \({l_1},\) \({l_2},\) the auxiliary equation takes the form

\[\left( {{m_1} + {m_2}} \right){g^2} - {\omega ^2}\left( {{m_1} + {m_2}} \right)\left( {{l_1} + {l_2}} \right)g + {\omega ^4}{m_1}{l_1}{l_2} = 0.\]

Thus, we have a biquadratic equation for the frequencies \(\omega.\) The general solution for this equation is somewhat cumbersome. Therefore we consider the case when the lengths of the rods of both pendulums are equal: \({l_1} = {l_2} = l.\) Then the normal frequencies will be determined by a compact formula

\[\omega _{1,2}^2 = \frac{g}{l}\left[ {1 + \mu \pm \sqrt {\left( {1 + \mu } \right)\mu } } \right],\;\; \text{where}\;\;\mu = \frac{{{m_2}}}{{{m_1}}}.\]

As it can be seen, the eigenfrequencies \({\omega _{1,2}}\) depend only on the mass ratio \(\mu = {\frac{{{m_2}}}{{{m_1}}}}.\) The dependencies of the frequencies \({\omega _1},\) \({\omega _2}\) on the parameter \(\mu\) (provided \({\frac{g}{l}} = 1\)) are shown in Figure \(2\).

Dependencies of the eigenfrequencies on the mass ratio — Figure 2.

For equal masses \({m_1} =\) \( {m_2} = m,\) that is when \(\mu = 1,\) the frequencies are given by

\[{\omega _{1,2}} = \sqrt {\frac{g}{l}} \sqrt {2 \pm \sqrt 2 } .\]

Now, after the eigenfrequencies \({\omega _{1,2}}\) are known, we still have to determine the eigenvectors \({\mathbf{H}_{1,2}}\) to describe the normal modes. They can be found by solving the vector-matrix equation

\[\left( {K - {\omega ^2}M} \right)\mathbf{H} = \mathbf{0}.\]

Let the eigenvector \({\mathbf{H}_1} = {\left( {{H_{11}},{H_{21}}} \right)^T}\) (the superscript \(^T\) denotes transposition) be corresponded to the normal frequency \({\omega _1}.\) Then we have the following equation for \({\mathbf{H}_1}:\)

\[\left( {K - \omega _1^2M} \right){\mathbf{H}_1} = \mathbf{0},\;\; \text{where}\;\; \omega _1^2 = \frac{g}{l}\left[ {1 + \mu + \sqrt {\left( {1 + \mu } \right)\mu } } \right],\;\;\Rightarrow {m_1}l \left[ {\begin{array}{*{20}{c}} {\left( {1 \text{+} \mu } \right)\left( {g - \omega _1^2l} \right)}&{ - \omega _1^2\mu l}\\ { - \omega _1^2\mu l}&{\mu \left( {g - \omega _1^2l} \right)} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{H_{11}}}\\ {{H_{21}}} \end{array}} \right] = \mathbf{0}.\]

The coordinates of the eigenvector \({\mathbf{H}_1}\) satisfy the equation

\[\left( {1 + \mu } \right)\left( {g - \omega _1^2l} \right){H_{11}} - \omega _1^2\mu l{H_{21}} = 0,\]

\[\Rightarrow \frac{{{H_{21}}}}{{{H_{11}}}} = - \frac{{\left( {1 + \mu } \right)\left[ {\mu + \sqrt {\left( {1 + \mu } \right)\mu } } \right]}}{{\mu \left[ {1 + \mu + \sqrt {\left( {1 + \mu } \right)\mu } } \right]}} = - \sqrt {\frac{{1 + \mu }}{\mu }} .\]

Thus, the eigenvector \({\mathbf{H}_1}\) is given by

\[{\mathbf{H}_1} = \left[ {\begin{array}{*{20}{c}} {{H_{11}}}\\ {{H_{21}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1\\ { - \sqrt {\frac{{1 + \mu }}{\mu }} } \end{array}} \right].\]

Similarly, we find the coordinates of the second eigenvector \({\mathbf{H}_2} = {\left( {{H_{12}},{H_{22}}} \right)^T}:\)

\[{\mathbf{H}_2} = \left[ {\begin{array}{*{20}{c}} {{H_{12}}}\\ {{H_{22}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1\\ { \sqrt {\frac{{1 + \mu }}{\mu }}} \end{array}} \right].\]

The general solution of the matrix equation can be written as

\[\boldsymbol{\alpha} \left( t \right) = \left[ {\begin{array}{*{20}{c}} {{\alpha _1}\left( t \right)}\\ {{\alpha _2}\left( t \right)} \end{array}} \right] = \text{Re}\left( {\left[ {\begin{array}{*{20}{c}} {{\mathbf{H}_1}}\\ {{\mathbf{H}_2}} \end{array}} \right]{e^{i\omega t}}} \right) = {C_1}\left[ {\begin{array}{*{20}{c}} 1\\ { - \sqrt {\frac{{1 + \mu }}{\mu }} } \end{array}} \right]\cos \left( {{\omega _1}t + {\varphi _1}} \right) + {C_2}\left[ {\begin{array}{*{20}{c}} 1\\ {\sqrt {\frac{{1 + \mu }}{\mu }} } \end{array}} \right]\cos \left( {{\omega _2}t + {\varphi _2}} \right),\]

where the constants \({C_1},\) \({C_2},\) \({\varphi _1},\) \({\varphi _2}\) depend on the initial positions and velocities of the pendulums.

Consider the character of small oscillations for a specific set of initial data. Suppose, for example, that the initial positions and velocities of the pendulums have the following values:

\[{\alpha _1}\left( {t = 0} \right) = 0,\;\; {\alpha _2}\left( {t = 0} \right) = \frac{\pi }{6},\;\; {\dot \alpha _1}\left( {t = 0} \right) = 0,\;\; {\dot \alpha _2}\left( {t = 0} \right) = 0.\]

In this case, the initial phases are zero: \({\varphi _1} = {\varphi _2} = 0.\) Determine the constants \({C_1}\) and \({C_2}:\)

\[\left\{ \begin{array}{l} {\alpha _1}\left( 0 \right) = {C_1} + {C_2} = 0\\ {{\alpha _2}\left( 0 \right) = - {C_1}\sqrt {\frac{{1 + \mu }}{\mu }} + {C_2}\sqrt {\frac{{1 + \mu }}{\mu }} = \frac{\pi }{6}} \end{array} \right.,\;\; \Rightarrow {C_1} = - {C_2},\;\; \Rightarrow 2{C_2}\sqrt {\frac{{1 + \mu }}{\mu }} = \frac{\pi }{6},\;\; \Rightarrow {C_2} = \frac{\pi }{{12}}\sqrt {\frac{\mu }{{1 + \mu }}} ,\;\; \Rightarrow {C_1} = - \frac{\pi }{{12}}\sqrt {\frac{\mu }{{1 + \mu }}} .\]

Then the law of oscillations of the pendulums is expressed by the formulas

\[{\alpha _1}\left( t \right) = - \frac{\pi }{{12}}\sqrt {\frac{\mu }{{1 + \mu }}} \cos \left( {{\omega _1}t} \right) + \frac{\pi }{{12}}\sqrt {\frac{\mu }{{1 + \mu }}} \cos \left( {{\omega _2}t} \right),\]

\[{\alpha _2}\left( t \right) = \frac{\pi }{{12}}\cos \left( {{\omega _1}t} \right) + \frac{\pi }{{12}}\cos \left( {{\omega _2}t} \right),\]

where the angular frequencies \({\omega _{1,2}}\) are given by

\[{\omega _{1,2}} = \sqrt {\frac{g}{l}} \sqrt {1 + \mu \pm \sqrt {\left( {1 + \mu } \right)\mu } } .\]

Here the angles \({\alpha _1}\left( t \right),\) \({\alpha _2}\left( t \right)\) are expressed in radians, and the time \(t\) in seconds. Figures \(3-5\) show plots of small oscillations for three values of \(\mu:\) \({\mu _1} = 0.2,\) \({\mu _2} = 1,\) \({\mu _3} = 5,\) provided \(l = {l_1} = {l_2} = 0.25\,\text{m},\) \(g = 9.8{\frac{\text{m}}{{{\text{s}^2}}}}.\) For convenience, the deflection angles of the pendulums are given in degrees.

Small oscillations of the double pendulum for mu=0.2 — Figure 3.

Small oscillations of the double pendulum for mu=1 — Figure 4.

Small oscillations of the double pendulum for mu=5 — Figure 5.

As you can see from the graphs, the energy is cyclically transferred between the two bobs. When one of the pendulums almost stops, the other swings with maximum amplitude. After some time, the bobs "switch roles" and so on.

The resulting oscillations are represented as oscillations at the higher frequency \(\bar{\omega} = \frac{{{\omega_1} + {\omega_2}}}{2}\) with a periodic amplitude modulation with frequency \(\Delta\omega = \frac{{{\omega_1} - {\omega_2}}}{2},\) which is also known as the beat frequency.

Thus, the small oscillations of the double pendulum look as periodic changes and are described by the sum of two harmonics with frequencies \({\omega _1},\) \({\omega _2}\) depending on the system parameters.

Strictly speaking, the small oscillations of the double pendulum will be periodic if only the ratio of the eigenfrequencies \({\omega _1},\) \({\omega _2}\) is equal to a rational number. If the ratio of the frequencies is an irrational number, the small oscillations cannot possibly be periodic. For more details see the topic Sum of two periodic functions is periodic?

Legendre Transformation and Hamilton Equations

We now turn back to the original nonlinear system of equations and examine the character of oscillations of arbitrary amplitude. This system of equations can not be solved analytically. Therefore, we consider a numerical model of the double pendulum.

The Lagrange equations given above are second order differential equations. It is more conveniently to convert them into the form of Hamilton's canonical equations. As a result, instead of the two second-order equations, we obtain a system of four differential equations of the first order.

In Hamiltonian mechanics, the state of a system is determined by the generalized coordinates and generalized momenta. In our case, we can use again as in the Lagrange equations the angles \({\alpha _1},{\alpha _2}\) as the generalized coordinates. Instead of the generalized velocities \({\dot \alpha _1},{\dot \alpha _2}\) (in the Lagrangian), we now introduce the generalized momenta \({p_1},{p_2},\) related to the velocities by the formulas

\[{p_1} = \frac{{\partial L}}{{\partial {{\dot \alpha }_1}}},\;\; {p_2} = \frac{{\partial L}}{{\partial {{\dot \alpha }_2}}},\]

or in the brief description:

\[{p_i} = \frac{{\partial L}}{{\partial {{\dot \alpha }_i}}},\;\; i = 1,2.\]

The transition from the Lagrangian to the Hamiltonian form of the equations is performed using the Legendre transformation, which is defined as follows.

Suppose that \(f\left( x \right)\) is a smooth convex downward function (Figure \(6\)).

Consider the line \(y = px\) passing through the origin. The distance between the line \(y = px\) and the function \(y = f\left( x \right)\) along the \(y\)-axis depends on the coordinate \(x.\) This distance will be maximal at a certain value of \(x.\) Clearly, it depends on the slope of the line, that is on the parameter \(p.\) Thus, we introduce a new function \(g\left( p \right):\)

\[g\left( p \right) = \max\limits_x \left( {px - f\left( x \right)} \right).\]

Such transformation of the function \(f\left( x \right)\) into the conjugate function \(g\left( p \right)\) is called the Legendre transformation. Note that the function \(g\left( p \right)\) reaches a maximum value with respect to the variable \(x\) when \(p = \frac{{df}}{{dx}}.\) Indeed,

\[\frac{d}{{dx}}\left( {px - f\left( x \right)} \right) = p - \frac{{df\left( x \right)}}{{dx}} = 0,\;\; \Rightarrow p = \frac{{df\left( x \right)}}{{dx}} = p\left( x \right).\]

Knowing the dependence \(p\left( x \right),\) we can find the inverse function \(x\left( p \right).\) Then the Legendre transform is expressed by the relationship

\[g\left( p \right) = px\left( p \right) - f\left( {x\left( p \right)} \right),\;\; \text{where}\;\;p = \frac{{df}}{{dx}}.\]

The Legendre transformation is easily generalized to the case of functions of several variables. In the model of the double pendulum, the transition from the Lagrangian to the Hamiltonian is described by the Legendre transformation of the form:

\[H\left( {{\alpha _1},{\alpha _2},{p_1},{p_2}} \right) = \sum\limits_{i = 1}^2 {{{\dot \alpha }_i}{p_i}} - L\left( {{\alpha _1},{\alpha _2},{{\dot \alpha }_1},{{\dot \alpha }_2}} \right) = {\dot \alpha _1}{p_1} + {\dot \alpha _2}{p_2} - L\left( {{\alpha _1},{\alpha _2},{{\dot \alpha }_1},{{\dot \alpha }_2}} \right),\;\; \text{where}\;\;{p_1} = \frac{{\partial L}}{{\partial {{\dot \alpha }_1}}},\; {p_2} = \frac{{\partial L}}{{\partial {{\dot \alpha }_2}}}.\]

In this expression, \(L\) is the Lagrangian, and the function \(H\) is the Hamiltonian of the system, which depends on the generalized coordinates \({\alpha _1},{\alpha _2}\) and the generalized momenta \({p_1},{p_2}.\)

As a result of this transformation, each Lagrange equation becomes a system of two Hamilton's canonical equations of the form:

\[\frac{d}{{dt}}\frac{{\partial L}}{{\partial {{\dot \alpha }_i}}} = \frac{{\partial L}}{{\partial {\alpha _i}}},\;\; \Rightarrow \left\{ {\begin{array}{*{20}{l}} {{{\dot p}_i} = - \frac{{\partial H}}{{\partial {\alpha _i}}}}\\ {{{\dot \alpha }_i} = \frac{{\partial H}}{{\partial {p_i}}}} \end{array}} \right..\]

Omitting the formal (and cumbersome) algebraic transformations, we write the Hamilton's canonical equations for the double pendulum in the final form:

\[\left\{ \begin{array}{l} {{\dot \alpha }_1} = \frac{{{p_1}{l_2} - {p_2}{l_1}\cos \left( {{\alpha _1} - {\alpha _2}} \right)}}{{l_1^2{l_2}\left[ {{m_1} + {m_2}{\sin^2}\left( {{\alpha _1} - {\alpha _2}} \right)} \right]}}\\[0.4em] {{\dot \alpha }_2} = \frac{{{p_2}\left( {{m_1} + {m_2}} \right){l_1} - {p_1}{m_2}{l_2}\cos \left( {{\alpha _1} - {\alpha _2}} \right)}}{{{m_2}{l_1}l_2^2\left[ {{m_1} + {m_2}\,{\sin^2}\left( {{\alpha _1} - {\alpha _2}} \right)} \right]}}\\[0.4em] {{{\dot p}_1} = - \left( {{m_1} + {m_2}} \right)g{l_1}\sin {\alpha _1} }-{ {A_1} + {A_2}} \\[0.4em] {{{\dot p}_2} = - {m_2}g{l_2}\sin {\alpha _2} }+{ {A_1} - {A_2}} \end{array} \right.,\]

where

\[{A_1} = \frac{{{p_1}{p_2}\sin\left( {{\alpha _1} - {\alpha _2}} \right)}}{{{l_1}{l_2}\left[ {{m_1} + {m_2}\,{{\sin }^2}\left( {{\alpha _1} - {\alpha _2}} \right)} \right]}},\]

\[{A_2} = \frac{1}{{2l_1^2l_2^2{{\left[ {{m_1} + {m_2}\,{{\sin }^2}\left( {{\alpha _1} - {\alpha _2}} \right)} \right]}^2}}} \left[ {p_1^2{m_2}l_2^2 - 2{p_1}{p_2}{m_2}{l_1}{l_2}\cos \left( {{\alpha _1} - {\alpha _2}} \right) + p_2^2\left( {{m_1} + {m_2}} \right)l_1^2} \right] \sin \left[ {2\left( {{\alpha _1} - {\alpha _2}} \right)} \right].\]

We can now proceed to the numerical analysis of the equations.

Numerical Simulation of Chaotic Oscillations

The most common method of numerical solution of differential equations is the \(4\)th or \(5\)th order Runge-Kutta method. Different variations of this method are used in most mathematical packages (usually with an automatic error control and adaptive time-stepping).

To model the motion of the double pendulum, we also use the classical \(4\)th order Runge-Kutta method \(\left(RK4\right).\) We somewhat simplify the differential equations assuming that the lengths of the pendulums are the same: \({l_1} =\) \( {l_2} = l.\) By introducing the parameter \(\mu\) equal to the mass ratio: \(\mu = \frac{{{m_2}}}{{{m_1}}},\) we can write the system in the following form:

\[\left\{ \begin{array}{l} {{\dot \alpha }_1} = \frac{{{p_1} - {p_2}\cos \left( {{\alpha _1} - {\alpha _2}} \right)}}{{{m_1}{l^2}\left[ {1 + \mu \,{\sin^2}\left( {{\alpha _1} - {\alpha _2}} \right)} \right]}}\\[0.4em] {{\dot \alpha }_2} = \frac{{{p_2}\left( {1 + \mu } \right) - {p_1}\mu \cos \left( {{\alpha _1} - {\alpha _2}} \right)}}{{{m_1}{l^2}\left[ {1 + \mu \,{\sin^2}\left( {{\alpha _1} - {\alpha _2}} \right)} \right]}}\\[0.4em] {{{\dot p}_1} = - {m_1}\left( {1 + \mu } \right)gl\sin {\alpha _1} - {A_1} + {A_2}}\\[0.4em] {{{\dot p}_2} = - {m_1}\mu gl\sin {\alpha _2} + {A_1} - {A_2}} \end{array} \right.,\]

where

\[{A_1} = \frac{{{p_1}{p_2}\sin \left( {{\alpha _1} - {\alpha _2}} \right)}}{{{m_1}{l^2}\left[ {1 + \mu \,{{\sin }^2}\left( {{\alpha _1} - {\alpha _2}} \right)} \right]}},\]

\[{A_2} = \frac{1}{{2{m_1}{l^2}{{\left[ {1 + \mu \,{{\sin }^2}\left( {{\alpha _1} - {\alpha _2}} \right)} \right]}^2}}} \left[ {p_1^2\mu - 2{p_1}{p_2}\mu \cos \left( {{\alpha _1} - {\alpha _2}} \right) + p_2^2\left( {1 + \mu } \right)} \right] \sin \left[ {2\left( {{\alpha _1} - {\alpha _2}} \right)} \right] .\]

This system can be rewritten in vector form:

\[\mathbf{Z'} = \mathbf{f}\left( \mathbf{Z} \right),\;\; \text{where}\;\; \mathbf{Z} = {\left( {{\alpha _1},{\alpha _2},{p_1},{p_2}} \right)^T},\;\; \mathbf{f} = {\left( {{f_1},{f_2},{f_3},{f_4}} \right)^T}.\]

The vector \(\mathbf{Z}\) is composed of \(4\) canonical variables of the system, and components of the vector \(\mathbf{f}\) correspond to the right hand sides of the differential equations.

The Runge-Kutta \(\left(RK4\right)\) method requires at each step sequential evaluation of the four intermediate vectors:

\[{\mathbf{Y}_1} = \tau \mathbf{f}\left( {{\mathbf{Z}_n}} \right),\;\; {\mathbf{Y}_2} = \tau \mathbf{f}\left( {{\mathbf{Z}_n} + \frac{1}{2}{\mathbf{Y}_1}} \right),\;\; {\mathbf{Y}_3} = \tau \mathbf{f}\left( {{\mathbf{Z}_n} + \frac{1}{2}{\mathbf{Y}_2}} \right),\;\; {\mathbf{Y}_4} = \tau \mathbf{f}\left( {{\mathbf{Z}_n} + {\mathbf{Y}_3}} \right).\]

The value of the vector \({\mathbf{Z}_{n + 1}}\) in the next node is given by

\[{\mathbf{Z}_{n + 1}} = {\mathbf{Z}_n} + \frac{1}{6}\left( {{\mathbf{Y}_1} + 2{\mathbf{Y}_2} + 2{\mathbf{Y}_3} + {\mathbf{Y}_4}} \right).\]

The total error of the algorithm on a finite interval has the order \(O\left( {{\tau ^4}} \right),\) i.e. the computational accuracy increases by \(16\) times while reducing the time step \(\tau\) twice.

The described model is the basis for the animation shown at the beginning of this page. For simplicity, we chose the same initial angles of deflection of the pendulum bobs: \({\alpha _1} = {\alpha _2} = \alpha .\) Some combinations of the parameters \(\alpha\) and \(\mu\) produce strange attractors like those shown in Figure 7. As you can see, the double pendulum has amazing nonlinear dynamics.

An attractor of the double pendulum oscillations — Figure 7 (μ = 2.75, α = 171°).