Deterministic chaos: A pedagogical review of the double pendulum case

Abstract

This review aims to provide a comprehensive evaluation of the dynamics of the double pendulum, with a particular emphasis on its chaotic behavior. It examines the complicated and unpredictable behavior of the double pendulum. It highlights the pedagogical value of the double pendulum as it bridges concepts across physics, mathematics, and computer science, providing a tangible demonstration of chaotic dynamics. Numerical methods, such as the 4th-order Runge-Kutta method, estimate solutions to the autonomous Hamiltonian equations that govern a system. We simulate the motion of a double pendulum, enabling the visualization of intricate trajectories and the analysis of emerging patterns with the Python and Fortran programming languages. We discuss the importance of Poincaré’s maps and the Lyapunov exponents in characterizing and quantifying the rate at which trajectories diverge in phase space to elucidate the chaotic nature of the system.

The educational significance of the double pendulum is emphasized in teaching key concepts in Classical Mechanics, Differential Calculus, Linear Algebra, and Numerical Methods for solving ordinary differential equations (ODEs). We also explore the interdisciplinary teaching opportunities presented by the double pendulum.

Keywords:
Double pendulum dynamics; chaotic behavior; pedagogical value; interdisciplinary approach; numerical simulation

1. Introduction

The pendular motion played more than a critical scientific role in forming the modern world. The pendulum was crucial in the horological revolution, closely linked to the scientific revolution [¹]. The double pendulum, a seemingly simple mechanical system, is a paradigmatic example of deterministic chaos [²] in classical mechanics [³]. This intriguing system, consisting of two pendulums attached end to end may exhibit fundamentally unpredictable behavior over long time scales despite being deterministic [⁴]. As such, the double pendulum is an invaluable pedagogical tool, bridging concepts in physics, mathematics, and computer science, providing a tangible demonstration of chaotic dynamics [⁵].

Deterministic chaos (see chapter 2 of [⁶] or page 102 of [⁷]) refers to the phenomenon where a system’s initial conditions determine its future behavior entirely yet remain inherently unpredictable due to its sensitive dependence on those initial conditions. In the realm of the double pendulum, slight variations in initial states can lead to vastly divergent trajectories, a hallmark of chaotic systems. This sensitivity poses unique challenges and opportunities for understanding the underlying mechanics of such systems.

The study of the double pendulum’s chaotic behavior often involves numerical methods [⁸], with the 4th-order Runge-Kutta method being a prominent tool. This method, a cornerstone in the numerical integration of ordinary differential equations, provides approximating solutions to the autonomous Hamiltonian equations that govern the pendulum’s motion. One can capture the essence of the pendulum’s chaotic dynamics through these computational approaches.

Python, the programming language, with its accessibility and widespread adoption, has become a popular tool for simulating the double pendulum’s motion. With an internet browser, one can implement Python scripts to visualize the pendulum’s intricate trajectories and analyze the emergent patterns. These simulations demonstrate the unpredictable nature of chaotic systems and serve as a practical application of computational physics. Further, the construction of Poincaré maps allows for the visualization of the system’s phase space behavior, unveiling features of chaos [⁹, ¹⁰].

An essential aspect of studying chaotic systems is the Lyapunov exponent [⁹]. This metric quantifies how rapidly trajectories diverge in phase space, with a positive Lyapunov exponent indicating a system’s sensitive dependence on initial conditions and its inherent chaotic nature. Calculating the Lyapunov exponent for the double pendulum provides a quantitative measure of its chaotic behavior, offering more profound insights into the dynamics of such complicated systems [¹¹].

This review highlights the double pendulum as an outstanding example of a deterministic chaotic system. We will delve into the exploration of its chaotic dynamics through the application of the 4th-order Runge-Kutta method to its governing autonomous Hamiltonian equations.

Moreover, we will highlight the educational value of the double pendulum, underscoring its utility in teaching key concepts in Classical Mechanics, Differential Calculus, Linear Algebra, and Numerical Methods for solving ordinary differential equations (ODEs). Additionally, we will discuss the role of Computational Coding in obtaining and analyzing numerical solutions, demonstrating the convergence of theoretical and practical aspects in the study of chaotic systems.

An excellent and concise review [¹²] emphasizes the analysis of the linear motion in the neighborhood of the equilibrium points of the double pendulum. Our review advances to the nonlinear chaotic regime with numerical and computational tools.

We introduce the simple pendulum in ^{Section 2}, which serves as a foundation for understanding the basic concepts and analytical tools needed to tackle the double-pendulum problem. In ^{Section 3}, we discuss the existence and uniqueness theorem and the Picard iteration process for solving initial value problems (IVPs) of ordinary differential equations (ODEs). We outline the steps involved in the Picard iteration process, including the definition of a sequence of functions, the Lipschitz condition, and the convergence of the sequence to a continuous function. We then apply the Runge-Kutta method to the simple pendulum problem, demonstrating how to estimate the computer time, time step, and numerical errors to keep under a tolerance limit. ^{Section 4} is dedicated to the double pendulum per se. We discuss the Hamiltonian formalism and its relevance to studying the system, splitting the second-order equation into two first-order equations. We also emphasize the importance of computational efficiency and memory usage in numerical simulations of ODEs, mentioning the use of optimized numerical methods and the need to monitor CPU time and memory usage during extended simulations. We present some tools [¹³] to discover exciting properties of the double pendulum, such as the small perturbations from stable equilibrium points, numerical solutions of its equations, the Poincaré’s maps, and Lyapunov characteristic exponents. In ^{Section 5}, we highlight the ample opportunities for teaching STEM subjects in an interdisciplinary manner, and we invite educators to explore these opportunities. Studying deterministic chaos in systems like the double pendulum is not merely an academic exercise. It provides valuable insights into real-world phenomena where chaos and unpredictability play crucial roles, from meteorological patterns to stock market fluctuations. The principles of chaos theory find applications in diverse fields, making the study of systems like the double pendulum both relevant and far-reaching in itsimpact.

2. Simple Pendulum

The simple pendulum, a fundamental concept in physics, is a gateway to understanding complex systems. It consists of a mass m (the pendulum bob) suspended from a pivot by a rod of length L, swinging back and forth under gravitational force mg. We consider the motion of the pendulum in the x - z plane, with the position of the mass given by:

(1) $$x=L\sin \theta ,z=-L\cos \theta$$

where θ is the angle between the rod and the vertical line (See Figure 1).

Figure 1
The mass point m can move on the circumference of radius L. θ is the angle of the vertical downwards direction and the Pendulum fixed length rod.

2.1. Equation of motion

The movement of the pendulum is driven by the libration angle θ, and the following differential equation can describe its motion on the circumference:

(2) $$mL\frac{d^2\theta }{d{t}^2}=-mg\sin \theta$$

Interestingly, in this ideal model (ignoring air resistance and friction), the mass m cancels out, reflecting Galileo’s principle that all bodies fall at the same rate under gravity g.

See Figure 1.

When we simplify the equation by removing the mass m and introducing a time scaling factor $\sqrt{L/g}$, we obtain a more manageable form:

(3) $$\frac{d^2\theta }{d{t}^2}+\sin \theta =0$$

a homogenous second-order ordinary differential equation non-linear in θ. Here, t represents time scaled in $\sqrt{L/g}$ units. For example, with L = 108cm and g = 981 cm/s², $\sqrt{L/g}=\text{1/3 s}$ making the mathematical time t = 6 equivalent to 2 physical seconds.

2.2. Energy conservation and phaseplane analysis

A pivotal concept in pendulum dynamics is energy conservation. By manipulating the motion equation (by multiplying equation (3) by dθ/dt ≠ 0), we can derive an expression for the system’s energy:

(4) $$\frac{d\theta }{dt}\left(\frac{d^2\theta }{d{t}^2}+\sin \theta \right)=0\Rightarrow \frac{d}{dt}\left(\frac{1}{2}{\left(\frac{d\theta }{dt}\right)}^2-\cos \theta \right)=0,$$

which leads us to the energy equation in terms of angular velocity w and angle θ

