Reappraising Double Pendulum Dynamics across Multiple Computational Platforms

1. Introduction

The double pendulum forms an interesting physics system consisting of two pendululum hanging in a rigidly connected line. It does not only rely on the observation of immediate effects but also the study of the underlying structures that make things change. Though at first glance this machine may look mundane, its behavior will prove to be exciting and the underlying idea of balance between order and disorder that has intrigued scientists and mathematicians for centuries will be inspirationally brought to life. It is an example of the fascinating ways simple systems can exhibit complex behavior. From people like Henri Poincaré, who first took a look at non-linear systems and laid the foundation for chaos theory, to today’s powerful computational models, the double pendulum remains a symbol that the secret of nonlinear dynamics is a profound one [1].

A double pendulum consists of two massless rods, each with a concentrated point mass at its end, connected by a frictionless hinge. This seemingly simple setup actually displays a remarkable sensitivity to initial conditions, a characteristic of chaotic systems that was first discovered by British mathematician and physicist Mary Cartwright in the early 20th century [1,2]. Even tiny changes in the starting position can cause significant divergences in the paths the system takes - a concept that challenges the traditional idea put forth by French scholar Pierre-Simon Laplace known as Laplace’s demon, which suggests that knowing all initial conditions of a system guarantees the ability to predict its future evolution with total certainty [3].

The double pendulum is a prominently studied dynamical system in climate science due to its exquisite sensitivity to initial conditions, making it a valuable model for studying chaos theory and its applications in understanding Earth’s complex climate systems [4]. Gaining insights into chaotic systems like the double pendulum could prove vital for tackling one of our most pressing global issues - anthropogenic climate change driven by human activities. Several studies have used simplified models like the double pendulum to gain insights into the non-linear dynamics underlying atmospheric-oceanic flows and long-range climate predictions [5,6,7,8,9].

The double pendulum exemplifies how deterministic systems can exhibit unpredictable, chaotic behavior, bridging the gap between the simple, ordered world of classical physics and the apparent randomness we observe in complex phenomena like weather and climate patterns [10]. This chaotic complexity represents a major frontier in our scientific understanding - we have robust theories for simple systems and stochastic models for randomness, but lack a unified framework to explain the rich dynamical behaviors that emerge in between the two extremes [11]. Unraveling the chaos inherent in multi-scale systems like the double pendulum may unlock deeper insights into the fundamental laws governing our climate and the universe at large.

In this paper, our main goal is to thoroughly investigate and analyze the complex movements of a double pendulum system using detailed numerical simulations. We utilized open-source computing platforms to develop engaging visuals that can be used as effective educational resources. This graphical representation will help students, especially those with limited abstract mathematics knowledge, better understand the concepts of calculus of variations. The calculus of variations is a branch of mathematical analysis concerned with optimizing functionals, which are mappings from a space of functions to the real numbers. It addresses problems where the goal is to find a function that extremizes a given functional, either minimizing or maximizing it [12]. This field is particularly important in areas such as physics, engineering, and economics, where one seeks to optimize quantities that depend on functions, such as energy, action, or cost [13]. By presenting this information visually and in a practical context, we hope to make it easier for students to connect theoretical concepts with real-world applications.

Furthermore, allowing students to freely access the simulation code enables them to explore the computational side of the project, enhancing their comprehension of how theory translates into practice. This also gives them the opportunity to improve and perfect the models, promoting a practical approach to learning physics and applied mathematics through coding. To ensure optimal performance and efficiency, we thoroughly evaluated the free and open-source computing environments utilized in this study. This study will offer helpful information for teachers, scholars, and professionals, helping them choose the best tools for their individual computing requirements and improving resource management and promoting excellence in scientific computing.

2. Methods

In order to predict the positional paths of two point masses in a double pendulum system, we need a classical mechanics framework represented by the Euler-Lagrange equation. This is because the Euler-Lagrange equation provides a systematic approach to deriving the equations of motion for complex systems like the double pendulum, taking into account the constraints and forces involved. On the other hand, using Newtonian mechanics alone for the double pendulum can be challenging due to the system’s nonlinear nature and the presence of constraints such as the lengths of the pendulum arms. While it’s possible to analyze simpler pendulum systems using Newton’s laws directly, the double pendulum’s motion involves coupled, nonlinear differential equations that are more conveniently handled using the Euler-Lagrange formalism [14].

The Euler-Lagrange equation is essential in classical mechanics and field theory, describing the motion of particles or fields by minimizing a functional known as the action. To derive this equation from scratch, we started with the action S, defined as the integral of the Lagrangian $\left(\mathcal{L}\right)$ over time:

$$S={\int }_{t_1}^{t_2}\mathcal{L}(q,\dot{q},t)dt$$

(1)

Here, q represents the generalized coordinates, $\dot{q}$ is the derivative of q with respect to time, and t is time. The Lagrangian $\left(\mathcal{L}\right)$ has a simple and concise definition:

$$\mathcal{L}\equiv T-V$$

(2)

The kinetic energy (T) and potential energy (V) together make up the classical Lagrangian, which is just the difference between these energies in the system. This applies to classical mechanics with conservative systems, where the total energy is the sum of kinetic and potential energies. Our next step was to find the path $q\left(t\right)$ that keeps the action constant, even if the path changes a little. This basic idea is called the principle of least action.

To find this stationary path, we used the calculus of variations. Let $\delta q\left(t\right)$ be a small variation in the path $q\left(t\right)$, such that $q\left(t\right)$ becomes $q\left(t\right)+\delta q\left(t\right)$. The variation in the action is then:

$$\delta S={\int }_{t_1}^{t_2}\left[\mathcal{L}(q+\delta q,\dot{q}+\delta \dot{q},t)-\mathcal{L}(q,\dot{q},t)\right]dt$$

(3)

We expanded $L(q+\delta q,\dot{q}+\delta \dot{q},t)$ in a Taylor series around q and $\dot{q}$:

$$\mathcal{L}(q+\delta q,\dot{q}+\delta \dot{q},t)=\mathcal{L}(q,\dot{q},t)+\frac{\partial \mathcal{L}}{\partial q}\delta q+\frac{\partial \mathcal{L}}{\partial \dot{q}}\delta \dot{q}+\mathcal{O}(\delta q^2,\delta {\dot{q}}^2)$$

(4)

Substituting the expression back into equation 3 we would consider terms up to second order or higher in $\delta q$ and $\delta \dot{q}$ in the variation of the action S, $\mathcal{O}(\delta q^2,\delta {\dot{q}}^2)$:

$$\delta S={\int }_{t_1}^{t_2}\left[\frac{\partial \mathcal{L}}{\partial q}\delta q+\frac{\partial \mathcal{L}}{\partial \dot{q}}\delta \dot{q}\right]dt+\mathcal{O}(\delta q^2,\delta {\dot{q}}^2)$$

(5)

Integrating the first term by parts with respect to t and assuming that variations $\delta q\left(t_1\right)$ and $\delta q\left(t_2\right)$ vanish (boundary conditions), we got:

$$\delta S={\int }_{t_1}^{t_2}\left[\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{q}}\right)-\frac{\partial \mathcal{L}}{\partial q}\right]\delta qdt$$

(6)

For the action to be stationary, $\delta S$ must be zero for all possible variations $\delta q$. This leads to the Euler-Lagrange equation:

$$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{q}}\right)-\frac{\partial \mathcal{L}}{\partial q}=0$$

(7)

This equation governs the dynamics of the system and provides the equations of motion for the generalized coordinates $q\left(t\right)$.

To apply the Euler-Lagrange equation to the double pendulum system, we must first establish its Lagrangian, which encapsulates both the system’s kinetic and potential energies. This process started by delineating the geometric relationships governing the vertical and horizontal positions within the double pendulum system, as illustrated in Figure 1. These relationships were formalized through the following equations, denoted as equations 8 and 9:

$$\begin{array}{c}{x}_1=L_{1}sin\left({\theta }_1\right)\hfill \\ x_2=x_1+L_{2}sin\left({\theta }_2\right)\hfill \end{array}$$

(8)

$$\begin{array}{c}{y}_1=L_{1}cos\left({\theta }_1\right)\hfill \\ y_2=y_1+L_{2}cos\left({\theta }_2\right)\hfill \end{array}$$

(9)

Here, $x_1$ and $x_2$ represent the horizontal positions, while $y_1$ and $y_2$ represent the vertical positions. These positions are determined by the lengths of the pendulum arms ($L_1$ and $L_2$) and the angles (${\theta }_1$ and ${\theta }_2$).

Since the positions x and y are functions of the angles ${\theta }_1$ and ${\theta }_2$, respectively, obtaining their first derivatives requires applying the chain rule. This application results in the following expressions:

$$\begin{array}{c}{\dot{x}}_1=L_1{\dot{\theta }}_{1}cos\left({\theta }_1\right)\hfill \\ {\dot{x}}_2=L_1{\dot{\theta }}_{1}cos\left({\theta }_1\right)+L_2{\dot{\theta }}_{2}cos\left({\theta }_2\right)\hfill \end{array}$$

(10)

$$\begin{array}{c}{\dot{y}}_1=-L_1{\dot{\theta }}_{1}sin\left({\theta }_1\right)\hfill \\ {\dot{y}}_2=-L_1{\dot{\theta }}_{1}sin\left({\theta }_1\right)-L_2{\dot{\theta }}_{2}cos\left({\theta }_2\right)\hfill \end{array}$$

(11)

These derivatives ${\dot{x}}_1,{\dot{x}}_2,{\dot{y}}_1,$ and ${\dot{y}}_2$ represent the rates of change of the horizontal and vertical positions with respect to time, taking into account the angular velocities ${\dot{\theta }}_1$ and ${\dot{\theta }}_2$.

Starting from the general formula for kinetic energy $T=\frac{1}{2}m{v}^2$ (where m is the total of two-point masses of $m_1+m_2$ and v is velocity), we substituted the velocities ${\dot{x}}_1,{\dot{x}}_2,{\dot{y}}_1,$ and ${\dot{y}}_2$ from the systems of equations 10 and 11 into equation 12. This substitution yields the expanded form:

$$\begin{array}{cc}\hfill T& =\frac{1}{2}\left(m_1({\dot{x}}_1^2+{\dot{y}}_1^2)+m_2({\dot{x}}_2^2+{\dot{y}}_2^2)\right)\hfill \\ \hfill & =\frac{1}{2}\left(m_1{\left(L_1{\dot{\theta }}_{1}cos\left({\theta }_1\right)\right)}^2+m_1{\left(-L_1{\dot{\theta }}_{1}sin\left({\theta }_1\right)\right)}^2\right.\hfill \\ \hfill & \left.+m_2{\left(L_1{\dot{\theta }}_{1}cos\left({\theta }_1\right)+L_2{\dot{\theta }}_{2}cos\left({\theta }_2\right)\right)}^2+m_2{\left(-L_1{\dot{\theta }}_{1}sin\left({\theta }_1\right)-L_2{\dot{\theta }}_{2}cos\left({\theta }_2\right)\right)}^2\right)\hfill \end{array}$$

(12)

Expanding and simplifying each term step by step, we obtained expressions for the kinetic energy components. These components involve terms related to the masses $m_1$ and $m_2$, the lengths $L_1$ and $L_2$ of the pendulum arms, and the angular velocities ${\dot{\theta }}_1$ and ${\dot{\theta }}_2$. After combining and simplifying the terms, we arrived at the final form of the kinetic energy:

$$T=\frac{1}{2}(m_1+m_2)L_1^2{\dot{\theta }}_1^2+\frac{1}{2}{m}_2{L}_2^2{\dot{\theta }}_2^2+m_2{L}_1{L}_2{\dot{\theta }}_1{\dot{\theta }}_{2}cos({\theta }_1-{\theta }_2)$$

(13)

This expression captures the kinetic energy T of the double pendulum system comprehensively, incorporating the masses $m_1$ and $m_2$, the lengths $L_1$ and $L_2$ of the pendulum arms, and the angular velocities ${\dot{\theta }}_1$ and ${\dot{\theta }}_2$ in a way that reflects the system’s dynamics and interactions.

