# Effects of Symmetric and Asymmetric Nonlinearity on the Dynamics of a Third-Order Autonomous Duffing–Holmes Oscillator

Authors: Isaac Sami Doubla; Jacques Kengne; Raoul Blaise Wafo Tekam; Zeric Tabekoueng Njitacke; Clotaire Thierry Sanjong Dagang
Journal: Complexity
 (2020)

---

## Abstract

A generalized third-order autonomous Duffing–Holmes system is proposed and deeply investigated. The proposed system is obtained by adding a parametric quadratic termmx2 to the cubic nonlinear term −x3 of an existing third-order autonomous Duffing–Holmes system. This modification allows the system to feature smoothly adjustable nonlinearity, symmetry, and nontrivial equilibria. A particular attention is given to the effects of symmetric and asymmetric nonlinearity on the dynamics of the system. For the specific case of m=0, the system is symmetric and interesting phenomena are observed, namely, coexistence of symmetric bifurcations, presence of parallel branches, and the coexistence of four (periodic-chaotic) and six (periodic) symmetric attractors. For m≠0, the system loses its symmetry. This favors the emergence of other behaviors, such as the coexistence of asymmetric bifurcations, involving the coexistence of several asymmetric attractors (periodic-periodic or periodic-chaotic). All these phenomena have been numerically highlighted using nonlinear dynamic tools (bifurcation diagrams, Lyapunov exponents, phase portraits, time series, frequency spectra, Poincaré section, cross sections of the attraction basins, etc.) and an analog computer of the system. In fact, PSpice simulations of the latter confirm numerical results. Moreover, amplitude control and synchronization strategies are also provided in order to promote the exploitation of the proposed system in engineering.


---

## Body

## 1. Introduction

The study of nonlinear dynamic systems has been attracting researchers from various fields of science. Thanks to analytical, numerical, and experimental tools, striking nonlinear phenomena involved for instance in the evolution of stock market [1], neural networks [2–6], and population dynamics [7] have been uncovered. However, systems exhibiting chaotic dynamics are of particular interest owing to their potential applications in various fields of science and engineering such as image encryption [8], random bits generation [9], and secure communications [10]. Yet, transient chaos [11, 12], antimonotonicity [13–15], multistability [13–16], and many other dynamical features of chaotic systems are relevant for their applications [15, 17–21].Multistability, namely, the occurrence of two or more attractors for a given set of system parameter values due to initial states (or/and noise sensitivity) [22–24], has explained a large number of events in biology [25, 26] and chemistry [27], just to mention a few. Some important works based on the construction and study of multistable chaotic systems or chaotic systems with the coexistence of multiple attractors have recently been published in [28–33]. The work of Lai et al. [28] was based on the study of an extended Lü system; their results showed that the system has a pair of strange attractors, a pair of limit cycles, and a pair of point attractors for different initial conditions. The authors [29] have constructed an extremely simple chaotic system with an infinite number of coexisting chaotic attractors. Lai et al. [30] reported a new chaotic system they generated from the simplest chaotic memristor circuit by introducing a simple nonlinear feedback control input. The main characteristic of the new system is that it has an infinite number of coexisting equilibria and abundant attractors. The results of [31] report a new chaotic system without equilibrium with hidden attractors and coexisting attractors. This special feature (multistability) can sometimes be justified by the fact that the system is symmetric [34]. In this line, Kengne and collaborators have highlighted the coexistence of four, six, and eight attractors from a cubic, hyperbolic tangent, and hyperbolic tangent-cubic nonlinearity-based chaotic jerk system, respectively [35–37]. The authors of [38] showed the presence of nine coexisting attractors in a hyperjerk chaotic system with hyperbolic sine nonlinearity. On the other hand, Negou and Kengne [34] recently introduced a jerk system with adjustable symmetry and nonlinearity. The system was obtained by replacing a hyperbolic sine nonlinearity in an existing jerk system with a smoothly adjustable nonlinearity. This modification demonstrated the relevance of adjustable symmetries in the emergence of striking multistability. The authors of [39–41] did similar works by considering adjustable symmetry and nonlinearities in the dynamics of a simple jerk system, snap system, and jerk circuit, respectively. As far as symmetry breaking is concerned, the results obtained highlighted issues, such as the presence of parallel bifurcation branches, hysteresis, and coexisting multiple asymmetric attractors in the mentioned systems and circuits. These results are of particular practical interest given that perfect symmetry is not a practical reality.Motivated by these results, we propose in this article a third-order chaotic system with adjustable symmetry, nonlinearity, and nontrivial equilibria as well. The novel system is a generalized form of the third-order autonomous Duffing–Holmes oscillator introduced by Lindberg et al. [42]. The novelty is brought by adding a parametric quadratic term to the cubic nonlinear term of the former system, resulting in a nonlinear function of the form φmx=x−mx2−x3. The quadratic term has been considered here because it models a possible imperfection that breaks the perfect symmetry of the original system. This modification can also be regarded as mathematical techniques to discover new nonlinear patterns. Historically, forced Duffing–Holmes system is one of the well-known simple nonautonomous second-order chaotic system [43, 44]. In 2009, Lindberg et al. [42] suggested an autonomous third-order Duffing–Holmes type system, with a smoother spectrum. In spite of the numerous studies recently devoted to uncovering multistability in nonlinear systems, no report has shown that the autonomous third-order Duffing–Holmes system suggested by Lindberg et al. has such striking behavior (to the best of the authors’ knowledge). With the generalized third-order autonomous Duffing–Holmes system proposed in this work, our aim is twofold. Firstly, we show that the autonomous Duffing–Holmes system introduced by Lindberg et al. (which is a particular case of the generalized one) exhibits multistability. Secondly, we investigate the effects of adjustable symmetry, nonlinearity, and nontrivial equilibria on the dynamics previously obtained (with the particular case). Moreover, amplitude control and adaptive synchronization of the system are performed for promoting possible applications in engineering.The remainder of this article is structured as follows. Section2 is devoted to the description of the novel system and analytical examination of its basic properties. In Section 3, the system is numerically analyzed using appropriate tools of nonlinear dynamics. In Section 4, an electronic analog calculator of the system is designed and PSpice simulation results are also presented. Sections 5 and 6 are about amplitude control and adaptive synchronization, respectively. Finally, conclusions are drawn in Section 7.


## 2. System Description and Basic Properties

### 2.1. System Description

The autonomous Duffing–Holmes system derived by Lindberg et al. [42] is as follows:(1)x˙=y,y˙=x−x3+by−kz,z˙=wfy−z,where b represents the damping coefficient, k is the feedback coefficient, and wf is the cutoff frequency of the filter (see Ref. [42] for more details). The system under investigation is obtained by adding the parametric quadratic term mx2 to the second equation of the autonomous Duffing–Holmes system (1), that is:(2)x˙=y,y˙=x−x3+by−kz−mx2,z˙=wfy−z.The quadratic nonlinearity coefficientm is introduced to modify system’s nonlinearity, symmetry, and equilibria as it will be shown later. Note that system (1) can always be obtained by setting m=0. Though this has been studied in [42], some complex behaviors such as coexistence of attractors were not revealed. With the generalized model (system (2)), we are going to revisit system (1) and investigate the case with m≠0.


### 2.2. Dissipation and Existence of Attractors

The presence of attractors in a system can be determined by calculating its divergence. The divergence of system (2) is given by(3)Λ=∂x˙∂x+∂y˙∂y+∂z˙∂z=−wf−b.Given thatwf and b are strictly positive, the divergence of system (2) is strictly negative whenever wf>b (case considered in this work). This proves the existence of attractors because the volume elements contract after a unit of time. This contraction reduces a volume V0 by a factor e−wf−bt. This means that each volume containing the trajectory of system (2) converges to zero when t tends to infinity at an exponential rate of wf−b. Therefore, all the orbits of system (2) are finally limited to a specific subset having a zero volume. Thus, the asymptotic movement settles on an attractor [45–49].


