Chaotic Dynamics of a Mixed Rayleigh–Liénard Oscillator Driven by Parametric Periodic Damping and External Excitations

Department of Industrial and Technical Sciences, LARPET, ENSET-Lokossa, Université Nationale des Sciences, Technologies, Ingénierie et Mathématiques-Abomey, Abomey, Benin Département de Génie Mécanique et Productique (GMP), Institut National Supérieur de Technologie Industrielle (INSTI) Lokossa/UNSTIM-Abomey, Abomey, Benin Department of Physics, FAST-Natitingou, Université Nationale des Sciences, Technologies, Ingénierie et Mathématiques-Abomey, Abomey, Benin Department of Physics, ENS-Natitingou, Université Nationale des Sciences, Technologies, Ingénierie et Mathématiques-Abomey, Abomey, Benin


Introduction
Many dynamical systems in various scientific fields are represented in terms of nonlinear ordinary differential equations which appear in a vast range of applications [1][2][3][4][5][6][7]. One of the classes of equations widely used is the class of mixed Liénard-type equations that has the following form: where a dot represents the derivative with respect to time t. f, g, and h are arbitrary functions of x. It is important to notice that, from equation (1), when f(x) � 0, the classical Liénard-type equation is obtained. is latter class of equations plays an important role in science and engineering [8][9][10][11], and its mathematical properties have been intensively studied in the literature from mathematical and physical points of view [12][13][14][15].
Recently, the nonlinear dynamics of resonant RLC series circuit has been described using a general class of mixed Rayleigh-Liénard oscillators [16]. is new class of oscillators contains the class of oscillators given by equation (1). e class of mixed Rayleigh-Liénard oscillators is widely used in physics, electronics, mathematics, biology, engineering, and many other disciplines. In recent years, the nonlinear dynamics of this class of oscillators has been intensively studied, and interesting results such as perioddoubling leading to chaotic motion, strange attractors, reverse period-doubling bifurcation, symmetry breaking, antimonotonicity, existence of horseshoe chaos, and so on have been obtained [17][18][19][20][21][22][23][24][25]. In most of these studies, the Melnikov perturbation method [18,24,26,27] has been widely used to detect chaotic dynamics and to analyze nearhomoclinic motion with deterministic or random perturbation. is method is today considered as a powerful analytical tool to provide an approximate criterion for the occurrence of hetero/homoclinic chaos in a wide class of dynamical systems. On the other hand, in these nonlinear dynamics works, the considered nonlinear oscillators are usually subjected to linear and nonlinear damping with constant coefficients. is damping combination is necessary since experiments showed a strong nonlinear dependence of damping with the vibration amplitude of nonlinear vibrating systems. To this end, several phenomenological nonlinear damping models have been introduced into existing fundamental nonlinear oscillators to understand the large amplitude, quasi-periodic and chaotic vibrations of nonlinear vibrating systems [24,[28][29][30][31][32]. In most cases, the phenomenological nonlinear damping does not reproduce with accuracy experimental results. us, the problem of finding nonlinear dissipative force capable to reproduce with accuracy experimental results has been recently studied. Amabili [33] used a fractional standard linear solid model and first-order harmonic balance method, and derived for the first time a nonlinear dissipative force of the type x 2 _ x for the systems with only cubic stiffness. Recently again, Amabili [34] derived for the first time, from classical viscoelasticity and using the same method, a nonlinear dissipative force containing linear damping ( _ x), quadratic damping (x _ x), and cubic damping (x 2 _ x) terms. In addition, this author showed that in case of small damping, quadratic damping term disappears and the nonlinear damping is of the type x. In [35], Amabili showed that a mathematical model with linear damping and nonlinear damping of the type x 2 _ x for a continuous system is capable to reproduce with accuracy experimental results for geometrically nonlinear vibrations. ese results originally introduced by Amabili [34,35] have been confirmed by Balasubramanian et al. [36]. erefore, for a better understanding on the dynamical properties of the real systems, the dissipative force to choose must incorporate at least linear, quadratic, and cubic damping terms.
Although the problem of identification of nonlinear damping is known, the effect of parametric damping on the dynamical behavior of nonlinear nonautonomous systems constitutes another interesting research topic. In this perspective, Mokni et al. [37] investigated the effect of cubic nonlinear parametric damping on vibration isolation of a sdof system, and the obtained results have revealed that the periodic nonlinear viscous damping enhances vibration isolation more than nonlinear time-dependent damping. On the other hand, the effect of linear parametric damping has been numerically studied by Rajasekar and Lakshmanan [38] in chaos control of an Bonhoeffer-Van der Pol oscillator. Parametric excitation with periodic damping has been examined for energy-harvesting applications [39,40]. Kwuimy et al. [41] studied the effect of tilted harmonic excitation and parametric damping on the chaotic dynamics in an asymmetric magnetic pendulum. e obtained results have shown that the presence of parametric damping without a periodic fluctuation can enhance or suppress chaos, while a parametric damping with a periodic fluctuation can increase the region of regular motions significantly.
It is then important to underline that the study of nonlinear oscillators under the influence of applied periodic force with the presence of the parametric periodic damping has attracted attention of very few investigators. Moreover, the effects of the parametric periodic damping and periodic excitation on the chaotic behavior of a mixed Rayleigh-Liénard oscillator with symmetric fourth-order polynomial double-well potential have not been discussed yet.
erefore, studying chaotic dynamics of such an oscillator must be performed for a better understanding of the dynamical characteristics of this system. In this perspective, we consider the following mixed Rayleigh-Liénard oscillator driven by parametric periodic damping and periodic external excitations: where ω 0 , c, α 0,1 , β 0,1 , F 0,1 , η, ], and ω are constant parameters. Physically, ω 0 represents the natural frequency, c denotes the cubic stiffness parameter, α 0 means the linear damping coefficient, and α 1 and β 0,1 designate the nonlinear damping coefficients. F 0 denotes the constant excitation amplitude. η and F 1 represent the amplitude of the parametric damping and periodic external excitations, respectively. ] and ω represent, respectively, the frequency of the parametric and external periodic excitations. From equation (2), one can note that the nonlinear dissipative force takes into account at least linear, quadratic, and cubic damping terms required for a good description of the experimental results.
Another motivation for our interest in this system is that it has a wide range of applications in various scientific fields. When η � α 1 � β 0,1 � F 0 � 0, equation (2) serves as a basic model for self-excited oscillations in various disciplines. For α 1 � β 0,1 � F 0,1 � 0, equation (2) arises in energy-harvesting applications [39,40]. On the other hand, the system given by equation (2) belongs to a general class of mixed Rayleigh-Liénard oscillators recently proposed for describing the nonlinear dynamics of resonant RLC series circuit [16]. In nonlinear electronic circuits, note that the memristor represents the fourth basic circuit element. It is described by flux and charge [42,43]. As a new device with memory, the memristor is used due to its potential applications, in many fields including nonvolatile memories on the nanoscale, neuromorphic systems and chaotic circuits [44,45]. e study on memristive nonlinear electronic circuits is one of the active topics of research in the past few years, since 2 Complexity memristors with versatile nonlinearities can be used in memristor-based oscillators, information encryption, memory, and so on. Given its potential applications, memristors are conveniently introduced into some existing linear or nonlinear electronic circuits to build various novel memristive chaotic circuits [46]. In this way, numerous mathematical models with different complexities have been proposed in the literature [47][48][49][50][51][52]. On the other hand, the two-dimensional FitzHugh-Nagumo (FHN) model [53,54] is simplified from the four-dimensional Hodgkin-Huxley model for geometrical explanation of neuronal excitability and spiking under excitation. Several studies have been performed on the significant and complex dynamical aspects of the FHN model [51,[55][56][57]. rough all these studies, the results showed that the nonautonomous memristive chaotic circuits exhibit rich dynamical behaviors, such as chaos and hyperchaos, coexisting multiple attractors, hidden attractors, and hidden extreme multistability [47][48][49][50][51][52][58][59][60][61]. erefore, dynamical behaviors of nonlinear electronic circuits must be investigated from theoretical and experimental points of view. us, the main objective of this paper is to investigate analytically and numerically the effects of the parametric periodic damping and external excitations on the transitions from regular to chaotic motions. In order to attain our objective, we investigate the equilibrium stability and its evolutions (Section 2). We derive the condition of existence of horseshoes chaos using the Melnikov method, and we analyze the effect of certain parameters on chaotic behavior of the system (Section 3). Afterward, we analyze the bifurcation mechanisms and eventual routes to chaos in our system using the fourth-order Runge-Kutta ODE45 algorithm (Section 4). Finally, we end with a conclusion (Section 5).

Equilibrium Stability and Its Evolutions
e equilibrium points of system (2) can be obtained by solving _ x � _ y � 0 [46]. So, these points are expressed as where x * verifies the following equation: Equation (4) describes the dynamics of equilibrium state E of system (2) with time.
e Jacobian matrix at the equilibrium state E and for ω � ] is derived from (2) as e characteristic polynomial of matrix (5) is where m x * � ηc From equation (6), two eigenvalues can be expressed as where According to the sign of m (x * ), n(x * ), and l (x * ), the equilibrium stabilities can be separately analyzed by the eigenvalues s 1,2 .
us, the following cases can be enumerated.
In this case, the eigenvalues are both pure imaginary. erefore, a Hopf bifurcation occurs when |x * | ≺ ω 0 / �� 3c and x * satisfies the following equation: . e roots of equation (9) can be derived as [46,62] )/2. According to Cardan discriminant [46,62], when Δ ≻ 0, one positive real equilibrium point is obtained from (11), since the equilibrium point cannot be a complex number. us, the system possesses in this case, one Hopf bifurcation point (HBP). However, two equilibrium points appear in equation (9) when Δ � 0. erefore, two Hopf bifurcation points appear in this case. Now, the system considered possesses three real equilibrium points from (10)-(12) when Δ ≺ 0. en, three Hopf bifurcation points appear in this condition.
Case 2. m(x * ) ≻ 0 and n(x * ) � 0. In this case, one obtains s 1 � − m(x * ) and s 2 � 0. en, the system presents two fold bifurcation points (FBPs) at E � (±ω 0 / �� 3c , 0). Case 3. m(x * ) ≻ 0 and n(x * ) ≺ 0. ere are two real roots with opposite signs, implying that the equilibrium state E is an unstable saddle point. Complexity In this case, there are two positive real roots, indicating that the equilibrium state E is an unstable node point.
ere is a pair of complex conjugate roots with positive real parts, which means that the equilibrium state E is an unstable node focus.
Case 6. m(x * ) ≻ 0, n(x * ) ≻ 0, and l(x * ) ≺ 0. One obtains, in this case, a pair of complex conjugate roots with negative real parts, indicating that the equilibrium state E is a stable node focus.
Case 7. m(x * ) ≻ 0, n(x * ) ≻ 0, and l(x * ) ≻ 0. Equation (6) has two negative real roots, which implies that the equilibrium state E is a stable node point. For and ω � ] � 0.6, we get numerically the values of the equilibrium state E as well as the Hopf and fold bifurcation points with the time evolution, as shown in Figure 1.
e above analytical and numerical results show that with the evolution of time, the equilibrium point is variable, and its stability has complex transitions between stable points and unstable points via three Hopf bifurcations and two fold bifurcations and then initiates the next period, leading to the generation of complex dynamical behaviors.

