# A Novel Megastable Hamiltonian System with Infinite Hyperbolic and Nonhyperbolic Equilibria

Authors: Gervais Dolvis Leutcho; Theophile Fonzin Fozin; Alexis Nguomkam Negou; Zeric Tabekoueng Njitacke; Viet-Thanh Pham; Jacques Kengne; Sajad Jafari
Journal: Complexity
 (2020)

---

## Abstract

The dimension of the conservative chaotic systems is an integer and equals the system dimension, which brings about a better ergodic property and thus have potentials in engineering application than the dissipative systems. This paper investigates the phenomenon of megastability in a unique and simple conservative oscillator with infinite of hyperbolic and nonhyperbolic equilibria. Using traditional nonlinear analysis tools, we found that the introduced oscillator possesses an invariable energy and displays either self-excited or hidden dynamics depending on the stability of its equilibria. Besides, the conservative nature of the new system is validated using theoretical measurement. Furthermore, an analog simulator of the oscillator is built and simulated in the PSpice environment to confirm that the previous results were not artifacts.


---

## Body

## 1. Introduction

Considerable efforts have been made recently by several authors in designing chaotic systems with striking features linked to their equilibria [1–7]. These dynamical systems include some with no equilibrium point [8], stable equilibria [1, 9–11], infinite equilibria [12–17], and unstable equilibria [2, 18, 19]. The latter belongs to self-excited nonlinear dynamical systems, while the formers are new and classified as systems with hidden attractors [20–22]. Several systems with hidden attractors are linked to multistability which refers to the coexistence of multiple attractors for the same group of parameters [23–25]. Indeed, a system may approach different stable states starting from different initial values. When the system basin of initial conditions is associated with unstable equilibrium points, one may find countable coexisting attractors. However, when the system has a great deal of the coexisting attractors, the phenomenon is defined as extreme multistability. Despite the intensive efforts in reporting systems with hidden attractors, most of the attention so far was oriented only on structure and characteristic of their equilibria [26], while other important specificities including shape and topology of the strange attractors were rarely considered. Few examples are chaotic systems with varying symmetry [27], chaotic systems with conditional symmetry [28–31], multiscroll attractors [32–34], attractor growing [35, 36], algebraically simple equations [37], and nested megastable attractors [38–45].Now in forefront of investigations, megastable systems are a special class of multistable systems with nested attractors [26]. First reported by Sprott and coauthors in [26], such systems are now widely investigated [38–45]. Although the phenomenon of megastability in chaotic systems was up to date linked to additional external excitation [46, 47], the recent work by Kathikeyan and colleagues demonstrates that the above condition is not mandatory [48]. However, very few megastable systems reported so far are conservative. Conservative chaotic flow is a category of incompressible chaos usually generated from dynamical systems without dissipation. Based on a long history and formalism from Euler, Lagrange Hamilton, and Jacobi [49], the main characteristic of such systems is that they do not have attractors because their phase space volumes are conserved [50]. However, such systems can exhibit complicated dynamical behaviors including chaos [37]. The dimension of the conservative chaotic systems is an integer and equals the system dimension, which brings about a better ergodic property than the dissipative system. Consequently, the probability distribution associated with conservative chaos is relatively flat like the uniform distribution of white noise, whereas dissipative chaos has either a single peak or multiple peaks, from which the main frequency characteristics can be identified. Clearly, the information encrypted based on the former can be decrypted or attacked based on frequency reconstruction [51] while the latter (conservative system) is adapted [37]. Therefore, with both having the same bandwidths, the conservative chaotic system is more suitable as a pseudorandom number generator than the dissipative chaotic system. Owing to the importance of conservative/Hamiltonian systems in the mechanical domain, we propose in this work an extremely simple megastable Hamiltonian system with its circuitry implementation.It is observed from Table1 that no 2D conservative system with infinite hyperbolic and nonhyperbolic equilibria and few terms (i.e., three) was reported in the literature. Also, Table 2 highlights the megastability (i.e., coexistence of infinite countable attractors) property of some simple conservative systems proposed within this work. To the best of the authors’ knowledge, no such simple conservative 2D system with both megastability and multistability was yet reported in the relevant literature so far and thus deserves dissemination.Table 1
Categorization of 2D chaotic megastable systems with number of terms and many equilibria.
No.Types of 2D megastable systemNumber of termsShape of equilibriaReference of papers12D dissipative systemThreeInfinite of equilibria[41–43, 45]22D dissipative systemThreeOne equilibrium point[26]32D dissipative systemThreeOne equilibrium point[38]42D dissipative systemFourInfinite of equilibria[44]52D dissipative systemFiveOne equilibrium point[39]62D dissipative & conservativeSixNone[40]72D conservative systemThreeInfinite of hyperbolic and nonhyperbolic equilibriaThis workTable 2
Categorization of reported conservative chaotic systems with various characteristics.
No.Types of conservative systemsNature/name of systemMultistabilityMegastabilityReference of papers13D systemDissipative and conservativeYes—coexistence of two different dynamical behaviorsNone[52–54]23D systemNose-hooverNoneNone[37]33D systemOther conservative systemsNoneNone[55, 56]44D systemOther conservative systemNoneNone[57]54D systemNonhyperbolic equilibriaYes—coexistence of two different attractorsNone[58]62D systemInfinite of hyperbolic and nonhyperbolic equilibriaYes—infinite attractorsYes—coexistence of infinite number of nested attractorsThis workThe organization of this paper is as follows. Section2 deals with the introduction and description of the new conservative megastable system. Preliminary analysis including investigations of steady states and their nature is presented. Then, Section 3 is dedicated to the conservative feature/property of the new system with two different methods, namely, divergence and Hamiltonian. In Section 4, numerical investigations of the forced derivative of the newly introduced system are performed based on well-known nonlinear dynamical tools including bifurcation diagrams, Lyapunov exponents, and phase portraits. In Section 5, an experimental study is carried out. The results of the corresponding electronic circuit implemented in the PSpice simulation environment show a good agreement with obtained numerical investigations. The last section concludes the paper.