### 2.3. Symmetry

Symmetry is an interesting feature of dynamic systems. It tells us if the solutions are unique or even. System (2) is invariant to the transformation x,y,z⇔−x,−y,−z when m=0. In this case, if sx,y,z is a solution of the system for a set of parameter values, then s−x,−y,−z is also a solution for the same values of parameters. This symmetry will be responsible for the appearance of the solutions in pair by polarity inversion of the initial conditions [50] and thus leads to the coexistence of attractors in the state space. However, the aforementioned symmetry disappears for m≠0.


### 2.4. Analysis of Fixed Points and Hopf Bifurcation

When studying the dynamics of a system, it is always interesting to analyze the stability of its fixed points. By setting the left-hand side of equation (2) equal to zero x˙=y˙=z˙=0, we find that system (2) has three fixed points: one trivial equilibrium p1x01,y01,z01=0,0,0 and two nontrivial equilibria p2,3x02,03,y02,03,z02,03=−m±m2+4/2,0,0. It should be noted that coordinates of nontrivial equilibria depend on the parameter m; it is therefore possible to adjust their positions by acting on the parameter m. The Jacobian matrix of system (2) associated with any of these fixed points is given by the following expression:(4)MJ=0101−3x2−2mxb−k0wf−wf.The nature of equilibria can be found from the following characteristic equationdetMJ−λI3=0:(5)λ3+wf−bλ2+3x2+2mx−bwf+kwf−1λ+wfx3x+2m−wf=0,where I3 denotes the 3 × 3 identity matrix. By applying Routh–Hurwitz stability criteria to characteristic equation (5), we obtain the following generalized stability conditions for the equilibrium points of system (2):(6)wf−b>0,wf−b3x2+2mx−bwf+kwf−1−wfx3x+2m+wf>0,wfx3x−2m−1>0.For the trivial fixed pointp1x01,y01,z01=0,0,0, we obtain the following stability conditions:(7)wf>b,k>1−1wfwf−bb,wf<0,while for the equilibrium points P2,3, we have the following stability conditions:(8)wf>b,k>b−wf2−3x02,032−2mx02,03+bwf+1wf−wf+b,wf>0,m<±21515,with(9)x02,03=−m±m2+42.Table1 confirms relations (7) and (8), for different values of the parameters wf,betm.Table 1
Stability of the equilibrium pointsP1,P2, and P3 obtained from equations (7) and (8) for different values of wf,b,k, and m.
Equilibrium pointsParameter valuesEigenvaluesStabilitywfbkmλ1λ2λ3P1−1−3−8—−0.2848 + 0.0000i−0.8576 + 1.6661i−0.8576 − 1.6661iStable−2−4−7—−0.4668 + 0.0000i−0.7666 + 1.9227i−0.7666 − 1.9227iStable1310—0.1761 + 0.0000i0.9120 + 2.2016i0.9120 − 2.2016iUnstable249—0.2329 + 0.0000i0.8836 + 2.7942i0.8836 − 2.7942iUnstableP2131.66960.6−0.7163 + 0.0000i1.3582 + 0.5695i1.3582 − 0.5695iUnstable312.25890.6−1.5882 + 0.0000i−0.2059 + 3.3692i−0.2059 − 3.3692iStableP3242.1936−0.6−1.4073 + 0.0000i1.7036 + 1.0422i1.7036 − 1.0422iUnstable423.7016−0.6−1.2907 + 0.0000i−0.3546 + 2.9277i−0.3546 − 2.9277iStableConsidering the change of variableλ=iw0w0>0 and substituting into (5), after separation of real part from imaginary part, we obtain relation (10) for the trivial point of equilibrium P1 when m=0.(10)wHopf=w0=wfb−wf,kc=1−1wfwf−bb.Differentiation of equation (5) with respect to k gives the following relation:(11)Real∂λ∂k|k=kc=−12wf1b−wf2−12+1≠0,∀wf>0.Equation (10) defines the frequency of stable oscillations as well as the critical value of k corresponding to the Hopf bifurcation of system (2). From relation (11), note that the transversality condition is always satisfied provided that wf is a strictly positive parameter.


## 2.1. System Description

The autonomous Duffing–Holmes system derived by Lindberg et al. [42] is as follows:(1)x˙=y,y˙=x−x3+by−kz,z˙=wfy−z,where b represents the damping coefficient, k is the feedback coefficient, and wf is the cutoff frequency of the filter (see Ref. [42] for more details). The system under investigation is obtained by adding the parametric quadratic term mx2 to the second equation of the autonomous Duffing–Holmes system (1), that is:(2)x˙=y,y˙=x−x3+by−kz−mx2,z˙=wfy−z.The quadratic nonlinearity coefficientm is introduced to modify system’s nonlinearity, symmetry, and equilibria as it will be shown later. Note that system (1) can always be obtained by setting m=0. Though this has been studied in [42], some complex behaviors such as coexistence of attractors were not revealed. With the generalized model (system (2)), we are going to revisit system (1) and investigate the case with m≠0.


## 2.2. Dissipation and Existence of Attractors

The presence of attractors in a system can be determined by calculating its divergence. The divergence of system (2) is given by(3)Λ=∂x˙∂x+∂y˙∂y+∂z˙∂z=−wf−b.Given thatwf and b are strictly positive, the divergence of system (2) is strictly negative whenever wf>b (case considered in this work). This proves the existence of attractors because the volume elements contract after a unit of time. This contraction reduces a volume V0 by a factor e−wf−bt. This means that each volume containing the trajectory of system (2) converges to zero when t tends to infinity at an exponential rate of wf−b. Therefore, all the orbits of system (2) are finally limited to a specific subset having a zero volume. Thus, the asymptotic movement settles on an attractor [45–49].


## 2.3. Symmetry

Symmetry is an interesting feature of dynamic systems. It tells us if the solutions are unique or even. System (2) is invariant to the transformation x,y,z⇔−x,−y,−z when m=0. In this case, if sx,y,z is a solution of the system for a set of parameter values, then s−x,−y,−z is also a solution for the same values of parameters. This symmetry will be responsible for the appearance of the solutions in pair by polarity inversion of the initial conditions [50] and thus leads to the coexistence of attractors in the state space. However, the aforementioned symmetry disappears for m≠0.


## 2.4. Analysis of Fixed Points and Hopf Bifurcation