Fixed Points and Melnikov Criterion for Homoclinic
Chaos. In this section, we first derive the fixed points and the phase portrait corresponding to the system of equation (2) when it is unperturbed. us, if we let α 0,1 � β 0,1 � F 0,1 � η � 0, the system of equation (2) becomes an unperturbed system and can be rewritten as follows: e system of equations (13) corresponds to an integrable Hamiltonian system with a potential function and the associated Hamiltonian function is given by From equation (13), we can determine the fixed points and analyze their stabilities.
us, if c ≺ 0, the system of equations (13) admits one fixed point (0, 0) which is a saddle. It has only a single-well potential. If c ≻ 0, the system of equation (12) has three fixed points (0, 0), (x 1 � − (ω 0 / � c √ ), 0), and (x 2 � (ω 0 / � c √ ), 0): ere are one saddle and two centers. In this case, it has a double-well symmetric potential. We here are interested in double-well potential. erefore, we fix c � 1 throughout this paper. e potential function and the phase portrait of the system of equations (13) are plotted in Figures 2(a) and 2(b), respectively.
Second, we now suppose that the unperturbed system discussed above is perturbed by a combination of dissipative and external excitations forces. In this section, we investigate analytically the condition for the onset of horseshoes chaos by applying the Melnikov method to equation (2). To this end, scaling α 0,1 � εα 0,1 , β 0,1 � εβ 0,1 , F 0,1 � εF 0,1 , and η � εη, the mixed Rayleigh-Liénard oscillator equation (2) can be rewritten as a first-order differential system in the form     Complexity where ε is a small parameter. It is shown that when ε � 0, the system of equations (16) is a Hamiltonian system with a pair of homoclinic trajectories connecting the saddle to itself given by with τ � t − t 0 . t 0 is the cross-section time of the Poincaré map and can be interpreted as the initial time of the forcing time. According to equation (17), the signs refer to the right and left half planes. Both solutions determine the separatrix orbit, since it separates two types of orbits in phase space. ''+'' refers to the homoclinic trajectory with x ≻ 0, and ''-'' refers to the homoclinic trajectory with x ≺ 0. Turning our attention to the perturbed system, that is, ε ≠ 0, we can compute the Melnikov integrals in order to find if the twodimensional stable and unstable manifolds of the periodic orbits intersect. In this case, the Melnikov function M (t 0 ) is a scalar, and we have where . Substituting equation (17) into (18) and evaluating the integral by using the standard integral tables [63], we obtain the Melnikov function: with If M(t 0 ) has simple zeros, then a homoclinic bifurcation occurs.
us, in the case where the amplitude of the parametric periodic damping is equal to zero (η � 0), the transverse crossings of stable and unstable manifolds exist if Melnikov function (19) has simple zeros, that is to say, erefore, the necessary condition for which the invariant manifolds intersect is given by If η ≠ 0, a necessary condition for the onset of Melnikov chaos in our system can be expressed from [64][65][66][67] as follows: with R � (A ± /C) In order to have a visual information on the range where the horseshoe chaos is possible, we have plotted in Figures 3(a)-3(d) the dependence of the amplitude η of the parametric periodic damping excitation on the frequency ω � ] for different values of F 1 , β 0 , β 1 , and F 0 . From Figure 3(a), we notice that as the amplitude F 1 of the external periodic excitation increases, the region where the horseshoes chaos appears decreases. e same observation can be made when F 0 decreases (see Figure 3(d)). Furthermore, when the nonlinear damping coefficients β 0,1 increase, horseshoes chaos domain decreases (see Figures 3(b) and  3(c)).