## 2. The New Simple Conservative Megastable Oscillator

Consider the following two-dimensional (2D) dissipative nonlinear oscillator which is investigated in [59]:(1)x˙=y,u˙=−x+ycosx.Inspired by system (1), we introduce the following 2D new and simple conservative nonlinear system:(2)x˙=−y,y˙=sinax+ksinbx.In this paper, the most following elegant [37] set of parameters (a=2/2, b=1/20, and k=2) are used in the detail investigation of the new conservative megastable system in order to uncover exciting and interesting dynamical behaviors. The equilibrium point or fixed points of system (2) can be calculated by equating to zero the right hand side of system (2). Therefore, we obtained the following transcendental equation:(3)ϕx=sinax+ksinbx=0.The graphical representation of functionϕx in the interval [−140, 140] is shown in Figure 1. According to this graph, one can see that there are infinite values of x which satisfy equation (3) when the range of the graphical representation is increased. Therefore, the new conservative system (2) has infinite fixed points, that is, Pixi,0 with i=1,2,3,…,∞.Figure 1
Graphical representation of the functionϕx showing some views of the different fixed points of system (1) in the range x∈−140,140.The Jacobian matrix for this system in any equilibrium points is(4)J=0−1−acosax˜+kbcosbx˜0x,y=x˜,0⟶0−1−acosax˜+kbcosbx˜0.Then, the eigenvalues can be calculated accordingly as(5)λI−J=0⟶λ−1−acosax˜+kbcosbx˜λ=0⟶λ2+acosax˜−kbcosbx˜=0,λ1,2=∓kbcosbx˜−acosax˜.Forkbcosbx˜−acosax˜>0, system (2) has infinite hyperbolic equilibria:(5a)λ1,2=∓kbcosbx˜−acosax˜.Forkbcosbx˜−acosax˜<0, system (2) has infinite nonhyperbolic equilibria:(5b)λ1,2=∓ikbcosbx˜−acosax˜.One can observe that the eigenvalues of the characteristic equation strongly depend on the fixed pointsPix˜,0. Table 3 shows an example of the two eigenvalues of equation (5) for some equilibrium points, that is, P1, P2, P3, P4, P5, and P6. Based on the analysis provided in Table 3, we concluded that the points P1,P2, and P4 are examples of hyperbolic fixed points, while P3,P5, and P6 are examples of nonhyperbolic fixed points. Using the previous parameter set, we carried out the numerical calculation of system (2). Some interesting and rare trajectories in this system are detected (see Figure 2). In particular, Figure 2 represents a plot of fourteen coexisting nested limit cycles in system (2). Remember that they are just an example from the infinite set of nested limit cycles around the x-axis and can be obtained under other proper initial conditions.Table 3
Six examples of equilibrium pointsP and their corresponding eigenvalues for parameters a=2/2, b=1/20, and k=2.0.
Equilibrium pointsPiEigenvaluesλ1,2Hyperbolic or nonhyperbolicP10,0λ1,2=±0.7791HyperbolicP2−7.6866,0λ1,2=±0.6124HyperbolicP358.4009,0λ1,2=±0.7332iNonhyperbolicP472.9614,0λ1,2=±0.5086HyperbolicP5−5.2078,0λ1,2=±0.8383iNonhyperbolicP6−198.2219,0λ1,2=±0.4035iNonhyperbolicFigure 2
Fourteen examples of limit cycles of system (2) for initial conditions 1,0, 2,0, 3,0, 4,0, 5,0, ±5.5,0, ±6,0, ±7,0, 10,0, 11,0, and 12,0.


