# Mixed-Mode Oscillation in a Class of Delayed Feedback System and Multistability Dynamic Response

Authors: Youhua Qian; Wenjing Meng
Journal: Complexity
 (2020)

---

## Abstract

In this paper, a class of two-parameter mixed-mode oscillation with time delay under the action of amplitude modulation is studied. The investigation is from four aspects. Firstly, a parametric equation is considered as a slow variable. By the time-history diagram and phase diagram, we can find that the system generates a cluster discovery image. Secondly, the Euler method is used to discrete the system and obtain the discrete equation. Thirdly, the dynamic characteristics of the system at different time scales are discussed when the ratio of the natural frequency and the excitation frequency of the system is integer and noninteger. Fourthly, we discuss the influence of time delay on the discovery of clusters of this kind of system. The research shows that the time lag does not interfere with the influence of the cluster image, but the dynamics of the upper and lower parts of the oscillation in each period will be delayed. So, we can improve peak performance by adjusting the time lag and obtain the desired peak. Finally, we explore the multistate dynamic response of a two-dimensional nonautonomous Duffing system with higher order. According to bifurcation diagram and time-history curve, bistable state will appear in the system within the critical range. With the gradual increase of parameters, the chaotic attractor will suddenly disappear which will lead to the destruction of the bistable state.


---

## Body

## 1. Introduction

In recent years, with the rapid development of science and technology, the nonlinear problems of the actual power systems in various fields have become more and more prominent. In the national economy, national defense industry, and engineering technology, a large number of practical problems urgently need to be processed by nonlinear dynamics theory and methods, which will promote the development of nonlinear dynamics into a more comprehensive and in-depth development period.The researcher can conduct a comprehensive analysis and discussion from two important aspects of amplitude and frequency. For example, when we focus on low-frequency forces, i.e., the external excitation frequency is much smaller than the natural frequency of the original system, and the oscillator can exhibit a typical fast-slow dynamic called mixed-mode oscillation (MMO). Sadhu [1] researched the canards and mixed-mode oscillations in a singularly perturbed two predators-one prey model. Upadhyay et al. [2] studied mixed-mode oscillations and the synchronous activity in the noise-induced modified Morris–Lecar neural system. Kingston and Thamilmaran [3] discussed the bursting oscillations and mixed-mode oscillations in the driven Lienard system. Shimizu et al. [4] made a thorough exploration of mixed-mode oscillations and chaos from a simple second-order oscillator under weak periodic perturbation. The oscillating behavior appeared in the above four articles is generally expressed as a periodic state characterized by a combination of a relatively large amplitude (spike state) close to the harmonics and a small-amplitude oscillation (stationary state). Due to its complexity and diversity, the system with delayed feedback has a wide range of practical backgrounds which are always an inevitable hysteresis when studying the laws of the motion of the objective. Therefore, it has important theoretical significance and practical value in the research of delayed feedback system. Inaba et al. [5] studied the feedback control problem of network systems with discrete delay and distributed delay. Thus, unlimited distributed delays were first introduced in discrete network domains. Weicker et al. [6] focused on the rapid transition layers among the plateaus and demonstrated their contribution to the total cycle. Porte et al. [7] experimentally characterized the mechanism of strong chaos in semiconductor lasers with delayed feedback. Sun et al. [8] proposed a method for time delay identification in a multidegree-of-freedom (MDOF) linear system with multiple feedbacks. Dmitrishin et al. [9] studied the robust stability problem of the linear delayed feedback control (DFC) mechanism.As one of the research directions of nonlinear dynamics, the multitime scale plays the nonlinear essential characteristics in the perspective of dynamics. Its theoretical method has been widely applied to neuroscience, chemistry, physics, bioscience, and other fields. Therefore, nonlinear systems with multiple time scales have attracted attention of many scholars at home and abroad. The effects of slow variables on the rupture of pancreatic cells were investigated [10] on the basis of the Chay–Keizer model with three time scales. Izumiet al. [11] discussed the relationship between fast scale bifurcation and slow scale bifurcation in the discontinuous circuit. Yu et al. [12] studied the delayed feedback control problem of bursting synchronization in the small-world neural network presented by the neural network in some areas of the cerebral cortex. Yu et al. [13] studied the generation of complex cluster patterns in the Duffing oscillator with delayed feedback and proposed the symmetric fold-fold and symmetric Hopf-Hopf bursting patterns. Cornforth and Lipson [14] introduced the fast and slow analysis method and applied it to the study of multitime scale problems in nonlinear systems. Han et al. [15] proposed a general method for analyzing the mixed-mode vibration of a system with two excitation frequencies. The validity of this method was verified by the Duffing and van der Pol equations. Yang et al. [16] discussed the influence of delay coupling on bursting synchronous differential feedback control in the modularized neural network. Meng et al. [17] presented and analyzed two different types of bursting in a two-compartment neuron model with the current feedback control due to totally different generation mechanisms. Li et al. [18] investigated the Brusselator with different time scales, which behave in the classical slow-fast effect. Zhou et al. [19] established a 3D discrete system featuring a new series of complex fast-slow behaviors caused by different bursters. Ding and Li [20] studied the Rulkov model with self-inhibiting synapses and time delays and compared them with the Rulkov model without self-inhibiting synapses. Fan and Wang [21] studied the effects of different time delays and coupling intensities on the synchronization and cluster transition of Hindmarsh–Rose neuron system. Bertram and Rubin [22] described the fast-slow analysis technique and applied it to relaxation oscillations, neuronal bursting oscillations, canard oscillations, and mixed-mode oscillations. Han et al. [23] studied the dynamics of the bursting by Duffing system with multifrequency excitation. Qian and Yan [24] studied a two-degree-of-freedom nonlinear-coupled Duffing system with an external excitation and two external excitations by the fast-slow analysis method. Han et al. [25] proposed two new bursting modes, fork-shaped delay, and multifrequency excitation of the Duffing system. Different patterns of electrical bursting were proposed, and the types and generation mechanisms of these bursting oscillations were analyzed by using fast-slow dynamics. For instance, Shen et al. [26] introduced the fast and slow analysis method and applied it to the study of multitime scale problems in nonlinear systems. Yu et al. [27] studied the generation of some new cluster modes in multidelay-controlled oscillators. The bifurcation condition of the fast subsystem and its stability related to time delay were calculated. Zhang et al. [28] analyzed the effects of time scales on the dynamic behavior of the system. McKenna and Bertram [29] explained the mechanism behind the oscillation in cells by using the fast-slow analysis method. Han et al. [30] proposed an approximate frequency-truncation fast-slow analysis method to analyze the dynamics of a fast-slow system with two incommensurate excitation frequencies. Yu and Wang [31] analyzed the dynamics involving different waves in a double-well potential oscillator coupling amplitude modulation control of low frequency. Zhou et al. [32] investigated the bursting in Sprott B system with a single excitation and showed that Hopf bifurcation delay may exhibit due to the effect of slow passage through the supercritical Hopf bifurcation. Wang et al. [33] dealt with transitions through Melnikov thresholds and the corresponding fast-slow dynamics in a family of biparametric mechanical oscillators.From nature to humanity society, the phenomenon of time delay is everywhere. In natural and social phenomena, the changes and development of many systems are not only related to the current state of the system but also depend on some past state of the system. In other words, time delay is inevitable in the system. Plaut and Hsich [34] discussed parametric excitation systems with time delay, through numerical simulation, and they found that the system had very complicated dynamics. Then, through the method of multiple scales, they studied the weak nonlinear time-delay system which is only in damping and discussed the main resonance, the harmonic resonance, superharmonic resonance, and time-delayed effect on the steady-state motion frequency amplitude curve. Raghothama and Narayanan [35] used the incremental harmonic equilibrium (IHB) method to analyze the dynamic response of systems with quadratic and cubic nonlinear time-delay parameters, studied the stability of the periodic solution of the system by Floquet theory, and obtained the bifurcation diagram of the system by combining the stability analysis with the path-following algorithm with arc-length parametric continuation. Maccari [36] studied dynamic response with a time-delay state feedback of van der Pol by using the asymptotic perturbation (asymptotic perturbation) method to get the system amplitude and frequency equation of two groups of slowly varying. Ji and Leung [37] considered a parametric excitation of Duffing time-delay feedback problems through the multiscale method to study the main parameters of the resonance system and analyzed the stability of the steady-state solution. It is found that saddle bifurcation and subcritical fork bifurcation exist in the system equilibrium point.It is also an important part in the field of nonlinear dynamics to study the multistability problem of systems. Multistability means that a system is neither stable nor completely unstable but switches between two or more mutually exclusive states over time. Multistable systems are also susceptible to noise, initial conditions, or system parameters. The methods of solving the problem include the analytical method, numerical analysis method, and experimental method. In the fields of chemistry, electricity, ecology, neuroscience and so on, the characteristic of the multistable state has been widely applied and has produced the vital influence to the research and development of these disciplines. Therefore, the system containing the multistable state has been paid attention by many scholars at home and abroad. Loukaides et al. [38] verified that the multistable structure could be produced by single additive manufacturing operation through analyzing examples, numerical simulation, and physical prototype of selective laser sintering production of titanium alloy. Yang and Ma [39] systematically studied the mechanical responses of two new two-dimensional (2D) mechanical metamaterials and realized phase transition/shape reconstruction and zero Poisson’s ratio on the basis of the multistable mechanism, achieving great morphological changes. Huang and Xü [40] obtained the mathematical model by introducing time-delay feedback to a plane autonomous nonlinear system, and the results showed that time delay can not only make the system Hopf bifurcation and produce periodic vibration but also make the system appear multistable periodic motion or periodic attractor. Schmitz et al. [41] discussed an example of the application of the multistable state to chemistry, showing that the phenomenon of the multistable state is usually described by discussing the steady-state solution of a nonlinear process, which is given by an abstract mathematical model of single variable x and evolves according to the differential equation. Lai et al. [42] studied the dynamic behaviors such as multistability and bifurcation of a class of neural network systems with time delay. The results show that the system has 16 kinds of stable states and its own attractive region. Huang et al. [43] discussed the phenomenon of multistable synchronization in the synchronous region of the Kuramoto phase oscillator on a one-dimensional closed loop under the action of asymmetric coupling and further theoretically analyzed its steady-state law and steady-state stability.Here, we describe analytical and numerical studies of a class of two-parameter mechanical systems with delayed feedback:(1)x¨+x˙−αx˙−axt−τ+bx3=f1+f2cosω1tcosω2t,where a is the linear restoring parameter and α>0 is the nonlinear damping coefficient. f1>0 is the unmodulated amplitude, f2 is the degree of forcing modulation, ω1 is the modulation frequency, and ω2 is the forcing frequency.First of all, we analyze the case forf2=0. By the singular perturbation methods, equation (1) can be given by(2)x¨+x˙−αx˙−axt−τ+bx3=f1cosω2t.Letγ=f1cosω2t, and use the Euler method obtained the following discrete systems:(3)xn+1=kn,kn+1=a1xn+c1kn−b1xn3−d1xn−e1kn+γ,where a1=Δt2+τΔt, b1=bΔt6, c1=2+αΔt, d1=1+αΔt2, and e1=τΔt. Set a=b=1,α=1,ω2=0.01, and f1=1. Figure 1 shows the time-history curve and phase portraits of the system when τ=0.3. As shown in the figure, we can found that the trajectory of system (2) undergoes symmetrical folding bifurcation, and it is a typical fast and slow oscillation system.Figure 1
Phase portraits (a) and time series (b, c) of fast-slow system (2) for a=b=α=1, ω2=0.01, f1=1, and τ=0.3.
(a)
(b)
(c)From the previous developments, this paper focuses on a class of two-parameter mixed-mode vibrations with time delay under the action of amplitude modulation. The structure of this paper is as follows. InSection 2, the bifurcation of the undisturbed model will be investigated. In Section 3, we will discuss the oscillating dynamics when the natural frequency is equal to the excitation frequency. In Section 4, we will study the mixed-mode oscillation dynamics reflected by the system when the natural frequency and the resonant frequency are not equal. In Section 5, the influence of the time delay is discussed. In Section 6, we explore the multistate dynamic response of a two-dimensional nonautonomous Duffing system with higher order. In Section 7, further conclusions are presented.