When studying the dynamics of a system, it is always interesting to analyze the stability of its fixed points. By setting the left-hand side of equation (2) equal to zero x˙=y˙=z˙=0, we find that system (2) has three fixed points: one trivial equilibrium p1x01,y01,z01=0,0,0 and two nontrivial equilibria p2,3x02,03,y02,03,z02,03=−m±m2+4/2,0,0. It should be noted that coordinates of nontrivial equilibria depend on the parameter m; it is therefore possible to adjust their positions by acting on the parameter m. The Jacobian matrix of system (2) associated with any of these fixed points is given by the following expression:(4)MJ=0101−3x2−2mxb−k0wf−wf.The nature of equilibria can be found from the following characteristic equationdetMJ−λI3=0:(5)λ3+wf−bλ2+3x2+2mx−bwf+kwf−1λ+wfx3x+2m−wf=0,where I3 denotes the 3 × 3 identity matrix. By applying Routh–Hurwitz stability criteria to characteristic equation (5), we obtain the following generalized stability conditions for the equilibrium points of system (2):(6)wf−b>0,wf−b3x2+2mx−bwf+kwf−1−wfx3x+2m+wf>0,wfx3x−2m−1>0.For the trivial fixed pointp1x01,y01,z01=0,0,0, we obtain the following stability conditions:(7)wf>b,k>1−1wfwf−bb,wf<0,while for the equilibrium points P2,3, we have the following stability conditions:(8)wf>b,k>b−wf2−3x02,032−2mx02,03+bwf+1wf−wf+b,wf>0,m<±21515,with(9)x02,03=−m±m2+42.Table1 confirms relations (7) and (8), for different values of the parameters wf,betm.Table 1
Stability of the equilibrium pointsP1,P2, and P3 obtained from equations (7) and (8) for different values of wf,b,k, and m.
Equilibrium pointsParameter valuesEigenvaluesStabilitywfbkmλ1λ2λ3P1−1−3−8—−0.2848 + 0.0000i−0.8576 + 1.6661i−0.8576 − 1.6661iStable−2−4−7—−0.4668 + 0.0000i−0.7666 + 1.9227i−0.7666 − 1.9227iStable1310—0.1761 + 0.0000i0.9120 + 2.2016i0.9120 − 2.2016iUnstable249—0.2329 + 0.0000i0.8836 + 2.7942i0.8836 − 2.7942iUnstableP2131.66960.6−0.7163 + 0.0000i1.3582 + 0.5695i1.3582 − 0.5695iUnstable312.25890.6−1.5882 + 0.0000i−0.2059 + 3.3692i−0.2059 − 3.3692iStableP3242.1936−0.6−1.4073 + 0.0000i1.7036 + 1.0422i1.7036 − 1.0422iUnstable423.7016−0.6−1.2907 + 0.0000i−0.3546 + 2.9277i−0.3546 − 2.9277iStableConsidering the change of variableλ=iw0w0>0 and substituting into (5), after separation of real part from imaginary part, we obtain relation (10) for the trivial point of equilibrium P1 when m=0.(10)wHopf=w0=wfb−wf,kc=1−1wfwf−bb.Differentiation of equation (5) with respect to k gives the following relation:(11)Real∂λ∂k|k=kc=−12wf1b−wf2−12+1≠0,∀wf>0.Equation (10) defines the frequency of stable oscillations as well as the critical value of k corresponding to the Hopf bifurcation of system (2). From relation (11), note that the transversality condition is always satisfied provided that wf is a strictly positive parameter.


## 3. Numerical Investigations

### 3.1. Methodology

In this section, we do a thorough analysis of the dynamics of system (2) using numerical tools such as bifurcation diagrams, maximum Lyapunov exponent, Poincaré section, phase portraits, time series, cross section of the basin of attraction, and so on. Bifurcation diagrams are obtained from the local maxima of one of the state variables for different values of the control parameters. Maximum Lyapunov exponents are computed using Wolf‘s algorithm [51]. In fact, the Maximum Lyapunov exponent makes it possible to estimate the rate of convergence or divergence between two close trajectories. Thus, if the maximum Lyapunov exponent is (i) positive, the system is considered as chaotic because for a small perturbation, its trajectory diverges (this implies the occurrence of a strange attractor in phase space); (ii) negative, the system is considered stable because for all perturbations, it returns to its equilibrium point; (iii) null, the system is considered periodic or quasiperiodic.All the numerical integrations are based on the fourth-order Runge–Kutta algorithm with an integration steph=5.10−3, for better precision.


### 3.2. Coexistence of Bifurcations and Multiple Attractors

#### 3.2.1. Consideringb as a Control Parameter with m=0

In this section, we considerb as a control parameter and m is equal to zero. The dynamics obtained is that represented by the bifurcation diagrams of Figure 1(a), carried out in the ascending direction for 0.22≤b≤0.32,wf=0.5,andk=1.6. This figure presents a very interesting phenomenon which is the coexistence of two symmetric bifurcation diagrams. This coexistence is highlighted through the superposition of two diagrams (red and blue) obtained by symmetrically changing the initial conditions 0.25,0,0⇔−0.25,0,0. They are in good agreement with the corresponding graphs of Lyapunov exponents as shown in Figure 2(a). Figure 3 shows some 3D phase portraits obtained for different values of the control parameter with initial conditions as indicated. Figures 4(a)–4(d) show a symmetric chaotic attractor, its time series, its frequency spectrum, and its Poincaré section (these representations are consistent with the chaotic behavior), respectively, for b=0.32. All these confirm the symmetric nature of the dynamic of system (2).Figure 1
Bifurcation diagrams as a function of the parameter b showing destruction of the symmetry of the system for some values ofm: (a) the symmetric system for m=0.00; (b–d) destruction of the system’s symmetry for m=0.05,m=0.10, and m=0.25, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the traces in red and blue, respectively.
(a)(b)(c)(d)Figure 2
Diagrams of the maximum Lyapunov exponents, corresponding to the bifurcation diagrams of Figure1. This figure shows the impact of the value of m on the symmetry of system (2) considering b as a control parameter.
(a)(b)(c)(d)Figure 3
Three-dimensional representationsx−y−z of trajectories in phase space, at different points of the bifurcation diagram of the symmetry system (Figure 1(a)) for different values of b: (a) cycle of period-1 for b=0.22, (b) cycle of period-2 for b=0.28, (c) cycle of period-4 for b=0.284, and (d, e) single and double chaotic bands for b=0.288 and b=0.3, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the curves in red and blue, respectively.
(a)(b)(c)(d)(e)Figure 4
Three-dimensional portrait (a), time series (b), frequency spectrum (c), and Poincaré section (d), confirming symmetric chaotic behavior of the system form=0.00 and b=0.32. The other parameters are specified in the text.
(a)(b)(c)(d)


#### 3.2.2. Consideringb as a Control Parameter with m≠0

In this section, we still considerb as a control parameter but m≠0.0. With 0.22≤b≤0.32,wf=0.5,andk=1.6, bifurcation diagrams shown in Figures 1(b)–1(d) are plotted for m equal to 0.05,0.10,and0.25, respectively. Let us recall that red color is used for x0,y0,z0=0.25,0,0 and blue color is used for x0,y0,z0=−0.25,0,0. Comparing Figures 1(b)–1(d) with Figure 1(a), we can observe that bifurcation diagrams in red are no more symmetric to their counterpart (in blue). There is a gap (discontinuity) in the blue bifurcation diagrams (Figures 1(b)–1(d)), which increases as m is increased. This observation is caused by a shift in equilibrium points (see Table 2) and implies the existence of coexisting attractors. The coexistence of a period-2 limit cycle with a period-3 limit cycle (Figure 5(c)) and a chaotic attractor with a periodic attractor (Figure 5(d)), for instance, is revealed in Figure 5. Figures 6(a)–6(d) represent an asymmetric chaotic attractor, its time series, its frequency spectrum, and its Poincaré section (these representations are specific to chaotic attractors), respectively, for b=0.32. The aforementioned results indicated an increase in system’s complexity when m≠0.Table 2
Modifications of nontrivial equilibrium points for different values ofm.
Values ofmNontrivial equilibrium pointsP2P30.00(1, 0, 0)(−1, 0, 0)0.05(0.9753, 0, 0)(−1.0253, 0, 0)0.10(0.9512, 0, 0)(−1.0512, 0, 0)0.125(0.9395, 0, 0)(−1.0645, 0, 0)0.25(0.8828, 0, 0)(−1.1328, 0, 0)Figure 5
Three-dimensional representationsx−y−z of trajectories in the phase space, at different points of the bifurcation diagram of the asymmetry system (Figure 2(a)), for different values of b with m = 0.05: (a, b) cycle of period-1 for b=0.22 and cycle of period-2 for b=0.25; (c) cycles of period-4 and period-2 for b=0.2824; (d) single-band chaotic attractor with cycle of period-4 for b=0.2856; (e, f) single and chaotic dual-band attractor for b=0.2892 and b=0.3, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the traces in red and blue, respectively.
(a)(b)(c)(d)(e)(f)Figure 6
Three-dimensional portrait (a), time series (b), frequency spectrum (c), and Poincaré section (d), confirming asymmetric chaotic behavior of the system form=0.25 and b=0.32. The other parameters are specified in the text.
(a)(b)(c)(d)


#### 3.2.3. Consideringk as a Control Parameter with m=0