## 3. Conservative Nature of the New Oscillator

### 3.1. Conservative Using Divergence Property

The nonlinear dynamical system (2) can be expressed in vector notation as(6)Y˙=fY=f1x,y,f2x,y,where Y∈R2 is a state vector and fY is the smooth function given by(7)f1x,y=−y,f2x,y=sinax+ksinbx.The divergence of the vector fielddivf=∇⋅f for nonlinear system (6) can be calculated as(8)divf=∇⋅f=∂f1∂x+∂f2∂y.However, system (2) is dissipative if ∇⋅f<0. Otherwise, system (2) is conservative.Thus, the divergence of system (2) can be easily obtained using (8):(9)divf=∇⋅f=∂y∂x+∂sinax+ksinbx∂y=0.As a result, the megastable system (2) has a conservative nature.Next, suppose that system (2) is expressed as X˙=hX. Let Ω be any space in R2 with a smooth boundary/surface and also ϕfΩ=Ωt, where ϕf is the flow of hX. Also, consider the volume element Vt of Ωt in the phase space. The volume element into a new smooth boundary Ωt+dt in an infinitesimal time dt is given by Vt+dt and can be obtained accordingly as [60](10)Vt+dt=Vt+∫h⋅ndtdS,where S represents the area of the smooth boundary and n is the outward normal apply on the surface Ω. It is clear from (10) that(11)dVdt=V˙t=Vt+dt−Vtdt=∫Ωh⋅ndS.Thus, by using the Ostrogradsky theorem (i.e., divergence theorem) in (11), we get(12)dVdt=V˙t=∰VΔ⋅hdV,where Δ⋅h is the divergence of h.Since the divergence of system (2) is zero, we can write (11) as(13)dVdt=V˙t=0.From (13), we get(14)Vt=V0.It is observed from (14) that the volume in the phase space of system (2) is a constant. So, system (2) has conservative behavior.


### 3.2. Conservative Nature Using Hamiltonian-Based Analysis