## 2. Representation and Local Bifurcation of the Unperturbed Model

We consider the left side of equation (1). When τ=0.3, the unperturbed form(4)x˙=y,y˙=αy+axt−τ−bx3−y3.Using Taylor’s expansion, we havext−τ≈xt−τx′t,(5)x˙=y,y˙=αy+ax−aτy−bx3−y3.By calculation, the system has three equilibrium points:E±±a/b,0 and E0=0,0. With α>0 and a>0, from the stability analysis, we know, E± are unstable and E0 is a saddle point. In order to improve the bifurcation analysis of system (5), we employ the Melnikov method. Using the following transformations,(6)x=ε1/3X,y=ε2/3Y,t=ε1/3t1,a=ε2/3ξ1,τ=ε2/3τ1,α=ε4/3ξ2,so system (5) turns into(7)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y.Settingξ1=1 and ε=0, we get the integrable Hamiltonian system as follows:(8)X˙=Y,Y˙=X−bX3.And the corresponding Hamiltonian function isHX,Y=1/2Y2−1/2X2+b/4Y4. As shown in Figure 2, the phase portraits can express the homoclinic trajectories for HX,Y=0.Figure 2
The phase portraits of Hamiltonian function.In order to discuss the value of parameters, we set(9)x0t=±2bsect,y0t=∓2bsecttanht.The Melnikov functions(10)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdtare given to verify the existence of Hamiltonian bifurcations. By calculating, we get the Melnikov function as(11)M±ξ2=4ξ23b−1635b2.From Melnikov function theory, whenM±ξ2=4ξ2/3b−16/35b2=0, that is, ξ2=12/35b, we get the approximate bifurcation of homoclinic orbits. Let ξ1=1, and then according to transformation equation (7), the Melnikov threshold α=12a2/35b. To illustrate the existence of saddle-node bifurcations in limit cycles at this time, we lead into a periodic orbit γe by the level energy of e and period of Te. In the case of e∈−1/4,0, there is a small-amplitude period; if e>0,γe, there becomes a large-amplitude period. Therefore, the Melnikov function for cyclic orbits is(12)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt.Through numerical simulation, it can be seen that the saddle-node bifurcations of the two periodic orbits are very close to the Melnikov threshold. Figure3 shows that when the parameter enters region 2 from region 1 through the saddle junction, the system transits from the equilibrium state of the two small-amplitude limit cycles to a larger amplitude periodic orbit. It also indicates that the Melnikov threshold curve leads to the disappearance of the two limit cycles, resulting in a stable large-amplitude limit cycle. In order to facilitate the expansion discussed below, we only consider the parameters of the region where there are two small-amplitude periods and the unstable saddle point of the ordinary equilibrium point. Therefore, the numerical simulation parameter values given in this paper are fixed at α=a=b=1.Figure 3
Melnikov threshold curve of the unperturbed model and the qualitative branching atb=1.