Beginning with equation 14, which defines potential energy as a function of the vertical positions $y_1$ and $y_2$, and the gravitational constant g, we can express it as:

$$V=m_{1}g{y}_1+m_{2}g{y}_2$$

(14)

Substituting the expressions for $y_1$ and $y_2$ from the system of equations 9 into the above equation yields:

$$V=m_{1}g(L_{1}cos\left({\theta }_1\right))+m_{2}g(L_{1}cos\left({\theta }_1\right)+L_{2}cos\left({\theta }_2\right))$$

(15)

We then simplified this expression to:

$$V=g((m_1+m_2)L_{1}cos\left({\theta }_1\right)+m_2{L}_{2}cos\left({\theta }_2\right))$$

(16)

Finally, to express potential energy V solely in terms of the angles ${\theta }_1$ and ${\theta }_2$, we derived:

$$V=-g((m_1+m_2)L_{1}cos\left({\theta }_1\right)+m_2{L}_{2}cos\left({\theta }_2\right))$$

(17)

Here, g represents the gravitational acceleration, and the potential energy V is derived from the vertical positions of the pendulum components and their respective masses, articulated in terms of the angles ${\theta }_1$ and ${\theta }_2$.

After deriving the expressions for kinetic energy T and potential energy V, we can now formulate the Lagrangian $\mathcal{L}$ for the double pendulum system. The Lagrangian is defined in the definition 2. Substituting the expressions we derived for T and V into this equation, we obtained:

$$\begin{array}{cc}\hfill \mathcal{L}& =\frac{1}{2}(m_1+m_2)L_1^2{\dot{\theta }}_1^2+\frac{1}{2}{m}_2{L}_2^2{\dot{\theta }}_2^2+m_2{L}_1{L}_2{\dot{\theta }}_1{\dot{\theta }}_{2}cos({\theta }_1-{\theta }_2)\hfill \\ \hfill & -\left(-g((m_1+m_2)L_{1}cos\left({\theta }_1\right)+m_2{L}_{2}cos\left({\theta }_2\right))\right)\hfill \\ \hfill & =\frac{1}{2}(m_1+m_2)L_1^2{\dot{\theta }}_1^2+\frac{1}{2}{m}_2{L}_2^2{\dot{\theta }}_2^2+m_2{L}_1{L}_2{\dot{\theta }}_1{\dot{\theta }}_{2}cos({\theta }_1-{\theta }_2)\hfill \\ \hfill & +g((m_1+m_2)L_{1}cos\left({\theta }_1\right)+m_2{L}_{2}cos\left({\theta }_2\right))\hfill \end{array}$$

(18)

This expression for the Lagrangian $\mathcal{L}$ encapsulates the dynamic behavior of the double pendulum system. It incorporates the masses $m_1$ and $m_2$, the lengths $L_1$ and $L_2$ of the pendulum arms, the angular velocities ${\dot{\theta }}_1$ and ${\dot{\theta }}_2$, and the gravitational constant g, as well as the angles ${\theta }_1$ and ${\theta }_2$ that describe the positions of the pendulum components. The Lagrangian $\mathcal{L}$ serves as a fundamental quantity in the analysis of the system’s motion and dynamics, providing a comprehensive representation of its energy and interactions.

Starting with the Lagrangian $\mathcal{L}$ derived earlier (equation 18), we applied the Euler-Lagrange equation (equation 7) to derive the equations of motion for a system with two point masses (the double pendulum in this case). The Euler-Lagrange equation for a variable $q_i$ is given by:

$$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial {\dot{q}}_i}\right)-\frac{\partial \mathcal{L}}{\partial q_i}=0$$

(19)

Applying this equation to the variables ${\theta }_1$ and ${\theta }_2$ separately, we obtained two set of equations:

$$\begin{array}{c}\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{{\theta }_1}}\right)-\frac{\partial \mathcal{L}}{\partial {\theta }_1}=0\hfill \\ \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{{\theta }_2}}\right)-\frac{\partial \mathcal{L}}{\partial {\theta }_2}=0\hfill \end{array}$$

(20)

For ${\theta }_1$, we first calculated the partial derivative of $\mathcal{L}$ with respect to ${\dot{\theta }}_1$:

$$\frac{\partial \mathcal{L}}{\partial {\dot{\theta }}_1}=(m_1+m_2)L_1^2{\dot{\theta }}_1+m_2{L}_1{L}_2{\dot{\theta }}_{2}cos({\theta }_1-{\theta }_2)$$

(21)

Then, we took the derivative of this with respect to time t:

$$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial {\dot{\theta }}_1}\right)=(m_1+m_2)L_1^2{\ddot{\theta }}_1+m_2{L}_1{L}_2{\ddot{\theta }}_{2}cos({\theta }_1-{\theta }_2)$$

(22)

Next, we calculated the partial derivative of $\mathcal{L}$ with respect to ${\theta }_1$:

$$\frac{\partial \mathcal{L}}{\partial {\theta }_1}=-m_2{L}_1{L}_2{\dot{\theta }}_1{\dot{\theta }}_{2}sin({\theta }_1-{\theta }_2)-g(m_1+m_2)L_{1}sin\left({\theta }_1\right)$$

(23)

Finally, substituting these derivatives into the Euler-Lagrange equation (equation 20) for ${\theta }_1$:

$$(m_1+m_2)L_1^2{\ddot{\theta }}_1+m_2{L}_1{L}_2{\ddot{\theta }}_{2}cos({\theta }_1-{\theta }_2)=0$$

(24)

To derive the equation of motion for ${\theta }_2$ using the Euler-Lagrange equation (equation 20), we started by calculating the partial derivative of the Lagrangian $\mathcal{L}$ with respect to the derivative of ${\theta }_2$, denoted as ${\dot{\theta }}_2$. This yields:

$$\frac{\partial \mathcal{L}}{\partial {\dot{\theta }}_2}=m_2{L}_2^2{\dot{\theta }}_2+m_2{L}_1{L}_2{\dot{\theta }}_{1}cos({\theta }_1-{\theta }_2)$$