The electrical activity of the systems is related to their energy modification [44]. Recently, the Hamiltonian energy function-based method is employed for the study of chaotic systems. Some technics to find the Hamilton energy function in chaotic systems have been presented [61]. In particular, a general procedure for investigation of generalized Hamiltonian realization of an n-dimensional autonomous dynamical system is reported [62].Consider the following dynamical system:(15)x˙=Fx,where x∈Rn and F:U⟶Rn is a smooth function so that U⊆Rn. The dynamical system (15) can be written in a generalized Hamiltonian form [63, 64] as(16)x˙=Mx∇H,where Mx is the local structure matrix and ∇H is the gradient vector of a smooth energy function Hx. For Hamiltonian system, Mx is a skew-symmetric matrix which satisfies the Jacobian identity. Otherwise, where Mx is not skew-symmetric, the system is generalized Hamiltonian. However, for a generalized Hamiltonian system, Mx can be divided as a sum of one skew-symmetric matrix Jx plus one symmetric matrix Rx. Therefore, equation (16) can be written accordingly as(17)x˙=Jx+Rx∇H.As a result,Jx and Rx can be obtained as follows:(18)H˙=dHdt=∇HTRx∇H,∇HTJx∇H=0.According to Helmholtz’s theorem, system (15) can be regarded as the vector field to discuss the energy problems and can be decomposed into the rotational tensor Fcx and gradient vector Fdx, that is,(19)Fx=Fcx+Fdx,thus(20)H˙=∇HTFdx,Jx∇HT=Fcx.ForJx to be a skew-symmetric matrix,(21)∇HTFcx=0,where Hx is the Hamiltonian function and ∇H is the gradient vector of Hx defined by(22)∇H=∂H∂x,∂H∂yT.The new 2D conservative autonomous system (2) can be rewritten as(23)x˙y˙=−ysinax+ksinbx+00,where(24)Fcx=−ysinax+ksinbx,Fdx=00.SinceFdx=0,0T, it means that there is no dissipation term in system (2) and trajectories are corresponding to the unique conservative field Fcx. That is, when the trajectories begin from an isosurface of the constant Hamilton energy Hx, they cannot escape from the isosurface [57, 58].According to (22), the Hamilton energy function Hx,y will satisfy the partial differential equation given by(25)y∂H∂x−sinax+ksinbx∂H∂y=0.So, a general solution of (25) can be approached by(26)Hx,y,a,b,k=12y2−1acosax−kbcosbx.The time derivative of (26) can be calculated by(27)H˙x,y,a,b,k=y.y˙+sinax+ksinbx⋅x˙=ysinax+ksinbx−ysinax+ksinbx=0=∇HTfdx.Since the time derivative of Hamiltonian energy of system (2) is zero, this Hamilton energy function is invariable. This justifies the conservative nature of the new system (2). Thus, the conservative property of system (2) is validated using theoretical measures.


## 3.1. Conservative Using Divergence Property

The nonlinear dynamical system (2) can be expressed in vector notation as(6)Y˙=fY=f1x,y,f2x,y,where Y∈R2 is a state vector and fY is the smooth function given by(7)f1x,y=−y,f2x,y=sinax+ksinbx.The divergence of the vector fielddivf=∇⋅f for nonlinear system (6) can be calculated as(8)divf=∇⋅f=∂f1∂x+∂f2∂y.However, system (2) is dissipative if ∇⋅f<0. Otherwise, system (2) is conservative.Thus, the divergence of system (2) can be easily obtained using (8):(9)divf=∇⋅f=∂y∂x+∂sinax+ksinbx∂y=0.As a result, the megastable system (2) has a conservative nature.Next, suppose that system (2) is expressed as X˙=hX. Let Ω be any space in R2 with a smooth boundary/surface and also ϕfΩ=Ωt, where ϕf is the flow of hX. Also, consider the volume element Vt of Ωt in the phase space. The volume element into a new smooth boundary Ωt+dt in an infinitesimal time dt is given by Vt+dt and can be obtained accordingly as [60](10)Vt+dt=Vt+∫h⋅ndtdS,where S represents the area of the smooth boundary and n is the outward normal apply on the surface Ω. It is clear from (10) that(11)dVdt=V˙t=Vt+dt−Vtdt=∫Ωh⋅ndS.Thus, by using the Ostrogradsky theorem (i.e., divergence theorem) in (11), we get(12)dVdt=V˙t=∰VΔ⋅hdV,where Δ⋅h is the divergence of h.Since the divergence of system (2) is zero, we can write (11) as(13)dVdt=V˙t=0.From (13), we get(14)Vt=V0.It is observed from (14) that the volume in the phase space of system (2) is a constant. So, system (2) has conservative behavior.