## 3. Oscillatory Dynamics with Equal Frequencies

This section focuses on the oscillating dynamics when the natural frequency is equal to the excitation frequency. We discuss the case whenf2=0 and f1 is a variable. Then, the case where both f2 and f1 are variables is discussed.


### 3.1. Oscillation Mechanism of MMOs forf2=0

Forf2=0, we set β=cosω2t, and the fast system is driven by(13)x¨+x˙3−αx˙−axt−τ+bx3=f1β.We obtain its perturbed form by scale transformation as follows:(14)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y+f1β.By further calculating the Melnikov function, we have(15)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdt+f1β∫−∞+∞y0tdt.Thus, the Melnikov function of the local periodic orbit can be expressed as(16)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt+f1β∫0Tey0tdt.In equation (16), these basic integrals are constants. When the value of y0t is given by a determined value, we can get the critical threshold of the function. With the change of β, f1β will periodically affect the saddle-junction bifurcation, which leads to the phase trajectory transition from a small-amplitude periodic orbit to a large-amplitude periodic orbit, which further determines the oscillation phenomenon.Given the parameter valuea=b=α=1,ω2=0.1,τ=0.3, and f1=1, Figure 4 shows the trajectory phase diagram and time-history curve of the system. The transition of the system between two small-amplitude periodic orbits can be seen from Figure 4. This excitation oscillation can be expressed as a closed singular orbit formed by two fast and slow motion trajectories.Figure 4
The phase trajectory (a) and time series (b) for the parametersa=b=α=1,ω2=0.1,f1=1,f2=0, and τ=0.3.
(a)
(b)


### 3.2. Oscillation Mechanism of MMOs as Varyingf1 and f2

Whenf1≠0,f2≠0, and ω1=ω2, set β=cosω1t=cosω2t, and we assume(17)x¨+x˙3−αx˙−axt−τ+bx3=f1+f2ββ.Use the scale change to obtain its perturbed form as(18)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y+f1+f2ββ.Then, we get the Melnikov function(19)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdt+f1+f2ββ∫−∞+∞y0tdt.And the Melnikov function of the local periodic orbit can be written as(20)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt+f1+f2ββ∫0Tey0tdt.Once the value ofy0t is determined in equation (20), we can obtain the critical threshold of the function. With the change of β, f1+f2ββ will periodically affect the saddle-junction bifurcation, which leads to the phase trajectory transit from a small-amplitude periodic orbit to a large-amplitude periodic orbit, which determines the oscillation phenomenon.Take the parameter valuea=b=α=1,ω2=ω1=0.01,τ=0.3,f1=1, and f2=0.5. The time-history curve of the system is seen in Figure 5. From the figures, we can observe that the mixed-mode oscillation at this time has two asymmetric small-amplitude periodic orbits and a set of independent large-amplitude periodic orbits in each period.Figure 5
The time series for the parametersa=b=α=1,ω2=ω1=0.01,f2=0.5,f1=1, and τ=0.3.


## 3.1. Oscillation Mechanism of MMOs forf2=0

Forf2=0, we set β=cosω2t, and the fast system is driven by(13)x¨+x˙3−αx˙−axt−τ+bx3=f1β.We obtain its perturbed form by scale transformation as follows:(14)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y+f1β.By further calculating the Melnikov function, we have(15)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdt+f1β∫−∞+∞y0tdt.Thus, the Melnikov function of the local periodic orbit can be expressed as(16)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt+f1β∫0Tey0tdt.In equation (16), these basic integrals are constants. When the value of y0t is given by a determined value, we can get the critical threshold of the function. With the change of β, f1β will periodically affect the saddle-junction bifurcation, which leads to the phase trajectory transition from a small-amplitude periodic orbit to a large-amplitude periodic orbit, which further determines the oscillation phenomenon.Given the parameter valuea=b=α=1,ω2=0.1,τ=0.3, and f1=1, Figure 4 shows the trajectory phase diagram and time-history curve of the system. The transition of the system between two small-amplitude periodic orbits can be seen from Figure 4. This excitation oscillation can be expressed as a closed singular orbit formed by two fast and slow motion trajectories.Figure 4
The phase trajectory (a) and time series (b) for the parametersa=b=α=1,ω2=0.1,f1=1,f2=0, and τ=0.3.
(a)
(b)


## 3.2. Oscillation Mechanism of MMOs as Varyingf1 and f2

Whenf1≠0,f2≠0, and ω1=ω2, set β=cosω1t=cosω2t, and we assume(17)x¨+x˙3−αx˙−axt−τ+bx3=f1+f2ββ.Use the scale change to obtain its perturbed form as(18)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y+f1+f2ββ.Then, we get the Melnikov function(19)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdt+f1+f2ββ∫−∞+∞y0tdt.And the Melnikov function of the local periodic orbit can be written as(20)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt+f1+f2ββ∫0Tey0tdt.Once the value ofy0t is determined in equation (20), we can obtain the critical threshold of the function. With the change of β, f1+f2ββ will periodically affect the saddle-junction bifurcation, which leads to the phase trajectory transit from a small-amplitude periodic orbit to a large-amplitude periodic orbit, which determines the oscillation phenomenon.Take the parameter valuea=b=α=1,ω2=ω1=0.01,τ=0.3,f1=1, and f2=0.5. The time-history curve of the system is seen in Figure 5. From the figures, we can observe that the mixed-mode oscillation at this time has two asymmetric small-amplitude periodic orbits and a set of independent large-amplitude periodic orbits in each period.Figure 5
The time series for the parametersa=b=α=1,ω2=ω1=0.01,f2=0.5,f1=1, and τ=0.3.


## 4. Oscillatory Dynamics with Two Resonant Frequencies

This section mainly discusses the mixed-mode oscillation dynamics reflected by the system when the natural frequency and the resonant frequency are not equal. Without loss of generality, assume that the two frequencies are proportional, andω1 is a positive integer multiple of ω2, i.e., ω1=nω2, n is an integer greater than 1, or ω1 is a noninteger multiple of ω2, i.e., n is an irrational number, and then system (1) can be turned into(21)x¨+x˙3−αx˙−axt−τ+bx3=f1+f2cosnω2cosω2.Then, we will discuss three cases: (1) the dynamics of the system withn being a small positive integer; (2) the oscillation behavior with n being a large positive integer; and (3) the coupling of time-delay systems with n being an irrational number.


### 4.1. Oscillation Mechanism of MMOs for a Relatively Smalln