Fractal Basin Boundaries.
e aim of this section is to verify the validity of the analytical predictions obtained in the previous section and to study the effect of η, F 0 , and F 1 on basin of attraction of the mixed Rayleigh-Liénard oscillator under study. For this purpose, we numerically solve the equation of motion (2) by using the fourth-order Runge-Kutta ODE45 algorithm with time step 10 -3 for generating basins of attraction in a plane element .005, and ω � ] � 0.6, the Melnikov threshold for the right plane is given analytically by 1.1397. erefore, below this critical value, the system under consideration may display a chaotic behavior. However, regular behavior may appear above this critical value. Now, for η � 2 chosen in the regular behavior domain predicted by the Melnikov criterion, the basin of attraction shows regular solutions (see Figure 4(d)). When the amplitude η of the parametric periodic damping excitation is chosen in the chaotic behavior region, a fractal structure of the basin of attraction of the initial conditions appears and becomes more and more visible as η decreases (see Figures 4(a)-4(c)). We can conclude that the analytical and numerical predictions are in good agreement. e effects of F 1 and F 0 on basin of attraction are investigated through Figures 5 and 6 for η � 0.8. From Figure 5, we notice that the basin of attraction becomes less eroding as the amplitude of the external periodic excitation decreases. However, the basin of attraction is destroyed and the fractal behavior becomes more and more visible when the constant term F 0 increases (see Figure 6). It is important to point out that the geometric shape of chaotic attractors has completely changed as F 1 and F 0 increase. When ω ≠ ] and η � 0.8, the fractal behavior increases with increasing F 0 and F 1 (see Figure 7). Moreover, the geometric shape of chaotic attractors also changes when F 1 and F 0 increase.