(5) $$\frac{1}{2}{w}^2-\cos \theta =\epsilon$$

where ε is the total energy (scaled in units of mgL).

Using this energy invariant, we can explore the phase [⁵] plane of the pendulum, a graphical representation of its state over time. Each point on the phase plane represents a state of the system, defined by θ and w. The pendulum’s behavior is predictable and follows distinct patterns in this plane, as shown in Figure 2.

Figure 2
The simple pendulum phase curves for ε = −0.71, 0.0, 1.0, 2.0, 3.0 from inside out. The portrait is periodic of period 2π along the horizontal θ axis. The phase plane also has the symmetry of the vertical axis w. ↔ −w.

2.3. Predictable yet complicated motion

The simple pendulum exhibits predictable motion, but its analytical solution can be complicated. The equation of motion becomes a first-order separable differential equation when expressed in terms of angular velocity:

(6) $$\frac{d\theta }{dt}=\pm \sqrt{2(\epsilon +\cos \theta )}$$

The integration of this equation involves special functions (Jacobi elliptic functions) and respects the symmetries observed in the phase plane.

That is, one gets the energy (in units of mgL):

(7) $$\frac{1}{2}{w}^2-\cos \theta =\epsilon ,$$

in which $w\equiv \frac{d\theta }{dt}$ is the angular velocity. For example, if m = 1 kg, and the same values used above, we have mgL ≈ 10,7 Joules.

From the mathematical point of view, the invariant ε, given by the initial state of θ(t₀ = θ₀) and w(t₀) = w₀ at initial time t₀ implicitly defines the motion:

(8) $$\epsilon =\frac{1}{2}{\left(w_0\right)}^2-\cos {\theta }_0\ge -1$$

Thus, one gets the so-called phase plane for the pendulum. See Figure 2. For ε = −1, there is no motion, θ = 0 = w. The green curve of ε = 1 is a separatrix between bound and unbound motions. In other words, the pendulum at bound motion, −1 ＜ ε ＜ 1, swings between the θ_max = ±arcoss (−ε) and maximum angular velocity at $w_{max}=\pm \sqrt{2(1+\epsilon )}$. For ε ＞ 1 the pendulum rotates the full circumference having the angular velocity between $w_{\pm }\text{=}\pm \sqrt{2(\epsilon \pm 1)}$ – the pendulum becomes a “rotator”.

All the curves in the phase plane are closed, either the librations between fixed values for |θ| ＜ π or for |θ| ＞ π, because the physical position of the pendulum, θ ↔ θ + 2π. Thus, the 2D phase space is cylindrical.

Therefore, the simple pendulum is highly predictable from the beginning; there is no chaos. Nevertheless, in general, the analytical solution is complicated.

2.4. Computational approach to the pendulum’s motion

While the analytical approach provides deep insights, studying the pendulum’s motion using numerical methods is often more practical. Especially for understanding the total period of the pendulum’s motion, whether in a bound or unbound state.

For bound motions, one gets the period integrating over the range of motion, which involves handling the singularity at the limits of motion using advanced numerical techniques like the Radau quadrature. The period formula for bound motion is:

(9) $$T_{\text{bound}}=4{{\displaystyle \int }}_0^{{\theta }_m}\frac{d\theta }{\sqrt{2(\epsilon +\cos \theta )}}$$

where θ_m = |arccos(−ε)|.

In contrast, the period for unbound motions is more straightforward to compute numerically, as the integrand does not encounter singularities within the integration range. This period is given by

From equation (7) one gets w² = 2(ε + cos θ). Therefore

(10) $$w=\frac{d\theta }{dt}=\pm \sqrt{2(\epsilon +\cos \theta )}$$

which is a first-order separable differential equation. One has to chose possibles intervals in which θ and t are one-to-one so that

(11) $$dt=\pm \frac{d\theta }{\sqrt{2(\epsilon +\cos \theta )}}$$

and then integrate. The unusual Jacobi elliptic functions can express the results. The symmetries we pointed out at the phase plane could give the remaining motion.

In practice, we do not follow this equation approach (11) to study the motion itself, but one can use it to compute the total period of the motion. For the bound motions

(12) $$T_{\text{bound}}=4{{\displaystyle \int }}_0^{{\theta }_m}\frac{d\theta }{\sqrt{2(\epsilon +\cos \theta )}}$$

in which θ_m = |arccos (−ε)|. The integral is not feasible to compute by the Cotes method because ε+cos θ_m = 0. So, one has to use a quadrature formula or a method of Radau that avoids the integrand’s singularity.

For unbound motions

(13) $$T_{\text{unbound}}=2{{\displaystyle \int }}_0^{\pi }\frac{d\theta }{\sqrt{2(\epsilon +\cos \theta )}}$$

is easily computed numerically. Remember that for unbound motions ε ＞ 1 ⇒ ε + cos θ ＞ 0, ∀θ.

Thus, instead of equation (11), we recommend to solve numerically the following initial valued problem (IVP) from the equation (3):

(14) $$\begin{array}{l}\frac{d\theta }{dt}=w\\ \frac{dw}{dt}=-\sin \theta \\ \theta (t_0)={\theta }_0\\ w(t_0)=w_0\end{array}$$

Let us set, as most of numerical methods texts do, the arrays

(15) $$y(t)=\begin{array}{l}\left(\begin{array}{l}\theta (t)\\ w(t)\end{array}\right)\hfill \end{array},\mathrm{}{y}_0=\begin{array}{l}\left(\begin{array}{l}{\theta }_0\\ w_0\end{array}\right)\hfill \end{array},\mathrm{}f(y)=\begin{array}{c}\left(\begin{array}{c}w\\ -\sin \theta \end{array}\right)\hfill \end{array}$$

so that our IVP is

(16) $$\frac{dy}{dt}\text{=}f\text{(}y\text{)}$$

(17) $$y(t_0)\text{=}{y}_{\text{0}}$$

a so-called autonomous dynamical system because f does not depend explicitly on t.

Let us show next that the equation (16), together with its initial condition (17), has a unique solution with simple requests on the function f, and thus, a numerical computational approach is worth it.

3. Initial Valued Problems, Existence, Uniqueness Solutions, and Runge-Kutta Methods

A fundamental result in the theory of differential equations [¹⁴], particularly for autonomous systems where the rate of change of the dependent variable is a function of the variable itself and not explicitly of the independent variable, is the Existence and Uniqueness Theorem for first-order ordinary differential equations (ODEs):

Theorem 1.

Let f: D ⊂ ℝ² → ℝ be a function such that f(y) is Lipschitz continuous on an open interval I containing y₀. Then, for any point (t₀, y₀) in the domain D, there exists an interval J containing t₀ such that there is a unique function y(t), which is a solution to the differential equation y′ = f(y) on J, satisfying the initial condition y(t₀) = y₀.

There are several proofs [¹⁵] of ^{Theorem 1.} One proof [¹⁶] applies the Picard-Lindelöf iteration method and the contraction mapping principle. Here is a sketch of this proof:

Proof.

• 1.

  Convert the IVP into an equivalent integral equation:

$$y(t)=y_0+{{\displaystyle \int }}_{t_0}^{t}f(y(s))ds$$

• 2.

  Define a sequence of functions {y_n(t)} by the Picard iteration process:

$$y_{n+1}(t)=y_0+{{\displaystyle \int }}_{t_0}^{t}f(y_n(s))ds$$

where y₁(t) = y₀ is the initial approximation.

• 3.

  Show that the Picard iteration is a contraction mapping under the given conditions. This involves demonstrating that there exists a constant L ＜ 1 such that for all t in J:

$$|y_{n+1}(t)-y_n(t)|\le L|y_n(t)-y_{n-1}(t)|$$

This step relies on the Lipschitz condition of f.

• 4.

  Prove that the sequence {y_n(t)} is a Cauchy sequence in the space of continuous functions on J. The completeness of this function space ensures that {y_n(t)} converges uniformly to a continuous function y(t).

• 5.

  Verify that the limit function y(t) satisfies the integral equation, and hence the original differential equation, bypassing the limit inside the integral.