Sincen is a small real number, then the natural frequency ω1 is still a low frequency. Thus, ω1 and ω2 can be treated as a fast and slow form with a single slow variable. Borrowing De Moivre’s formula, we arrive at(22)cosnx=Cn0cosnx+Cn2cosn−2xisinx2+⋯+Cnmcosn−mxisinxm,where m is the maximum value not greater than n. Let us see some simple examples. When n takes 2 and 3, respectively, we can get cos2x=2cos2x−1 and cos3x=4cos3x−3cosx. Let n=2, the oscillation behavior of the system is shown in Figure 6, and the values of the parameters are a=b=1,α=0.1,ω1=0.02,ω2=0.01, τ=0.5,f1=0.5, and f2=2, respectively.Figure 6
The time history of fast-slow system (1) at ω1=2ω2, where the parameters are a=b=1,α=0.1,ω1=0.02,ω2=0.01, τ=0.3,f1=0.1, and f2=2.From Figure6, the oscillation mode can be clearly divided into different parts, and the upper and lower oscillations are connected by the jump connection. And compared with previous Figure 6, the number of different oscillations is significantly increased in the gentle region and the peak region in each cycle.Now, considering the situation of divergence whenn=2 and setting β=cosω1t and cosω2t=2β2−1A, we have(23)x¨+x˙3−αx˙−axt−τ+bx3=f1β2+f2β3−f2β.Perturbation form by reference to the scale change method is(24)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y+f1β2+f2β3−f2β.Then, we obtain the Melnikov function(25)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdt+f1β2+f2β3−f2β∫−∞+∞y0tdt.The Melnikov function of the local periodic orbit can be written as(26)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt+f1β2+f2β3−f2β∫0Tey0tdt.Ify0t in equation (26) is given, the critical threshold of the function can be determined. With the change of β, f1β2+f2β3−f2β periodically affects the saddle-junction bifurcation, causing the phase trajectory to change from a small-amplitude periodic orbit to a large-amplitude periodic orbit and generating oscillation behavior.


### 4.2. Oscillation Mechanism of MMOs for a Relatively Largen

Whenn is a sufficiently large integer, there is a large step gap between the natural frequency and the excitation frequency. Since there are many oscillatory components in each oscillation mode, the De Moivre’s formula cannot be used for analysis. Now, in this example, we discuss the following.Lettingβ=cosω2t, then system (1) can be described as(27)x¨+x˙−αx˙−axt−τ+bx3=f1+f2cosω1tβ.Accordingly, its perturbation form is(28)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y+f1+f2cosω1tβ.So, the Melnikov function of the system becomes(29)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdt+f1+f2cosω1t∫−∞+∞y0tdt.Thus, we obtain the Melnikov function of the periodic orbit as follows:(30)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt+f1+f2cosω1t∫0Tey0tdt.If we can determine the value ofy0t in equation (30), then we also can get the critical threshold of the function. The value of the natural frequency ω also affects the function M±ξ2,e. In this case, as β changes, the system periodically crosses the saddle-junction of the limit cycle, resulting in more complex oscillatory behavior.Through numerical simulation, we use Figure7 to illustrate the dynamic behavior of the mixed-mode oscillation when there is a gap between the natural frequency and the excitation frequency. Here, we set the parameters of the system a=b=1,α=0.1,ω1=2,ω2=0.01, τ=0.3, and f1=f2=0.5. From Figure 7, we can find the transition between two local limit cycles.Figure 7
The time history of fast-slow system (1) at ω1=200ω2, where the parameters are a=b=1,α=0.1,ω1=2,ω2=0.01, τ=0.3, and f1=f2=0.5.


### 4.3. Oscillation Mechanism of MMOs for an Irrational Numbern

In this section, we consider the coupled Duffing equation with multiple-frequency external forces and delayed feedbacks as follows:(31)x″+δ1x′−yt−τ+x3=β1cosω1t+β3cosω2t,y″+δ2y′−xt−τ+y3=β2cosω1t+β4cosω2t,where δii=1,2δ>0 are the dampings, βii=1,2,3,4 are the excitation amplitudes, and ωii=1,2 are the corresponding excitation frequencies. ττ≥0 is the time delay and ω1=Οε, where ε≪1.


#### 4.3.1. General Method

We study the dynamic behavior of fast-slow system (31) with two slow variables by using the fast-slow analysis method. The two slow variables can be represented by functions containing γt, that is, cosω1t=f1γt and cosω2t=f2γt, and the system can be transformed into a fast-slow system with only one slow variable γt. Therefore, the system can be studied by the traditional fast-slow analysis method.According to Taylor series expansion, we obtain(32)cosnt=fn∗cost,where(33)fn∗x=Cn0xn−Cn2xn−21−x2+Cn4xn−41−x22−⋯+imCnmxn−m1−x2m/2.Setω2=nω1 (n is a positive integer), and then equation (32) becomes(34)x″+δ1x′−yt−τ+x3=β1γt+β3fn∗γt,y″+δ2y′−xt−τ+y3=β2γt+β4fn∗γt,where γt=cosω1t. Furthermore, the fast subsystem can be regarded as(35)x″+δ1x′−yt−τ+x3=β1γ+β3fn∗γ,y″+δ2y′−xt−τ+y3=β2γ+β4fn∗γ.


#### 4.3.2. The Duffing System with Commensurate Excitation Frequencies

We begin our analysis by considering the case whenω1=0.01 and ω2=0.03, name ω2=3ω1. Setting δ1=2, δ2=1, β1=β2=1.2, β3=0.5, and β4=0.555. Figure 8 shows the time-history curve of the system when τ=0, τ=1, and τ=2. As shown in the figure, periodic oscillations occur in these states. Each periodic oscillation is composed of large-amplitude oscillations and small-amplitude oscillations. When appearing small-amplitude oscillation, the system is said to be in a resting state. When the system trajectory passes through the bifurcation point, it loses its equilibrium state and shows a large oscillation. At this point, the system moves from the resting state to the excited state. When the orbit of the system jumps back and forth between the resting state and the excited state, it is called bursting. The two frequencies are engaged in the bursting of the system. It is seen that each bursting pattern can be divided into two parts, i.e., the upper and lower oscillations which are connected by catastrophic jumps.Figure 8
Time series of the bursting forδ1=2, δ2=1, β1=β2=1.2, β3=0.5, β4=0.555, ω1=0.01, ω2=0.03. a1,a2 τ=0, b1,b2 τ=1, and c1,c2 τ=2.
(a)
(b)
(c)To further study the dynamic behavior in Figure8, the equilibrium curve and the transformation phase diagram of fast subsystem (32) were drawn by setting ω1=0.01, ω2=0.03, δ1=2, δ2=1, β1=β2=1.2, β3=0.5 and β4=0.555 (see Figure 8). So, we have the fast subsystem(36)x″+δ1x′−yt−τ+x3=β1γ+β3f3∗γ,y″+δ2y′−xt−τ+y3=β2γ+β4f3∗γ,where γt=cos0.01t is the control parameter. Figure 9 shows the equilibrium curve and the transformation phase diagram τ=0. It leads to an S-shaped equilibrium curve. The solid line represents the stable equilibrium point, and the dotted line represents the unstable equilibrium point; FB indicates the fold bifurcation point, and it can be seen from the figure that there are two stable equilibrium points and two fold bifurcation points in the fast subsystem. With the change of γ, the system trajectory moves to the right along the stable equilibrium and loses the balance after meeting the fold bifurcation point and then jumps to the lower part of the balance curve. The system is excited from a static state. Due to the attraction of the lower balance, the rail line tends to be stable and exits the excited state into the static state until γ attains its maximum value, and then the rail line moves to the left. By this way, two jumps are completed in a cycle.Figure 9
Fast-slow analysis of the burstinga1,b1,a2,b2, c1,d1,c2,d2, and e1,f1,e2,f2 is related to the bursting patterns in Figure 8 a1,a2, b1,b2, and c1,c2.
(a)
(b)
(c)
(d)
(e)
(f)Figures9(c)–9(f) are the equilibrium curve and the transformation phase diagram when τ=1 and τ=2. The solid line represents the stable equilibrium point, while the dotted line represents the unstable equilibrium point; FB indicates the fold bifurcation point. The similar fold bifurcation occurred in Figures 9(a) and 9(b), and almost the same rail line was drawn. Combining the equilibrium curve and the transformation phase diagram, we found that the occurrence of bursting is not significantly disturbed even if the time delay increased, and the dynamic behavior of the system do not lose its typical characteristics with the change of time delay.To further study the influence of time delay on the bursting, Figure10 shows the plot of the peak parts of Figures 8(a)–8(c). We find that the time delay causes the dynamic changes in the upper and lower parts of the oscillation within each period.Figure 10
The plot of the peak part of Figures8(a)–8(c).
(a)
(b)
(c)


