Comment on “Hyperchaos in constrained Hamiltonian system and its control” by J. Li, H. Wu and F. Mei

Abstract

The aim of this comment is to show that discovery of hyperchaos in three systems investigated in Li et al. (Nonlinear Dyn 94(3):1703–1720, 2018) is not correct. It is justified both theoretically and numerically. Corrected calculations of Lyapunov exponents and corresponding bifurcation diagram are given. Examples of hyperchaotic Hamiltonian multiple pendulum systems are presented.

1 Introduction

Chaos and hyperchaos in nonlinear dynamical systems are well-known but still fascinating phenomena. We can find many interesting applications of chaotic and hyperchaotic systems in various fields of science. However, it seems that the hyperchaos in Hamiltonian systems still does not attract sufficient attention. The aim of the authors in [1] was to propose a Hamiltonian system with two degrees of freedom and its restricted versions with a holonomic or a non-holonomic constraint exhibiting hyperchaotic phenomena.

We show in this comment that hyperchaos cannot appear neither in Hamiltonian systems with two degrees of freedom nor in its restricted versions with a holonomic or a non-holonomic constraint. This fact follows from the well-known properties of Lyapunov exponents for Hamiltonian systems. Incorrect conclusions of [1] were caused probably by insufficient accuracy of numerical calculations, or by too small number of iterations. Moreover, the initial conditions chosen in [1] for restricted Hamiltonian systems do not satisfy equations of constraints.

Having in mind the importance of the subject and in order to be constructive, we show hyperchaotic phenomena in Hamiltonian systems with more than two degrees of freedom through examples with triple and quartic pendula.

This comment is organized as follows. In Sect. 2, we recall basic properties of Lyapunov exponents that will be used in explanations of results of our corrections made for numerical calculations in [1]. In Sect. 3 will be presented the basic Hamiltonian system considered in this article and with preliminary explanations and description of its dynamics on the base of Poincaré sections. Section 4 contains corrected Lyapunov exponents spectra with the corresponding bifurcations diagrams. Section 5 concerns the restricted version of the considered Hamiltonian system with a holonomic or a non-holonomic constraint. It starts from the general explanations concerning the applied restriction procedure. Then in two subsections is considered the analysed Hamiltonian system with a holonomic or a non-holonomic constraint, respectively. Finally, in Sect. 6 we show the true hyperchaotic dynamics in Hamiltonian multiple pendulum systems.

2 Properties of Lyapunov exponents

Let us consider the general autonomous system of ordinary differential equations

$$\begin{aligned} \dot{\varvec{x}}=\varvec{v}(\varvec{x}), \qquad \varvec{x}\in \mathbb {R}^n. \end{aligned}$$

(1)

Let $\varvec{y}=R(\varvec{x})$ will be an invertible transformation. The transformed system takes the form

$$\begin{aligned} \dot{\varvec{y}}=R'(R^{-1}(y))\varvec{v}\left( R^{-1}(\varvec{y})\right) , \end{aligned}$$

(2)

where $R'$ denotes the Jacobi matrix of the transformation with respect to its argument. If the vector field $\varvec{v}$ has the property

$$\begin{aligned} R'(R^{-1}(y))\varvec{v}\left( R^{-1}(\varvec{y})\right) =-\varvec{v}(\varvec{y}) \end{aligned}$$

(3)

or equivalently

$$\begin{aligned} R'(\varvec{x})\varvec{v}(\varvec{x})=-\varvec{v}\left( R(\varvec{x})\right) , \end{aligned}$$

(4)

then we say that system is time reversal because transformation $(\varvec{x},t)\rightarrow \left( R(\varvec{x}),-t\right)$ preserves the form of system (1). For time-reversal systems, if $\lambda$ is a Lyapunov exponent, then also $-\lambda$ is a Lyapunov exponent [2,3,4].

The most known symmetry of this type is $R(\varvec{q},\varvec{p})=(\varvec{q},-\varvec{p})$ for classical Hamiltonian systems. Benettin et al. in [5] proved that for a Hamiltonian system the Lyapunov exponents obey the pairing rule

$$\begin{aligned} \lambda _i=-\lambda _{N+i},\qquad i=1,2,\ldots N, \end{aligned}$$

(5)

where N denotes the number of degrees of freedom.

Symmetries and conservation laws have also influence on a spectrum of Lyapunov’s exponents. Namely, their presence manifests itself by vanishing of exponents, see Sec. 2.5.6 in [6]. For an autonomous dynamical system, its invariance with respect to time shift, which gives one zero Lyapunov exponent, is the example of such symmetry. Preservation of volume in the phase space implies that sum of all Lyapunov exponents must be equal zero, see, for example, Sec. 2.3 in [7]. Finally, zero exponents may also (non-generically) occur at bifurcation points, where some direction is (linearly) marginally stable and for periodic or quasi-periodic orbits [6, 7]. For dissipative systems the sum of all Lyapunov’s exponents must be negative. Moreover, the Lyapunov exponents characterize the attractor, see e.g., Sec. 5.3 in [8].