Bifurcation Mechanisms and
Transition to Chaos e aim of this section is to investigate some bifurcation mechanisms in our system when certain system parameters evolve. For this purpose, we numerically solve the equation of motion (2) under consideration using the fourth-order Runge-Kutta ODE45 algorithm with time step 10 -2 . e initial conditions used to perform the numerical simulations are X(0) � (0.5, 0.5) and Y(0) � (− 0.5, − 0.5).
e bifurcation diagram and its Lyapunov exponent as function of the amplitude η are plotted in Figure 8. rough this figure, various bifurcations such as period doubling, quasiperiodic, periodic windows, intermittency, and chaos are exhibited by system (2). By comparing the two sets of data, one can notice from Figure 8 that system (2) displays the monostability and bistability phenomena as well as the phenomenon of coexistence of attractors. ese phenomena are illustrated in Figure 9 for several different values of η.
When the amplitude F 1 is used as control parameter, the dynamics of system (2) becomes rich. One can see that the system presents periodic windows, reverse periodic windows, periodic bubble, intermittency, and chaotic oscillations (see Figure 10). In addition, we also notice that system (2) exhibits the coexistence of attractors as well as the monostability and bistability phenomena. In order to have an idea about these attractors, various phase portraits of system (2) are plotted in Figure 11 for different values of F 1 . We have also investigated through Figure 12 the effect of the constant term F 0 on the bifurcation diagram when F 1 � 0.1 and ω � ] � 0.2. We clearly see that system (2) exhibits periodic windows, reverse periodic windows, intermittency, and chaos. We also notice that the chaotic behavior is more dominant than periodic behavior when F 0 evolves. On the other hand, the system under consideration exhibits monostability and bistability phenomena as well as the phenomenon of coexistence of attractors.
When ω � 0.2, ] � (( � 5 √ − 1)/2), and F 1 � 0.1, the bifurcation diagram and its corresponding Lyapunov exponents as F 0 evolves are shown in Figure 13. We observe that system (2) transits from a chaotic motion to a period-3 motion followed by a chaotic motion. Afterward, the chaotic behavior disappears and a period-2 orbit takes place in the system. We also notice that the periodic behavior is more dominant than chaotic behavior. By comparing the two sets of data, one can notice from Figure 13 that system (2) shows the monostability and bistability phenomena as well as the phenomenon of coexistence of attractors.