Consideringk as a control parameter with m equal to zero, the dynamics of system (2) is governed by the superimposed bifurcation diagrams and Lyapunov exponents of Figures 7(a) and 8(a), respectively. These diagrams are obtained by varying the control parameter (downward) in the range 2≤k≤3 for initial conditions of opposite signs ±0.75,0,0 with wf=0.5andb=0.41. Figure 9 shows an enlargement of the diagram of Figure 7(a) (in red) in the intervals 2.6575≤k≤2.68 (Figure 9(a)) and 2.671≤k≤2.672 (Figure 9(b)) illustrating parallel branches. They show the coexistence of four (two periodic and two chaotic) and six attractors occurring in symmetric pairs. Figures 10 and 11 illustrate the 2D representations of these different attractors in the phase space. These attractors are obtained for the values of k=2.667and2.6714, respectively. The different values of the initial conditions x0,y0,z0 are as follows: for the coexistence of the four attractors—period-1 ±0.12,0,0 and chaotic spiral ±0.39,0,0; for the coexistence of the six attractors—n-finite periods ±0.28,0,0, period-6±0.06,0,0, and period-1 ±0.108,0,0. These multistabilities are well illustrated on the cross sections of the basins represented in Figures 12(a)–12(d), where we can observe several regions distinguished by different colors that correspond to each attractor.Figure 7
Bifurcation diagrams as function of the parameterk, showing destruction of the symmetry of the system for some values of m: (a) the symmetric system for m = 0; (b, c) destruction of the system symmetry for m = 0.10 and m=0.125, respectively. The initial conditions are x0,y0,z0=0.75,0,0 and −0.75,0,0 for the traces in blue and red, respectively.
(a)(b)(c)Figure 8
Diagrams of the maximum Lyapunov exponents corresponding to the bifurcation diagrams of Figure7 showing the impact of the value of m on the symmetry of system (2) considering k as a control parameter.
(a)(b)(c)Figure 9
Enlargement of the bifurcation diagram of Figure7(a), in the intervals 2.6575 ≤ k ≤ 2.68 (a) and 2.671 ≤ k ≤ 2.672 (b), revealing the coexistence of attractors in the form of parallel branches.
(a)(b)Figure 10
Coexistence of four attractors at the pointk = 2.66 of Figure 9(a), for different values of the initial conditions (x(0), y(0), z(0)): (a) a period-1 attractor and its symmetry for (±0.12, 0, 0); (b) a chaotic spiral attractor and its symmetry for (±0.39, 0, 0).
(a)(b)Figure 11
Coexistence of six periodic attractors at pointk = 2.6714 of Figure 9(b), for different values of the initial conditions (x(0), y(0), z(0)): (a) a finite n-period attractor and its symmetry for (±0.028, 0, 0); (b) a period-6 attractor and its symmetry for (±0.06, 0, 0); (c) a period-1 attractor and its symmetry for (±0.108, 0, 0).
(a)(b)(c)Figure 12
Cross sections of the basins of attraction in thex0−y0 and y0−z0 plane: on (a, b) blue and yellow colors indicate regions of period-1 while green and magenta indicate regions of chaotic attractors; on (c, d) magenta and yellow colors indicate the regions of n-finite periods, cyan and black indicate the regions of period-6 while green and blue the regions indicate of period-1. Regions of unbounded behavior are in red.
(a)(b)(c)(d)


#### 3.2.4. Consideringk as a Control Parameter with m≠0

Keepingk2≤k≤3 as the control parameter, with m≠0m=0.10,0.125, wf=0.5,andb=0.41, we obtain the asymmetric bifurcation diagrams of Figures 7(b) and 7(c), which govern the evolution of system (2) for m=0.10and0.125, respectively. These diagrams are obtained for the initial conditions ±0.75,0,0. The asymmetry observed in these diagrams is due to the destruction of the system’s symmetry. This is caused by a nonsymmetric displacement of nontrivial equilibrium points with respect to the values of m (see Table 2). This can also be observed on the corresponding Lyapunov exponents diagrams of Figures 8(b) and 8(c).


### 3.3. Two Parameters and Lyapunov Exponents

The changes that occur in chaotic dynamic systems may not be only related to small differences of initial conditions but also rely on the variation of parameters to which the system is sensitive. Figure13 shows an interesting way to represent the behavior of system (2) according to its parameters b,k,andwf through the maximum Lyapunov exponents λmax [52, 53]. In this figure, we can appreciate behavioral changes marked by transitions, from periodic zones (where λmax<0) in blue to chaotic zones (where λmax>0) in red (see Figures 13(a) and 13(b)). This way of presenting the dynamics of the system (2) thus makes it possible to prove the existence of chaotic solution for a combination of the values of its parameters.Figure 13
Behaviors of system (2) with respect to b−k,wf−k and maximum Lyapunov exponents λmax, for m=0: the periodic zones (where λmax<0) are in blue and the other zones (where λmax>0) are chaotic.
(a)(b)


## 3.1. Methodology

In this section, we do a thorough analysis of the dynamics of system (2) using numerical tools such as bifurcation diagrams, maximum Lyapunov exponent, Poincaré section, phase portraits, time series, cross section of the basin of attraction, and so on. Bifurcation diagrams are obtained from the local maxima of one of the state variables for different values of the control parameters. Maximum Lyapunov exponents are computed using Wolf‘s algorithm [51]. In fact, the Maximum Lyapunov exponent makes it possible to estimate the rate of convergence or divergence between two close trajectories. Thus, if the maximum Lyapunov exponent is (i) positive, the system is considered as chaotic because for a small perturbation, its trajectory diverges (this implies the occurrence of a strange attractor in phase space); (ii) negative, the system is considered stable because for all perturbations, it returns to its equilibrium point; (iii) null, the system is considered periodic or quasiperiodic.All the numerical integrations are based on the fourth-order Runge–Kutta algorithm with an integration steph=5.10−3, for better precision.


## 3.2. Coexistence of Bifurcations and Multiple Attractors

### 3.2.1. Consideringb as a Control Parameter with m=0

In this section, we considerb as a control parameter and m is equal to zero. The dynamics obtained is that represented by the bifurcation diagrams of Figure 1(a), carried out in the ascending direction for 0.22≤b≤0.32,wf=0.5,andk=1.6. This figure presents a very interesting phenomenon which is the coexistence of two symmetric bifurcation diagrams. This coexistence is highlighted through the superposition of two diagrams (red and blue) obtained by symmetrically changing the initial conditions 0.25,0,0⇔−0.25,0,0. They are in good agreement with the corresponding graphs of Lyapunov exponents as shown in Figure 2(a). Figure 3 shows some 3D phase portraits obtained for different values of the control parameter with initial conditions as indicated. Figures 4(a)–4(d) show a symmetric chaotic attractor, its time series, its frequency spectrum, and its Poincaré section (these representations are consistent with the chaotic behavior), respectively, for b=0.32. All these confirm the symmetric nature of the dynamic of system (2).Figure 1
Bifurcation diagrams as a function of the parameter b showing destruction of the symmetry of the system for some values ofm: (a) the symmetric system for m=0.00; (b–d) destruction of the system’s symmetry for m=0.05,m=0.10, and m=0.25, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the traces in red and blue, respectively.
(a)(b)(c)(d)Figure 2
Diagrams of the maximum Lyapunov exponents, corresponding to the bifurcation diagrams of Figure1. This figure shows the impact of the value of m on the symmetry of system (2) considering b as a control parameter.
(a)(b)(c)(d)Figure 3
Three-dimensional representationsx−y−z of trajectories in phase space, at different points of the bifurcation diagram of the symmetry system (Figure 1(a)) for different values of b: (a) cycle of period-1 for b=0.22, (b) cycle of period-2 for b=0.28, (c) cycle of period-4 for b=0.284, and (d, e) single and double chaotic bands for b=0.288 and b=0.3, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the curves in red and blue, respectively.
(a)(b)(c)(d)(e)Figure 4
Three-dimensional portrait (a), time series (b), frequency spectrum (c), and Poincaré section (d), confirming symmetric chaotic behavior of the system form=0.00 and b=0.32. The other parameters are specified in the text.
(a)(b)(c)(d)