## 3.2. Conservative Nature Using Hamiltonian-Based Analysis

The electrical activity of the systems is related to their energy modification [44]. Recently, the Hamiltonian energy function-based method is employed for the study of chaotic systems. Some technics to find the Hamilton energy function in chaotic systems have been presented [61]. In particular, a general procedure for investigation of generalized Hamiltonian realization of an n-dimensional autonomous dynamical system is reported [62].Consider the following dynamical system:(15)x˙=Fx,where x∈Rn and F:U⟶Rn is a smooth function so that U⊆Rn. The dynamical system (15) can be written in a generalized Hamiltonian form [63, 64] as(16)x˙=Mx∇H,where Mx is the local structure matrix and ∇H is the gradient vector of a smooth energy function Hx. For Hamiltonian system, Mx is a skew-symmetric matrix which satisfies the Jacobian identity. Otherwise, where Mx is not skew-symmetric, the system is generalized Hamiltonian. However, for a generalized Hamiltonian system, Mx can be divided as a sum of one skew-symmetric matrix Jx plus one symmetric matrix Rx. Therefore, equation (16) can be written accordingly as(17)x˙=Jx+Rx∇H.As a result,Jx and Rx can be obtained as follows:(18)H˙=dHdt=∇HTRx∇H,∇HTJx∇H=0.According to Helmholtz’s theorem, system (15) can be regarded as the vector field to discuss the energy problems and can be decomposed into the rotational tensor Fcx and gradient vector Fdx, that is,(19)Fx=Fcx+Fdx,thus(20)H˙=∇HTFdx,Jx∇HT=Fcx.ForJx to be a skew-symmetric matrix,(21)∇HTFcx=0,where Hx is the Hamiltonian function and ∇H is the gradient vector of Hx defined by(22)∇H=∂H∂x,∂H∂yT.The new 2D conservative autonomous system (2) can be rewritten as(23)x˙y˙=−ysinax+ksinbx+00,where(24)Fcx=−ysinax+ksinbx,Fdx=00.SinceFdx=0,0T, it means that there is no dissipation term in system (2) and trajectories are corresponding to the unique conservative field Fcx. That is, when the trajectories begin from an isosurface of the constant Hamilton energy Hx, they cannot escape from the isosurface [57, 58].According to (22), the Hamilton energy function Hx,y will satisfy the partial differential equation given by(25)y∂H∂x−sinax+ksinbx∂H∂y=0.So, a general solution of (25) can be approached by(26)Hx,y,a,b,k=12y2−1acosax−kbcosbx.The time derivative of (26) can be calculated by(27)H˙x,y,a,b,k=y.y˙+sinax+ksinbx⋅x˙=ysinax+ksinbx−ysinax+ksinbx=0=∇HTfdx.Since the time derivative of Hamiltonian energy of system (2) is zero, this Hamilton energy function is invariable. This justifies the conservative nature of the new system (2). Thus, the conservative property of system (2) is validated using theoretical measures.


## 4. The Forced Chaotic Oscillator

Like in the van der Pol equation, a temporally periodic forcing termAsinωt can be added to nonlinear system (2), and thus we introduced the following system to describe the forced version of (2):(28)x˙=−y,y˙=sinax+ksinbx+Asinwt.In the following sections, the dynamical analysis of system (28) is discussed using various numerical tools including two-parameter Lyapunov exponents, Lyapunov spectrum plots, and phase trajectories. Here, the main objective is finding possible domains of chaotic attractors in (28). In this regard, several combinations of A,w,k can be used to obtain chaos. For more simplicity, we choose w=0.7, and we consider the same set of system parameter a,b as in the previous section. Therefore, the dynamics of the new system (28) is investigated in terms of A,k.