Conclusions
is paper deals with chaotic dynamics in mixed Rayleigh-Liénard oscillator driven by parametric periodic damping and external excitations. e equilibrium point and its evolutions with the time are theoretically investigated. It is found that the equilibrium stability has complex dynamical transitions between stable and unstable points via Hopf/fold bifurcations. Moreover, e Melnikov perturbation method is used to derive the analytical criterion for the appearance of chaos in the system under consideration. A convenient demonstration of the use and accuracy of the method is obtained from the basin of attraction. e effects of the parametric periodic damping and external excitations are investigated. It is found that for ω � ], the horseshoe chaos decreases and disappears as the amplitude of the parametric periodic damping increases. For η � 0.8, the basin of attraction is destroyed and the fractal behavior becomes more and more visible as the constant term F 0 increases. However, the fractality of the basin of attraction decreases when the amplitude F 1 of the external periodic excitation increases. Moreover, the geometric shape of chaotic attractors is modified when F 0 and F 1 increase. For ω ≠ ] and η � 0.8, the fractality of the basin of attraction becomes visible when the amplitude of the external periodic excitation and the constant term increase. Bifurcation diagrams, Lyapunov exponents, and phase portraits are plotted  using the fourth-order Runge-Kutta ode 45 algorithm, and period-doubling, periodic windows, reverse periodic windows, periodic bubble, intermittency, and chaotic behavior are obtained. It is important to underline that the system under consideration can vibrate from a chaotic motion to a period-3 motion followed by a chaotic motion. In addition, we also notice that the system presents monostable and bistable oscillations as well as the coexistence of attractors.

Data Availability
No data were used in this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.