For autonomous Hamiltonian flows, one pair of Lyapunov’s exponents is always zero due to time shift continuous symmetry and conservation in time of Hamiltonian. Moreover, the presence of any additional first integral, functionally independent with the Hamiltonian, leads to vanishing of another pair of Lyapunov’s exponents. Thus, the determination of a full spectrum of Lyapunov’s exponents can be used as an indicator of a number of independent integrals of motion that a considered dynamical system may possess. A detailed explanation can be found in the classical books [6, 9,10,11], while for the rigorous proofs please consult [2, 5].

According to the generally accepted terminology for autonomous differential systems, chaos appears when at least one Lyapunov exponent is positive, while hyperchaos is manifested by at least two positive Lyapunov exponents. Thus, the minimal dimension for a dissipative non-time-reversal hyperchaotic system is 4 as, for example, in the first hyperchaotic system due to Rösler [12], and for time-reversal ones, the minimal dimension is 5 [13]. For Hamiltonian systems, hyperchaos can appear when the number of degrees of freedom is at least 3.

3 Description of the system its dynamics and motivation

General three-parameter family of Hamiltonian system analysed in [1] is given by the Hamiltonian function

$$\begin{aligned}&H=T+V,\nonumber \\&T=\frac{1}{2}\left( p_1^2+p_2^2\right) ,\nonumber \\&V=-\frac{\theta }{2}\left( q_1^2+q_2^2\right) +\eta q_1\left( q_2-\frac{\eta }{2}q_1\right) +\frac{\beta }{4}\left( q_1^4+q_2^4\right) ,\nonumber \\ \end{aligned}$$

(6)

where $(q_1,q_2)$ are the generalized coordinates and $(p_1,p_2)$ their corresponding momenta, and $\theta ,\eta ,\beta \in \mathbb {R}^+$. Equations of motion generated by Hamiltonian (6) are the following

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q}_1=p_1,\\ \dot{q}_2=p_2,\\ \dot{p}_1=(\theta +\eta ^2)q_1-\eta q_2-\beta q_1^3,\\ \dot{p}_2=-\eta q_1+\theta q_2-\beta q_2^3. \end{array}\right. } \end{aligned}$$

(7)

Fig. 1

Poincaré cross sections of system (7) made for the parameters $\beta =1, \theta =1, \eta =0.01,$ with cross-plain $q_2=0$ and $p_2>0$

Since evolution of system (7) takes place in a four-dimensional phase space, we can get a quick insight into its dynamics by making several Poincaré cross sections for fixed values of the parameters, for instance, $\beta =1, \theta =1, \eta =0.01$. Points creating patterns presented in Fig. 1 were obtained as traces of intersections of orbits calculated numerically with a properly generic surface of section $q_2=0$ with $p_2>0$, in a three-dimensional hypersurface defined by a constant energy level $H=E$. Figures are ordered according to increasing energy of the system. They show that, for the chosen values of the parameters, the system possesses chaotic behaviour and it is in general not integrable. As typical in Hamiltonian systems, we can observe a coexistence of three types of motion: periodic (finite numbers of points of orbit), quasi-periodic (closed loops) and chaotic (scattered points). We do not observe an appearance of strange attractors because Hamiltonian systems are divergence-free, and thus, there are no contractions of volume in the whole phase space (see Liouville theorem). Thus, it is unclear why the authors of [1] have used the notion of “hyperchaotic attractor” in the captions of Poincaré maps visible in Figs. 10, 11, 16 of paper [1], which is misleading.

Although the periodic, quasi-periodic and chaotic behaviour of the system can be easily detected in the Poincaré sections, its hyperchaotic behaviour can be recognized only in their Lyapunov exponents spectrum. Indeed, authors of the paper [1] have been searching for hyperchaotic behaviour of system (7) by computing the Lyapunov spectrum for certain values of the parameters. Let us mention just three of them. Namely, for the chosen initial condition

$$\begin{aligned} (q_{10},q_{20},p_{10},p_{20})=(1,0.1,0.01,0.001), \end{aligned}$$

(8)

and for specified values of the parameters

$$\begin{aligned}&\theta =1,\quad \eta =0.01,\quad \beta \in [0.1,10], \end{aligned}$$

(9)

$$\begin{aligned}&\beta =1, \quad \eta =0.01,\quad \theta \in [0.1,10], \end{aligned}$$

(10)

$$\begin{aligned}&\theta =1, \quad \beta =1,\qquad \eta \in [0.001,10], \end{aligned}$$

(11)

they calculated Lyapunov’s exponents spectra using MATLAB, without specifying by which method/algorithm these exponents were calculated. The results are depicted in Figs. 2, 4 and 5 in [1]. Surprisingly enough, the obtained spectra show that system (7) is hyperchaotic for a wide range values of the parameters. Moreover, for $\theta =1, \eta =0.01$ over the range $\beta \in [1.131,1.735]$, the system has three positive Lyapunov exponents.

Fig. 2

Lyapunov exponents spectra of system (7) versus the control parameters with initial condition (8)

As it is explained in Sect. 2 for a two-degree-of-freedom Hamiltonian system, such as system (7), the Lyapunov exponents spectrum is symmetric about zero, i.e. $\Lambda =\{\lambda ,0,0,-\lambda \}$. Depending on chosen initial conditions, periodic and quasi-periodic motions manifest through all zero Lyapunov exponents, while for chaotic orbit we have $\lambda >0$. From the fundamental reasons, these systems cannot possess a hyperchaotic behaviour as the authors of [1] claimed. In fact, Figs. 2, 4 and 5 in [1] show that accuracy of numerical calculations is not sufficient because two Lyapunov exponents, which should vanish, differ significantly from zero.