### 4.1. Two-Parameter Lyapunov Exponent Analysis and Multistability

The dynamical regions of the simple Hamiltonian megastable system are drawn in this part by means of two-parameter Lyapunov exponent analysis when varying the control parameterA and k in the appropriate domain. However, the global dynamic behaviors of the system can be well quantified exploiting two-parameter Lyapunov exponents by adjusting simultaneously these two last parameters of the 2D system (28) through the establishment of appropriate colorful diagrams. The obtained colorful diagrams are produced by the Lyapunov exponent spectrum on a grid of 500×500 values of the two control space parameters (i.e., A and k). The main advantage of the two-parameter Lyapunov exponent is that the complexity of the system can be directly analyzed on the variation of two parameters, which offers a major advantage over the traditional exponent which is defined on a single parameter. To properly obtain this result, system (28) is numerically solved using the standard fourth-order Runge–Kutta integration method with a fixed step size equal to 4×10−3. We also used 26×104 steps to compute the Lyapunov exponent (LE). Figure 3 shows the two-parameter Lyapunov exponents in the parameter space A,k using two different initial conditions. Periodic and quasiperiodic regions are indicated by blue and cyan colors, respectively, while chaotic oscillations are materialized by yellowish-reddish color. The graphs are obtained by increasing both control parameters. From these graphs, one can observe the regions of coexisting behaviors (i.e., hysteretic dynamics) characterized by the presence of different colors in the map. To better highlight the presence of multiple coexisting attractors, we have plotted the trajectories showing different coexisting strange attractors of system (28) for A=0.58 and k=2 as shown in Figure 4.Figure 3
Coexisting of two-parameter sweep diagrams plotted in the same range of parameters with just two different initial conditions30,0.2 and 1,0.2. The coexisting solution regions are synthesized by the simultaneous appearance of different colors on figures for the same range of parameters. The rest of the system parameters are fixed as in the text.
(a)(b)Figure 4
Trajectories showing different coexisting strange attractors of system (28) for A=0.58. Each trajectory is for a duration of 1200 seconds.


### 4.2. Infinite Number of Coexisting Nested Strange Attractors