#### 4.3.3. The Duffing System with Incommensurate Excitation Frequencies

Now, we consider that the ratio of two frequencies is not the ratio of two integers, i.e., there is at least one of the two excitation frequencies which is irrational. We assume that, in system (31), ω1 is a rational frequency, while ω2 is an irrational frequency, e.g., ω2=π/100. Then, the irrational frequency ω2 leads to a rational sequence Ωn=π/100, in which πn denotes the 10−n-grade truncated π, i.e., an approximation of π. For example, π3=3.141 and π6=3.141592. Set δ1=2, δ2=1, β1=β2=1.2, β3=0.5, β4=0.555, and τ=1. Figure 11 shows the time-history curve of the system when ω1=0.01, ω2=π/100, ω1=0.01, ω2=3/100, ω1=π/100, ω2=0.01, and ω1=3/100, ω2=0.01. The complex bursting behavior in the system with delayed feedback is universal. Compared with Figure 8, the upper and lower vibration of each bursting curve is more intense.Figure 11
Time series of the bursting forδ1=2, δ1=1, β1=β2=1.2, β3=0.5, β4=0.555, and τ=1. a1,a2 ω1=0.01, ω2=π/100 and b1,b2 ω1=0.01, ω2=3/100. c1,c2 ω1=π/100, ω2=0.01 and d1,d2 ω1=3/100, ω2=0.01.
(a)
(b)
(c)
(d)Since the ratio of two excitation frequencies is not the ratio of two integers, we cannot directly use the traditional fast and slow analysis method to analyze the behavior of the system with incommensurate excitation frequencies. So, we try to use the numerical simulation to draw the overlap of the bursting with incommensurate excitation frequencies and the one with truncated and commensurate excitation frequency (see Figure12). It is shown that the yellow curve agrees well with the red curve, i.e., the bursting pattern in Figure 12(a) with 10−2-grade truncated frequency agrees well with the one in Figure 11(a), and the bursting pattern in Figure 12(b) with 10−3-grade truncated frequency agrees well with the one in Figure 11(b). Therefore, we can conclude that the bursting mode related to the incommensurate excitation frequency can be well approximated to the 10−n-grade truncated excitation frequency.Figure 12
Bursting pattern with incommensurate excitation frequencies (red curve) agrees well with the one with truncated, commensurate excitation frequencies (yellow curve).
(a)
(b)Then, the bursting in Figure10(a) can be used to analyze the following system:(37)x″+δ1x′−yt−τ+x3=β1f50∗γ+β3f157∗γ,y″+δ2y′−xt−τ+y3=β2f50∗γ+β4f157∗γ,where γt=cos0.0002t is the control parameter so that we can continue to study the fast subsystem by the fast-slow analysis method.


## 4.1. Oscillation Mechanism of MMOs for a Relatively Smalln

Sincen is a small real number, then the natural frequency ω1 is still a low frequency. Thus, ω1 and ω2 can be treated as a fast and slow form with a single slow variable. Borrowing De Moivre’s formula, we arrive at(22)cosnx=Cn0cosnx+Cn2cosn−2xisinx2+⋯+Cnmcosn−mxisinxm,where m is the maximum value not greater than n. Let us see some simple examples. When n takes 2 and 3, respectively, we can get cos2x=2cos2x−1 and cos3x=4cos3x−3cosx. Let n=2, the oscillation behavior of the system is shown in Figure 6, and the values of the parameters are a=b=1,α=0.1,ω1=0.02,ω2=0.01, τ=0.5,f1=0.5, and f2=2, respectively.Figure 6
The time history of fast-slow system (1) at ω1=2ω2, where the parameters are a=b=1,α=0.1,ω1=0.02,ω2=0.01, τ=0.3,f1=0.1, and f2=2.From Figure6, the oscillation mode can be clearly divided into different parts, and the upper and lower oscillations are connected by the jump connection. And compared with previous Figure 6, the number of different oscillations is significantly increased in the gentle region and the peak region in each cycle.Now, considering the situation of divergence whenn=2 and setting β=cosω1t and cosω2t=2β2−1A, we have(23)x¨+x˙3−αx˙−axt−τ+bx3=f1β2+f2β3−f2β.Perturbation form by reference to the scale change method is(24)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y+f1β2+f2β3−f2β.Then, we obtain the Melnikov function(25)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdt+f1β2+f2β3−f2β∫−∞+∞y0tdt.The Melnikov function of the local periodic orbit can be written as(26)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt+f1β2+f2β3−f2β∫0Tey0tdt.Ify0t in equation (26) is given, the critical threshold of the function can be determined. With the change of β, f1β2+f2β3−f2β periodically affects the saddle-junction bifurcation, causing the phase trajectory to change from a small-amplitude periodic orbit to a large-amplitude periodic orbit and generating oscillation behavior.


## 4.2. Oscillation Mechanism of MMOs for a Relatively Largen

Whenn is a sufficiently large integer, there is a large step gap between the natural frequency and the excitation frequency. Since there are many oscillatory components in each oscillation mode, the De Moivre’s formula cannot be used for analysis. Now, in this example, we discuss the following.Lettingβ=cosω2t, then system (1) can be described as(27)x¨+x˙−αx˙−axt−τ+bx3=f1+f2cosω1tβ.Accordingly, its perturbation form is(28)X˙=Y,Y˙=ξ1X−bX3+εξ2Y−Y3−τξ1Y+f1+f2cosω1tβ.So, the Melnikov function of the system becomes(29)M±ξ2=ξ2∫−∞+∞y02tdt−∫−∞+∞y04tdt+f1+f2cosω1t∫−∞+∞y0tdt.Thus, we obtain the Melnikov function of the periodic orbit as follows:(30)M±ξ2,e=ξ2∫0Tey02tdt−∫0Tey04tdt+f1+f2cosω1t∫0Tey0tdt.If we can determine the value ofy0t in equation (30), then we also can get the critical threshold of the function. The value of the natural frequency ω also affects the function M±ξ2,e. In this case, as β changes, the system periodically crosses the saddle-junction of the limit cycle, resulting in more complex oscillatory behavior.Through numerical simulation, we use Figure7 to illustrate the dynamic behavior of the mixed-mode oscillation when there is a gap between the natural frequency and the excitation frequency. Here, we set the parameters of the system a=b=1,α=0.1,ω1=2,ω2=0.01, τ=0.3, and f1=f2=0.5. From Figure 7, we can find the transition between two local limit cycles.Figure 7
The time history of fast-slow system (1) at ω1=200ω2, where the parameters are a=b=1,α=0.1,ω1=2,ω2=0.01, τ=0.3, and f1=f2=0.5.


## 4.3. Oscillation Mechanism of MMOs for an Irrational Numbern