In the next section, we give properly calculated Lyapunov exponents spectra of system (7) for the same values of the parameters and initial conditions as the authors of the work [1] have done.

4 Lyapunov exponents spectra and bifurcations diagrams for system (7)

Algorithms for computing the numerical values of Lyapunov’s spectra of continuous dynamical systems have been well established over the years [5, 14,15,16,17,18,19,20]. These different methods are mainly focused on the so-called standard algorithm introduced by Benettin et al. [5, 14]. This method is based on the integration of variational equations for n initial conditions with successive applications of the Gram–Schmidt orthonormalization procedure. This procedure can be easily implemented in Mathematica [21]. In this work, we use the standard algorithm with adopting a sufficient amount of k steps, so that the convergence of Lyapunov exponents is assured. We keep the relative and absolute error up to $10^{-11}$.

Figure 2 presents the Lyapunov exponent’s spectra of system (7) made for parameters (9) under initial condition (8). As we can notice, these plots are completely different comparing to the ones calculated in [1]. In contrary to [1], computed exponents obey the paring rule (5), as it should be for Hamiltonian system. Of course, there is no hyperchaotic behaviour. Instead of that, we can observe the coexistence of regular and chaotic behaviour depending on values of the control parameter.

From the Lyapunov exponents spectrum, we cannot deduce whether the regular pattern corresponds to periodic or whether to quasi-periodic motion. In order to detect periodic orbits, regular regions (windows) between chaotic layers in chaos, etc., a bifurcation diagram can be constructed. A bifurcation diagram gives insight to the dynamics of a considered system by plotting dependence a state variable with respect to a change of a control parameter [22].

In Fig. 3 of paper [1], the authors present the bifurcation diagram, by plotting the time series of variable $p_1(t)$ as a function of a certain control parameter. However, from this diagram we cannot deduce any conclusions since chaotic, quasi-periodic and periodic motions are not distinguishable. In order to get more useful information, we make the following. For given values of constant parameters $\beta =1,\theta =1$, and for initial condition (8), we effectively constituted the bifurcation diagram by sampling points periodically in the same way as the Poincaré sections with surface of sections chosen by $q_2 = -1$, with gradually increasing the control parameter $\eta \in (0,001,10)$. Figure 3a presents the comparison of the largest Lyapunov exponent with the bifurcation diagram. As the first sight, we can observe their very good agreement. From the bifurcation diagram, we are able to select these values of the control parameter for which system possesses chaotic, quasi-periodic and periodic behaviours. For instance, Fig. 3b shows the enlargement of Fig. 3a, which is depicted in colours values of the control parameter $\eta$ for which system possesses three different types of motions. For better understanding, particular orbits corresponding to these types of motion with marked points of intersections with Poincaré cross-sectional plane are presented in Fig. 4.

Fig. 3

The largest Lyapunov exponent and the bifurcation diagram versus $\eta$. The remaining parameters were chosen by $\beta =1, \theta =1$ with initial condition (8). Cross-sectional plain was defined by $q_2=-1$ with $p_2>0$

Fig. 4

Chaotic, quasi-periodic and periodic phase curves of system (7) for parameters: $\beta =1, \theta =1,$ with varying $\eta$. Points of intersections with Poincaré cross-sectional plane are depicted by small bullets. Cross-plain was specified as $q_2=-1$ with direction $p_2>0$, and initial condition was chosen by (8)

It is worth mentioning that in the case of Hamiltonian systems results for bifurcation diagrams and Lyapunov exponents depend strongly on the chosen initial condition. To show that, we present in Fig. 5 the largest Lyapunov exponent with its corresponding bifurcation diagram, computed for the same values of the parameters as Fig. 3 but for different initial conditions. As we can notice, this bifurcation diagram presents completely different structure. First of all, it is symmetric about zero and, moreover, in this diagram can observe the periodic “windows” between completely chaotic layers.

Fig. 5

The largest Lyapunov exponent and the bifurcation diagram versus $\eta$. The remaining parameters were chosen as: $\beta =1, \theta =1,$ with the initial condition: $(q_{10},q_{20} ,p_{10},p_{20})=(0,0,2,0.5)$. Cross-plain was specified by $q_2=0$ with $p_2>0$

5 Dynamics of restricted systems

In Section 3 of [1], authors consider Hamiltonian system (7) subjected to constraints. There are many approaches to describe such systems depending on types of constraints and other conditions. Authors of [1] have chosen approach described in [23]. Here, we present its simplified version when constrains are given just by one equation.

The dynamics of a natural Hamiltonian

$$\begin{aligned} \dot{q}=\frac{\partial H}{\partial p},\qquad \dot{p}=-\frac{\partial H}{\partial q}, \qquad q, p\in \mathbb {R}^n \end{aligned}$$

(12)

is restricted by a single constraint $\varphi (q,p)=0$ satisfying $\nabla _p \varphi (q,p)\ne 0$. The constrained Hamilton’s equations have the form

