# A numerical analysis of chaos in the double pendulum

Tomasz Stachowiak, Toshio Okada

► To cite this version:

Tomasz Stachowiak, Toshio Okada. A numerical analysis of chaos in the double pendulum. Chaos, Solitons & Fractals, 2006, 29 (2), pp.417 - 422. 10.1016/j.chaos.2005.08.032. hal-01389907

HAL Id: hal-01389907

https://hal.science/hal-01389907

Submitted on 30 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A numerical analysis of chaos in the double pendulum

Tomasz Stachowiak<sup>a</sup>, Toshio Okada<sup>b</sup>

We analyse the double pendulum system numerically, using a modified mid point integrator. Poincaré sections and bifurcation diagrams are constructed for certain, characteristic values of energy. The largest Lyapunov characteristic exponents are also calculated. All three methods confirm the passing of the system from the regular low energy limit into chaos as energy is increased.

## 1. Introduction

The planar double pendulum, as shown in Fig. 1, is one of the mechanical systems that exhibits chaotic behaviour, despite their simple construction. As it is a Hamiltonian system, its energy is a conserved quantity determining the global properties of motion at the two limits of zero and infinite energy, the system's motion is regular. (And in the infinite limit, the total angular momentum might be considered as the second conserved quantity.) It is the interval between these limits that presents typical chaotic features. Although there has been no rigorous proof of the system being non integrable, the analysis of chaos has given many results indicating such a possibility (see [1] for further references).

In this paper, we concentrate on three chaotic aspects. In Section 2 we set the notation and derive the equations used, together with the parameter values chosen for simulations. The subsequent sections contain a presentation of Poincaré sections given for increasing values of energy, bifurcation diagrams for regular, quasi regular and chaotic regions, and finally, calculation of the largest Lyapunov characteristic exponents.

## 2. Main equations

With the notation used in Fig. 1, we can immediately obtain the Lagrange function for the system:

$$ \mathcal{L} = \frac{1}{2} [l_1^2(m_1 + m_2)\dot{\varphi}_1^2 + l_2^2 m_2 \dot{\varphi}_2^2 + 2l_1 l_2 m_2 \cos(\varphi_1 - \varphi_2)\dot{\varphi}_1\dot{\varphi}_2] + gl_1(m_1 + m_2)\cos\varphi_1 + gl_2 m_2 \cos\varphi_2, \quad (1) $$

<footnote><sup>a</sup> Astronomical Observatory, Jagiellonian University, Orla 171, 30 244 Kraków, Poland</footnote>

<footnote><sup>b</sup> Tokyo University of Agriculture and Technology, Naka machi, Koganei shi, Tokyo 184 8588, Japan</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. The physical setup of the system.</figcaption></figure>

where the dot denotes the derivative with respect to time. We choose not to employ the Hamilton's formalism, and use as two additional variables the angular velocities instead of momenta. Thus, the second order equations of motion:

$$ 
\begin{aligned}
l_1(m_1 + m_2)\ddot{\varphi}_1 + l_2 m_2 [\cos(\varphi_1 \quad \varphi_2)\ddot{\varphi}_2 + \sin(\varphi_1 \quad \varphi_2)\dot{\varphi}_2^2] + g(m_1 + m_2) \sin \varphi_1 &= 0, \\
l_2 \ddot{\varphi}_2 + l_1 [\cos(\varphi_1 \quad \varphi_2)\dot{\varphi}_1 \sin(\varphi_1 \quad \varphi_2)\dot{\varphi}_1^2] + g \sin \varphi_2 &= 0
\end{aligned} 
\quad (2) $$

are reduced to a system of four equations of the first order:

$$ \dot{\varphi}_1 = \omega_1, $$

$$ \dot{\varphi}_2 = \omega_2, $$

$$ \dot{\omega}_1 = \frac{\sin(\varphi_1 \quad \varphi_2)[l_1 \cos(\varphi_1 \quad \varphi_2)\omega_1^2 + \omega_2^2]}{2l_1[1 + m_1 \cos^2(\varphi_1 \quad \varphi_2)]} \quad \frac{(1 + 2m_1) \sin \varphi_1 + \sin(\varphi_1 \quad 2\varphi_2)}{l_1[1 + m_1 \cos^2(\varphi_1 \quad \varphi_2)]}, \quad (3) $$

$$ \dot{\omega}_2 = \sin(\varphi_1 \quad \varphi_2) \frac{(1 + m_1)(\cos \varphi_1 + l_1 \omega_1^2) + \cos(\varphi_1 \quad \varphi_2) \omega_2^2}{1 + m_1 \cos^2(\varphi_1 \quad \varphi_2)}, $$

where $m_2 = 1$ is the unit of mass, $l_2 = 1$ of length and $g = 1$, which effectively sets a unit of time. For practical application, we choose particular values of the remaining parameters as: $m_1 = 3$ and $l_1 = 2$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Poincaré section for $E = E_1 = E_{01} + 0.01$; regular behaviour.</figcaption></figure>
There are four critical (or stationary) points:

$$ (\varphi_1, \varphi_2, \omega_1, \omega_2) = \begin{cases} (0, 0, 0, 0), & E_{01} = 1, \quad l_1(1+m_1) = 9, \\ (0, \pi, 0, 0), & E_{02} = +1, \quad l_1(1+m_1) = 7, \\ (\pi, 0, 0, 0), & E_{03} = 1+l_1(1+m_1) = 7, \\ (\pi, \pi, 0, 0), & E_{04} = +1+l_1(1+m_1) = 9, \end{cases} \quad (4) $$