In this section, we consider the coupled Duffing equation with multiple-frequency external forces and delayed feedbacks as follows:(31)x″+δ1x′−yt−τ+x3=β1cosω1t+β3cosω2t,y″+δ2y′−xt−τ+y3=β2cosω1t+β4cosω2t,where δii=1,2δ>0 are the dampings, βii=1,2,3,4 are the excitation amplitudes, and ωii=1,2 are the corresponding excitation frequencies. ττ≥0 is the time delay and ω1=Οε, where ε≪1.


### 4.3.1. General Method

We study the dynamic behavior of fast-slow system (31) with two slow variables by using the fast-slow analysis method. The two slow variables can be represented by functions containing γt, that is, cosω1t=f1γt and cosω2t=f2γt, and the system can be transformed into a fast-slow system with only one slow variable γt. Therefore, the system can be studied by the traditional fast-slow analysis method.According to Taylor series expansion, we obtain(32)cosnt=fn∗cost,where(33)fn∗x=Cn0xn−Cn2xn−21−x2+Cn4xn−41−x22−⋯+imCnmxn−m1−x2m/2.Setω2=nω1 (n is a positive integer), and then equation (32) becomes(34)x″+δ1x′−yt−τ+x3=β1γt+β3fn∗γt,y″+δ2y′−xt−τ+y3=β2γt+β4fn∗γt,where γt=cosω1t. Furthermore, the fast subsystem can be regarded as(35)x″+δ1x′−yt−τ+x3=β1γ+β3fn∗γ,y″+δ2y′−xt−τ+y3=β2γ+β4fn∗γ.


### 4.3.2. The Duffing System with Commensurate Excitation Frequencies

We begin our analysis by considering the case whenω1=0.01 and ω2=0.03, name ω2=3ω1. Setting δ1=2, δ2=1, β1=β2=1.2, β3=0.5, and β4=0.555. Figure 8 shows the time-history curve of the system when τ=0, τ=1, and τ=2. As shown in the figure, periodic oscillations occur in these states. Each periodic oscillation is composed of large-amplitude oscillations and small-amplitude oscillations. When appearing small-amplitude oscillation, the system is said to be in a resting state. When the system trajectory passes through the bifurcation point, it loses its equilibrium state and shows a large oscillation. At this point, the system moves from the resting state to the excited state. When the orbit of the system jumps back and forth between the resting state and the excited state, it is called bursting. The two frequencies are engaged in the bursting of the system. It is seen that each bursting pattern can be divided into two parts, i.e., the upper and lower oscillations which are connected by catastrophic jumps.Figure 8
Time series of the bursting forδ1=2, δ2=1, β1=β2=1.2, β3=0.5, β4=0.555, ω1=0.01, ω2=0.03. a1,a2 τ=0, b1,b2 τ=1, and c1,c2 τ=2.
(a)
(b)
(c)To further study the dynamic behavior in Figure8, the equilibrium curve and the transformation phase diagram of fast subsystem (32) were drawn by setting ω1=0.01, ω2=0.03, δ1=2, δ2=1, β1=β2=1.2, β3=0.5 and β4=0.555 (see Figure 8). So, we have the fast subsystem(36)x″+δ1x′−yt−τ+x3=β1γ+β3f3∗γ,y″+δ2y′−xt−τ+y3=β2γ+β4f3∗γ,where γt=cos0.01t is the control parameter. Figure 9 shows the equilibrium curve and the transformation phase diagram τ=0. It leads to an S-shaped equilibrium curve. The solid line represents the stable equilibrium point, and the dotted line represents the unstable equilibrium point; FB indicates the fold bifurcation point, and it can be seen from the figure that there are two stable equilibrium points and two fold bifurcation points in the fast subsystem. With the change of γ, the system trajectory moves to the right along the stable equilibrium and loses the balance after meeting the fold bifurcation point and then jumps to the lower part of the balance curve. The system is excited from a static state. Due to the attraction of the lower balance, the rail line tends to be stable and exits the excited state into the static state until γ attains its maximum value, and then the rail line moves to the left. By this way, two jumps are completed in a cycle.Figure 9
Fast-slow analysis of the burstinga1,b1,a2,b2, c1,d1,c2,d2, and e1,f1,e2,f2 is related to the bursting patterns in Figure 8 a1,a2, b1,b2, and c1,c2.
(a)
(b)
(c)
(d)
(e)
(f)Figures9(c)–9(f) are the equilibrium curve and the transformation phase diagram when τ=1 and τ=2. The solid line represents the stable equilibrium point, while the dotted line represents the unstable equilibrium point; FB indicates the fold bifurcation point. The similar fold bifurcation occurred in Figures 9(a) and 9(b), and almost the same rail line was drawn. Combining the equilibrium curve and the transformation phase diagram, we found that the occurrence of bursting is not significantly disturbed even if the time delay increased, and the dynamic behavior of the system do not lose its typical characteristics with the change of time delay.To further study the influence of time delay on the bursting, Figure10 shows the plot of the peak parts of Figures 8(a)–8(c). We find that the time delay causes the dynamic changes in the upper and lower parts of the oscillation within each period.Figure 10
The plot of the peak part of Figures8(a)–8(c).
(a)
(b)
(c)


### 4.3.3. The Duffing System with Incommensurate Excitation Frequencies

Now, we consider that the ratio of two frequencies is not the ratio of two integers, i.e., there is at least one of the two excitation frequencies which is irrational. We assume that, in system (31), ω1 is a rational frequency, while ω2 is an irrational frequency, e.g., ω2=π/100. Then, the irrational frequency ω2 leads to a rational sequence Ωn=π/100, in which πn denotes the 10−n-grade truncated π, i.e., an approximation of π. For example, π3=3.141 and π6=3.141592. Set δ1=2, δ2=1, β1=β2=1.2, β3=0.5, β4=0.555, and τ=1. Figure 11 shows the time-history curve of the system when ω1=0.01, ω2=π/100, ω1=0.01, ω2=3/100, ω1=π/100, ω2=0.01, and ω1=3/100, ω2=0.01. The complex bursting behavior in the system with delayed feedback is universal. Compared with Figure 8, the upper and lower vibration of each bursting curve is more intense.Figure 11
Time series of the bursting forδ1=2, δ1=1, β1=β2=1.2, β3=0.5, β4=0.555, and τ=1. a1,a2 ω1=0.01, ω2=π/100 and b1,b2 ω1=0.01, ω2=3/100. c1,c2 ω1=π/100, ω2=0.01 and d1,d2 ω1=3/100, ω2=0.01.
(a)
(b)
(c)
(d)Since the ratio of two excitation frequencies is not the ratio of two integers, we cannot directly use the traditional fast and slow analysis method to analyze the behavior of the system with incommensurate excitation frequencies. So, we try to use the numerical simulation to draw the overlap of the bursting with incommensurate excitation frequencies and the one with truncated and commensurate excitation frequency (see Figure12). It is shown that the yellow curve agrees well with the red curve, i.e., the bursting pattern in Figure 12(a) with 10−2-grade truncated frequency agrees well with the one in Figure 11(a), and the bursting pattern in Figure 12(b) with 10−3-grade truncated frequency agrees well with the one in Figure 11(b). Therefore, we can conclude that the bursting mode related to the incommensurate excitation frequency can be well approximated to the 10−n-grade truncated excitation frequency.Figure 12
Bursting pattern with incommensurate excitation frequencies (red curve) agrees well with the one with truncated, commensurate excitation frequencies (yellow curve).
(a)
(b)Then, the bursting in Figure10(a) can be used to analyze the following system:(37)x″+δ1x′−yt−τ+x3=β1f50∗γ+β3f157∗γ,y″+δ2y′−xt−τ+y3=β2f50∗γ+β4f157∗γ,where γt=cos0.0002t is the control parameter so that we can continue to study the fast subsystem by the fast-slow analysis method.