$$\begin{aligned} \dot{q}=\frac{\partial H}{\partial p},\qquad \dot{p}=-\frac{\partial H}{\partial q}+Q(q,p), \end{aligned}$$

(13)

where term Q(q, p) describes additional ‘forces’ caused by the constraint. Its explicit form derived in [23] is

$$\begin{aligned} Q(q,p,t) = H_{pp}^{-1}A_p^T (A_p H_{pp} A_p^T)^{-1}\left( A_p\cdot \frac{\partial H}{\partial q} +b_p\right) , \end{aligned}$$

(14)

where

$$\begin{aligned} A_p=\nabla _p\varphi ,\qquad b_p=-\nabla _q \varphi \cdot \frac{\partial H}{\partial p},\qquad H_{pp}=\frac{\partial ^2 H}{\partial p^2}. \end{aligned}$$

(15)

In the above formula, the gradient of a function is considered as a row vector, while the partial derivative as a column vector. For details, see [23], where also the more general cases with several constraints were considered.

As it is explicitly stated in [23], the constrained dynamics of the system takes place on the manifold

$$\begin{aligned} \varSigma =\left\{ (q,p)\in \mathbb {R}^{2n}\,|\, \varphi (q,p)=0\right\} . \end{aligned}$$

(16)

Thus, to study properties of the restricted system we have to consider only those solutions (q(t), p(t)) of constrained Hamilton’s equations (13) which lie on $\varSigma$. To ensure this, it is enough to require that $(q(t_0), p(t_0)\in \varSigma$ for a certain $t_0$. This is because we have the following.

Lemma 1

Function $\varphi (q,p)$ is a first integral of restricted system (13).

Proof

We have

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d} t }\varphi (q,p)= & {} \nabla _q\varphi \cdot \dot{q}+ \nabla _p\varphi \cdot \dot{p} \\= & {} \nabla _q\varphi \cdot \frac{\partial H}{\partial p} + \nabla _p\varphi \cdot \left( -\frac{\partial H}{\partial q} +Q \right) \\= & {} -b_p -A_p\cdot \frac{\partial H}{\partial q} + A_p\cdot Q. \end{aligned}$$

But using definition Q given by (14), we easily obtain that

$$\begin{aligned} A_p\cdot Q= b_p + A_p\cdot \frac{\partial H}{\partial q}. \end{aligned}$$

The above shows that the total time derivative of $\varphi (q,p)$ vanishes identically, so $\varphi (q,p)$ is a first integral. $\square$

Without an additional justification, a study of solutions of constrained Hamilton’s equations (13) with initial conditions $(q(t_0), p(t_0)\not \in \varSigma$, as it was made in [1], does not make big sense because these solutions are related to neither the constrained nor unconstrained dynamics of the system.

For Hamiltonian system (12) restricted by a holonomic constraint $\varphi _0(q)=0$, we define

$$\begin{aligned} \varphi (q,p):=\nabla _q\varphi _0(q)\cdot \dot{q}=\nabla _q\varphi _0(q)\cdot \frac{\partial H}{\partial p}. \end{aligned}$$

(17)

Then, by Lemma 1, we know that $\varphi (q,p)$ is a first integral of constrained system (13). Moreover, the total time derivative of $\varphi _0(q)$ is

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d} t }\varphi _0(q)=\nabla _q\varphi _0(q)\cdot \frac{\partial H}{\partial p}=\varphi (q,p). \end{aligned}$$

The above facts show that $2(n-1)$-dimensional manifold

$$\begin{aligned} \varPi :=\left\{ (q,p)\in \mathbb {R}^{2n}\,|\, \varphi _0(q)=0 \ \mathrm {and} \ \varphi (q,p)=0\right\} \subset \mathbb {R}^{2n} \end{aligned}$$

(18)

is invariant with respect to the phase flow generated by system (13). That is, if an initial condition $(q(t_0), p(t_0))$ belongs to $\varPi$, then the whole solution (q(t), p(t)) of the constrained Hamilton’s equations (13) belongs to $\varPi$. Moreover, only these solutions describe the dynamics of constrained system. If we introduce local parameterization of $\varPi$, then the constrained Hamilton’s equations are canonically mapped into Hamilton’s equations in $\mathbb {R}^{2n-2}$ with respect to a certain symplectic structure; for details, see Proposition 2.1 in [24].

5.1 Holonomic constraint

Following [1], we first consider Hamiltonian system (7) with the following holonomic constraint

$$\begin{aligned} \varphi _0(q_1,q_2)=q_1^2+q_2^2-L^2=0, \end{aligned}$$

(19)

for which, according to (17), we have

$$\begin{aligned} \varphi (q,p)=2q_1p_1+2q_2 p_2. \end{aligned}$$

(20)

Then, taking into account formulae (15), we obtain

$$\begin{aligned} H_{pp}=\mathbb {I}_2, \quad A_p=(q_1,q_2),\quad b_p=-(p_1^2+p_2^2), \end{aligned}$$

(21)