• 6.

  Prove that if y₁(t) and y₂(t) are two solutions to the IVP, then they must be identical by showing that the difference |y₁(t) − y₂(t)| is zero, using the Gronwall’s inequality.

The Lipschitz condition is crucial for the uniqueness part of the theorem. It requires the function f(y) to be continuous and does not change too much with y. If d f(y)/dy is also continuous, then f satisfies the Lipschitz condition, but the Lipschitz condition is not enough to guarantee the continuity of d f(y)/dy.

^{Theorem 1} is a cornerstone in the study of differential equations, as it guarantees that under reasonable conditions on f, the behavior of a dynamical system can be predicted uniquely from its initial state. Furthermore, it guarantees that a numerical approach will not be in vain.

The Runge-Kutta methods are a family of iterative methods used for the approximate numerical solutions of ODEs [⁸]. These methods are among the most widely used and efficient techniques that provide a systematic procedure to generate accurate solutions to ODEs without the need for analytical derivations or integrations of the ODE.

The basic idea of the Runge-Kutta method is to approximate the solution of an ODE at successive points, with each step involving one or more evaluations of the derivative specified by the ODE. The method combines these evaluations to produce an approximation that is more accurate than the Euler method, which is a first-order method. The ’order’ of a Runge-Kutta method refers to the degree of accuracy of the approximation. Higher-order methods yield more accurate results but require more computational effort. The fourth-order Runge-Kutta (RK4) method is prevalent due to its balance between accuracy and computational efficiency.

The RK4 method is a fourth-order method, meaning that the error per step is of the order of h⁵, while the total error of several steps of the same size h is of the order of h⁴ [⁸].

The algorithm for our autonomous system is the following. Given an initial value problem specified by dy/dt = f(y), y(t₀) = y₀, the RK4 method calculates the value y_n+1 at t_n+1 = t_n + h using the following steps:

• Compute k₁ = f(y_n), the slope at the beginning of the interval.

• Compute k₂ = f(y_n + hk₁/2), the slope at the midpoint of the interval, using k1 to estimate the value of y at t_n + h/2.

• Compute k₃ = f(y_n + hk₂/2), another slope estimate at the midpoint, but now using k₂ to estimate y at t_n + h/2.

• Compute k₄ = f(y_n + hk₃), the slope at the end of the interval, with y estimated using k₃.

• Weight average these slopes to calculate the next value of y, i.e.,

$$y_{n+1}\text{=}{y}_n\text{+}\frac{h}{6}\left(k_1+2\mathrm{}\left(k_2+k_3\right)+k_4\right)$$

See the Figure 3.

Figure 3
The fourth-order Runge-Kutta method uses four slope estimates to compute a good approximation for y(t₀ + h) using y(t₀) and the ODE function f(y).

See the piece of the Python code in the Appendix 1.

While RK4 is efficient for many problems, it is not adaptive; the step size remains constant. For problems where the solution varies rapidly in certain regions, adjusting the step size based on the solution’s behavior may be more efficient. Furthermore, the error in the approximated solution is not included in the numerical evolution by RK4. Of course, we may use other methods with numerical error estimates and adjustable step sizes.

For the autonomous system, one can use the constant of motion to check whether or not the error is within a given tolerance.

One could explore a step as small as our float system computer system allows, but it would require too many steps, long computer times, and large computer memory to evolute the system from the initial time t₀ to the final t_f. So, it is prudent to evaluate the time RK4 takes for one step. Moreover, this time highly affects the computational cost of computing f(y).