## 4.3.1. General Method

We study the dynamic behavior of fast-slow system (31) with two slow variables by using the fast-slow analysis method. The two slow variables can be represented by functions containing γt, that is, cosω1t=f1γt and cosω2t=f2γt, and the system can be transformed into a fast-slow system with only one slow variable γt. Therefore, the system can be studied by the traditional fast-slow analysis method.According to Taylor series expansion, we obtain(32)cosnt=fn∗cost,where(33)fn∗x=Cn0xn−Cn2xn−21−x2+Cn4xn−41−x22−⋯+imCnmxn−m1−x2m/2.Setω2=nω1 (n is a positive integer), and then equation (32) becomes(34)x″+δ1x′−yt−τ+x3=β1γt+β3fn∗γt,y″+δ2y′−xt−τ+y3=β2γt+β4fn∗γt,where γt=cosω1t. Furthermore, the fast subsystem can be regarded as(35)x″+δ1x′−yt−τ+x3=β1γ+β3fn∗γ,y″+δ2y′−xt−τ+y3=β2γ+β4fn∗γ.


## 4.3.2. The Duffing System with Commensurate Excitation Frequencies

We begin our analysis by considering the case whenω1=0.01 and ω2=0.03, name ω2=3ω1. Setting δ1=2, δ2=1, β1=β2=1.2, β3=0.5, and β4=0.555. Figure 8 shows the time-history curve of the system when τ=0, τ=1, and τ=2. As shown in the figure, periodic oscillations occur in these states. Each periodic oscillation is composed of large-amplitude oscillations and small-amplitude oscillations. When appearing small-amplitude oscillation, the system is said to be in a resting state. When the system trajectory passes through the bifurcation point, it loses its equilibrium state and shows a large oscillation. At this point, the system moves from the resting state to the excited state. When the orbit of the system jumps back and forth between the resting state and the excited state, it is called bursting. The two frequencies are engaged in the bursting of the system. It is seen that each bursting pattern can be divided into two parts, i.e., the upper and lower oscillations which are connected by catastrophic jumps.Figure 8
Time series of the bursting forδ1=2, δ2=1, β1=β2=1.2, β3=0.5, β4=0.555, ω1=0.01, ω2=0.03. a1,a2 τ=0, b1,b2 τ=1, and c1,c2 τ=2.
(a)
(b)
(c)To further study the dynamic behavior in Figure8, the equilibrium curve and the transformation phase diagram of fast subsystem (32) were drawn by setting ω1=0.01, ω2=0.03, δ1=2, δ2=1, β1=β2=1.2, β3=0.5 and β4=0.555 (see Figure 8). So, we have the fast subsystem(36)x″+δ1x′−yt−τ+x3=β1γ+β3f3∗γ,y″+δ2y′−xt−τ+y3=β2γ+β4f3∗γ,where γt=cos0.01t is the control parameter. Figure 9 shows the equilibrium curve and the transformation phase diagram τ=0. It leads to an S-shaped equilibrium curve. The solid line represents the stable equilibrium point, and the dotted line represents the unstable equilibrium point; FB indicates the fold bifurcation point, and it can be seen from the figure that there are two stable equilibrium points and two fold bifurcation points in the fast subsystem. With the change of γ, the system trajectory moves to the right along the stable equilibrium and loses the balance after meeting the fold bifurcation point and then jumps to the lower part of the balance curve. The system is excited from a static state. Due to the attraction of the lower balance, the rail line tends to be stable and exits the excited state into the static state until γ attains its maximum value, and then the rail line moves to the left. By this way, two jumps are completed in a cycle.Figure 9
Fast-slow analysis of the burstinga1,b1,a2,b2, c1,d1,c2,d2, and e1,f1,e2,f2 is related to the bursting patterns in Figure 8 a1,a2, b1,b2, and c1,c2.
(a)
(b)
(c)
(d)
(e)
(f)Figures9(c)–9(f) are the equilibrium curve and the transformation phase diagram when τ=1 and τ=2. The solid line represents the stable equilibrium point, while the dotted line represents the unstable equilibrium point; FB indicates the fold bifurcation point. The similar fold bifurcation occurred in Figures 9(a) and 9(b), and almost the same rail line was drawn. Combining the equilibrium curve and the transformation phase diagram, we found that the occurrence of bursting is not significantly disturbed even if the time delay increased, and the dynamic behavior of the system do not lose its typical characteristics with the change of time delay.To further study the influence of time delay on the bursting, Figure10 shows the plot of the peak parts of Figures 8(a)–8(c). We find that the time delay causes the dynamic changes in the upper and lower parts of the oscillation within each period.Figure 10
The plot of the peak part of Figures8(a)–8(c).
(a)
(b)
(c)


## 4.3.3. The Duffing System with Incommensurate Excitation Frequencies

Now, we consider that the ratio of two frequencies is not the ratio of two integers, i.e., there is at least one of the two excitation frequencies which is irrational. We assume that, in system (31), ω1 is a rational frequency, while ω2 is an irrational frequency, e.g., ω2=π/100. Then, the irrational frequency ω2 leads to a rational sequence Ωn=π/100, in which πn denotes the 10−n-grade truncated π, i.e., an approximation of π. For example, π3=3.141 and π6=3.141592. Set δ1=2, δ2=1, β1=β2=1.2, β3=0.5, β4=0.555, and τ=1. Figure 11 shows the time-history curve of the system when ω1=0.01, ω2=π/100, ω1=0.01, ω2=3/100, ω1=π/100, ω2=0.01, and ω1=3/100, ω2=0.01. The complex bursting behavior in the system with delayed feedback is universal. Compared with Figure 8, the upper and lower vibration of each bursting curve is more intense.Figure 11
Time series of the bursting forδ1=2, δ1=1, β1=β2=1.2, β3=0.5, β4=0.555, and τ=1. a1,a2 ω1=0.01, ω2=π/100 and b1,b2 ω1=0.01, ω2=3/100. c1,c2 ω1=π/100, ω2=0.01 and d1,d2 ω1=3/100, ω2=0.01.
(a)
(b)
(c)
(d)Since the ratio of two excitation frequencies is not the ratio of two integers, we cannot directly use the traditional fast and slow analysis method to analyze the behavior of the system with incommensurate excitation frequencies. So, we try to use the numerical simulation to draw the overlap of the bursting with incommensurate excitation frequencies and the one with truncated and commensurate excitation frequency (see Figure12). It is shown that the yellow curve agrees well with the red curve, i.e., the bursting pattern in Figure 12(a) with 10−2-grade truncated frequency agrees well with the one in Figure 11(a), and the bursting pattern in Figure 12(b) with 10−3-grade truncated frequency agrees well with the one in Figure 11(b). Therefore, we can conclude that the bursting mode related to the incommensurate excitation frequency can be well approximated to the 10−n-grade truncated excitation frequency.Figure 12
Bursting pattern with incommensurate excitation frequencies (red curve) agrees well with the one with truncated, commensurate excitation frequencies (yellow curve).
(a)
(b)Then, the bursting in Figure10(a) can be used to analyze the following system:(37)x″+δ1x′−yt−τ+x3=β1f50∗γ+β3f157∗γ,y″+δ2y′−xt−τ+y3=β2f50∗γ+β4f157∗γ,where γt=cos0.0002t is the control parameter so that we can continue to study the fast subsystem by the fast-slow analysis method.