and the constrained Hamilton’s equations (13) take the form

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q}_1=p_1,\\ \dot{q}_2=p_2,\\ \dot{p}_1=-\beta q_1^3 - \eta q_2 + \eta ^2 q_1-\dfrac{q_1 P(q_1,q_2,p_1,p_2)}{q_1^2+q_2^2},\\ \dot{p}_2=-\beta q_2^3 - \eta q_1 -\dfrac{q_2 P(q_1,q_2,p_1,p_2)}{q_1^2+q_2^2}, \end{array}\right. } \end{aligned}$$

(22)

where

$$\begin{aligned} P=p_1^2+p_2^2-\beta (q_1^4+q_2^4)-2\eta q_1q_2+\eta ^2 q_1^2. \end{aligned}$$

The constrained dynamics takes place on the two-dimensional manifold $\varPi$, which is given by

$$\begin{aligned} \varPi =\left\{ (q_1,q_2,p_1,p_2)\in \mathbb {R}^{4}\,|\, q_1^2+q_2^2=L^2, \, q_1p_1 + q_2p_2=0 \right\} . \end{aligned}$$

On this manifold, the dynamics is described by solutions of an autonomous Hamiltonian system with one degree of freedom so it is integrable and this precludes any chaotic behaviour.

One can ask whether we can find a chaotic behaviour in the whole phase space of the constrained Hamilton’s equations (13) investigating its solutions with arbitrary initial conditions $(q_1,q_2,p_1,p_2)\in \mathbb {R}^4$. Numerical experiments show that for arbitrary solution of system (22) all Lyapunov exponents vanish. To understand why it is so, we introduce polar coordinates $(r,\vartheta )$ and the corresponding momenta $(p_r,p_{\vartheta })$ in the standard way

$$\begin{aligned}&q_1=r\cos \vartheta ,\quad q_2=r\sin \vartheta ,\\&p_1=\cos \vartheta p_r-\sin \vartheta \frac{p_{\vartheta }}{r},\quad p_2=\sin \vartheta p_r+\cos \vartheta \frac{p_{\vartheta }}{r}. \end{aligned}$$

In these coordinates, the restricted system (22) takes the form

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{r}=p_r,\\ \dot{\vartheta }=\dfrac{p_{\vartheta }}{r^2},\\ \dot{p}_r=-\dfrac{p_r^2}{r},\\ \dot{p}_{\vartheta }=\dfrac{r^2}{4}\left[ \beta r^2\sin (4\vartheta )-2\eta (2\cos (2\vartheta )+\eta \sin (2\vartheta ))\right] . \end{array}\right. } \end{aligned}$$

(23)

Equations for r and $p_r$ form the closed subsystem and constraints take the forms

$$\begin{aligned} \varphi _0=r^2-L^2,\qquad \varphi =2r p_r. \end{aligned}$$

(24)

Thus, the invariant manifold $\varPi$ is given by

$$\begin{aligned} \varPi = \{ (r, \vartheta , p_r, p_{\vartheta }) \, |\, p_r=0, \ r = L \}. \end{aligned}$$

(25)

On this manifold, system (23) reduces to the system with one degree of freedom

$$\begin{aligned}&\dot{\vartheta }=\dfrac{p_{\vartheta }}{L^2},\nonumber \\&\dot{p}_{\vartheta }= \frac{L^2}{4}\left[ \beta L^2\sin (4\vartheta )-2\eta (2\cos (2\vartheta )+ \eta \sin (2\vartheta ))\right] ,\nonumber \\ \end{aligned}$$

(26)

which always has the energy first integral. In this case, all Lyapunov exponents therefore must be equal zero, so the motion is periodic or quasi-periodic but cannot be chaotic.

System (26) is Hamiltonian with the following Hamiltonian function

$$\begin{aligned} \widetilde{H}= \frac{p_{\vartheta }^2}{2L^2}- \frac{\eta L^2}{4}\left[ \eta \cos (2\vartheta )-2\sin (2\vartheta )\right] + \frac{L^4\beta }{16}\cos (4\vartheta ). \end{aligned}$$

(27)

We can obtain the same result using the standard procedure. In fact, first we take the Lagrange function corresponding to the Hamiltonian (6). Then, as local coordinates on the constraint, we can choose the polar angle $\vartheta$ by setting $q_1=L\cos \vartheta$ and $q_2=L\sin \vartheta$. As a result, we obtain the Lagrange function corresponding to system (6), which is

$$\begin{aligned} \mathscr {L}= & {} \frac{L^2}{2}\dot{\vartheta }^2+\frac{\eta L^2}{4}\left[ \eta \cos (2\vartheta )-2\sin (2\vartheta )\right] \\&-\frac{L^4\beta }{16}\cos (4\vartheta )+C_0, \end{aligned}$$

where $C_0=\tfrac{L^2}{16}(4\eta ^2+8\theta -3\beta L^2)$ is irrelevant constant. The canonical momentum conjugated with $\vartheta$ is $p_{\vartheta }=\partial L/\partial \dot{\vartheta }=L^2\dot{\vartheta }$. Neglecting $C_0$ the corresponding Hamiltonian is exactly (27).

Fact that we obtained completely different results than those in [1] needs an explanation. First of all in the cited paper, the authors do not investigate the constrained Hamilton’s equations (22) but their modification

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q}_1=p_1,\\ \dot{q}_2=p_2,\\ \dot{p}_1=-\beta q_1^3 - \eta q_2 + \eta ^2 q_1-\dfrac{q_1 P_L(q_1,q_2,p_1,p_2)}{L^2},\\ \dot{p}_2=-\beta q_2^3 - \eta q_1 -\dfrac{q_2P_L(q_1,q_2,p_1,p_2)}{L^2}, \end{array}\right. } \end{aligned}$$