For more simplicity, we setk=2 and consider A as a control parameter of system (28). Figure 5 shows the Lyapunov spectrum plots when parameter A is varied from 0 to 1.2 with initial conditions 1,0. One can note from this diagram that the forced version of new system (2) presents a very large range of chaotic dynamic provided that in the region 0.12<A≤1.2, λ1>0, λ2=0, and λ3<0. However, the conservative nature of the new two-dimensional megastable system is also appreciated because ∑i=13λi=0 [55, 56, 65]. That is, the maximum value of LE λ1 and the minimum value λ3 are symmetric with respect to λ2. When ICs are 1,0 and A=0.58, the corresponding LEs are λ1=0.0943, λ2=0.000, and λ3=−0.0943.Figure 5
Lyapunov exponent spectrum for system (28) versus forcing amplitude A when w=0.7. This graph is plotted with initial condition (1, 0).However,∑i=12λi=λ1>0 and ∑i=13λi=λ1+λ2=0. The Kaplan–Yorke dimension (Lyapunov dimension) is defined as(29)DKY=D+∑i=1DλDλD+1,where D is the largest integer satisfying ∑i=1DλD≥0 and ∑i=1D+1λD<0. Thus, DKY=3 provided that ∑i=13λi=∑i=1Dλi=0.Furthermore, using four different initial conditions, the nested coexisting bifurcation diagrams of system (28) are plotted and superimposed as shown in Figure 6. The graphs represent the plot of local maximal values of y in terms of the control parameter A that is varied (i.e., increased) in the range 0<A≤1.2. The presence of four different coexisting data in Figure 6 further confirms the complicated dynamical behaviors observed in this work. As a result, the presence of infinite number of nested attractors for A=1.2 is depicted in Figure 7(a) and the zoom (i.e., attractors in the red box of Figure 7(a)) is shown in Figure 7(b). Note that these plots are just examples in an infinite number of stable states around the x-axis, which could be found under other proper initial conditions.Figure 6
Nested bifurcation diagrams of local maxima of “y” versus forcing amplitude A in system (28) with w=0.7. Four sets of data are superimposed with initial conditions blue (1, 0.2), black (25, 0.2), red (30, 0.2), and green (43, 0.2) without reinitiating.Figure 7
(a) Trajectories in system (28) for A=1.2, with 75 initial conditions located on the x-axis (from x=−150 to y=150 with steps equal to 4). Each trajectory is plotted for 1200 seconds. (b) The enlarged version of the selected window in (a) (for interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
(a)(b)


## 4.1. Two-Parameter Lyapunov Exponent Analysis and Multistability

The dynamical regions of the simple Hamiltonian megastable system are drawn in this part by means of two-parameter Lyapunov exponent analysis when varying the control parameterA and k in the appropriate domain. However, the global dynamic behaviors of the system can be well quantified exploiting two-parameter Lyapunov exponents by adjusting simultaneously these two last parameters of the 2D system (28) through the establishment of appropriate colorful diagrams. The obtained colorful diagrams are produced by the Lyapunov exponent spectrum on a grid of 500×500 values of the two control space parameters (i.e., A and k). The main advantage of the two-parameter Lyapunov exponent is that the complexity of the system can be directly analyzed on the variation of two parameters, which offers a major advantage over the traditional exponent which is defined on a single parameter. To properly obtain this result, system (28) is numerically solved using the standard fourth-order Runge–Kutta integration method with a fixed step size equal to 4×10−3. We also used 26×104 steps to compute the Lyapunov exponent (LE). Figure 3 shows the two-parameter Lyapunov exponents in the parameter space A,k using two different initial conditions. Periodic and quasiperiodic regions are indicated by blue and cyan colors, respectively, while chaotic oscillations are materialized by yellowish-reddish color. The graphs are obtained by increasing both control parameters. From these graphs, one can observe the regions of coexisting behaviors (i.e., hysteretic dynamics) characterized by the presence of different colors in the map. To better highlight the presence of multiple coexisting attractors, we have plotted the trajectories showing different coexisting strange attractors of system (28) for A=0.58 and k=2 as shown in Figure 4.Figure 3
Coexisting of two-parameter sweep diagrams plotted in the same range of parameters with just two different initial conditions30,0.2 and 1,0.2. The coexisting solution regions are synthesized by the simultaneous appearance of different colors on figures for the same range of parameters. The rest of the system parameters are fixed as in the text.
(a)(b)Figure 4
Trajectories showing different coexisting strange attractors of system (28) for A=0.58. Each trajectory is for a duration of 1200 seconds.


## 4.2. Infinite Number of Coexisting Nested Strange Attractors

For more simplicity, we setk=2 and consider A as a control parameter of system (28). Figure 5 shows the Lyapunov spectrum plots when parameter A is varied from 0 to 1.2 with initial conditions 1,0. One can note from this diagram that the forced version of new system (2) presents a very large range of chaotic dynamic provided that in the region 0.12<A≤1.2, λ1>0, λ2=0, and λ3<0. However, the conservative nature of the new two-dimensional megastable system is also appreciated because ∑i=13λi=0 [55, 56, 65]. That is, the maximum value of LE λ1 and the minimum value λ3 are symmetric with respect to λ2. When ICs are 1,0 and A=0.58, the corresponding LEs are λ1=0.0943, λ2=0.000, and λ3=−0.0943.Figure 5
Lyapunov exponent spectrum for system (28) versus forcing amplitude A when w=0.7. This graph is plotted with initial condition (1, 0).However,∑i=12λi=λ1>0 and ∑i=13λi=λ1+λ2=0. The Kaplan–Yorke dimension (Lyapunov dimension) is defined as(29)DKY=D+∑i=1DλDλD+1,where D is the largest integer satisfying ∑i=1DλD≥0 and ∑i=1D+1λD<0. Thus, DKY=3 provided that ∑i=13λi=∑i=1Dλi=0.Furthermore, using four different initial conditions, the nested coexisting bifurcation diagrams of system (28) are plotted and superimposed as shown in Figure 6. The graphs represent the plot of local maximal values of y in terms of the control parameter A that is varied (i.e., increased) in the range 0<A≤1.2. The presence of four different coexisting data in Figure 6 further confirms the complicated dynamical behaviors observed in this work. As a result, the presence of infinite number of nested attractors for A=1.2 is depicted in Figure 7(a) and the zoom (i.e., attractors in the red box of Figure 7(a)) is shown in Figure 7(b). Note that these plots are just examples in an infinite number of stable states around the x-axis, which could be found under other proper initial conditions.Figure 6
Nested bifurcation diagrams of local maxima of “y” versus forcing amplitude A in system (28) with w=0.7. Four sets of data are superimposed with initial conditions blue (1, 0.2), black (25, 0.2), red (30, 0.2), and green (43, 0.2) without reinitiating.Figure 7
(a) Trajectories in system (28) for A=1.2, with 75 initial conditions located on the x-axis (from x=−150 to y=150 with steps equal to 4). Each trajectory is plotted for 1200 seconds. (b) The enlarged version of the selected window in (a) (for interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
(a)(b)


## 5. Circuit Implementation

In this part of our work, an analog circuit is constructed and used to make a comparison between the theoretical/numerical results obtained previously and the experimental results. The circuit implementation defines the analog circuit of our model, which confirms a possibility of laboratory measurement. Unlike other simulation software, PSpice remains the best software by its capacity to integrate the transient time. The circuit diagram that allows us to perform the various simulations in the PSpice software is presented in Figure8. This electronic implementation of the model is very singular because the model is built based only on linear and trigonometric terms. The circuit of the new system is designed using two capacitors C1,C2 and ten resistors R1,…,R10 including five op-amps TL082CD. In addition, it uses two sine blocs, one AC source, and a symmetric power supply. The circuit equation using Kirchhoff’s electrical circuit laws can be obtained as(30)C1dXdt=−1R1Y,C2dYdt=1R6sinR3R2X+1R7sinR5R4X+VemaxR8sinWt.Figure 8
Analog computer of the new Hamiltonian megastable oscillator.SetC1=C2=C=10nF, Vemax=1V, and R=Ri=10KΩ except R2, R5, R10, R12, and R13; adopt the rescale of time t=τRC and variables X=10V×x and Y=1V×y.System (30) is the same as the introduced system (28) with the following expression of parameters:R1=10R=100KΩ, R2=R/10a=1.414KΩ, R4=R/10b=20KΩ, R7=R/K=5KΩ, W=w0w=7000, and R8=VemaxR/A=50KΩ with w0=1/RC.The dynamics of the circuit provided in Figure8 is simulated in PSpice as provided in Figure 9. This obtained result agrees with the one obtained based on Pascal simulation. From this figure, it can be seen that some of the coexisting attractors of the new megastable oscillator have been captured and support the fact that the obtained results were not artifacts.Figure 9
Some coexisting attractors obtained from the new megastable oscillator forR8=50KΩ (A = 0.2) with initial conditions 5V,1V, 10V,1V, and 8V,1V, respectively.
(a)(b)(c)


## 6. Concluding Remark

This paper was focused on the dynamical analysis of the simple 2D autonomous system with infinite of hyperbolic and nonhyperbolic equilibria reported to date. The investigations of the Hamiltonian and dissipation properties of the model reveal that it is conservative and displays the coexistence of a large number of stable states. Equally, it is demonstrated under the presence of a sinusoidal excitation that the proposed model exhibits the striking phenomenon of megastability characterized by the coexistence of an infinite number of countable stable states. Finally, the electronic implementation of the introduced model has been provided using PSPice simulations to further support our results. In the future, it would be interesting to carry out the analytical methods which can predict megastability behaviors in a nonlinear system and address the megastable system-based secure communication.



---

*Source: 9260823-2020-09-28.xml*