## 5. The Effect of Time Delay on Cluster Discovery

Time delay is an important parameter of time-delay systems, which affects the dynamic bifurcation and stability of the system. This section mainly discusses the influence of time lag on the occurrence and development of cluster discovery before and after approximation and then confirms the effective time lag range by numerical simulation to reduce the error size and obtain the desired peak dynamics. Here, we mainly consider the system(38)x¨+x˙−αx˙−axt−τ+bx3=f1cosω2t.Setxt−τ≈xt−τx′t in the numerical simulation of the two cases before and after the approximation of the time lag which is performed.Now, we take the parametersa=b=α=1,ω2=0.01,f1=1, τ=1,τ=0.6, τ=0.5, and τ=0.3. Figures 13–16 compare the time-history curves of the system before and after approximation. From the comparison of these figures, the time lag does not affect the generation of cluster images, respectively. By comparing the preapproximation and the approximation of the system delay, it is found that the approximation of the system in Figures 13 and 14 causes the dynamics of the upper and lower portions of the oscillation to change in each cycle. Figures 15 and 16 show the approximation of the system is almost consistent with the dynamic behavior before the approximation. Therefore, the effective time range of the system should be 0∼0.3. In the above sections, we generally take τ=0.3.Figure 13
Time history of system (36). The red line represents the curve before time-delay approximation, and the blue line represents the curve after time-delay approximation when the parameters are a=b=α=1,ω2=0.01,f1=1, and τ=1.Figure 14
Time history of system (36). The red line represents the curve before time-delay approximation, and the blue line represents the curve after time-delay approximation when the parameters are a=b=α=1,ω2=0.01,f1=1, and τ=0.6.Figure 15
Time history of system (36). The red line represents the curve before time-delay approximation, and the blue line represents the curve after time-delay approximation when the parameters are a=b=α=1,ω2=0.01,f1=1, and τ=0.5.Figure 16
Time history of system (36). The red line represents the curve before time-delay approximation, and the blue line represents the curve after time-delay approximation when the parameters are a=b=α=1,ω2=0.01,f1=1, and τ=0.3.


## 6. The Multistate Dynamic Response of the System

In dynamic systems, multistability is the property that there are multiple stable equilibrium points in the vector space spanned by the state of the system. According to the mathematical, there must also be points of instability between the stable points. Near the unstable equilibrium point, any system will be sensitive to noise, initial conditions, and system parameters, which will lead to the development of the system in many different directions. Here, we mainly consider the system(39)x˙=y,y˙=−a1x+a2y+y5−fcosΩt,where xt and yt are the real function. f>0 is the amplitude, and ω is the frequencies. a1 and a2 are the physical parameters.Use the Euler method to obtain the following discrete systems:(40)xn+1=yn,yn+1=ayn−bxn−cyn5+β,where a=2−a2Δt,b=a1Δt2+a2Δt+1,c=Δt5, and β=fcosΩt. Next, we will study the multistate dynamic response of the above discrete system. Set a and b as the real parameter and β as the control parameter. Controllable slow periodic motion is expressed as Zn=fcosΩt, in which the external excitation frequency is taken as 0.001. For convenience of calculation, we set c=1.Now, we takeβ as the bifurcation parameter to explore the coexistence of chaotic attractors, periodic attractors, and chaotic attractors or the occurrence of numbers and sizes when parameters a and b take different values.Case 1.
Dynamic response without jumping phenomenon.
According to Figures17 and 18, when a=1.2 and β=0.24, the system is in a single periodic motion. When −0.6<β<−0.5 and 0.5<β<0.6, the system is in double periodic motion. And there is no jump occurring at β=0.Figure 17
Bifurcation diagram ofβ.Figure 18
Time series diagram.Case 2.
Whena=1.2 and b=0.23, dynamic occurs in the jumping phenomenon at β=0.
From Figures19 and 20, we know, when −0.53<β<−0.43 and 0.43<β<0.53, the system is in double periodic motion. When −0.43<β<0.43, the system is in a single periodic motion. But, there appear jumping phenomenon when β is near 0.Figure 19
Bifurcation diagram ofβ.Figure 20
Time series diagram.Case 3.
The transition of the system from the bistable state to the chaotic attractor to the monoperiodic attractor.
As shown in Figure21, when βc±=±0.91, the system is in a state of bistable, and chaotic attractor and haploid cycle attractor coexist. When β>0.91 and β<−0.91, the system enters into haploid periodic motion, and the bistable state disappears. And we also found that as β increases or decreases, namely, more than the critical value, the system of the chaotic attractor suddenly disappeared.Figure 21
Bifurcation diagram ofβ.Case 4.
The transition of the system from the bistable state to the chaotic attractor to the double periodic attractor.
Whena=1.9 and b=0.17, from Figure 22, we obtain βc±=±0.042, and the system is in a state of bistable, and chaotic attractor and double period attractor coexist. When β>0.042 and β<−0.042, the system enters into double periodic motion, and the bistable state disappears. And as β increases or decreases, namely, more than the critical value, the system of the chaotic attractor disappeared.Figure 22
Bifurcation diagram ofβ.Case 5.
The transition of the system from the bistable state to the chaotic attractor to the quadruple periodic attractor.
Whena=1.9 and b=0.23, from Figure 23, we obtain βc±=±0.026, and the system is in a state of bistable, and chaotic attractor and haploid cycle attractor coexist. When β>0.026 and β<−0.026, the system enters into quadruple periodic motion, and the bistable state disappears. And as β increases or decreases, namely, more than the critical value, the system of the chaotic attractor disappeared.Figure 23
Bifurcation diagram ofβ.Case 6.
The transition of the system from the bistable state to the chaotic attractor to the chaotic attractor.
From Figure24, we obtain βc±=±0.026, and when a=1.9 and b=0.14, the system is in a state of bistable, and chaotic attractor and haploid cycle attractor coexist. When β>0.026 and β<−0.026, the system enters into quadruple periodic motion, and the bistable state disappears. And as β increases or decreases, namely, more than the critical value, the system of the chaotic attractor disappeared.Figure 24
Bifurcation diagram ofβ.


## 7. Conclusion

In this paper, we study the mixed-mode dynamics of a class of oscillators with time-delay modulated amplitude. Combine theoretical analysis and numerical simulation to make time-history graphs and phase diagrams to explain when the natural frequency and excitation frequency are the same or proportional. We use the Melnikov method to explore the influence of threshold parameters on the bifurcation. We also discuss the coupled time-delay Duffing system to illustrate the dynamic behavior when the frequency ratio is irrational. Meanwhile, we analyze the effect of time-delay approximation before and after the delay approximation on the system clustering phenomenon. Finally, we explore the bifurcation behavior of the fast subsystem. The numerical simulation results show that, near the critical valueβc, the chaotic attractor will coexist with multiple periodic orbits or chaos. By studying the bifurcation of stability and multistability of these systems with parameter excitation, these results have certain practicability and interest for the mathematical modeling of the systems studied in the fields of physics, chemistry, and mechanics. We can also idealize the required systems in more fields.



---

*Source: 4871068-2020-02-06.xml*