(28)

where

$$\begin{aligned} P_L=p_1^2+p_2^2-\beta (L^4-2q_1^2q_2^2)-2\eta q_1q_2+\eta ^2 q_1^2, \end{aligned}$$

and L is considered as a parameter. Although this system coincides with (22) on manifold $\varPi$, it has completely different global properties. For instance, function (17) is a first integral of equations (22), but it is not a first integral of (28).

Fig. 6

Lyapunov exponents spectrum of system (28) versus the control parameter $\beta$ with initial condition (8). The remaining parameters were specifies as: $\eta =0.01$, and $L=1$

As system (28) has time-reversal symmetry $R(\varvec{q},\varvec{p})=(\varvec{q},-\varvec{p})$, the Lyapunov exponents obey the paring rule (5), and Lyapunov’s exponents spectrum is symmetric about zero: $\{\lambda ,0,0-\lambda \}$ with $\lambda >0$. Thus, in system (28) the hyperchaotic behaviours are precluded. To confirm this, we present in Fig. 6 the Lyapunov exponents spectrum of system (28) made for the parameters $\eta =0.01, L=1, \beta \in [0,10]$ and initial condition (8), which clearly not lie on $\varPi$. This spectrum was made for the same values of the parameters and the initial condition as the spectrum depicted in Fig. 7 in [1], where the authors have found hyperchaos for the entire range of parameter $\beta$. As we can notice, Fig. 6 is symmetric about zero, and the Lyapunov exponents obey the parting rule, as it should be.

5.2 Non-holonomic constraint

As the second example, we follow [1], where the authors considered Hamiltonian system (7) subjected to a non-holonomic constraint of the form

$$\begin{aligned} \varphi (q,p)=p_1^2+p_2^2-l^2=0. \end{aligned}$$

(29)

For this constraint, formulae (15) give

$$\begin{aligned} H_{pp}=\mathbb {I}_2,\quad A_{p}=(p_1,p_2),\quad b_p=0, \end{aligned}$$

(30)

where $\mathbb {I}_2$ is two-dimensional identity matrix. Thus, constrained system (13) looks as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q}_1=p_1,\\ \dot{q}_2=p_2,\\ \dot{p}_1=-\beta q_1^3 - \eta q_2 + (\eta ^2+\theta ) q_1+\dfrac{p_1 S(q_1,q_2,p_1,p_2)}{p_1^2+p_2^2},\\ \dot{p}_2=-\beta q_2^3 - \eta q_1 +\theta q_2+\dfrac{p_2 S(q_1,q_2,p_1,p_2)}{p_1^2+p_2^2}, \end{array}\right. } \end{aligned}$$

(31)

where

$$\begin{aligned} S=\left( \beta q_1^3 + \eta q_2- (\eta ^2 + \theta )q_1\right) p_1+ \left( \beta q_2^3 + \eta q_1 - \theta q_2\right) p_2. \end{aligned}$$

This system has a quadratic first integral of the form

$$\begin{aligned} \varphi (q,p)=p_1^2+p_2^2. \end{aligned}$$

(32)

Moreover, the divergence of system (31) is

$$\begin{aligned} {\text {div}}(w)=\dfrac{S}{p_1^2+p_2^2}. \end{aligned}$$

Fig. 7

Lyapunov exponents spectrum and the bifurcation diagram of system (33) versus the control parameter $\beta$ with cross-plain $p_2=0$. The remaining parameters were specified as: $\theta =1, \eta =1, l=10$, and the initial condition was chosen by (8)

In [1], the authors did not investigate the truly constrained Hamiltonian system (31) under constraint (29) but its modification

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q}_1=p_1,\\ \dot{q}_2=p_2,\\ \dot{p}_1=-\beta q_1^3 - \eta q_2 + (\eta ^2+\theta ) q_1+\dfrac{p_1 S(q_1,q_2,p_1,p_2)}{l^2},\\ \dot{p}_2=-\beta q_2^3 - \eta q_1 +\theta q_2+\dfrac{p_2 S(q_1,q_2,p_1,p_2)}{l^2}, \end{array}\right. } \end{aligned}$$

(33)

where l was fixed and considered as a control parameter. Again, although this system agrees with the one given in Eq. (31) on constraint $p_1^2+p_2^2=l^2$, it has very different global properties. For instance, the Lie derivative of function (32) along the vector field (33) is given by

$$\begin{aligned} \frac{\mathrm {d} }{\mathrm {d}t}\varphi (q,p)=2\left( p_1^2+p_2^2-l^2\right) \frac{S}{l^2}, \end{aligned}$$

(34)