(25)

Taking the time derivative of this expression gave us:

$$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial {\dot{\theta }}_2}\right)=m_2{L}_2^2{\ddot{\theta }}_2+m_2{L}_1{L}_2{\ddot{\theta }}_{1}cos({\theta }_1-{\theta }_2)-m_2{L}_1{L}_2{\dot{\theta }}_1{\dot{\theta }}_{2}sin({\theta }_1-{\theta }_2)$$

(26)

Next, we calculated the partial derivative of $\mathcal{L}$ with respect to ${\theta }_2$, denoted as $\frac{\partial \mathcal{L}}{\partial {\theta }_2}$, which was given by:

$$\frac{\partial \mathcal{L}}{\partial {\theta }_2}=-m_2{L}_1{L}_2{\dot{\theta }}_1{\dot{\theta }}_{2}sin({\theta }_1-{\theta }_2)-m_2{L}_{2}gsin\left({\theta }_2\right)$$

(27)

Substituting the expressions for $\frac{\partial \mathcal{L}}{\partial {\dot{\theta }}_2}$ and $\frac{\partial \mathcal{L}}{\partial {\theta }_2}$ into the Euler-Lagrange equation and simplifying this equation further results in the equation of motion for ${\theta }_2$, given as:

$$m_2{L}_2{\ddot{\theta }}_2+m_2{L}_1{\ddot{\theta }}_{1}cos({\theta }_1-{\theta }_2)-m_2{L}_1{\dot{\theta }}_1^{2}sin({\theta }_1-{\theta }_2)+m_{2}gsin\left({\theta }_2\right)=0$$

(28)

Combining the final forms of the Euler-Lagrange equations for ${\theta }_1$ (equation 24) and ${\theta }_2$ (equation 28) together form the complete system of equations of motion for the double pendulum system. They are coupled second-order ordinary differential equations (ODEs) because the acceleration terms ${\ddot{\theta }}_1$ and ${\ddot{\theta }}_2$ are dependent on each other due to the cosine and sine terms involving ${\theta }_1$ and ${\theta }_2$. This coupling reflects the interdependence of the pendulum’s motions and positions, making the system dynamically rich and challenging to analyze without numerical or advanced analytical techniques.

For the numerical analysis, we substituted the following expressions $({\dot{\theta }}_1={\omega }_1,{\dot{\theta }}_2={\omega }_2,\Delta \theta ={\theta }_1-{\theta }_2)$ to simplify the system of equations for computational purposes. Substituting these into the given equations (equations 24 and 28), we obtained the following system:

$$\begin{array}{c}(m_1+m_2)L_1{\dot{\omega }}_1+m_2{L}_2{\dot{\omega }}_{2}cos\left(\Delta \theta \right)+m_2{L}_2{\omega }_2^{2}sin\left(\Delta \theta \right)+(m_1+m_2)gsin{\theta }_1=0\hfill \\ m_2{L}_2{\dot{\omega }}_2+m_2{L}_1{\dot{\omega }}_{1}cos\left(\Delta \theta \right)-m_2{L}_1{\omega }_1^{2}sin\left(\Delta \theta \right)+m_{2}gsin\left({\theta }_2\right)=0\hfill \end{array}$$

(29)

This transformed system allows us to numerically solve for the angular frequencies ${\omega }_1$ and ${\omega }_2$ given initial conditions and system parameters.

To simplify the numerical analysis, we first rewrote the original system of equations in terms of the defined variables:

$$\left\{\begin{array}{c}\alpha =(m_1+m_2)L_1\hfill \\ \beta =m_2{L}_{2}cos\left(\Delta \theta \right)\hfill \\ \gamma =m_2{L}_{1}cos\left(\Delta \theta \right)\hfill \\ \delta =m_2{L}_2\hfill \\ \epsilon =-m_2{L}_2{\omega }_2^{2}sin\left(\Delta \theta \right)-m_{2}gsin\left({\theta }_2\right)\hfill \\ \zeta =m_2{L}_2{\omega }_1^{2}sin\left(\Delta \theta \right)-m_{2}gsin\left({\theta }_2\right)\hfill \end{array}\right.$$

(30)

We then rewrote the system of equations in matrix form:

$$\left(\begin{array}{cc}\alpha & \beta \\ \gamma & \delta \end{array}\right)\left(\begin{array}{c}{\dot{\omega }}_1\\ {\dot{\omega }}_2\end{array}\right)=\left(\begin{array}{c}\epsilon \\ \zeta \end{array}\right)$$

(31)

Finally, solving for ${\dot{\omega }}_1$ and ${\dot{\omega }}_2$, we obtained:

$$\begin{array}{c}{\dot{\omega }}_1=\frac{\epsilon \delta -\beta \zeta }{\alpha \delta -\beta \gamma }\hfill \\ {\dot{\omega }}_2=\frac{\alpha \zeta -\gamma \epsilon }{\alpha \delta -\beta \gamma }\hfill \end{array}$$

(32)

These derived equations were expressed in terms of the defined variables, providing a more organized and manageable representation for the numerical experiments.

For the numerical simulation needs in this study, we used parameters consistent with the physical properties of the system. These include an acceleration due to gravity (g) of $9.8{\mathrm{m}/\mathrm{s}}^2$, the lengths of the pendulum arms $L_1$ and $L_2$ set to 2 meters and 1 meter respectively, and the masses of the pendulums ($m_1$ and $m_2$) set at 1 kilogram and 2 kilograms respectively. The simulation time $\left(t\right)$ is chosen to be 10 seconds with 10,000 linearly time spacing, providing sufficient duration to observe the system’s behavior.

With these parameters established, we set the initial conditions for the simulation. The initial angles (${\theta }_1$ and ${\theta }_2$) are set to 3.14 radians (which is equivalent to ${180}^{\circ }$) for ${\theta }_1$ and 1.57 radians (equivalent to ${90}^{\circ }$) for ${\theta }_2$. Additionally, the initial angular velocities (${\omega }_1$ and ${\omega }_2$) are initialized to 0 radians per second, representing a starting point where the pendulums are at rest. For the subsequent simulation, we changed the both angular velocities to $0.001$ radians per seconds fo the sensitivity test with the initial conditions.

We used the Runge-Kutta-Fehlberg (RKF) method for approximating the solution of our problem (system of equations 32). RKF is a powerful numerical integration technique used to solve systems of ODEs with high accuracy and computational efficiency [13]. Its adaptability makes it particularly suitable for dynamic systems like the double pendulum, where the motion can be complex and highly nonlinear.

We started with the initial conditions and set the initial step size h. Then, we proceeded to define the predictor step by using the 4th order Runge-Kutta formulas to predict the solution at the next time step. For a general ODE of the form $\frac{dy}{dt}=f(t,y)$, the 4th order Runge-Kutta formulas are:

$$\begin{array}{c}{k}_1=hf(t_n,y_n)\hfill \\ k_2=hf(t_n+\frac{h}{2},y_n+\frac{k_1}{2})\hfill \\ k_3=hf(t_n+\frac{h}{2},y_n+\frac{k_2}{2})\hfill \\ k_4=hf(t_n+h,y_n+k_3)\hfill \end{array}$$

(33)

The predicted solution at $t_{n+1}$ is then given by:

$$y_{n+1}^{\left(4\right)}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$$

(34)

We used the 5th order Runge-Kutta formulas to compute a more accurate estimate of the solution at the next time step. The 5th order Runge-Kutta formulas are similar to the 4th order ones but include an additional evaluation point:

$$\begin{array}{c}{k}_1=hf(t_n,y_n)\hfill \\ k_2=hf(t_n+\frac{h}{4},y_n+\frac{k_1}{4})\hfill \\ k_3=hf(t_n+\frac{3h}{8},y_n+\frac{3k_1}{32}+\frac{9k_2}{32})\hfill \\ k_4=hf(t_n+\frac{12h}{13},y_n+\frac{1932k_1}{2197}-\frac{7200k_2}{2197}+\frac{7296k_3}{2197})\hfill \\ k_5=hf(t_n+h,y_n+\frac{439k_1}{216}-8k_2+\frac{3680k_3}{513}-\frac{845k_4}{4104})\hfill \\ k_6=hf(t_n+\frac{h}{2},y_n-\frac{8k_1}{27}+2k_2-\frac{3544k_3}{2565}+\frac{1859k_4}{4104}-\frac{11k_5}{40})\hfill \end{array}$$

(35)

The corrected solution at $t_{n+1}$ is then given by:

$$y_{n+1}^{\left(5\right)}=y_n+\frac{1}{90}(7k_1+32k_3+12k_4+32k_5+7k_6)$$

(36)

Then, we calculated the local truncation error by comparing the 4th and 5th order solutions. The error estimate ${ϵ}_n$ is given by:

$${ϵ}_n=|y_{n+1}^{\left(5\right)}-y_{n+1}^{\left(4\right)}|$$

(37)

We compared the error estimate ${ϵ}_n$ to a predefined tolerance. If ${ϵ}_n$ is within the tolerance, accept the step and update the solution. If ${ϵ}_n$ exceeds the tolerance, reduce the step size h and we repeated the process until the error is acceptable. Finally, we continued integrating until reaching the desired end time or number of steps.

The adaptive step-size control implemented in the RKF method ensures accurate integration while minimizing computational costs. By dynamically adjusting the step size based on error estimation, the RKF method provides accurate numerical approximations of the system’s behavior over time, rendering it an effective tool for analyzing complex dynamic systems such as the double pendulum. However, in the present study, we did not employ the RKF method directly from the beginning. Instead, we utilized several platform-specific tools, including SciPy in Python [15], deSolve in R [16], DifferentialEquations.jl in Julia [17], and the built-in function ode45 in GNU Octave [18].

The motivation behind using multiple programming languages and their respective computing environments was twofold. First, it allowed us to leverage the strengths and unique features of each language and environment, enabling a comprehensive exploration of the double pendulum problem from diverse perspectives. Second, it facilitated a comparative analysis of the performance and accuracy of different numerical solvers across these platforms.

In Python, the SciPy library provided a robust and well-established collection of scientific computing tools, including numerical solvers for ODEs. The flexibility and ease of use of Python, combined with the power of SciPy, made it a suitable choice for implementing and testing numerical methods for the double pendulum problem [19,20]. However, the interpreted nature of Python may introduce performance overhead compared to compiled languages.

The R programming language, with its deSolve package, offered a specialized environment for solving ODEs and differential algebraic equations (DAEs) [21]. R’s strong emphasis on statistical analysis and data visualization made it an attractive option for exploring the double pendulum problem, enabling efficient data analysis and visual representation of the results. Nonetheless, the high-level nature of R may lead to performance limitations for computationally intensive simulations.

Julia, a relatively new language designed for scientific computing, provided a high-performance computing environment through its DifferentialEquations.jl package. The combination of Julia’s dynamic programming capabilities and its efficient just-in-time (JIT) compilation made it a promising choice for solving the double pendulum problem with potentially improved computational performance [22,23]. However, the relatively young ecosystem of Julia may present challenges in terms of package maturity and community support.

Finally, GNU Octave, a high-level language primarily intended for numerical computations, offered the built-in ode45 function, which implements a versatile Runge-Kutta method for solving ODEs. GNU Octave’s compatibility with MATLAB^® syntax and its open-source nature made it an accessible option for researchers and students alike [18]. However, its performance may be limited compared to lower-level languages or specialized numerical libraries. By employing these various computing environments, we aimed to evaluate the trade-offs between performance, accuracy, and ease of use for each approach.

When conducting scientific research involving computational simulations, it is crucial to ensure reproducibility and transparency in the methods used. In this particular study, we aim to benchmark the performance of four different computing environments: Python, R, Octave/MATLAB, and Julia, by running a script that simulates the motion of a double pendulum. The script was executed 1,000 times in each computing environment, and the runtime and memory usage for each run will be measured and recorded. The benchmarking was performed on a Fedora Linux 39 (Budgie) x86_64 system with a 20LB0021US ThinkPad P52s laptop equipped with an Intel i7-8550U (8) @4.000GHz CPU.

The rationale behind this benchmarking approach is multifaceted. Firstly, it promotes reproducibility by providing the scripts and the benchmarking script, allowing other researchers to easily replicate the computational experiments and verify the results. Additionally, it facilitates a direct comparison of the runtime and memory usage across different computing environments for the same task, providing valuable insights into the relative performance of each environment. This information can guide researchers in choosing the most suitable tool for their specific computational needs.

Furthermore, running the scripts 1,000 times in each environment helps to account for potential variability and ensures that the performance measurements are consistent and reliable. This is particularly important when dealing with computationally intensive simulations like the double pendulum simulation, where minor fluctuations in hardware or software configurations can impact the results. With a large number of runs, it becomes possible to perform statistical analysis on the collected data, such as calculating confidence intervals, identifying outliers, and determining the statistical significance of any observed performance differences.

Measuring memory usage alongside runtime can provide insights into the scalability and resource requirements of each computing environment. This information is crucial when working with large-scale simulations or data-intensive applications, where efficient memory management is essential. By including the benchmarking methodology and results in the scientific paper, researchers demonstrate transparency and allow others to critically evaluate the computational approaches used in the study [24]. This aligns with the principles of open science and facilitates future replication and extension of the research.

After conducting the benchmarking process and collecting the runtime and memory usage data for each computing environment, we performed statistical analyses to determine if there were significant differences in performance among the platforms. Specifically, we employed the Kruskal-Wallis (K-W) test and Dunn’s test with Bonferroni adjustment, which are commonly used in various scientific fields for comparing multiple groups or treatments.

The K-W test is a non-parametric alternative to the one-way analysis of variance (ANOVA) and is used when the assumptions of normality and homogeneity of variances are violated [25]. It is a rank-based test that evaluates whether the populations from which the samples were drawn have the same distribution. The test statistic for the K-W test was calculated as:

$$H=\frac{12}{N(N+1)}\sum _{i=1}^k\frac{R_i^2}{n_i}-3(N+1),$$

(38)

where N is the total number of observations across all groups, k is the number of groups, $R_i$ is the sum of ranks for group i, and $n_i$ is the number of observations in group i. The null hypothesis for the K-W test is that the populations have the same distribution, while the alternative hypothesis is that at least one population has a different distribution from the others.

If the K-W test indicates significant differences among the groups, post-hoc tests are typically performed to determine which specific groups differ from each other. In our case, we used Dunn’s test with Bonferroni adjustment, which is a multiple comparison procedure that adjusts the significance level to control the family-wise error rate (FWER) [26]. The Dunn’s test statistic for comparing groups i and j was calculated as:

$$Z=\frac{R_i-R_j}{\sqrt{\frac{N(N+1)}{12}\left(\frac{1}{n_i}+\frac{1}{n_j}\right)}}$$

(39)

, where $R_i$ and $R_j$ are the average ranks for groups i and j, respectively, N is the total number of observations across all groups, and $n_i$ and $n_j$ are the number of observations in groups i and j, respectively. The Bonferroni adjustment was applied by dividing the desired significance level ($\alpha$) by the number of pairwise comparisons made, resulting in an adjusted significance level of $\alpha /\left(\genfrac{}{}{0pt}{}{k}{2}\right)$, where k is the number of groups.

The K-W test and Dunn’s test with Bonferroni adjustment are widely used in various scientific fields [27,28,29]. These tests are particularly useful when dealing with non-normal data or when the assumptions of parametric tests (such as ANOVA) are violated. They provide a robust and reliable way to compare multiple groups or treatments, ensuring that any observed differences are statistically significant and not due to chance alone. By applying these statistical tests to our benchmarking data, we aimed to determine if there were significant differences in performance among the four computing environments (Python, R, GNU Octave, and Julia) for the double pendulum simulation task. These non-parametric procedures were implemented in Python using the SciPy stats module [15] and the scikit-posthoc library [30].

After conducting the benchmarking process to measure the runtime and memory usage of the double pendulum simulation across different computing environments, we further analyzed the obtained time series data to gain insights into the underlying dynamics of the system. One of the analysis techniques we employed was the calculation of Shannon entropy [31], which is a measure derived from information theory that quantifies the amount of information or uncertainty present in a random variable or time series. The Shannon entropy was calculated using the following equation:

$$H\left(X\right)=-\sum _{i=1}^{n}p\left(x_i\right){log}_{2}p\left(x_i\right)$$

(40)

, where $H\left(X\right)$ is the Shannon entropy, n is the number of unique values in the time series, and $p\left(x_i\right)$ is the probability of the occurrence of the value $x_i$ in the time series. A higher Shannon entropy value indicates a higher degree of uncertainty or unpredictability in the time series, while a lower entropy value suggests a more predictable or regular pattern. This measure can provide valuable insights into the complexity and dynamics of the double pendulum system.

In our analysis, we calculated the Shannon entropy for each variable (e.g., $x_1$, $y_1$, ${\theta }_1$, ${\theta }_2$, ${\omega }_1$, ${\omega }_2$) in the time series data obtained from the double pendulum simulation. To interpret the entropy scores, we established thresholds based on the following criteria: low entropy (predictable) for entropy scores $\le 0.5$, medium entropy (some predictability) for entropy scores between 0.5 and 1.0, and high entropy (unpredictable) for entropy scores $>1.0$. These thresholds are commonly used in various applications and provide a convenient way to interpret the degree of predictability or uncertainty present in the time series data. By calculating and analyzing the Shannon entropy of the time series data, we can gain insights into the predictability and complexity of the double pendulum system’s dynamics. A high entropy score for a particular variable suggests that the corresponding time series is highly unpredictable or complex, while a low entropy score indicates a more regular or predictable pattern.

The choice of Shannon entropy as an analysis technique was motivated by its strong theoretical foundation in information theory and its widespread use in various scientific fields for quantifying the complexity and uncertainty of dynamical systems. Furthermore, the interpretation of entropy scores based on predefined thresholds provides a convenient and standardized way to categorize the time series data into different levels of predictability or complexity.

After conducting the entropy measurement, we employed the Kolmogorov-Smirnov (K-S) test to assess the sensitivity of the system to initial conditions, which is a characteristic feature of chaotic systems. The K-S test is a non-parametric statistical test that compares the cumulative distribution functions (CDFs) of two samples to determine if they are drawn from the same underlying distribution [32]. In the context of dynamical systems analysis, the K-S test can be used to compare the time series obtained from the original system with slightly perturbed initial conditions, allowing us to quantify the divergence between the trajectories. Let $F\left(x\right)$ and $G\left(x\right)$ be the empirical CDFs of the original time series and the perturbed time series, respectively. The K-S statistic is defined as the maximum absolute difference between these two CDFs:

$$D_{n,m}=\underset{x}{sup}|F\left(x\right)-G\left(x\right)|$$

(41)

, where ${sup}_x$ represents the supremum (least upper bound) of the set of absolute differences between the CDFs over all possible values of x.

The null hypothesis for the K-S test is that the two samples are drawn from the same continuous distribution, while the alternative hypothesis is that they are drawn from different distributions. The null hypothesis is rejected if the K-S statistic $D_{n,m}$ exceeds a critical value that depends on the chosen significance level and the sample sizes n and m. In our analysis, we compared the original time series $(x_1,y_1,{\theta }_1,{\theta }_2,{\omega }_1,{\omega }_2)$ obtained from the double pendulum simulation with slightly perturbed time series, where the initial conditions for ${\omega }_1$ and ${\omega }_2$ were perturbed by 0.001 rad/s. By applying the K-S test to each pair of original and perturbed time series, we can assess whether the small perturbation in the initial conditions leads to a significant divergence in the trajectories over time. If the K-S test rejects the null hypothesis, indicating that the original and perturbed time series are drawn from different distributions, it suggests that the double pendulum system is sensitive to initial conditions, which is a hallmark of chaotic behavior. Conversely, if the test fails to reject the null hypothesis, it implies that the system is less sensitive to small perturbations in the initial conditions, suggesting a more predictable or regular dynamics.

The choice of the K-S test was motivated by its nonparametric nature, which means that it does not make assumptions about the underlying distribution of the data, making it suitable for analyzing complex dynamical systems where the distribution is often unknown or difficult to model parametrically. Additionally, the K-S test is widely used in various scientific fields for comparing distributions and detecting deviations from a hypothesized distribution [33,34,35], further justifying its application in our analysis. We employed SciPy’s stats module [15] to conduct an automated K-S test.

3. Results and Discussion

The K-W test, a non-parametric test for comparing multiple groups, yielded a test statistic of 3523.203 and a p-value of 0.000 for the runtime data, indicating that at least one group’s median runtime significantly differs from the others. Dunn’s post-hoc test, with Bonferroni adjustment for multiple comparisons, showed that all pairs of groups had p-values below 0.05, indicating significant differences in their median runtimes.

The results (Figure 2) indicate that R and GNU Octave exhibited the fastest runtimes, with mean values of 1.914 seconds and 1.944 seconds, respectively. The fast performance of R can be attributed to the use of the deSolve package, which provides efficient solvers for ODEs. GNU Octave, on the other hand, likely leverages optimized numerical libraries or solvers for ODE systems. Python followed closely with a mean runtime of 2.503 seconds, benefiting from the SciPy and NumPy libraries, which provide efficient numerical computations and ODE solvers. Notably, Julia had the slowest runtime performance, with a mean of 35.701 seconds, which is significantly longer than the other environments. One potential factor contributing to Julia’s slower performance could be the overhead associated with its JIT compilation process. JIT compilation can introduce additional runtime overhead, particularly for computationally intensive tasks like solving ODEs, which may have impacted Julia’s overall runtime performance in this specific implementation.

The K-W test for memory usage data yielded a test statistic of 3375.563 and a p-value of 0.000, indicating that at least one group’s median memory usage significantly differs from the others. Dunn’s post-hoc test showed that the pairs of groups with p-values below 0.05, indicating significant differences in their median memory usage, were GNU Octave-Julia, GNU Octave-Python, Julia-Python, and Julia-R. The results (Figure 2) show that Python, R, and GNU Octave exhibited similar memory usage patterns, with a mean of approximately 142.84 MB. Julia, on the other hand, had significantly higher memory usage, with a mean of 1031.703 MB, which is about seven times higher than the other environments. The differences in memory usage across the computing environments can be attributed to various factors, such as the memory management strategies of the respective programming languages and the specific implementation of the double pendulum simulation in each environment. However, it is worth noting that the pure Julia implementation using DifferentialEquations.jl may have different memory requirements or optimization strategies compared to the packages or libraries used in the other environments.

When selecting the appropriate computing environment for the double pendulum simulation, the statistical significance of the runtime and memory usage differences should be considered alongside other factors, such as existing codebase, familiarity with the language, and potential optimizations. If runtime performance is the primary concern, R and GNU Octave would be the recommended choices based on the provided results, with R leveraging the deSolve package and GNU Octave likely utilizing optimized numerical libraries or solvers for ODE systems. Python, with the SciPy and NumPy libraries, also exhibited a relatively good runtime performance and could be a viable option, especially if existing Python code or familiarity with the language is a consideration. If memory usage is a critical factor and the higher memory footprint of Julia is not a concern, Julia could be considered, potentially with further optimization efforts or alternative implementations. It is important to note that the pure Julia implementation using DifferentialEquations.jl may have different memory requirements or optimization strategies compared to the packages or libraries used in the other environments.

Figure 3 and Figure 4 present time series plots of a double pendulum system with slightly different initial angular velocities for the inner and outer pendulum. Each subplot displays the time evolution of the angular positions of both pendulums over 10 seconds with 1000 time steps.

The sensitivity to initial conditions is evident as the trajectories diverge significantly, despite only a small difference of 0.001 rad/s in the initial angular velocities (${\omega }_1$ and ${\omega }_2$) of the inner and outer pendulums. This phenomenon is characteristic of chaotic systems, where minuscule changes in initial conditions can lead to vastly different long-term behaviors, making precise predictions challenging [1]. The time series exhibit intricate patterns with recurring oscillations, indicative of the underlying nonlinear dynamics. However, the trajectories quickly diverge, showcasing the system’s sensitivity to initial conditions, a hallmark of chaos theory [36]. The time series initially following similar paths but rapidly diverging due to the infinitesimal differences in the initial angular velocities, further demonstrating the sensitive dependence on initial conditions in chaotic systems.

These results highlight the importance of accurately measuring and accounting for initial conditions in chaotic systems, as even minute uncertainties can amplify over time, leading to substantial deviations in the system’s behavior. The double pendulum serves as a compelling example of the intricate dynamics and unpredictability that can arise in nonlinear systems due to their extreme sensitivity to initial conditions, a fundamental concept in chaos theory [37].

To further quantify and statistically validate this divergence, we performed the K-S test on the time series data for various variables, including positions (x, y, $\theta$) and angular velocities $({\omega }_1,{\omega }_2)$, before and after the initial perturbation. The K-S test results revealed that for each of these variables, the distributions of the time series data before and after the perturbation were significantly different, with p-values less than 0.05, leading to the rejection of the null hypothesis. This statistically confirms that the slight change in the initial angular velocity has a profound impact on the system’s dynamics, causing the distributions of the position and velocity variables to diverge substantially.

This divergence can be attributed to the chaotic nature of the double pendulum system, where the nonlinear equations governing its motion exhibit extreme sensitivity to initial conditions. Even minuscule differences in the starting values can rapidly amplify over time, leading to vastly different trajectories in the phase space [36]. The K-S test results corroborate this fundamental characteristic of chaos, demonstrating that the distributions of the system’s variables become statistically distinct due to the exponential divergence of initially close trajectories, a phenomenon known as the "butterfly effect" [38].

Figure 5 and Figure 5 show the trajectories of the inner and outer pendulums in the $x-y$ plane, providing a visualization of their motion over time. However, Figure 5 and Figure 5 represent the phase space diagrams of the coupled pendulum system. A phase space diagram is a powerful tool for studying dynamical systems, as it provides a geometric representation of the system’s state at any given time [13]. In the case of the coupled pendulum system, the phase space diagrams plot the angular position $\left(\theta \right)$ of each pendulum against its angular velocity $\left(\omega \right)$, capturing the evolution of the system’s state over time.

The complex patterns observed in Figure 5 and Figure 5 are indicative of the chaotic nature of the coupled pendulum system. The phase space trajectories exhibit intricate, aperiodic behavior, never repeating the same pattern or revisiting the same state. The presence of chaos in the coupled pendulum system has significant consequences for its predictability and controllability [36]. Even with precise knowledge of the initial conditions and governing equations, the system’s behavior becomes increasingly difficult to predict over long time scales due to the exponential divergence of nearby trajectories in the phase space.

Furthermore, chaotic systems often exhibit strange attractors in the phase space, which are geometrically complex structures that the system’s trajectories are confined. The presence of strange attractors in the phase space diagrams of the coupled pendulum system suggests that the system’s dynamics are governed by an underlying deterministic process, despite the apparent randomness and unpredictability of its behavior [1].

The high Shannon entropy values (greater than 1.0) mentioned for the time series data of the coupled pendulum system are consistent with the observed chaotic dynamics. For chaotic systems, the Shannon entropy is expected to be high due to the inherent unpredictability and lack of regularity in the system’s behavior [39].

4. Conclusions

The double pendulum system serves as a compelling example of the intricate dynamics and unpredictability that can arise in nonlinear systems due to their extreme sensitivity to initial conditions, a fundamental concept in chaos theory. Through numerical simulations and extensive analysis, we have demonstrated the chaotic behavior of the double pendulum, characterized by the exponential divergence of trajectories, the presence of strange attractors in the phase space, and high Shannon entropy values. These findings underscore the importance of accurately measuring and accounting for initial conditions in chaotic systems, as even minute uncertainties can amplify over time, leading to substantial deviations in the system’s behavior. Furthermore, our comparative study of different computing environments (Python, R, GNU Octave, and Julia) has revealed significant differences in their runtime performance and memory usage for the double pendulum simulation task. This information can guide researchers in selecting the most appropriate computing environment based on their specific needs and resource constraints.

To further enhance the computational performance of the double pendulum simulations, future work could explore the use of accelerated computing techniques. Numba, a JIT compiler for Python [40], can be leveraged to optimize numerical computations by compiling Python code to efficient machine instructions, potentially improving the runtime performance of the Python-based simulations [41,42,43]. Additionally, the utilization of GPU-accelerated computing libraries like CuPy (CUDA for Python) [44] could significantly accelerate the simulations by offloading computationally intensive tasks to the highly parallel architecture of modern graphics processing units (GPUs) [45,46,47]. The massive parallelism provided by GPUs can lead to substantial speedups, especially for large-scale simulations or ensemble runs. While the double pendulum system studied in this work was modeled as a deterministic system, future studies could explore the effects of incorporating stochastic elements. Real-world systems often exhibit random fluctuations or noise, which can significantly impact the system’s dynamics and introduce additional complexity. By incorporating stochastic components into the double pendulum model, researchers could investigate the interplay between deterministic chaos and random noise, potentially revealing new insights into the behavior of complex dynamical systems.

Future investigations could also focus on studying chaotic synchronization, a phenomenon where two or more chaotic systems can become synchronized, exhibiting correlated behavior despite their inherent unpredictability. Exploring chaotic synchronization in coupled double pendulum systems or the potential applications of such synchronization in various fields, such as secure communication, signal processing, and control systems, could yield valuable insights. As mentioned in the introduction, the double pendulum system serves as a valuable model for studying chaos theory and its applications in understanding Earth’s complex climate systems. Future work could explore the potential use of the double pendulum as a simplified model for investigating the nonlinear dynamics underlying atmospheric-oceanic flows and developing improved long-range climate predictions. By pursuing these future directions, researchers can continue to deepen our understanding of chaotic systems, leverage advanced computational techniques for efficient simulations, and explore the potential applications of chaos theory in various scientific and engineering domains.