In the case of the simple pendulum equation (16), one can use the IPython magic function %timeit to time a particular piece of code (a single execution statement, namely

(18) $$\%\text{timeit rk4(f, y0, h, 2)}$$

resulting¹ in μ s ± 818 ns, in seven runs with 10,000 loops each, the time per loop has a mean of 23.2 μs ± 818 ns of standard deviation. With this estimate, one can decide how many steps n we can afford to run.

The RK4 method is compatible with the Taylor expansion up to the fourth time derivative. Thus, the local error, according to the Taylor theorem, is

(19) $$erro{r}_{\text{RK}4}=\frac{d^{(5)}y}{d{t}^5}(\xi )\frac{h^5}{5!},\mathrm{}\xi \in [t_n,t_n+h]$$

This error is an analytical, theoretical result that, in general, it is difficult to take advantage of because it is complicated to obtain d⁽⁵⁾y/dt⁵ at a hypothetical ξ ∈ [t_n, t_n + h].

For the simple pendulum, we certainly can compute

(20) $$\frac{d^{(5)}\theta }{d{t}^5}=\left(4{\mathrm{}\cos }^2\mathrm{}\theta -3+w^2\cos \theta \right)w$$

using successively the equations (14). One can further use the constant (7) to get

(21) $$\frac{d^{(5)}\theta }{d{t}^5}=\sqrt{2\mathrm{}\left(\epsilon +\cos \theta \right)}\mathrm{}\left(2\epsilon \cos \theta +6{\mathrm{}\cos }^2\mathrm{}\theta -3\right)$$

Thus, one gets an upper bound, say

(22) $${\theta }_{\text{max}}^{(5)}\equiv \text{max}\mathrm{}\left|\frac{d^{(5)}\theta }{d{t}^5}\right|$$

Therefore, one has the upper bound on the error, and from equation (19), we get an upper bound on h for a given tolerance tol:

(23) $$h^5\le \frac{5!}{{\theta }_{\max }^{(5)}}\text{tol}$$

For the example we run below, with tol= 10⁻¹⁴, θ₀ = 5π/6, w = 0 ⇒ ε ≈ 0.8660, we get h ≤ 2.65 10⁻³. For these values, the pendulum swings with the period given by (12) T ≈ 6.392568. Thus, for three full swings, the final time t_f ≈ 33.22, which requires 12,524 steps. From our estimate for our computer, the evolution would run in 12524 × 23.2μs ≈ 0.29 s. Very fast indeed. See Figure 4.

Figure 4
The simple pendulum evolution for θ₀ = 5π/6, w₀ = 0 ⇒ ε ≈ 0.8660 in three full period of swing.

We also checked the conservation of the energy ε along the evolution: Maximum Relative and Absolute Error ≤ 6.7 10⁻¹³ and 5.8 10⁻¹³, respectively. High precision, indeed.

import time

start_time = time . process_time ()

y = rk4(f, y0 , h, n_steps )

a = y[: , 0]; w = y[: , 1]

end_time = time . process_time ()

elapsed_time = end_time - start_time

et = elapsed_time

print (" Elapsed CPU time : %3.2 f " % et , " seconds ")

As expected, we got Elapsed CPU time: 0.29 seconds.

One has to achieve the highest possible accuracy with affordable computer run-time to study chaotic motion or for extended times run. High-order methods are generally more accurate than low-order methods if we use sufficiently small step sizes at the expense of the more expensive work of completing each step [¹⁷]. Thus, for a given error tolerance, the time step for an RK5 method could be larger than that for the RK4 one; that is, one needs fewer steps to get to the final t_f. Nevertheless, the computer time spent in each step is longer since more evaluations are needed.

In our single pendulum, we can compute the error for the RK5, as done above for the RK4:

(24) $$erro{r}_{\text{RK}5}=\frac{d^{(6)}y}{d{t}^6}(\xi )\frac{h^6}{6!},\mathrm{}\xi \in [t_n,t_n+h]$$

We certainly can compute, using successively the equations (14) and the constant (7) to get

(25) $$\frac{d^{(6)}\theta }{d{t}^6}=\sin \theta \mathrm{}\left(3-11w^2\cos \theta -w^4-4{\mathrm{}\cos }^2\theta \right)$$

and

(26) $$\frac{d^{(6)}\theta }{d{t}^6}=\sin \mathrm{}\theta \mathrm{}\left(3-30\cos \theta (\epsilon +\cos \theta )-4{\epsilon }^2\right)$$

Therefore, one has the upper bound on the error, and from equation (24), we get an upper bound on h for a given tolerance tol. Let

(27) $${\theta }_{\text{max}}^{(6)}\equiv max\mathrm{}\left|\frac{d^{(6)}\theta }{d{t}^6}\right|,$$

from which one has the upper bound on the error and from equation (24). Thus, we get an upper bound on the step $\widehat{h}$ for a given tolerance tol:

(28) $${\widehat{h}}^6\le \frac{6!}{{\theta }_{\text{max}}^{(6)}}\text{tol}$$

For the same run above, with tol= 10⁻¹⁴, we get $\widehat{h}\le \text{8}{\text{.27 10}}^{-\text{3}}$, which means that RK5 would require a bit over one-third of the number of steps of RK4. Nevertheless, RK5 requires at least six stages, while RK4 requires four. So, RK4 requires about one-third less computer time than RK5; therefore, for long runs, there is not a big difference between these two methods.

Splitting the second-order equation (3) into two first-order equations was not by chance. They are the equations one gets from Hamiltonian formalism instead of Lagrangian formalism.

The following section will deal with a double pendulum system, starting with the Hamiltonian formalism.

4. The Double Pendulum

4.1. A brief overview

A double pendulum is a system that consists of two pendulums attached end-to-end. The coordinates of each pendulum can be defined in terms of the angles θ₁ and θ₂, and their corresponding Cartesian coordinates x₁, z₁, x₂, and z₂.

(29) $$x_1=L\text{sin}\mathrm{}{\theta }_{\text{1}}\text{,}{z}_{\text{1}}\text{=}-L\text{cos}\mathrm{}{\theta }_{\text{1}}$$

(30) $$x_2=x_{\text{1}}+L\text{sin}\mathrm{}{\theta }_{\text{2}}\text{,}{z}_{\text{2}}\text{=}{z}_{\text{1}}-L\text{cos}\mathrm{}{\theta }_{\text{2}}$$

See Figure 5. One can derive the coupled second-order differential equations that describe the system from Newton’s laws or the Lagrangian formalism. The Lagrangian is the difference between the potential energy U and the kinetic energy T [³]:

(31) $$\mathcal{L}=T-U=\frac{m}{2}\left({\left(\frac{d{z}_1}{dt}\right)}^2+{\left(\frac{d{z}_2}{dt}\right)}^2\right)-mg\mathrm{}\left(z_1+z_2\right),$$

Figure 5
A double pendulum with equal lengths L and masses m.

which, in the variables θ_i, using the notation $\stackrel{.}{()}\equiv d\text{()/}dt$, is

(32) $$\begin{array}{l}\mathcal{L}\text{=}\frac{m{L}^2}{2}\left(2{\dot{\theta }}_1{}^2\text{+ 2}{\dot{\theta }}_1{\dot{\theta }}_2\text{cos}\left({\theta }_1-{\theta }_2\right)\right)\\ \text{+}mgL\mathrm{}\text{(2 cos}\mathrm{}{\theta }_1\text{+ cos}\mathrm{}{\theta }_2)\end{array}$$

One can rescale the Lagrangian by mgL and use the time variable t in units of $\sqrt{g/L}$, as we did in the simple pendulum case, with no lack of generality. In other words, we use time, length, and energy units such as m = L = g = 1. Let us call them math units.

From the coordinates θ_i one gets their conjugate momenta

(33) $$p_i\equiv p_{{\theta }_i}=\frac{\partial \mathcal{L}}{\partial \dot{{\theta }_i}}$$

namely

(34) $$p_1\text{= 2}{\dot{\theta }}_1\text{+}{\dot{\theta }}_2\text{cos (}{\theta }_1-{\theta }_2\text{)}$$

(35) $$p_2\text{=}{\dot{\theta }}_2\text{+}{\dot{\theta }}_1\text{cos (}{\theta }_1-{\theta }_2\text{)}$$

Let us set q_i ≡ θ_i to ease the notation. The Hamiltonian is obtained as $H\text{=}{{\displaystyle \sum }}_i{p}_i{\dot{q}}_i-\mathcal{L}$, with the angular velocities ${\dot{\theta }}_i$ written in terms of the conjugate momenta p_i from equations (28, 28). After a few algebraic and trigonometric transformations, we find the Hamiltonian

(36) $$H=\frac{p_1^2+2p_2^2-2p_1{p}_2\cos \left(q_1-q_2\right)}{2\mathrm{}\left(1+{\sin }^2\mathrm{}\left(q_1-q_2\right)\right)}-2\cos q_1-\cos q_2$$

From this, one gets Hamilton’s equations, which describe the system’s motion in phase space [³].

$$\frac{d{q}_i}{dt}\text{=}\frac{\partial H}{\partial p_i},\frac{d{p}_i}{dt}\text{=}-\frac{\partial H}{\partial q_i},\frac{dH}{dt}\text{=}\frac{\partial H}{\partial t}$$

Since H does not depend explicitly on t, one concludes from the last equation above that the Hamiltonian (the total energy) is constant (conserved) along the motion. Let us define the auxiliary functions:

(37) $$A_1\text{≡}\frac{p_1{p}_2{\text{sin(q}}_{\text{1}}-{\text{q}}_{\text{2}}\text{)}}{1+{\sin }^2{\mathrm{}\text{(q}}_{\text{1}}-{\text{q}}_{\text{2}}\text{)}}$$

(38) $$A_2\text{≡}\frac{\begin{array}{l}(p_1{}^2\text{+ 2}{p}_2{}^{\text{2}}-\text{2}{p}_{\text{1}}{p}_{\text{2}}\text{cos(}{q}_{\text{1}}-q_{\text{2}}\text{))}\\ \text{sin(}{q}_{\text{1}}-q_{\text{2}}\text{) cos(}{q}_{\text{1}}-q_{\text{2}}\text{)}\end{array}}{{(1+{\sin }^2{\mathrm{}\text{(q}}_{\text{1}}-{\text{q}}_{\text{2}}\text{))}}^{\text{2}}}$$

Then, the equations of motion are:

(39) $$\frac{d{q}_1}{dt}\text{=}\frac{p_1-p_{\text{2}}\cos (q_1-q_2)}{2-{\cos }^2\mathrm{}(q_1-q_2)}$$

(40) $$\frac{d{q}_2}{dt}\text{=}\frac{2p_2-p_{\text{1}}\cos (q_1-q_2)}{2-{\cos }^2\mathrm{}(q_1-q_2)}$$

(41) $$\frac{d{p}_1}{dt}\text{=}-A_{\text{1}}\text{+}{A}_{\text{2}}-\text{2sin}{q}_{\text{1}}$$

(42) $$\frac{d{p}_2}{dt}\text{= +}{A}_{\text{1}}-A_{\text{2}}-\text{sin}{q}_2$$

The motion of the system in phase space is the curve, the orbit of the state y; the solution of the following equations:

(43) $$y(t)=\begin{array}{l}\left(\begin{array}{l}{q}_1\\ q_2\\ p_1\\ p_2\end{array}\right)\hfill \end{array},\frac{dy}{dt}=f(y),f(y)=\begin{array}{l}\hfill \left(\begin{array}{c}\frac{\partial H}{\partial p_1}\\ \frac{\partial H}{\partial p_2}\\ -\frac{\partial H}{\partial q_1}\\ -\frac{\partial H}{\partial q_2}\end{array}\right)\hfill \end{array}$$

y ∈ ℝ⁴, which we call the phase space where the curve y(t) is the orbit of the motion in the phase space.

4.2. Small perturbations from stable equilibrium points

The equilibrium points of (43) are the fixed points denoted by q₁ = jπ, q₂ = kπ with p₁ = 0 = p₂, ∀j, k ∈ ℤ. All the points are unstable except for the case of q₁ = q₂ = p₁ = p₂ = 0.

To find the normal modes of the system, we can study small perturbations (SP) around the stable equilibrium [¹²], the fixed point (q₁, q₂, p₁, p₂) = (0, 0, 0, 0).

By keeping only the first order in q_i, p_i, we get the following equations:

(44) $$\frac{d{q}_1}{dt}\approx p_{\text{1}}-p_{\text{2}}\text{,}\frac{d{q}_1}{dt}\approx -2q_1,$$

(45) $$\frac{d{q}_2}{dt}\approx \text{2}{p}_{\text{2}}-p_{\text{1}}\text{,}\frac{d{p}_2}{dt}\approx -q_2.$$

We can then use the eigenvalue, eigenvector SP equations

(46) $$\begin{array}{l}\frac{dy}{dt}\text{=}\text{λ}y\iff \left(\begin{array}{llll}\hfill 0\hfill & 0\hfill & \hfill 1\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill \\ \hfill -2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right)\left(\begin{array}{c}\hfill q_1\hfill \\ \hfill q_2\hfill \\ \hfill p_1\hfill \\ \hfill p_2\hfill \end{array}\right)\\ \text{=}\text{λ}\left(\begin{array}{c}\hfill q_1\hfill \\ \hfill q_2\hfill \\ \hfill p_1\hfill \\ \hfill p_2\hfill \end{array}\right)\end{array}$$

to find the eigenvalues and eigenvectors.

The eigenvalue equation is given by:

(47) $${\text{λ}}^4\text{+ 4}{\text{λ}}^2\text{+ 2 = 0}\Rightarrow {\text{λ}}^2\text{=}-\left(2\pm \sqrt{2}\right)$$

which yields four imaginary solutions λ = ±iw±. Therefore, the SP states around the rest point are

(48) $$y_{sp}\propto e^{\pm iw\pm t},$$

which are trigonometric functions sin and cos of two fundamental frequencies, w₊ and w₋ given [¹⁸] by:

(49) $$w_{+}\text{=}\sqrt{2+\sqrt{2}}\approx 1.84776,\mathrm{}{\text{w}}_{-}\text{=}\sqrt{2-\sqrt{2}}\approx 0.76537$$

For comparison, recall that SP around the rest stable point for a single pendulum has the equations $\dot{\theta }\text{=}w$, $\dot{w}\text{=}-\theta$, which has the imaginary eigenvalues ±i with the fundamental mode’s frequency w₀ = 1 and period T₀ = 2π, in the units we are employing. See (3) and (14) for small θ. Accordingly, the general solution in this case is θ = a cos t + b sin t, for a, b constants.

One can find the eigenvectors of (46) to get the most general SP solution for the double pendulum:

(50) $$\begin{array}{l}{y}_{sp}\text{=}{C}_{+}\left(\begin{array}{c}\hfill 1\hfill \\ \hfill -\sqrt{2}\hfill \\ \hfill i{w}_{+}(2-\sqrt{2})\hfill \\ \hfill -i{w}_{+}(\sqrt{2}-1)\hfill \end{array}\right)e^{i{w}_{+}t}\\ \text{+}{C}_{-}\left(\begin{array}{c}\hfill 1\hfill \\ \hfill \sqrt{2}\hfill \\ \hfill i{w}_{-}(2+\sqrt{2})\hfill \\ \hfill i{w}_{-}(\sqrt{2}+1)\hfill \end{array}\right)e^{i{w}_{-}t}+\text{compl}\text{. conj}\text{.}\end{array}$$

where C_± are complex constants, we added the complex conjugate partner to get real solutions for y_sp.

The general solution for the double pendulum is more complicated and nonlinear. However, we can visualize an explicit solution for each fundamental mode. Note that the motion is harmonic, and given the frequency, one computes the period of motion, namely

(51) $$T_{+}\text{=}\frac{2\pi }{w_{+}}\approx 3.40044,T_{-}\text{=}\frac{2\pi }{w_{-}}\approx 8.20938$$

These values are helpful for numerical evolution and plotting to decide how long we should run a motion and the limits of plotting.

One solution for the plus mode is

(52) $$y_{1+}\text{=}{\theta }_{\text{0}}\mathrm{}\left[\begin{array}{l}\text{cos(}{w}_{\text{+}}t\text{)}\hfill \\ -\sqrt{2}\text{cos(}{w}_{\text{+}}t\text{)}\hfill \\ -(2-\sqrt{2})w_{\text{+}}\text{sin(}{w}_{\text{+}}t\text{)}\hfill \\ w_{+}\text{(}\sqrt{2}-1\text{) sin(}{w}_{\text{+}}t\text{)}\hfill \end{array}\right]$$

which is a configuration starting at t = 0 at rest, with θ₀ and $-\sqrt{2}{\theta }_0$ for each pendulum.

One solution for the minus mode is

(53) $$y_{1-}\text{=}{\theta }_{\text{0}}\mathrm{}\left[\begin{array}{l}\text{cos(}{w}_{-}t\text{)}\hfill \\ \sqrt{2}\text{cos(}{w}_{-}t\text{)}\hfill \\ -(2+\sqrt{2})w_{-}\text{sin(}{w}_{-}t\text{)}\hfill \\ -\text{(1+}\sqrt{2}\text{)}{w}_{-}\text{sin(}{w}_{-}t\text{)}\hfill \end{array}\right]$$

which is a configuration starting at t = 0 at rest, with θ₀ and $\sqrt{2}{\theta }_0$ for each pendulum. See Figures 6 and 7. The simulated amplitudes are exaggerated for visual purposes – they break the small perturbations approximation, but the relative proportionalities are respected, namely ${\theta }_2\text{=}\mathrm{\pm }\sqrt{2}{\theta }_1$ for the plus or minus modes.

Figure 6
Plus SP mode with frequency w₊ ≈ 1.848 and ${\theta }_2\text{=}-\sqrt{2}{\theta }_1$.

Figure 7
Minus SP mode with frequency w₋ ≈ 0.765 and ${\theta }_2\text{= +}\sqrt{2}{\theta }_1$.

The SP motions are interesting, but there are no nonlinear behaviors. So, let us work with the whole double pendulum equations. Of course, one has to request numerical computations because analytical solutions are elusive.

4.3. Numerical solutions of the double pendulum system

Based on the simple pendulum discussion, we will get numerical solutions using the 4th-order Runge-Kutta Method. See the code in Appendix 2. See Figure 8 for the case for q₁(0) = −0.8224, q₂(0) = 1.4335, p₁(0) = 1.5422, p₂(0) = 0, ⇒ ε ≈ −0.7552 which looks essentially periodic. The Figure 9 display the motion of the pendula in the Cartesian plane x − z.

Figure 8
The double pendulum evolution for q₁(0) = −0.8224, q₂(0) = 1.4335, p₁(0) = 1.5414, p₂(0) = 0, ⇒ ε ≈ −0.7552 is essentially periodic.

Figure 9
The double pendulum motion for q₁(0) = −0.8224, q₂(0) = 1.4335, p₁(0) = 1.5414, p₂(0) = 0 is essentially speriodic. We show the pendula bars’ initial states and the masses’ curves.

The motion is quite different for another set of initial conditions. See Figure 10 for q₁(0) = 0.5, q₂(0) = π, p₁(0) = 0, p₂(0) = 0, which has the same energy ε ≈ −0.7552 but it is not periodic – it seems chaotic. The Figure 11 display the motion of the pendula in the Cartesian plane x − z.

Figure 10
The double pendulum evolution for q₁(0) = 0.5, q₂(0) = π, p₁(0) = 0 = p₂(0), ⇒ ε ≈ −0.7552 seems chaotic.

Figure 11
The double pendulum motion for q₁(0) = 0.5, q1₂(0) = π, p₁(0) = 0 = p₂(0), seems chaotic. We show the pendula bars’ initial state and the masses’ curves on the plane.

We aim to find what is called Hamiltonian deterministic chaos in the motion of the double pendulum. Moreover, the Hamiltonian, the conserved energy, is a diagnosis of the accuracy of the numerical method employed. Since we want to distinguish whether a motion is quasi-periodic or chaotic, we “must sharpen our observational tools” [⁴, ¹³].

The double pendulum have four canonical variables, we called q₁, q₂, p₁, p₁, $\in {\mathcal{S}}^4\subset {ℝ}^4$, the 4-dimensional phase space. Nevertheless, since the Hamiltonian (36) is constant for a given initial condition, the variables are confined to a 3-dimensional hypersurface of ${\mathcal{S}}^4$. If another time invariant existed besides the Hamiltonian, the motion would be confined to a 2D surface of the phase space. No other invariant for the generic double pendulum has been found so far.

It is difficult to see the generic orbits in the 4D phase space on a 2D display paper or computer monitor. Remember that, for the simple pendulum, the phase space is a plane, so for each energy value, the orbits are just 1D curves on the phase plane. See Figure 2 – the orbits do not cross each other. The orbits are constrained to the constant energy (Hamiltonian) 3D manifold for the double pendulum motion.

4.4. Poincaré’s maps

To further study the double pendulum motion, we can take advantage of Poincaré maps. These maps offer a sectional view of orbits in the Phase Space, providing an insightful graphical representation to understand the system’s behavior. For instance, we can plot the outer mass coordinates x₂ and z₂, at moments when the momentum p₂ is zero and increasing. This is achieved by plotting the pair (x(t), z(t)) at the instants t when p₂(t) ≤ 0 and p₂(t + h) ＞ 0 for a given time-step h.

In these maps, nonperiodic orbits densely populate the constant energy 3D sub-manifold of the Phase Space, tracing paths that never intersect. They fully occupy the chosen cross-section. In contrast, periodic orbits form closed loops and manifest as either distinct segments or isolated points.

Creating detailed Poincaré maps necessitates extensive system evolution, which can be computationally intensive. For instance, our Python simulations for Figures 10 and 11, running up to t_f = 7000, required approximately 542 seconds. To mitigate this computational expense, we switched to Fortran language and its compiler for longer runs, experiencing a performance improvement of nearly a hundredfold over Python. See the code in the Appendix 3. Additionally, it’s crucial to monitor memory usage, as data storage can become a significant challenge during extended simulations.

We executed simulations for the initial conditions listed in Table 1, extending them to t_f = 70, 000 to generate their respective Poincaré maps. Each set of initial conditions underwent 61 million time steps. We constrained the system to a 3D hypersurface in the Phase Space with an energy level of ε ≈ −0.755165, maintaining an energy precision of 10⁻¹⁰.

Table 1
Initial Conditions for the Poincaré maps, in which p₂ = 0 and p₁ is such that ε ≈ −0.755165.

IC label	q 1	q 2	p 1	type
1	0.5	π	0	chaotic
2	0.5236	2.618	−0.6229	chaotic
3	0.5236	−2.618	0.4709	chaotic
4	−0.5	1.4	2.106	2-periodic
5	−0.822	1.4335	1.5414	periodic
6	−0.65	1.4	1.8973	i-periodic
7	−0.7	1.4	1.8156	i-periodic
8	−0.85	1.4	1.536	i-periodic
9	0.9	1.4	−1.2722	chaotic

Figure 12 showcases the superimposed Poincaré maps for these initial conditions. Among the nine sets, four exhibit chaotic behavior, as evidenced by the scattered distribution of points across the map.

Figure 12
The Poincaré Maps for nine initial conditions of Table 1 with ε ≈ −0.7552. The origin (0, 0) is plotted to make the reference easy. Each dot in the figure represents the position of the mass 2, such that p₂ is zero and increasing.

One periodic case is represented by a single point encircled by three closed loops, each indicative of periodic motion with quasi-periodic incommensurate frequencies (i-periodic). Furthermore, there are pairs of small, closed curves symbolizing a 2-periodic orbit, where successive crossings alternate between two distinct paths.

While Poincaré maps are instrumental in illustrating the chaotic characteristics of a dynamical system, they are just one tool in a broader analytical arsenal. Developing and employing various analytical methods to deepen our understanding of such complex systems is essential. For that purpose, we will present [¹³] the Lyapunov characteristic exponents next.

4.5. Lyapunov characteristic exponents

Let us study how sensitive is an orbit to the small change in the initial conditions.

Lyapunov Characteristic Exponents (LCEs) are fundamental in understanding the behavior of dynamical systems. They provide a quantitative measure of the system’s sensitivity to initial conditions, a key aspect in the study of chaos. By analyzing the divergence or convergence of trajectories in phase space, LCEs allow us to classify the stability and predictability of orbits in these systems.

Let y(t) and y(t) + ξ(t) be two nearby orbits at some instant t. So, from the equation (43) we get the time evolution of ξ:

(54) $$\frac{d\xi }{dt}\text{=}f\text{(}y\text{+}\xi \text{)}-f\text{(}y\text{+}\xi \text{)}\approx \frac{\partial f}{\partial y}\xi +\text{O (||}\xi {\text{||}}^2\text{)}\text{.}$$

We remember that y and y+ξ are in the 4D phase space. $\frac{\partial f}{\partial y}$ is well defined because we assume the vector field f is Liphictizian in y. Thus $\frac{\partial f}{\partial y}$ which is evaluated the fidutial orbit y(t), i.e. at ξ = 0, is a 4×4 matriz. Therefore, one can compute the eigenvalues of the matrix equation

(55) $$\frac{\partial f}{\partial y}\xi \text{=}\text{λ}\xi \Rightarrow \xi \approx \exp (\text{λ}t)$$

The linear approximation of the dynamics, as shown in the equation, is valid for small perturbations ξ. This linearization is fundamental in calculating the LCEs but also limits the analysis to local behavior around the trajectory.

The eigenvalues of the Jacobian matrix ∂f/∂y provide critical insights into the local dynamics of the system. A complex pair of eigenvalues, for instance, indicates oscillatory behavior, while the sign of the real parts of the eigenvalues reveals whether the system experiences divergence (positive sign), a hallmark of chaos, or convergence (negative sign) in the vicinity of the trajectory. This helps us understand the local stability and predictability of the system.

The matrix $\frac{\partial f}{\partial y}$ for a conservative Hamiltonian system possesses a symplectic structure.We see in equation (43) J has four blocks of 2 × 2 sub-matrices, namely

(56) $$\frac{\partial f}{\partial y}\text{=}\left(\begin{array}{cc}\hfill \frac{{\partial }^{2}H}{\partial q\partial p}\hfill & \hfill \frac{{\partial }^{2}H}{{\partial }^{2}p}\hfill \\ \hfill -\frac{{\partial }^{2}H}{{\partial }^{2}q}\hfill & -\frac{{\partial }^{2}H}{\partial p\partial q}\hfill \end{array}\right)$$

which imposes specific constraints on its eigenvalues. Notably, for every eigenvalue, its negative is also an eigenvalue. This reflects the conservation of phase space volume in Hamiltonian dynamics, as per Liouville’s theorem, and directly influences the computation and interpretation of the LCEs.

A chaotic system exhibits sensitive dependence on initial conditions, by definition, meaning that minor differences in initial conditions lead to significantly different outcomes. The positive exponent’s magnitude (s) measures the rate at which predictability is lost.

Let us discuss the numerical studies to compute LCE.

4.6. Numerical computations of LCE

Numerically computing Lyapunov Characteristic Exponents (LCEs) can be challenging due to the handling of infinities and infinitesimals in digital computations. Practical methods, like the QR decomposition technique, are employed to approximate these exponents. However, it is important to acknowledge the limitations of these numerical methods, especially with regards to their sensitivity to the choice of time steps, duration of orbit integration and initial conditions. It is essential to note that LCEs can vary significantly across different regions of the phase space.

The complete evolution of ξ may not be just exp(λt) of equation (55) because $\frac{\partial f}{\partial y}$ is itself a time depend matrix.

Furthermore, we assumed a O(||ξ||²) approximation in the equation (54). From the “distance” d(t) ≡ ||ξ|| between the orbits, Lyapunov defined his exponents as the asymptotic growth rate of d(t), the mean rate of exponential divergence [⁷]:

(57) $$d(t)\approx e^{\sigma t}{d}_0\Rightarrow \sigma \equiv \underset{t\to \infty }{\lim }\left(\frac{1}{t}\right)\underset{d(0)\to 0}{\lim }\text{In}\left(\frac{d(t)}{d(0)}\right)$$

If an orbit is constrained to a compact set in S⁴ and f(y) is differentiable and, the Jacobian ∂f/∂y is uniformly bounded [⁷]. Therefore, the growth of any vector norm is at most exponential.

Of course, the definition (57) can not be implemented numerically on our computers. Infinity and infinitesimal are not numbers. One can show that for

(58) $$\frac{\partial f}{\partial y}\text{bounded}\Rightarrow \left(\frac{1}{t}\right)\underset{d(0)\to 0}{\lim }\text{In}\left(\frac{d(t)}{d(0)}\right)\text{bounded,}$$

both above and below for t ＞ 0. Since any bounded sequence has limit points, we may define the LCE spectrum as the set of limit points [⁷]:

(59) $$LCE(y,\xi )\text{=}\left\{\sigma \text{=}\underset{j\to \infty }{\lim }\left(\frac{1}{t_j}\right)\underset{d(0)\to 0}{\lim }\text{In}\left(\frac{d(t_j)}{d(0)}\right)\text{for some sequence of}{t}_j\right\}$$

It is important to estimate how long we need to evolve the system and what sequence of t_j we could take. We can calculate the eigenvalues from (55) to get some guidance. For example, for our first IC of Table 1, the four eigenvalues are:

(60) $$\text{λ}\text{=}\left[\begin{array}{c}\hfill -\text{1}\text{.1394}\hfill \\ \hfill \text{1}\text{.0484}i\hfill \\ \hfill -\text{1}\text{.0484}i\hfill \\ \hfill \text{1}\text{.1394}\hfill \end{array}\right],$$

which indicates two directions in which the orbits are compelled to oscillate, one direction to diverge and another to converge. The sum of these eigenvalues is zero.

Therefore, our run must last at least δt = 1/1.1394 ≈ 0.87 to feel an exponential divergence in one direction. We can also compute the average time between a hit and another at the Poincaré section. For this IC, we ran the orbit for t_f = 70,000 with 8692 hits, so the average time is Δt ≈ 8.1.

We can use the fact that the orbits are confined in the energy constant 3D hypersurface, which constrains the possible values of the canonical variables. In particular, the momenta p₁ and p₂ are bounded by an ellipse in the plane p₁ × p₂. In other words, the distance ||ξ|| cannot grow without limits. If we start two nearby orbits with d₀ = 10⁻⁹, the distance may increase by ten orders of magnitude. We estimate the maximum distance to be of magnitude 10 for this IC. We get this estimate by paying attention to the possible values of p₁, p₂, and q₁:

(61) $$p_1^2\le 4\left(3+\epsilon \right),p_2^2\le 2\left(3+\epsilon \right),-\frac{1}{2}\left(1+\epsilon \right)\le \cos q_1\le 1$$

and based on our numerical evolution |q₂| ≤ 5π. Therefore, we have an estimate of the norm

(62) $$\left|\right|y|{|}^2=q_1^2\text{+}{q}_2^2\text{+}{p}_1^2+p_2^2\le \text{263}\text{.1}\Rightarrow \text{||}y\text{||}\le \text{16}\text{.2}$$

Thus, if we start two nearby orbits with d₀ = 10⁻⁹, the distance may increase by ten orders of magnitude. Therefore, if one uses only the largest eigenvalue, it would be necessary to run the system by

(63) $$e^{\text{λ}{T}_{max}}{\text{= 10}}^{\text{10}}\Rightarrow T_{max}{\text{= log (10}}^{\text{10}}\text{) /1}\text{.1394}\approx \text{20}\text{.2}$$

We gather together, from the discussion above, three relevant time scales one has to pay attention to compute the Lyapunov exponents for the fiducial orbit defined by the IC-1:

(64) $$\delta t\approx 0.87,\Delta t\approx \text{8}\text{.1,}{T}_{max}\approx 20$$

We say [⁵] that ξ, as a vector in the 4D vector space, is a linear combination of the basis vectors ${\widehat{e}}_i\mathrm{}(i=1,\mathrm{.\; .\; .},4)$ [¹⁹]. The stretching, contracting, or no change in each direction ${\widehat{e}}_i$ gives us the spectra of the Lyapunov characteristics exponents (LCE) σ_i (i = 1, . . . , 4). For the 4D phase space, there will be four such quantities that can be ordered as

(65) $${\sigma }_1\ge {\sigma }_2\ge {\sigma }_3\ge {\sigma }_4$$

For our conservative Hamiltonian autonomous 4D dynamical system, there is a special symmetry, namely

(66) $${\sigma }_1=-{\sigma }_4,{\sigma }_2=-{\sigma }_3$$

Thus, a stretching in one direction is compensated by a contraction in the other, as it should be from Liouville’s theorem. In other words, if at the moment t we have a small 4D spherical ball around the point y(t) of a fiducial orbit, the ball will be transformed into a scratched 4D ellipsoidal. This result is also useful for checking our numerical solutions for the LCE.

The computation of the complete spectrum of LCE requires tools [²⁰, ²¹] from linear algebra and vector space.

The spectrum of Lyapunov exponents generally depends on the fiducial or referential orbit. The Lyapunov exponents describe the behavior of vectors ξ, [²²] which formally are in the tangent space of the phase space and are defined from the Jacobian matrix

(67) $$J_{ij}(t)\text{=}{\left.\frac{\partial f_i(y)}{\partial y_j}\right|}_{y(t)}$$

this Jacobian defines the evolution of the tangent vectors, given by the matrix Ξ, via the equations

(68) $$\frac{d\text{Ξ}}{dt}\text{=}J\text{Ξ}$$

(69) $$\text{Ξ}(0)\text{= Identit}y$$

There are 16 equations in (68). The time derivative d/dt acts on each component of the matrix Ξ, which describes how a small change around the point y(0) propagates to around the final point y(t).

Consider the fundamental matrix

(70) $$\text{Ξ}(t)\text{= (}{\xi }_1\text{(}t\text{),}.\; .\; .,{\xi }_4\text{(}t\text{)),}$$

in which each column consists of the linearly independent solutions ξ_i(t) of linearized system with respect to the fiducial orbit y(). Thus, Ξ(t) is not singular. “For time-varying linearization of nonlinear systems, Lyapunov introduced the so-called Lyapunov characteristic exponents (LCEs) as the upper bounds of the exponential growth rates of solutions” as: [¹¹, ²²]

(71) $$LC{E}_i(y_0)\text{=}\underset{t\to \infty }{\lim }\text{sup}\frac{1}{t}\log \text{||}{\xi }_i(t)|,i\text{= 1,}.\; .\; .\text{, 4}\text{.}$$

The details of computing the Lyapunov exponents go beyond the scope of this review. Let us point out that the chosen method to compute the Lyapunov exponent requires solving (16) together with (68). For that purpose, we “align” 20 functions, evolving the system with our standard RK4 method (or equivalent) and factoring the matrix Ξ by the famous QR decomposition. We use Python’s libraries to get the code in the Appendix 3. We used Δt = 2³ and built a sequence with iter = 2¹⁴ to obtain:

(72) $${\text{IC}}_1:{\sigma }_{\text{1}}\text{= 0}\text{.143,}{\sigma }_{\text{2}}\text{= 5}\times {\text{10}}^{-\text{5}},{\displaystyle \sum }_i{\sigma }_i\text{=}-\text{5}\times {\text{10}}^{-\text{9}}.$$

See (66). This computation took about 37 minutes of CPU time.

We similarly run for the other IC, for example:

(73) $${\text{IC}}_5:{\sigma }_{\text{1}}\text{= 9}\times {\text{10}}^{-\text{5}}\text{,}{\sigma }_{\text{2}}\text{= 2}\times {\text{10}}^{-\text{5}},{\displaystyle \sum _i{\sigma }_i\text{= 7}}\times {\text{10}}^{-\text{9}}.$$

The computed LCEs for various initial conditions reveal the nuanced nature of chaotic behavior in dynamical systems. For instance, a significantly positive Lyapunov exponent, as observed for IC1, indicates chaotic behavior, where slight differences in initial conditions lead to vastly different trajectories over time. Conversely, smaller LCE values, like those computed for IC5, suggest more predictable and less sensitive dynamics. These distinctions are crucial in characterizing the nature of orbits in the phase space and understanding the system’s overall behavior.

5. Educational Implications

The double pendulum is an exemplary teaching aid for educators in STEM (Science, Technology, Engineering, and Mathematics) disciplines, offering a rich blend of theory and practical application. Its complex dynamics provide an ideal platform for introducing students to advanced physics and engineering concepts while honing their computational skills, crucial for contemporary scientific inquiry.

We advocate using online markdown notebooks to facilitate an interactive and immersive learning experience. These digital notebooks eliminate the need for downloads or installations, streamlining the learning process. They are particularly advantageous as they enable students and educators to engage in a dynamic, collaborative educational environment. Integrating computations, narrative text, vivid imagery, and elegantly formatted equations using LaTeX, all within a single, accessible document, makes these notebooks an invaluable resource.

Moreover, these notebooks offer the flexibility to tailor educational content. Educators can explore beyond the provided material and delve into more complex scenarios. For instance, they can delve into the system of double-pendulum with unequal masses and lengths, add a damping acceleration proportional to each angular velocity, consider an extra acceleration at the connecting pendula, impose physical constraints on the amplitude of some or both libration angles (simulate an articulate leg or arm, for example) or allow motion beyond the planar x − z. The possibilities are endless. This adaptability enhances the learning experience and encourages students to engage in creative problem-solving and critical thinking.

Incorporating these methodologies in the classroom enriches the curriculum. Educators can significantly elevate the learning experience through such innovative teaching approaches, making it more engaging, comprehensive, and reflective of real-world scientific practices, preparing students for the rigorous and diverse demands of modern scientific research, and fostering a deeper understanding of the underlying principles of STEM fields.

6. Conclusion

The double pendulum is a fascinating system that offers profound insights into the nature of deterministic chaos. Its study provides a comprehensive learning experience, encompassing Classical Mechanics, Hamiltonian formalism, and Numerical Methods for solving ODEs. With the computational tools available today, educators have an unprecedented opportunity to provide a holistic educational experience that prepares students for the challenges of contemporary research in physics and mathematics.

By delving into the intricacies of the double pendulum, students and educators alike can explore a wide array of scientific phenomena, from the deterministic yet unpredictable nature of chaotic systems to the practical applications of numerical methods. Therefore, the double pendulum serves as an invaluable resource for pedagogical and research-oriented endeavors in physics and mathematics.

Acknowledgments

This review was inspired by the stimulating magazine Computers in Physics article [⁴] and the compiled reprint selection book Hamiltonian Dynamical Systems [²³, ²⁵].

Visual Studio Code [²⁶] was used as the framework for editing and running serial (no parallel) codes. For grammar and style, Grammarly [²⁷] was used, and for identifying interesting articles, scite.ai [²⁸] platform and Google Scholar [²⁹] were employed. Wikipedia [³⁰] has been consulted for various entries. Finally, Zotero [³¹] and Google Scholar [²⁹] played an instrumental role in building the references.

Supplementary Material

The following online material is available for this article: Appendices 1, 2 and 3.

References

• ^[1]
  M.R. Matthews, Time for Science Education: How Teaching the History and Philosophy of Pendulum Motion Can Contribute to Science Literacy (Springer Science & Business Media, Berlin, 2000).
• ^[2]
  K.T. Alligood, T.D. Sauer and J.A. Yorke, Chaos: an introduction to dynamical systems (Springer, New York, 2010).
• ^[3]
  H. Goldstein, Classical Mechanics (Addison-Wesley, Boston, 1980), 2 ed.
• ^[4]
  N. Srivastava, C. Kaufman and G. Müller, Computer in Physics 4, 549 (1990).
• ^[5]
  M. Tabor, Chaos and Integrability in Nonlinear Dynamics: An Introduction (Wiley-Interscience, New Jersey, 1989).
• ^[6]
  J.L. McCauley, Chaos, dynamics, and fractals: an algorithmic approach to deterministic chaos (Cambridge University Press, Cambridge, 1994), v. 2.
• ^[7]
  J.D. Meiss, Differential dynamical systems (SIAM, Pennsylvania, 2007).
• ^[8]
  W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge, 2007), 3 ed.