with E being the corresponding energy values. The phase space sections and Lyapunov exponents have been obtained for three values of energy, close to the above, namely: $E_1 = E_{01} + 0.01$, $E_2 = E_{02} + 0.01$ and $E_3 = E_{03} + 0.01$.

The numerical analysis was performed using the Bulirsch Stoer method with Richardson extrapolation as described in [2].

### 3. Poincaré sections

As the whole phase space is four dimensional, it is convenient to look at its sections, i.e., the values of the coordinates whenever one of them assumes a certain value. Here the sections were constructed with the surface $\varphi_1 = 0$, and $\omega_1 = \dot{\varphi}_1 > 0$, and restricted to the plane of $\varphi_2$ and $\omega_2$.

The first one (Fig. 2) at the energy level of $E = E_1$ presents a very regular image, as would be expected from the regular limit $E = E_{01}$. The situation becomes more complex in Fig. 3, when the energy is increased to $E = E_2$, although without chaotic regions, the invariant tori are visibly deformed. Finally, in Fig. 4 almost all of the regular orbits decay into a global chaotic region.

### 4. Bifurcation diagrams

In order to detect periodic orbits and their number, regular windows in chaos, etc., a bifurcation diagram can be constructed. It depicts the dependence of a suitably chosen quantity reflecting the chaotic behaviour on some parameter of the system. As we fixed the values for mass and length, we used energy for that purpose. The other quantity, which should be different for different orbits, if they exist, could be the amplitude of one of the variables. Choosing

<figure><img src="image_4.png" /><figcaption>Fig. 3. Poincaré section for $E = E_2 = E_{02} + 0.01$; regular trajectories begin to deform.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Poincaré section for $E = E_3 = E_{03}$ 0.01; regular regions decay with the appearance of a global, completely chaotic region.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. The first appearance of global chaos, and “windows” that are magnified in Fig. 6.</figcaption></figure>

$\phi_2$ as this variable, we effectively needed to construct Poincaré sections with respect to $\dot{\phi}_2 = \omega_2$, with $\omega_2 < 0$, so that the maximum values of $\phi_2$ were obtained.

Fig. 5 shows the behaviour from the lowest energy $E_{01}$, up to the first appearance of global chaos. There are regions, shown magnified in Fig. 6, where the orbits' maxima tend to cluster at evenly spaced intervals of $\phi_2$, and “diverge” around such places, densely occupying the whole admissible region.

In Fig. 7 contains a window between two completely chaotic regions, showing smaller regions, of behaviour similar to Fig. 6, where there are only four possible values of $\phi_2$, which then chaotically increase.
<figure><img src="image_7.png" /><figcaption>Fig. 6. A “window” showing convergence of admissible values of $\phi_{2max}$; this is thus a regular behaviour.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. A gap between two globally chaotic regions, showing windows, and regular regions with bifurcations quickly leading to chaos.</figcaption></figure>

## 5. Lyapunov exponents

We limited our investigations to the largest Lyapunov exponents, as it effectively gives us the information on the divergence of two close trajectories. We used the simple method described, in [3]. The main formula is

$$ \chi = \frac{1}{tN} \sum_{i=1}^{N} \ln \frac{|d_i(t)|}{|d_i(0)|}, \qquad (5) $$

where $\chi$ denotes the Lyapunov exponents, the index $i$, consecutive initial positions, and $d$ is the separation between two close trajectories (Fig. 8).

The simulations were performed for the same values of energy as in the case of the Poincaré sections and the results are in very good agreement. For both $E_1$ and $E_2$, the exponents tend to zero, indicating practically regular behaviour.
<figure><img src="image_9.png" /><figcaption>Fig. 8. The largest Lyapunov exponents for the system. The regular behaviour of the two first cases, $E = E_1$ and $E = E_2$, makes the thin solid line, and the dashed line, respectively, almost coincide. The third line, corresponding to the chaotic case of $E_3$, indicates a positive value of the exponents.</figcaption></figure>

In the case of $E_3$, definite limit is not clearly visible, but it is obvious that the exponents stay positive in the vicinity of 0.179, thus corresponding to the chaos clearly visible in the third phase space section.

## 6. Summary

In this article we present further numerical results concerning the planar double pendulum system. Its chaotic features are fully exposed by all the tools used, namely, Poincaré sections, bifurcation diagrams and Lyapunov characteristic exponents. They reveal regular behaviour at the zero energy limit, which, as energy is increased, transforms into quasi regular and, finally, globally chaotic motion.

We note that the existence and decay of the invariant torii, is visible on the phase space sections, but the Lyapunov exponent analysis is not as evident. Together with the other aspects, though, it could be considered as a valid result.

The bifurcation diagrams are not standard, in the sense that energy was used as the parameter, instead of a quantity explicitly appearing in the equations of motion. However, it is with energy that the characteristics of the system change drastically, as the obtained diagrams proved. Thus, further evidence of chaotic behaviour such as windows and bifurcations were found, providing yet another image of chaos for the double pendulum.

## References

[1] Ohlhoff A, Richter PH. Forces in the double pendulum. Z Angew Math Mech 2000;80(8):517-34.

[2] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes in C. Cambridge: Cambridge University Press; 1992.

[3] Wolf A, Swift JB, Swinney HL, Vastano JA. Determining Lyapunov exponents from a time series. Physica D 1985;16:285-315.