### 3.2.2. Consideringb as a Control Parameter with m≠0

In this section, we still considerb as a control parameter but m≠0.0. With 0.22≤b≤0.32,wf=0.5,andk=1.6, bifurcation diagrams shown in Figures 1(b)–1(d) are plotted for m equal to 0.05,0.10,and0.25, respectively. Let us recall that red color is used for x0,y0,z0=0.25,0,0 and blue color is used for x0,y0,z0=−0.25,0,0. Comparing Figures 1(b)–1(d) with Figure 1(a), we can observe that bifurcation diagrams in red are no more symmetric to their counterpart (in blue). There is a gap (discontinuity) in the blue bifurcation diagrams (Figures 1(b)–1(d)), which increases as m is increased. This observation is caused by a shift in equilibrium points (see Table 2) and implies the existence of coexisting attractors. The coexistence of a period-2 limit cycle with a period-3 limit cycle (Figure 5(c)) and a chaotic attractor with a periodic attractor (Figure 5(d)), for instance, is revealed in Figure 5. Figures 6(a)–6(d) represent an asymmetric chaotic attractor, its time series, its frequency spectrum, and its Poincaré section (these representations are specific to chaotic attractors), respectively, for b=0.32. The aforementioned results indicated an increase in system’s complexity when m≠0.Table 2
Modifications of nontrivial equilibrium points for different values ofm.
Values ofmNontrivial equilibrium pointsP2P30.00(1, 0, 0)(−1, 0, 0)0.05(0.9753, 0, 0)(−1.0253, 0, 0)0.10(0.9512, 0, 0)(−1.0512, 0, 0)0.125(0.9395, 0, 0)(−1.0645, 0, 0)0.25(0.8828, 0, 0)(−1.1328, 0, 0)Figure 5
Three-dimensional representationsx−y−z of trajectories in the phase space, at different points of the bifurcation diagram of the asymmetry system (Figure 2(a)), for different values of b with m = 0.05: (a, b) cycle of period-1 for b=0.22 and cycle of period-2 for b=0.25; (c) cycles of period-4 and period-2 for b=0.2824; (d) single-band chaotic attractor with cycle of period-4 for b=0.2856; (e, f) single and chaotic dual-band attractor for b=0.2892 and b=0.3, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the traces in red and blue, respectively.
(a)(b)(c)(d)(e)(f)Figure 6
Three-dimensional portrait (a), time series (b), frequency spectrum (c), and Poincaré section (d), confirming asymmetric chaotic behavior of the system form=0.25 and b=0.32. The other parameters are specified in the text.
(a)(b)(c)(d)


### 3.2.3. Consideringk as a Control Parameter with m=0

Consideringk as a control parameter with m equal to zero, the dynamics of system (2) is governed by the superimposed bifurcation diagrams and Lyapunov exponents of Figures 7(a) and 8(a), respectively. These diagrams are obtained by varying the control parameter (downward) in the range 2≤k≤3 for initial conditions of opposite signs ±0.75,0,0 with wf=0.5andb=0.41. Figure 9 shows an enlargement of the diagram of Figure 7(a) (in red) in the intervals 2.6575≤k≤2.68 (Figure 9(a)) and 2.671≤k≤2.672 (Figure 9(b)) illustrating parallel branches. They show the coexistence of four (two periodic and two chaotic) and six attractors occurring in symmetric pairs. Figures 10 and 11 illustrate the 2D representations of these different attractors in the phase space. These attractors are obtained for the values of k=2.667and2.6714, respectively. The different values of the initial conditions x0,y0,z0 are as follows: for the coexistence of the four attractors—period-1 ±0.12,0,0 and chaotic spiral ±0.39,0,0; for the coexistence of the six attractors—n-finite periods ±0.28,0,0, period-6±0.06,0,0, and period-1 ±0.108,0,0. These multistabilities are well illustrated on the cross sections of the basins represented in Figures 12(a)–12(d), where we can observe several regions distinguished by different colors that correspond to each attractor.Figure 7
Bifurcation diagrams as function of the parameterk, showing destruction of the symmetry of the system for some values of m: (a) the symmetric system for m = 0; (b, c) destruction of the system symmetry for m = 0.10 and m=0.125, respectively. The initial conditions are x0,y0,z0=0.75,0,0 and −0.75,0,0 for the traces in blue and red, respectively.
(a)(b)(c)Figure 8
Diagrams of the maximum Lyapunov exponents corresponding to the bifurcation diagrams of Figure7 showing the impact of the value of m on the symmetry of system (2) considering k as a control parameter.
(a)(b)(c)Figure 9
Enlargement of the bifurcation diagram of Figure7(a), in the intervals 2.6575 ≤ k ≤ 2.68 (a) and 2.671 ≤ k ≤ 2.672 (b), revealing the coexistence of attractors in the form of parallel branches.
(a)(b)Figure 10
Coexistence of four attractors at the pointk = 2.66 of Figure 9(a), for different values of the initial conditions (x(0), y(0), z(0)): (a) a period-1 attractor and its symmetry for (±0.12, 0, 0); (b) a chaotic spiral attractor and its symmetry for (±0.39, 0, 0).
(a)(b)Figure 11
Coexistence of six periodic attractors at pointk = 2.6714 of Figure 9(b), for different values of the initial conditions (x(0), y(0), z(0)): (a) a finite n-period attractor and its symmetry for (±0.028, 0, 0); (b) a period-6 attractor and its symmetry for (±0.06, 0, 0); (c) a period-1 attractor and its symmetry for (±0.108, 0, 0).
(a)(b)(c)Figure 12
Cross sections of the basins of attraction in thex0−y0 and y0−z0 plane: on (a, b) blue and yellow colors indicate regions of period-1 while green and magenta indicate regions of chaotic attractors; on (c, d) magenta and yellow colors indicate the regions of n-finite periods, cyan and black indicate the regions of period-6 while green and blue the regions indicate of period-1. Regions of unbounded behavior are in red.
(a)(b)(c)(d)


### 3.2.4. Consideringk as a Control Parameter with m≠0

Keepingk2≤k≤3 as the control parameter, with m≠0m=0.10,0.125, wf=0.5,andb=0.41, we obtain the asymmetric bifurcation diagrams of Figures 7(b) and 7(c), which govern the evolution of system (2) for m=0.10and0.125, respectively. These diagrams are obtained for the initial conditions ±0.75,0,0. The asymmetry observed in these diagrams is due to the destruction of the system’s symmetry. This is caused by a nonsymmetric displacement of nontrivial equilibrium points with respect to the values of m (see Table 2). This can also be observed on the corresponding Lyapunov exponents diagrams of Figures 8(b) and 8(c).


## 3.2.1. Consideringb as a Control Parameter with m=0