so it vanishes only on constraint. In Figs. 12–17, of the work [1] the authors have performed the numerical analysis of their “constrained” system (33) with initial conditions not satisfying the constrained condition $p_1^2(t_0)+p_2^2(t_0)-l^2=0$. They have found certain values of the parameters for which system (33) is hyperchaotic and, moreover, it has three positive Lyapunov exponents. Of course, it is not true and their wrong conclusion was probably caused by not sufficient accuracy of numerical calculations. Since system (33) possesses time-reversal symmetry $R(\varvec{q},\varvec{p})=(\varvec{q},-\varvec{p})$, the Lyapunov exponents again obey the paring rule. Thus, in system (33) the hyperchaotic behaviours are precluded because only one Lyapunov exponent can be positive. To confirm this, we present in Fig. 7 the Lyapunov exponents spectrum and its corresponding bifurcation diagram of system (33). They were made for the same values of the parameters and the initial condition as Fig. 12 of in [1]. Since system (31) possesses first integral (32), we can reduce its dimension by one. For this purpose, we introduce new variables in the following way

$$\begin{aligned} q_1= q_1,\quad q_2=q_2,\quad p_1=p \cos \vartheta ,\quad p_2=p \sin \vartheta . \end{aligned}$$

(35)

In these variables, constrained system (31) reads

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q}_1=p\cos \vartheta ,\\ \dot{q}_2=p\sin \vartheta ,\\ \dot{\vartheta }=\dfrac{1}{p}\left[ \left( \eta q_2-(\eta ^2+\theta )q_1+\beta q_1^3\right) \sin \vartheta \right. \\ \left. -\left( \eta q_1-\theta q_2+\beta q_2^3\right) \cos \vartheta \right] ,\\ \dot{p}=0. \end{array}\right. } \end{aligned}$$

(36)

As we can see $p={\text {const}}$, and thus, we may consider it as a parameter $p=l$. Hence, the final form of the reduced, constrained system (31) is given by

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{q}_1=l\cos \vartheta ,\\ \dot{q}_2=l\sin \vartheta ,\\ \dot{\vartheta }=\dfrac{1}{l}\left[ \left( \eta q_2-(\eta ^2+\theta )q_1+\beta q_1^3\right) \sin \vartheta \right. \\ \left. -\left( \eta q_1-\theta q_2+\beta q_2^3\right) \cos \vartheta \right] . \end{array}\right. } \end{aligned}$$

(37)

It has time-reversal symmetry given by $R(q_1,q_2,\vartheta )=(-q_1,-q_2,\vartheta )$. In Figs. 8 and 9, the Poincaré cross section, the bifurcation diagram and its corresponding Lyapunov spectrum are visible. As we can notice, for the chosen values of the parameters reduced system (37) still possesses very rich and complex dynamics. However, only one Lyapunov exponent $\lambda >0$, the second by time-reversibility symmetry is $-\lambda$ and the third vanishes.

Fig. 8

Poincaré cross section of system (37) with cross-plain $\vartheta =0$. The parameters were specified as: $\beta =1,\theta =-1, \eta =1/100, l=1$

Fig. 9

Lyapunov exponents spectrum and the bifurcation diagram of system (37) versus the control parameter $\theta$ with cross-plain $\vartheta =0$. The remaining parameters and the initial condition were specified as: $\beta =1, \eta =1/100, l=1$, $(q_{10},q_{20},\vartheta _0)=(0,2,0)$

6 Hyperchaos in pendulums systems

As we have already mentioned, the hyperchaotic behaviours in Hamiltonian systems can appear only in systems with at least three degrees of freedom. Multiple pendulums seem to be natural candidates for such systems with many degrees of freedom. In this short section, we will search for the hyperchaotic behaviours of triple and quartic pendulums, which move under the influence of the constant gravity field. For this purpose, we will compute full Lyapunov’s exponents spectra and the bifurcations diagrams corresponding to them.

6.1 Simple triple pendulum

Simple triple pendulum is composed of three simple pendulums attached to each other. The system has three degrees of freedom, and it is fully described by three angles $\varphi _1,\varphi _2,\varphi _3$. Using the notations given in Fig. 10, we can write the corresponding Lagrange function in the form

$$\begin{aligned} \mathscr {L}&=\frac{1}{2}\sum _{i=1}^3\mu _{i3}l_i(l_i\dot{\phi }_i^2+2g \cos \phi _i)\nonumber \\&\quad +\sum _{\begin{array}{c} j=2, \\ j\ne i \end{array}}^3\sum _{i=1}^2\mu _{j3}l_il_j\dot{\phi }_i\dot{\phi }_j\cos (\phi _j-\phi _i), \end{aligned}$$

(38)

where for simplicity we have introduced

$$\begin{aligned} \mu _{nk}:=\left( \sum _{i=n}^k m_i\right) . \end{aligned}$$

Figure 11 presents the Lyapunov exponents spectrum and the bifurcation diagram for the simple triple pendulum system defined by Lagrange function (38). This diagram plots extremal values (amplitudes) of the system variable $\phi _{3}(t)=\phi _{3,\text {extr}}$, when $\dot{\phi }_3=p_3=0$. For values of the parameters $m_i=l_i=g=1$ with the initial condition: $(\phi _i=\phi _0,p_i=0)$ for $i=1,2,3$, we show how the change of initial angle of the pendulums deviation $\phi _{0}\in [0,\pi ]$ affects the dynamics of the system. As one would expect, for small angles the pendulums oscillate near the equilibrium. However, when the angular displacement amplitude of the pendulums is large enough, the regular pattern diverges and the motion finally becomes highly complex. The corresponding Lyapunov exponents spectrum indicates that motion is, in fact, hyperchaotic with two positive exponents.