• ^[9]
  A.M. De Almeida, Hamiltonian systems: Chaos and quantization (Cambridge University Press, Cambridge, 1988).
• ^[10]
  A.M.O. Almeida, Sistemas Hamiltonianos: caos e quantização (Editora da UNICAMP, São Paulo, 1995), 3 ed.
• ^[11]
  N.V. Kuznetsov, T.A. Alexeeva and G.A. Leonov, Nonlinear Dynamics 85, 195 (2016).
• ^[12]
  G.A. Monerat, E.V.C. Silva, G. Oliveira-Neto, A.R.P. Assumpção and A.R.R. Papa, Revista Brasileira de Ensino de Física 28, 177 (2006).
• ^[13]
  A.H. Nayfeh and B. Balachandran, Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods, (Wiley, New Jersey, 1995).
• ^[14]
  J. Palis Jr. and W. Melo, Introdução aos sistemas dinâmicos (IMPA, Rio de Janeiro, 2017).
• ^[15]
  A.C. King, J. Billingham and S.R. Otto, Differential Equations: Linear, Nonlinear, Ordinary, Partial (Cambridge University Press, Cambridge, 2003).
• ^[16]
  W.E. Boyce and R.C. DiPrima, Elementary Differential Equations and Boundary Value Problems (John Wiley & Sons, New York, 1997).