In this section, we considerb as a control parameter and m is equal to zero. The dynamics obtained is that represented by the bifurcation diagrams of Figure 1(a), carried out in the ascending direction for 0.22≤b≤0.32,wf=0.5,andk=1.6. This figure presents a very interesting phenomenon which is the coexistence of two symmetric bifurcation diagrams. This coexistence is highlighted through the superposition of two diagrams (red and blue) obtained by symmetrically changing the initial conditions 0.25,0,0⇔−0.25,0,0. They are in good agreement with the corresponding graphs of Lyapunov exponents as shown in Figure 2(a). Figure 3 shows some 3D phase portraits obtained for different values of the control parameter with initial conditions as indicated. Figures 4(a)–4(d) show a symmetric chaotic attractor, its time series, its frequency spectrum, and its Poincaré section (these representations are consistent with the chaotic behavior), respectively, for b=0.32. All these confirm the symmetric nature of the dynamic of system (2).Figure 1
Bifurcation diagrams as a function of the parameter b showing destruction of the symmetry of the system for some values ofm: (a) the symmetric system for m=0.00; (b–d) destruction of the system’s symmetry for m=0.05,m=0.10, and m=0.25, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the traces in red and blue, respectively.
(a)(b)(c)(d)Figure 2
Diagrams of the maximum Lyapunov exponents, corresponding to the bifurcation diagrams of Figure1. This figure shows the impact of the value of m on the symmetry of system (2) considering b as a control parameter.
(a)(b)(c)(d)Figure 3
Three-dimensional representationsx−y−z of trajectories in phase space, at different points of the bifurcation diagram of the symmetry system (Figure 1(a)) for different values of b: (a) cycle of period-1 for b=0.22, (b) cycle of period-2 for b=0.28, (c) cycle of period-4 for b=0.284, and (d, e) single and double chaotic bands for b=0.288 and b=0.3, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the curves in red and blue, respectively.
(a)(b)(c)(d)(e)Figure 4
Three-dimensional portrait (a), time series (b), frequency spectrum (c), and Poincaré section (d), confirming symmetric chaotic behavior of the system form=0.00 and b=0.32. The other parameters are specified in the text.
(a)(b)(c)(d)


## 3.2.2. Consideringb as a Control Parameter with m≠0

In this section, we still considerb as a control parameter but m≠0.0. With 0.22≤b≤0.32,wf=0.5,andk=1.6, bifurcation diagrams shown in Figures 1(b)–1(d) are plotted for m equal to 0.05,0.10,and0.25, respectively. Let us recall that red color is used for x0,y0,z0=0.25,0,0 and blue color is used for x0,y0,z0=−0.25,0,0. Comparing Figures 1(b)–1(d) with Figure 1(a), we can observe that bifurcation diagrams in red are no more symmetric to their counterpart (in blue). There is a gap (discontinuity) in the blue bifurcation diagrams (Figures 1(b)–1(d)), which increases as m is increased. This observation is caused by a shift in equilibrium points (see Table 2) and implies the existence of coexisting attractors. The coexistence of a period-2 limit cycle with a period-3 limit cycle (Figure 5(c)) and a chaotic attractor with a periodic attractor (Figure 5(d)), for instance, is revealed in Figure 5. Figures 6(a)–6(d) represent an asymmetric chaotic attractor, its time series, its frequency spectrum, and its Poincaré section (these representations are specific to chaotic attractors), respectively, for b=0.32. The aforementioned results indicated an increase in system’s complexity when m≠0.Table 2
Modifications of nontrivial equilibrium points for different values ofm.
Values ofmNontrivial equilibrium pointsP2P30.00(1, 0, 0)(−1, 0, 0)0.05(0.9753, 0, 0)(−1.0253, 0, 0)0.10(0.9512, 0, 0)(−1.0512, 0, 0)0.125(0.9395, 0, 0)(−1.0645, 0, 0)0.25(0.8828, 0, 0)(−1.1328, 0, 0)Figure 5
Three-dimensional representationsx−y−z of trajectories in the phase space, at different points of the bifurcation diagram of the asymmetry system (Figure 2(a)), for different values of b with m = 0.05: (a, b) cycle of period-1 for b=0.22 and cycle of period-2 for b=0.25; (c) cycles of period-4 and period-2 for b=0.2824; (d) single-band chaotic attractor with cycle of period-4 for b=0.2856; (e, f) single and chaotic dual-band attractor for b=0.2892 and b=0.3, respectively. The initial conditions are x0,y0,z0=0.25,0,0 and −0.25,0,0 for the traces in red and blue, respectively.
(a)(b)(c)(d)(e)(f)Figure 6
Three-dimensional portrait (a), time series (b), frequency spectrum (c), and Poincaré section (d), confirming asymmetric chaotic behavior of the system form=0.25 and b=0.32. The other parameters are specified in the text.
(a)(b)(c)(d)


## 3.2.3. Consideringk as a Control Parameter with m=0

Consideringk as a control parameter with m equal to zero, the dynamics of system (2) is governed by the superimposed bifurcation diagrams and Lyapunov exponents of Figures 7(a) and 8(a), respectively. These diagrams are obtained by varying the control parameter (downward) in the range 2≤k≤3 for initial conditions of opposite signs ±0.75,0,0 with wf=0.5andb=0.41. Figure 9 shows an enlargement of the diagram of Figure 7(a) (in red) in the intervals 2.6575≤k≤2.68 (Figure 9(a)) and 2.671≤k≤2.672 (Figure 9(b)) illustrating parallel branches. They show the coexistence of four (two periodic and two chaotic) and six attractors occurring in symmetric pairs. Figures 10 and 11 illustrate the 2D representations of these different attractors in the phase space. These attractors are obtained for the values of k=2.667and2.6714, respectively. The different values of the initial conditions x0,y0,z0 are as follows: for the coexistence of the four attractors—period-1 ±0.12,0,0 and chaotic spiral ±0.39,0,0; for the coexistence of the six attractors—n-finite periods ±0.28,0,0, period-6±0.06,0,0, and period-1 ±0.108,0,0. These multistabilities are well illustrated on the cross sections of the basins represented in Figures 12(a)–12(d), where we can observe several regions distinguished by different colors that correspond to each attractor.Figure 7
Bifurcation diagrams as function of the parameterk, showing destruction of the symmetry of the system for some values of m: (a) the symmetric system for m = 0; (b, c) destruction of the system symmetry for m = 0.10 and m=0.125, respectively. The initial conditions are x0,y0,z0=0.75,0,0 and −0.75,0,0 for the traces in blue and red, respectively.
(a)(b)(c)Figure 8
Diagrams of the maximum Lyapunov exponents corresponding to the bifurcation diagrams of Figure7 showing the impact of the value of m on the symmetry of system (2) considering k as a control parameter.
(a)(b)(c)Figure 9
Enlargement of the bifurcation diagram of Figure7(a), in the intervals 2.6575 ≤ k ≤ 2.68 (a) and 2.671 ≤ k ≤ 2.672 (b), revealing the coexistence of attractors in the form of parallel branches.
(a)(b)Figure 10
Coexistence of four attractors at the pointk = 2.66 of Figure 9(a), for different values of the initial conditions (x(0), y(0), z(0)): (a) a period-1 attractor and its symmetry for (±0.12, 0, 0); (b) a chaotic spiral attractor and its symmetry for (±0.39, 0, 0).
(a)(b)Figure 11
Coexistence of six periodic attractors at pointk = 2.6714 of Figure 9(b), for different values of the initial conditions (x(0), y(0), z(0)): (a) a finite n-period attractor and its symmetry for (±0.028, 0, 0); (b) a period-6 attractor and its symmetry for (±0.06, 0, 0); (c) a period-1 attractor and its symmetry for (±0.108, 0, 0).
(a)(b)(c)Figure 12
Cross sections of the basins of attraction in thex0−y0 and y0−z0 plane: on (a, b) blue and yellow colors indicate regions of period-1 while green and magenta indicate regions of chaotic attractors; on (c, d) magenta and yellow colors indicate the regions of n-finite periods, cyan and black indicate the regions of period-6 while green and blue the regions indicate of period-1. Regions of unbounded behavior are in red.
(a)(b)(c)(d)


## 3.2.4. Consideringk as a Control Parameter with m≠0

Keepingk2≤k≤3 as the control parameter, with m≠0m=0.10,0.125, wf=0.5,andb=0.41, we obtain the asymmetric bifurcation diagrams of Figures 7(b) and 7(c), which govern the evolution of system (2) for m=0.10and0.125, respectively. These diagrams are obtained for the initial conditions ±0.75,0,0. The asymmetry observed in these diagrams is due to the destruction of the system’s symmetry. This is caused by a nonsymmetric displacement of nontrivial equilibrium points with respect to the values of m (see Table 2). This can also be observed on the corresponding Lyapunov exponents diagrams of Figures 8(b) and 8(c).


## 3.3. Two Parameters and Lyapunov Exponents