Fig. 10

Geometry of simple triple pendulum

Fig. 11

Lyapunov exponents spectrum and the bifurcation diagram versus the initial amplitude $\phi _0$, for the simple triple pendulum. The parameters were chosen as: $m_i=l_i=g=1$ with the initial condition: $(\phi _i(0)=\phi _0,p_i(0)=0)$ for $i=1,2,3$. Here $\phi _{3,\text {extr}}$ denotes the extremal values of the system variable $\phi _3$, when $\dot{\phi }_3=p_3=0$

6.2 Flail triple pendulum

Let us now look for the hyperchaotic nature of the flail triple pendulum moving under the influence of the gravitational force. The geometry of the system is presented in Fig. 12. It consists of three simple pendulums with masses $m_1, m_2, m_3$ and lengths $l_1, l_2, l_3$. The first pendulum is attached to a fixed point, and to its end mass, the other two pendulums are joined. This system, which was in the absence of the gravity field intensively studied in [25], has the following Lagrange function:

$$\begin{aligned} \mathscr {L}= & {} \sum _{{\begin{array}{c} i=1,j=2, \\ i+j\in \mathbb {Z}_\text {even} \end{array}}}^3 \mu _{ij}l_i\left( l_i\dot{\phi }_i^2+2g\cos \phi _i\right) \nonumber \\&+\, l_1\dot{\phi }_1\sum _{i=2}^3m_il_i\dot{\phi }_i\cos (\phi _1-\phi _i). \end{aligned}$$

(39)

Figure 13 presents the Lyapunov exponents spectrum and its corresponding bifurcation diagram versus the initial angle of deviation $\phi _0\in [0,\pi ]$. These plots were made for fixed values of the parameters $m_i=l_i=g=1$ with the initial condition:$(\phi _i(0),p_i(0))=(0,\phi _0,-\phi _0,0,0,0)$ for $i=1,2,3$. As in the case of the simple triple pendulum, $\phi _{3,\text {extr}}$ states for the extremal values of the system variable $\phi _3$, when $\dot{\phi }_3=p_3=0$.

As expected, for small values of $\phi _0$, the pendulums oscillate near the equilibrium, and we can conclude that for such small values of $\phi _0$ the motion is almost integrable. Nonetheless, for $\phi _0\approx \pi /2$ the motion starts to be more complex and finally becomes chaotic and later hyperchaotic, which is evidenced in Fig. 13. Surprisingly enough, over the range $\phi _0\in [2.246, 2.293]$, we can still detect stable solutions localized in the gap between completely hyperchaotic regions. This corresponds to the case when the first pendulum is at the rest, and other two pendulums oscillate in the opposite directions, as shown in Fig. 14 presenting the time series of the angular positions of each mass of flail triple pendulum (39) with initial condition: $(\phi _i(0),p_i(0))=(0,2.25,-2.25,0,0,0)$.

Fig. 12

Geometry of flail-like triple pendulum

Fig. 13

Lyapunov exponents spectrum and the bifurcation diagram versus the initial amplitude $\phi _0$ for the flail pendulum. The parameters: $m_i=l_i=g=1$ with i.c.: $(\phi _i(0),p_i(0))=(0,\phi _0,-\phi _0,0,0,0)$. Here, $\phi _{3,\text {extr}}$ denotes the extremal values of $\phi _3$, when $\dot{\phi }_3=p_3=0$

Fig. 14

Time series with i.c.: $(\phi _i(0),p_i(0))=(0,2.25,-2.25,0,0,0)$

6.3 Quartic simple pendulum

Let us finally consider the simple quartic pendulum. This system is composed of four simple pendula of masses $m_i$ and lengths $l_i$ which are connected to each other and move under the influence of the gravity, as shown in Fig. 15. The Lagrange function of this model is as follows:

$$\begin{aligned} \mathscr {L}= & {} \frac{1}{2}\sum _{i=1}^4\mu _{i4}l_i(l_i\dot{\phi }_i^2+2g\cos \phi _i)\nonumber \\&+\sum _{\begin{array}{c} j=2, \\ j\ne i \end{array}}^4\sum _{i=1}^3\mu _{j 4}l_il_j\dot{\phi }_i\dot{\phi }_j\cos (\phi _j-\phi _i). \end{aligned}$$

(40)

This system is hyperchaotic with three positive Lyapunov exponents, provided the common initial amplitude $\phi _0$ is large enough. For instance, for $m_i=l_i=g=1$ with the initial condition: $(\phi _i(0)=\phi _0,p_i(0)=0)$, the motion starts to be hyperchaotic for $\phi _0>\pi /3$, as shown in Fig. 16.

Fig. 15

Geometry of simple quartic pendulum

Fig. 16

Lyapunov exponents spectrum and the bifurcation diagram versus the initial amplitude $\phi _0$ for the simple quartic pendulum. The parameters were chosen as: $m_i=l_i=g=1$ with the initial condition: $(\phi _i(0)=\phi _0,p_i(0)=0)$ for $i=1,2,3,4$. The cross-sectional plane was specified as: $\phi _1=0$ with $p_1>0$