• ^[17]
  J.C. Butcher, in: Numerical Methods for Ordinary Differential Equations, edited by J.C. Butcher (John Wiley & Sons, New Jersey, 2016).
• ^[18]
  H.J. Korsch and H.J. Jodl, in: Chaos: A Program Collection for the PC, edited by H.J. Korsch and H.J. Jodl (Springer Berlin Heidelberg, Berlin, 1994).
• ^[19]
  G.B. Arfken, H.J. Weber and F.E. Harris, Mathematical Methods for Physicists, (Academic Press, London, 2013), 7 ed.
• ^[20]
  G. Benettin, L. Galgani, A. Giorgilli and J.M. Strelcyn, Meccanica 15, 9 (1980).
• ^[21]
  G. Benettin, L. Galgani, A. Giorgilli and J.M. Strelcyn, Meccanica 15, 21 (1980).
• ^[22]
  N.V. Kuznetsov, Physics Letters A 380, 2142 (2016).
• ^[23]
  R.S MacKay and J.D Meiss, Hamiltonian Dynamical Systems: a reprint selection, (CRC Press, Florida, 2020).
• ^[24]
  M.V. Berry, AIP Conference Proceedings 46, 16 (1978).
• ^[25]
  M.V. Berry, in: Hamiltonian Dynamical Systems: a reprint selection edited by R.S. MacKay and J.D. Meiss (Adam Hilger, Bristol, 1987).
• ^[26]
  Visual Studio Code – Code Editing, available in: https://code.visualstudio.com/
  » https://code.visualstudio.com/
• ^[27]
  Grammarly, available in: https://app.grammarly.com/
  » https://app.grammarly.com/
• ^[28]
  scite, available in: https://scite.ai
  » https://scite.ai
• ^[29]
  Google Scholar, available in: https://scholar.google.com/
  » https://scholar.google.com/
• ^[30]
  Wikipedia, available in: https://www.wikipedia.org/
  » https://www.wikipedia.org/
• ^[31]
  Zotero – Your personal research assistant, available in: https://www.zotero.org/
  » https://www.zotero.org/

• 1
  We utilized Python 3.9.17 [GCC 11.2.0] Anaconda, Inc. on linux Ubuntu 22.04.3 LTS, 64 bits. Our processors were Intel Xeon(R) CPU E3-1225 v5 @ 3.30GHz x 4 with 8 GB memory.

Publication Dates

• Publication in this collection
  15 Apr 2024
• Date of issue
  2024

History

• Received
  19 Feb 2024
• Accepted
  28 Feb 2024