The changes that occur in chaotic dynamic systems may not be only related to small differences of initial conditions but also rely on the variation of parameters to which the system is sensitive. Figure13 shows an interesting way to represent the behavior of system (2) according to its parameters b,k,andwf through the maximum Lyapunov exponents λmax [52, 53]. In this figure, we can appreciate behavioral changes marked by transitions, from periodic zones (where λmax<0) in blue to chaotic zones (where λmax>0) in red (see Figures 13(a) and 13(b)). This way of presenting the dynamics of the system (2) thus makes it possible to prove the existence of chaotic solution for a combination of the values of its parameters.Figure 13
Behaviors of system (2) with respect to b−k,wf−k and maximum Lyapunov exponents λmax, for m=0: the periodic zones (where λmax<0) are in blue and the other zones (where λmax>0) are chaotic.
(a)(b)


## 4. Circuit Design and PSpice Simulations

The theoretical and numerical predictions previously obtained are verified in this section using an electronic analog of system (2). The circuit diagram of the proposed electronic analog computer is shown in Figure 14. It involves three integrators, an inverter, a subtractor, and two multipliers which are at the origin of the nonlinear terms.Figure 14
Circuit realization of the three-dimensional autonomous Duffing–Holmes like oscillator described by system (2) where Uii=1,2,3,4,5 are TL084CN and Ujj=6,7 are AD633JN.Using Kirchhoff’s laws, equations governing the circuit can be obtained as(12)1RCdxdτ=1RCy,1RCdydτ=1RCy−1RCx3+RbRCy−RkRCz−RmRCx2,1RCdzdτ=RwfRCy−z,with the following considerations:(13)x=X,y=Y,z=Z,τ=tRC.Equation (12) becomes(14)dXdt=1RC1Y,dYdt=1RC2X−1RC2X3+RbRC2Y−RkRC2Z−RmRC2X2,dZdt=RwfRC2Y−Z.Systems (14) and (2) are equivalent for(15)b=RbR,k=RkR,wf=RwfR,m=RmR.Let us carry out analog simulations of the circuit’s dynamics in PSpice. For these simulations, the following configurations of the transient analysis are adopted. Print step: 200 ns; final step: 500 ms; no-print delay: 480 ms; step ceiling: 4µs. For electronic components, C1=C2=C3=10nF and R=10kΩ while Rb,Rk,Rm, and Rwf are tunable. The different portraits presented in Figure 14 correspond to the transitions to chaos of the symmetric system obtained by fixing Rm=1MΩ,Rwf=20kΩ, and Rk=6.25kΩ and by varying Rb (see Table 3) for the initial conditions X0,Y0,Z0=±0.25V,0V,0V. We can observe similarities between phase portraits of Figure 15 with those of Figure 3 obtained from the bifurcation diagram of the symmetric system (Figure 1(a)). Figure 16 shows the transitions of the asymmetric system to chaos obtained for Rm=200kΩ,Rwf=20kΩ,Rk=6.25kΩ by varying Rb (see Table 3) with the initial conditions X0,Y0,Z0=±0.15V,0V,0V. These transitions are similar to those of Figure 5 obtained by numerical integrations at different points of the bifurcation diagram of the asymmetric system (Figure 1(b)). The coexistence of the four and six attractors (see Figures 17 and 18) was obtained for Rk=3.746kΩ and Rk=3.745kΩ, respectively (see Table 3), with Rm=1MΩ,Rwf=20kΩ,Rb=24.39kΩ. Concerning the values of the initial conditionsX0,Y0,Z0, they are ±0.13V,0V,0V and±0.15V,0V,0V for the four attractors and ±0.11V,0V,0V, ±0.25V,0V,0V, 0.17V,0V,0V, and −0.88V,0V,0V for the six attractors.Table 3
Comparison between the results of numerical integrations and PSpice.
SystemAttractorsNumerical integration resultsPSpice resultsbRb=10/bkΩkRk=10/kkΩRbkΩRkkΩSymmetricm=0.0,Rm=1MΩPeriod 10.2245.45451.66.25456.25Period 20.2835.714235.6Period 40.28435.211235.4Single band chaos0.28834.722235Double band chaos0.333.333333.33Coexistence of four0.4124.39022.663.759324.393.746Coexistence of six2.67143.743324.393.745Asymmetricm=0.05,Rm=200kΩPeriod 10.2245.45451.66.25376.25Period 20.254036Period 4-period 20.282435.410735.6Single band chaos-period 20.285635.0135.2Single band chaos0.289234.578134.8Double band chaos0.333.333334.631Figure 15
Phase portraits obtained from PSpice forRm=1MΩ (corresponding to the symmetric case) and different values of the control resistance Rb showing (a) cycle of period-1 for Rb=45kΩ, (b) cycle of period-2 for Rb=35.6kΩ, (c) cycle of period-4 for Rb=35.4kΩ, (d) single chaotic band for Rb=35kΩ, and (e) double-band chaotic attractor Rb=33.33kΩ. The pairs are obtained for X0,Y0,Z0=±0.25,0,0.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Figure 16
Phase portraits obtained from PSpice forRm=200kΩ (corresponding to the asymmetric case) and different values of the control resistance Rb showing (a, b) cycle of period-1 for Rb=37kΩ, (c, d) cycle of period-2 for Rb=36kΩ; (e) and (f) for cycle of period-4 and period-2 respectively with Rb=35.6kΩ; (g) and (h) for single-band chaotic attractor and cycle of period-4 respectively with Rb=35.2kΩ; (i, j) for single-band chaotic attractor with Rb=34.8kΩ; (k) and (l) for two asymmetric double-band chaotic attractor with Rb=34.631kΩ. The pairs are obtained for X0,Y0,Z0=±0.15,0,0.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 17
Representations of the four coexisting attractors obtained in PSpice withRm = 1MΩ, Rk = 3.746kΩ for different values of the initial conditions (X(0), Y(0), Z(0)): (a) cycle of period-1 and (b) its symmetry for (±0.13, 0, 0); (c) chaotic attractor and (d) its symmetry for (±0.15, 0, 0).
(a)(b)(c)(d)Figure 18
Representation of the six coexistence attractors obtained in PSpice withRm=1MΩ, Rk=3.745kΩ for different values of the initial conditions X0,Y0,Z0: (a) cycle of n-finite period and (b) its symmetry for ±0.25,0,0; (c) cycle of period-6 and (d) its symmetry for +0.17,0,0 and −0,88,0,0, respectively; (e) cycle of period-1 and (f) its symmetry for ±0.11,0,0.
(a)(b)(c)(d)(e)(f)


## 5. Total Amplitude Control (TAC)

Amplitude control of chaotic signals is important in engineering applications where the amplitude is desired for signal generation and transmission [54–58]. It makes it possible to get the attractor large or small by changing one or more variables in a range, without changing the dynamic and topological properties of the attractor [59]. The amplitude control can be total (TAC), in this case all variables are controlled in a linked way. When all the variables are controlled separately, we speak of composite amplitude control (CAC). In the case where only the amplitude of one or of n−1n≥3 variables is controlled, it is a partial amplitude control (PAC) [56, 59].Referring to [56, 59–61], it is possible to perform a total amplitude control (TAC) of system (2) by introducing an amplitude control function fr on the cubic term, which leads to the following system:(16)x˙=y,y˙=x−1frx3+by−kz−mx2,z˙=wfy−z.The functionfr makes it possible to control the amplitudes of the variables x, y, and z according to 1/frfr>0.Proof.
Let:u=x/fr, v=y/fr and w=z/fr. The resulting system (17) considering m=0 is identical to system (2):(17)u˙=v,v˙=u−u3+bv−kw−mfru2,w˙=wfv−w.
So, the functionfr controls all the amplitudes according to 1/fr. Consider fr=1/rr>0, with r being the total amplitude control parameter. By acting on the amplitude control parameter rr>0, the variables are modified while leaving the dynamics of the attractor unchanged. Figure 19(a) shows the reductions in the amplitudes of the x, y, and z variables of the dual-band chaotic attractor when increasing the control parameter r. Figure 19(b) shows the largest Lyapunov exponent diagram, which remains unchanged when the control parameter varies, which justifies the invariance of the properties of the attractor during amplitude control. The phase portraits of this dual-band chaotic attractor are shown in Figure 20 in the x−y, x−z, and y−z planes for the discrete values of r. These representations come to validate the total amplitude control (TAC) because the amplitudes of the three variables are reduced without altering the dynamic and topological properties of the attractor. This control can be achieved by introducing a potentiometer in the electronic circuit. This is convenient or beneficial in the context of secure communication as well as in various fields of information processing in engineering [54–58].Figure 19
Bifurcation diagrams (a) and largest Lyapunov exponent (b), illustrating the change in amplitude of the variablesx, y, and z (in red, blue, and black, respectively) as a function of the parameter r, without alteration of the properties of the system. The parameters b, wf, k, and m are set to 0.3, 0.5, 1.6, and 0.0, respectively.
(a)(b)Figure 20
Representation of the amplitude control of the variables, in the phase space, for three different values of the amplitude control parameter: the amplitudes of the curves in blue, red, and black are obtained for 0.5, 0.7, and 0.9, respectively. The other parametersb, wf, k, and m are set to 0.3, 0.5, 1.6, and 0.0, respectively.
(a)(b)(c)


## 6. Adaptive Synchronization

The importance of synchronizing chaotic and/or hyperchaotic systems could be its application in areas of engineering such as secure communication [62–66]. To date, there are several synchronization schemes that have already been developed with notable applications [67–71]. Among these, one of the best known methods is the adaptive method, which we will use in this section.Consider the master and slave subsystems described by expressions (18) and (19), respectively:(18)x˙1=y1,y˙1=x1−x13−mx12+by1−kz1,z˙1=wfy1−z1,(19)x˙2=y2+u1,y˙2=x2−x23−mx22+by2−kz2+u2,z˙2=wfy2−z2+u3,where x1,y1,z1 and x2,y2,z2 represent the state variables of the master and slave, respectively. uii=1,2,3 are the nonlinear controllers to choose so that the subsystems (18) and (19) synchronize.The synchronization errors between (18) and (19) are given by(20)e˙1=x˙2−x˙1=e2+u1,e˙2=y˙2−y˙1=e1−x23−x13−mx22−x12+be2−ke3+u2,e˙3=z˙2−z˙1=wfe2−e3+u3,where(21)e2=y2−y1,e3=z2−z1.For reasons of simplicity, consider the expressions of the controllers described by(22)u1=−e2−k1e1,u2=−e1+x23−x13+mx22−x12−b^e2+k^e3−k2e2,u3=−w^fe2−e3−k3e3,with b^, w^f, and k^ being the estimated parameters of b, wf, and k, respectively. kii=1,2,3 is the positive feedback gain.By substituting (22) in (19) and (20), we obtain the following expressions for the slave subsystem (23) and the synchronization errors (24):(23)x˙2=y1−k1e1,y˙2=x1−x13−mx12+y2eb+b^y1−z2ek−k^z1−k2e2,z˙2=y2−z2ewf+w^fy1−z1−k3e3,(24)e˙1=−k1e2,e˙2=e2eb−e3ek−k2e2,e˙3=e2−e3ewf−k3e3,where(25a)eb=b−b^,ewf=wf−w^f,ek=k−k^,such as(25b)e˙b=−b^˙,e˙wf=−w^˙f,e˙k=−k^˙.The synchronization problem amounts to impose that the synchronization erroreτ⟶0 (where e=e1,e2,e3T) when τ⟶∞ (asymptotically stable around its point of equilibrium) [72]. For this, let us consider the following positive definite function as Lyapunov candidate function:(26)Ve1,e2,e3,eb,ewf,ek=12e12+e22+e32+eb2+ewf2+ek2.By differentiating the Lyapunov candidate function (26) along the error trajectory (27), we have(27)dVdτ=e1e˙1+e2e˙2+e3e˙3+ebe˙b+ewfe˙wf+eke˙k,which yields(28)V˙=−k1e12−k2e22−k3e32+e˙b+e22eb+e˙wf+e3e2−e3ewf+ek−e2e3ek.The subsystems (18) and (19) will be synchronized if (29) is satisfied.(29)V˙<0.The condition for (29) to be satisfied is specified by the following expression:(30)e˙b=−e22,e˙wf=−e3e2−e3,e˙k=e2e3.Figure21 shows the evolutions of the estimated parameters b,wf, and k converging asymptotically to the values 0.32,0.5, and 1.6, respectively. On the other hand, Figure 22 shows evolutions of the errors e1,e2, and e3 which converge asymptotically towards zero with time. These evolutions prove that the subsystems (18) and (19) are well synchronized. For these numerical simulations, we have considered the initial conditions x10,y10,z10=1,−1,1 and x20,y20,z20=1,0,1 with the values of the gains k1,k2,k3=1,1,0. For parameter values, b=0.32,wf=0.5,m=0, and k=1.6. At the starting point, estimated parameters areb^0,w^f0,k^0=−0.15,−1,−1.Figure 21
Temporal evolutions of the estimated gains ofb,wf, and k converging asymptotically to (a) 0.32 (in blue), (b) 0.5 (in green), and (c) 1.6 (in red), respectively.
(a)(b)(c)Figure 22
Temporal evolutions of the errors between the two subsystems (11) and (12) (master and slave), starting from different initial conditions: convergence of the errors (a) e1 (in green), (b) e2 (in blue), and (c) e3in red towards zero.
(a)(b)(c)


## 7. Conclusion

This work focused on the effects of symmetric and asymmetric nonlinearity on the dynamics of a third-order chaotic system, namely, the Duffing–Holmes autonomous oscillator. This has been investigated by supplementing a parametric quadratic termmx2 to the cubic nonlinear term −x3 of the existing third-order autonomous Duffing–Holmes system [42]. An electronic circuit analog to the proposed system has also been designed. Using the Routh–Hurwitz criteria, stability conditions of the fixed points have been established and the existence of a Hopf bifurcation is obtained. The dynamics of the system has been studied by considering two control parameters (b and k) for discrete values of the parameter m. The parameter m made it possible to adjust the symmetry of the system by modifying nontrivial equilibrium points. With bas control parameter, the system exhibits a coexistence of symmetric bifurcation diagrams, confirmed by superposition of the Lyapunov exponent diagrams, phase portraits, and others when m=0. When m≠0, symmetry is destroyed. With k as control parameter, dynamics of the system is still symmetric in the case of m=0. This has been confirmed by the superposition of bifurcation diagrams obtained from different initial conditions. With this, the coexistence of four and six attractors in a set of symmetric pairs has been revealed at certain points of the bifurcation diagrams. The cross sections of the basins of attraction have been represented in these different points of coexistence. For the case of m≠0 (with k as control parameter), the system exhibits asymmetry behavior confirmed by bifurcation diagrams obtained from different initial conditions. These symmetric breakdowns (asymmetry) are related to the movements of the nontrivial equilibria, caused by the values of the parameter m. The representations of the dynamics of the system according to its parameters and the Lyapunov exponents revealed the presence of periodic and chaotic regions. The results from the analog simulations in PSpice show good agreement with the numerical results. The results obtained under total amplitude control (TAC) and synchronization analyses conducted using the adaptive method prove that this new system is suitable for applications in various fields of engineering such as encryption of images [73], secure communication [74], and random bit generator [75].The model reported in this work is a typical example of a 3D system with three rest points and a nonlinear function with a region of negative slope. Also, we conjecture that the results obtained in this work may also be found in the Chua circuit, Shinriki oscillator, and jerk system (with cubic, hyperbolic sine, or hyperbolic tangent) just to cite a few, hence the pertinence and relevance of our study.



---

*Source: 8891816-2020-12-14.xml*