# Threshold for Chaos of a Duffing Oscillator with Fractional-Order Derivative

Authors: Wuce Xing; Enli Chen; Yujian Chang; Meiqi Wang
Journal: Shock and Vibration
 (2019)

---

## Abstract

In this paper, the necessary condition for the chaotic motion of a Duffing oscillator with the fractional-order derivative under harmonic excitation is investigated. The necessary condition for the chaos in the sense of Smale horseshoes is established based on the Melnikov method, and then the chaotic threshold curve is obtained. The largest Lyapunov exponents are provided, and some other typical numerical simulation results, including the time histories, frequency spectrograms, phase portraits, and Poincare maps, are presented and compared. From the analysis of the numerical simulation results, it could be found that, near the chaotic threshold curve, the system generates chaos via the period-doubling bifurcation, from single periodic motion to period-2 motion and period-4 motion to chaotic motion. The effects of fractional-order parameters, the stiffness coefficient, and the damping coefficient on the threshold value of the chaotic motion are analytically discussed. The results show that the coefficient of the fractional-order derivative has great effect on the threshold value of the chaotic motion, while the order of the fractional-order derivative has less. The analysis results reveal some new phenomena, and it could be useful for designing or controlling dynamic systems with the fractional-order derivative.


---

## Body

## 1. Introduction

The development of fractional-order calculus has a history of more than 300 years, but until recently, the research in this field still focused on pure mathematical theory. In recent years, the research and application of fractional calculus in the mechanical engineering field has gradually increased [1], especially to accurately characterize the true constitutive relationship of viscoelastic materials [2, 3] and controlling engineering [4, 5].In many nonlinear dynamical systems, such as van der Pol, Mathieu, Hill’s oscillators, Duffing oscillators, and other systems, large number of results can be found in the literature [6–10]. The complex dynamic phenomena of typical nonlinear systems, such as periodic solutions [11, 12], bifurcation, and chaos [13], have been more and more concerned. Among these nonlinear dynamical phenomena, chaos is a unique dynamic behavior that involves nonperiodic motion originating from deterministic systems and depending on initial conditions. Its effect on mechanical systems cannot be ignored. Therefore, the study of chaotic threshold is of great significance to the understanding of the mechanism of chaos generation, which can lead to the control and anticontrol of chaos [14, 15]. A chaotic threshold is the boundary between the periodic motion and chaotic motion of the nonlinear dynamical system. A Duffing oscillator with the fractional-order derivative is a typical nonlinear dynamic system with complex dynamic behaviors. At present, the research on fractional-order Duffing oscillators mainly focuses on the analysis of dynamic response characteristics. Shen et al. [16] and Mobin and Saeed [17] studied Duffing oscillators with fractional order and obtained the amplitude-frequency response curve of the systems by using the average method. They also analyzed the effect of fractional-order parameters on the primary resonance amplitude value and frequency. Wen et al. [18, 19] derived the approximate analytical solution of the fractional-order Mathieu–Duffing equation based on the incremental harmonic balance method. Xu et al. [20] considered a complex Duffing system subjected to nonstationary random excitation. Chen et al. [21–23] studied the random jumping and bifurcation behaviors of fractional-order Duffing oscillators under combined harmonic and white noise excitation, and they designed a controller based on the fractional-order PID.There has been some research on chaotic motions and synchronization control of chaos to the fractional-order Duffing system. Liu et al. [24, 25] proposed a numerical simulation method called extended generalized unit mapping for a nonautonomous fractional-order Duffing system, and they studied the chaotic turbulence inside the system. Zufiria and Jiménez [26] described the chaotic behavior of fractional-order Duffing oscillators through the Lyapunov exponents and strange attractors in the effective phase space of the systems. Mohammad et al. [27] suggested a technique for finding unstable periodic orbits in chaotic fractional-order systems, using a fractional-order Duffing system as an example. Some scholars [28–30] analyzed the chaotic synchronization behavior of a modified van der Pol-Duffing system and revealed the effect of fractional order on chaotic synchronization controls. These examples show that applications of the fractional-order derivative in engineering fields are becoming increasingly important.The Melnikov method has certain advantages when it is used to analyze the chaotic motion threshold value of multiparameter nonlinear systems. For getting threshold values of the chaotic motion, many scholars have used the analytic Melnikov method to study the necessary condition for chaos of Duffing oscillators without a fractional-order derivative [31–35]. However, due to the complexity of the fractional-order derivative, it is rare to discuss its effect on the chaotic threshold of nonlinear fractional-order systems. Anague Tabejieu et al. [36] studied the Rayleigh beams with fractional order, but there was no relevant literature on the threshold value of the chaotic motion of Duffing oscillators with a fractional-order derivative.In this paper, the necessary condition for generating chaos in sense of Smale horseshoes in a typical Duffing oscillator with the fractional-order derivative under harmonic excitation is investigated based on the Melnikov method. The paper is organized as follows. In Section2, the Melnikov function is established based on Melnikov theory, and the necessary conditions for chaos are obtained. The fractional-order derivative is calculated by the numerical method, and the analytical expressions are obtained in the rest part. In Section 3, the necessary condition is verified by the maximum Lyapunov exponents, and the typical dynamical responses, including the time histories, frequency spectrograms, phase portraits, and Poincare maps of the system, are also investigated. The results show that this system has chaotic motions. The result obtained based on the numerical solution satisfies the chaos threshold curve, which fully verifies the correctness of the necessary conditions. In Section 4, the effects of fractional-order parameters, stiffness coefficient, and the damping coefficient of the system are analyzed. Finally, the detailed conclusions are summarized and analyzed.


## 2. Necessary Condition for Chaos of Fractional-Order Duffing Oscillators

Duffing oscillators are one of the most typical and important objects in nonlinear dynamics. In this paper, the ordinary Duffing oscillator with a fractional-order derivative under harmonic excitation is described with(1)mx¨−k1x+k3x3+c1x˙+hDtpxt=Fcosωt,where m, k1, k3, c1, F, and ω are the system mass, linear stiffness coefficient, nonlinear stiffness coefficient, linear viscous damping coefficient, excitation amplitude, and excitation frequency, respectively. The expression hDtpxt is the pth-order derivative from x(t) to t for the fractional coefficient h(h > 0), and the fractional order is restricted to 0 < p < 1. All the system parameters are positive and dimensionless.Using the transformations(2)b1=k1m,B=εb=c1m,b2=k3m,μ1=εμ=hm,F0=εf=Fm,where ε is a small real parameter, 0<ε≪1, equation (1) becomes(3)x¨−b1x+b2x3=εfcosωt−μDtpx−bx˙.


### 2.1. Homoclinic Orbits of Fractional-Order Duffing Oscillators

Equation (1) can be transformed into the form of a state equation:(4)x˙=y,y˙=b1x−b2x3+εfcosωt−μDtpx−by.Ifε = 0, meaning that neither system damping nor external excitation exists, then equation (3) can be transformed into(5)x˙=y=f1,y˙=b1x−b2x3=f2.Now the system is a Hamilton system whose Hamilton function satisfiesf1 = ∂H/∂y and f2 = −∂H/∂x.In our case,(6)Hx,y=12y2−b12x2+b24x4.Lettingx˙=0 and y˙=0 leads to three fixed points: 0,0, b1/b2,0, and −b1/b2,0.According to equation (5), the characteristic equation of the system is(7)0−λ1b1−3b2x20−λ=0.Using equation (7), one can obtain(8)λ=±b1−3b2x2.When the fixed point is (0, 0), the eigenvalue of the system isλ=±b1. When the fixed point is b1/b2,0 or −b1/b2,0, the eigenvalue of the system is λ=±2b1i. According to the properties of a fixed point [37], fixed point (0, 0) is the saddle point of the system and fixed points b1/b2,0 and −b1/b2,0 are the center of the system. Using equation (6), one can obtain the homoclinic orbit, and the result is shown in Figure 1. When the Hamilton function H(x, y) = 0, there are two homoclinic orbits connecting saddle points. The positions of a saddle point and two centers are indicated in Figure 1. The saddle point is unstable, and the centers are stable.Figure 1
Homoclinic orbits of the system (b1 = 1, b2 = 1).The homoclinic orbits of the system described in equation (1) can be obtained by solving the following expression:(9)H0,0=12y2−b12x2+b24x4=0,x˙=y.Equation (8) leads to(10)y=±b1x2−b22x4,dxdt=y.The variables in equation (9) can be separated and the integral done, and this leads to the parametric equations of the two homoclinic orbits:(11)xt=±2b1b2sechb1t,yt=∓b12b2sechb1ttanhb1t.An analysis of equation (10) leads to the conclusion that a plus sign of x(t) represents the positive axis of a homoclinic orbit and a minus sign represents the negative part. The plus sign of y(t) represents the upper half of a homoclinic orbit, and minus sign represents the lower half.


### 2.2. Melnikov Function Determines the Chaotic System’s Necessary Condition

Ifx→=xx˙, where fx→=yb1x−b2x3 and gx→,t=0εfcosωt−μDtpx−by, equation (3) can be transformed into(12)x→˙=fx→+εgx→,t.Melnikov theory [38], proposed by Melnikov, is a perturbation method originally for proving the existence of transverse homoclinic or heteroclinic orbits. Melnikov theory has been applied in many researches. The Melnikov theory was firstly applied to study a periodically driven Duffing oscillator by Holmes [39]. According to the Melnikov process, one can determine the Melnikov function with y(t) as the odd function:(13)Mt0=∫−∞∞fq0τ−t0∧gq0τ−t0,τdτ=ε∫−∞∞yτ−t0fcosωτ−μDtpx−byτ−t0dτ=t=τ−t0ε∫−∞∞ytfcosωt+t0−μDt+t0px−bytdt=εf∫−∞∞ytcosωtcosωt0−sinωtsinωt0dt−εμ∫−∞∞ytDt+t0pxdt−εb∫−∞∞y2tdt=−εfsinωt0∫−∞∞ytsinωtdt−εμ∫−∞∞ytDt+t0pxdt−εb∫−∞∞y2tdt.SubstitutingI1=∫−∞∞ytsinωtdt, I2=∫−∞∞ytDt+t0pxdt, and I3=∫−∞∞y2tdt into M(t0) yields(14)Mt0=−F0sinωt0I1−μ1I2−BI3.According to Melnikov theory [37], this leads to(15)Mt0=−F0I1sinωt0−μ1I2−BI3=0.Sincesinωt0≤1, equation (15) can be used to the get necessary condition:(16)1F≤1mI1μ1I2+BI3.The integral in equation (11) can be divided into three parts: I1, I2, and I3. These three parts were calculated in the following.


#### 2.2.1. CalculatingI1 and I3

We substituted the parametric equation (11) of the two homoclinic orbits into equation (14). Since I1 and I3 do not contain fractional-order derivatives, the residue theorem and contour integration can be used to calculate I1 and I3 [38, 40]:(17)I1=∓b12b2∫−∞∞sechb1ttanhb1tsinωtdt=±2b2πωsechπω2b1,(18)I3=2b12b2∫−∞∞sech2b1ttanh2b1tdt=4b123b2b1.


#### 2.2.2. CalculatingI2

SinceI2 in equation (14) contains a fractional-order derivative, it is impossible to integrate directly, so it can only be calculated using the definition of a fractional-order derivative. The Caputo definition and the Riemann–Liouville definition are as follows:(19)​t0CDtpxt=1Γ1-p∫t0tx′τt-τpdτ,​t0RLDtpxt=1Γ1-pddt∫t0tfτt-τpdτ,where Γ(x) is the gamma function satisfying Γ(x + 1) = xΓ(x).If the initial value ofx(t) is nonzero, it is relatively easy to use the Caputo definition [41] to calculate the fractional-order derivative. Thus, we used Caputo’s definition to calculate I2. Comparing the two definitions of the fractional-order derivative, the relationship between them can be obtained as follows:(20)​t0CDtpxt=​t0RLDtpxt−ft=​t0RLDtpxt−∑k=0m−1xktΓk−p+1t−t0k−p,where m = ΓpΊ and f(x) is a penalty function.The specific process of calculating the high-accuracy Caputo fractional-order derivative is below. First, the improved direct recursive algorithm [42] is used to calculate the Riemann–Liouville fractional-order derivative with high accuracy. Then, the penalty function f(t) is calculated, and finally, the penalty function is combined with the Riemann–Liouville fractional-order derivative, leading to the high-accuracy Caputo fractional-order derivative. In this way, the integral with the fractional-order derivative can be calculated, and the value of I2 can be obtained:(21)I2=∫−∞∞ytDt+t0pxdt=∓b12b2∫−∞∞sechb1ttanhb1tDt+t0p±2b1b2sechb1t+b1t0dt=2b1b1b2∫−∞∞sechb1ttanhb1tDt+t0psechb1t+b1t0dt=4b1b1b2H,where(22)H=∫0∞sechb1ttanhb1tDt+t0psechb1t+b1t0dt.


#### 2.2.3. Discussion of the Parametert0

Equation (20) contains the parameter t0, which is the initial moment when the stable manifold intersects the unstable manifold, that is to say, when the system generates chaos in the sense of the Smale horseshoe. The variable t0 can be regarded as a variable about time, t0∈0,T, where T is the period of the system motion. We analyzed the effect of t0 on the chaos threshold value (the system parameters were m = 1, k1 = –1, c1 = 0.5, k3 = 1, and h = 0.5). Figure 2 shows the results of this analysis.Figure 2
Effect oft0 on the chaos threshold.Figure2 shows that an increase of t0 would cause the threshold extremum of the chaotic motion to decrease, reducing the chaotic area and making it more difficult to generate the chaos, but the parameter t0 has little effect on the threshold extremum of chaotic motion (the chaos thresholds of point A and point B are 2.807 and 2.775, respectively). Therefore, we selected t0 = 1 when analyzing the chaotic threshold of system.


#### 2.2.4. Necessary Condition and Threshold Curves of a Chaotic System

Substituting the result of equations (17), (18), and (21) into equation (16), the necessary condition for generating chaos in the sense of Smale horseshoes can be obtained as follows:(23)1F≤1m2/b2πωsechπω/2b14μ1b1b1/b2H+4Bb12/3b2b1.Based on the relationship between the excitation amplitudeF and the excitation frequency ω obtained from equation (23), the necessary condition of the system for generating chaos can be determined. The chaotic threshold curve is shown in Figure 3 (the system parameters are m = 1, k1 = 1, c1 = 0.5, k3 = 1, h = 0.5, and p=0.5).Figure 3
Chaotic threshold curve.Figure3 shows that when the value of the parameter 1/F is located in the region below the curve (gray region), the system might generate chaos or not. Equation (23) shows a condition for chaos to occur, but it is a sufficient but not a necessary condition. When the value of the parameter 1/F is located in the region above the curve (white region), the system will have a stable periodic motion, which may be a single-cycle motion or a multicycle motion.From the analysis of Figure3, it could be concluded that whether chaotic motion occurs is affected by the force and frequency of external excitation. The following is an analysis of these two aspects:(1)
When the external excitation force of the system is small, the system will not generate chaotic motions, and the system will have regular periodic motions. When the external excitation force increases and exceeds the critical value for chaos, the chaotic motion may occur in the system, and the larger the excitation force, the higher the possibility of the chaotic motion.(2)
The chaotic threshold value first increases with the increasing external excitation frequency, then decreases, and finally approaches zero. At the resonance frequency, the system chaotic threshold value reaches a maximum. In the resonance region of the system, the vibration amplitude is larger, and chaos is more likely to occur.Therefore, in the practical engineering of the systems, we should try to avoid the resonance region of the system in order to avoid the irregular chaos phenomenon.


## 2.1. Homoclinic Orbits of Fractional-Order Duffing Oscillators

Equation (1) can be transformed into the form of a state equation:(4)x˙=y,y˙=b1x−b2x3+εfcosωt−μDtpx−by.Ifε = 0, meaning that neither system damping nor external excitation exists, then equation (3) can be transformed into(5)x˙=y=f1,y˙=b1x−b2x3=f2.Now the system is a Hamilton system whose Hamilton function satisfiesf1 = ∂H/∂y and f2 = −∂H/∂x.In our case,(6)Hx,y=12y2−b12x2+b24x4.Lettingx˙=0 and y˙=0 leads to three fixed points: 0,0, b1/b2,0, and −b1/b2,0.According to equation (5), the characteristic equation of the system is(7)0−λ1b1−3b2x20−λ=0.Using equation (7), one can obtain(8)λ=±b1−3b2x2.When the fixed point is (0, 0), the eigenvalue of the system isλ=±b1. When the fixed point is b1/b2,0 or −b1/b2,0, the eigenvalue of the system is λ=±2b1i. According to the properties of a fixed point [37], fixed point (0, 0) is the saddle point of the system and fixed points b1/b2,0 and −b1/b2,0 are the center of the system. Using equation (6), one can obtain the homoclinic orbit, and the result is shown in Figure 1. When the Hamilton function H(x, y) = 0, there are two homoclinic orbits connecting saddle points. The positions of a saddle point and two centers are indicated in Figure 1. The saddle point is unstable, and the centers are stable.Figure 1
Homoclinic orbits of the system (b1 = 1, b2 = 1).The homoclinic orbits of the system described in equation (1) can be obtained by solving the following expression:(9)H0,0=12y2−b12x2+b24x4=0,x˙=y.Equation (8) leads to(10)y=±b1x2−b22x4,dxdt=y.The variables in equation (9) can be separated and the integral done, and this leads to the parametric equations of the two homoclinic orbits:(11)xt=±2b1b2sechb1t,yt=∓b12b2sechb1ttanhb1t.An analysis of equation (10) leads to the conclusion that a plus sign of x(t) represents the positive axis of a homoclinic orbit and a minus sign represents the negative part. The plus sign of y(t) represents the upper half of a homoclinic orbit, and minus sign represents the lower half.


## 2.2. Melnikov Function Determines the Chaotic System’s Necessary Condition

Ifx→=xx˙, where fx→=yb1x−b2x3 and gx→,t=0εfcosωt−μDtpx−by, equation (3) can be transformed into(12)x→˙=fx→+εgx→,t.Melnikov theory [38], proposed by Melnikov, is a perturbation method originally for proving the existence of transverse homoclinic or heteroclinic orbits. Melnikov theory has been applied in many researches. The Melnikov theory was firstly applied to study a periodically driven Duffing oscillator by Holmes [39]. According to the Melnikov process, one can determine the Melnikov function with y(t) as the odd function:(13)Mt0=∫−∞∞fq0τ−t0∧gq0τ−t0,τdτ=ε∫−∞∞yτ−t0fcosωτ−μDtpx−byτ−t0dτ=t=τ−t0ε∫−∞∞ytfcosωt+t0−μDt+t0px−bytdt=εf∫−∞∞ytcosωtcosωt0−sinωtsinωt0dt−εμ∫−∞∞ytDt+t0pxdt−εb∫−∞∞y2tdt=−εfsinωt0∫−∞∞ytsinωtdt−εμ∫−∞∞ytDt+t0pxdt−εb∫−∞∞y2tdt.SubstitutingI1=∫−∞∞ytsinωtdt, I2=∫−∞∞ytDt+t0pxdt, and I3=∫−∞∞y2tdt into M(t0) yields(14)Mt0=−F0sinωt0I1−μ1I2−BI3.According to Melnikov theory [37], this leads to(15)Mt0=−F0I1sinωt0−μ1I2−BI3=0.Sincesinωt0≤1, equation (15) can be used to the get necessary condition:(16)1F≤1mI1μ1I2+BI3.The integral in equation (11) can be divided into three parts: I1, I2, and I3. These three parts were calculated in the following.


### 2.2.1. CalculatingI1 and I3

We substituted the parametric equation (11) of the two homoclinic orbits into equation (14). Since I1 and I3 do not contain fractional-order derivatives, the residue theorem and contour integration can be used to calculate I1 and I3 [38, 40]:(17)I1=∓b12b2∫−∞∞sechb1ttanhb1tsinωtdt=±2b2πωsechπω2b1,(18)I3=2b12b2∫−∞∞sech2b1ttanh2b1tdt=4b123b2b1.


### 2.2.2. CalculatingI2

SinceI2 in equation (14) contains a fractional-order derivative, it is impossible to integrate directly, so it can only be calculated using the definition of a fractional-order derivative. The Caputo definition and the Riemann–Liouville definition are as follows:(19)​t0CDtpxt=1Γ1-p∫t0tx′τt-τpdτ,​t0RLDtpxt=1Γ1-pddt∫t0tfτt-τpdτ,where Γ(x) is the gamma function satisfying Γ(x + 1) = xΓ(x).If the initial value ofx(t) is nonzero, it is relatively easy to use the Caputo definition [41] to calculate the fractional-order derivative. Thus, we used Caputo’s definition to calculate I2. Comparing the two definitions of the fractional-order derivative, the relationship between them can be obtained as follows:(20)​t0CDtpxt=​t0RLDtpxt−ft=​t0RLDtpxt−∑k=0m−1xktΓk−p+1t−t0k−p,where m = ΓpΊ and f(x) is a penalty function.The specific process of calculating the high-accuracy Caputo fractional-order derivative is below. First, the improved direct recursive algorithm [42] is used to calculate the Riemann–Liouville fractional-order derivative with high accuracy. Then, the penalty function f(t) is calculated, and finally, the penalty function is combined with the Riemann–Liouville fractional-order derivative, leading to the high-accuracy Caputo fractional-order derivative. In this way, the integral with the fractional-order derivative can be calculated, and the value of I2 can be obtained:(21)I2=∫−∞∞ytDt+t0pxdt=∓b12b2∫−∞∞sechb1ttanhb1tDt+t0p±2b1b2sechb1t+b1t0dt=2b1b1b2∫−∞∞sechb1ttanhb1tDt+t0psechb1t+b1t0dt=4b1b1b2H,where(22)H=∫0∞sechb1ttanhb1tDt+t0psechb1t+b1t0dt.


### 2.2.3. Discussion of the Parametert0

Equation (20) contains the parameter t0, which is the initial moment when the stable manifold intersects the unstable manifold, that is to say, when the system generates chaos in the sense of the Smale horseshoe. The variable t0 can be regarded as a variable about time, t0∈0,T, where T is the period of the system motion. We analyzed the effect of t0 on the chaos threshold value (the system parameters were m = 1, k1 = –1, c1 = 0.5, k3 = 1, and h = 0.5). Figure 2 shows the results of this analysis.Figure 2
Effect oft0 on the chaos threshold.Figure2 shows that an increase of t0 would cause the threshold extremum of the chaotic motion to decrease, reducing the chaotic area and making it more difficult to generate the chaos, but the parameter t0 has little effect on the threshold extremum of chaotic motion (the chaos thresholds of point A and point B are 2.807 and 2.775, respectively). Therefore, we selected t0 = 1 when analyzing the chaotic threshold of system.


### 2.2.4. Necessary Condition and Threshold Curves of a Chaotic System

Substituting the result of equations (17), (18), and (21) into equation (16), the necessary condition for generating chaos in the sense of Smale horseshoes can be obtained as follows:(23)1F≤1m2/b2πωsechπω/2b14μ1b1b1/b2H+4Bb12/3b2b1.Based on the relationship between the excitation amplitudeF and the excitation frequency ω obtained from equation (23), the necessary condition of the system for generating chaos can be determined. The chaotic threshold curve is shown in Figure 3 (the system parameters are m = 1, k1 = 1, c1 = 0.5, k3 = 1, h = 0.5, and p=0.5).Figure 3
Chaotic threshold curve.Figure3 shows that when the value of the parameter 1/F is located in the region below the curve (gray region), the system might generate chaos or not. Equation (23) shows a condition for chaos to occur, but it is a sufficient but not a necessary condition. When the value of the parameter 1/F is located in the region above the curve (white region), the system will have a stable periodic motion, which may be a single-cycle motion or a multicycle motion.From the analysis of Figure3, it could be concluded that whether chaotic motion occurs is affected by the force and frequency of external excitation. The following is an analysis of these two aspects:(1)
When the external excitation force of the system is small, the system will not generate chaotic motions, and the system will have regular periodic motions. When the external excitation force increases and exceeds the critical value for chaos, the chaotic motion may occur in the system, and the larger the excitation force, the higher the possibility of the chaotic motion.(2)
The chaotic threshold value first increases with the increasing external excitation frequency, then decreases, and finally approaches zero. At the resonance frequency, the system chaotic threshold value reaches a maximum. In the resonance region of the system, the vibration amplitude is larger, and chaos is more likely to occur.Therefore, in the practical engineering of the systems, we should try to avoid the resonance region of the system in order to avoid the irregular chaos phenomenon.


## 2.2.1. CalculatingI1 and I3

We substituted the parametric equation (11) of the two homoclinic orbits into equation (14). Since I1 and I3 do not contain fractional-order derivatives, the residue theorem and contour integration can be used to calculate I1 and I3 [38, 40]:(17)I1=∓b12b2∫−∞∞sechb1ttanhb1tsinωtdt=±2b2πωsechπω2b1,(18)I3=2b12b2∫−∞∞sech2b1ttanh2b1tdt=4b123b2b1.


## 2.2.2. CalculatingI2

SinceI2 in equation (14) contains a fractional-order derivative, it is impossible to integrate directly, so it can only be calculated using the definition of a fractional-order derivative. The Caputo definition and the Riemann–Liouville definition are as follows:(19)​t0CDtpxt=1Γ1-p∫t0tx′τt-τpdτ,​t0RLDtpxt=1Γ1-pddt∫t0tfτt-τpdτ,where Γ(x) is the gamma function satisfying Γ(x + 1) = xΓ(x).If the initial value ofx(t) is nonzero, it is relatively easy to use the Caputo definition [41] to calculate the fractional-order derivative. Thus, we used Caputo’s definition to calculate I2. Comparing the two definitions of the fractional-order derivative, the relationship between them can be obtained as follows:(20)​t0CDtpxt=​t0RLDtpxt−ft=​t0RLDtpxt−∑k=0m−1xktΓk−p+1t−t0k−p,where m = ΓpΊ and f(x) is a penalty function.The specific process of calculating the high-accuracy Caputo fractional-order derivative is below. First, the improved direct recursive algorithm [42] is used to calculate the Riemann–Liouville fractional-order derivative with high accuracy. Then, the penalty function f(t) is calculated, and finally, the penalty function is combined with the Riemann–Liouville fractional-order derivative, leading to the high-accuracy Caputo fractional-order derivative. In this way, the integral with the fractional-order derivative can be calculated, and the value of I2 can be obtained:(21)I2=∫−∞∞ytDt+t0pxdt=∓b12b2∫−∞∞sechb1ttanhb1tDt+t0p±2b1b2sechb1t+b1t0dt=2b1b1b2∫−∞∞sechb1ttanhb1tDt+t0psechb1t+b1t0dt=4b1b1b2H,where(22)H=∫0∞sechb1ttanhb1tDt+t0psechb1t+b1t0dt.


## 2.2.3. Discussion of the Parametert0

Equation (20) contains the parameter t0, which is the initial moment when the stable manifold intersects the unstable manifold, that is to say, when the system generates chaos in the sense of the Smale horseshoe. The variable t0 can be regarded as a variable about time, t0∈0,T, where T is the period of the system motion. We analyzed the effect of t0 on the chaos threshold value (the system parameters were m = 1, k1 = –1, c1 = 0.5, k3 = 1, and h = 0.5). Figure 2 shows the results of this analysis.Figure 2
Effect oft0 on the chaos threshold.Figure2 shows that an increase of t0 would cause the threshold extremum of the chaotic motion to decrease, reducing the chaotic area and making it more difficult to generate the chaos, but the parameter t0 has little effect on the threshold extremum of chaotic motion (the chaos thresholds of point A and point B are 2.807 and 2.775, respectively). Therefore, we selected t0 = 1 when analyzing the chaotic threshold of system.


## 2.2.4. Necessary Condition and Threshold Curves of a Chaotic System

Substituting the result of equations (17), (18), and (21) into equation (16), the necessary condition for generating chaos in the sense of Smale horseshoes can be obtained as follows:(23)1F≤1m2/b2πωsechπω/2b14μ1b1b1/b2H+4Bb12/3b2b1.Based on the relationship between the excitation amplitudeF and the excitation frequency ω obtained from equation (23), the necessary condition of the system for generating chaos can be determined. The chaotic threshold curve is shown in Figure 3 (the system parameters are m = 1, k1 = 1, c1 = 0.5, k3 = 1, h = 0.5, and p=0.5).Figure 3
Chaotic threshold curve.Figure3 shows that when the value of the parameter 1/F is located in the region below the curve (gray region), the system might generate chaos or not. Equation (23) shows a condition for chaos to occur, but it is a sufficient but not a necessary condition. When the value of the parameter 1/F is located in the region above the curve (white region), the system will have a stable periodic motion, which may be a single-cycle motion or a multicycle motion.From the analysis of Figure3, it could be concluded that whether chaotic motion occurs is affected by the force and frequency of external excitation. The following is an analysis of these two aspects:(1)
When the external excitation force of the system is small, the system will not generate chaotic motions, and the system will have regular periodic motions. When the external excitation force increases and exceeds the critical value for chaos, the chaotic motion may occur in the system, and the larger the excitation force, the higher the possibility of the chaotic motion.(2)
The chaotic threshold value first increases with the increasing external excitation frequency, then decreases, and finally approaches zero. At the resonance frequency, the system chaotic threshold value reaches a maximum. In the resonance region of the system, the vibration amplitude is larger, and chaos is more likely to occur.Therefore, in the practical engineering of the systems, we should try to avoid the resonance region of the system in order to avoid the irregular chaos phenomenon.


## 3. Numerical Simulation Analysis of a Fractional-Order Duffing System

In order to verify the validity of the necessary condition for generating chaos, a set of illustrative system parameters is chosen asm = 1, k1 = 1, k3 = 1, c1 = 0.5, h = 0.5, and p=0.5. The chaotic threshold curve of the external force F with respect to the excitation frequency ω is obtained. Above the curve is the region of the periodic motion, and below the curve is the region where chaos may occur. Different points in the two regions are selected, and the location of each point is shown in Figure 4. When F = 1, F = 1.5, and F = 1.75, the largest Lyapunov exponents are shown in Figures 5–7, respectively. The numerical simulation method for the largest Lyapunov exponents about the differential equation can be found in [43]. From the analysis of Figures 4–7, it can be found that the maximum Lyapunov exponents of points in the chaotic region is greater than zero, while the maximum Lyapunov exponents of points in the periodic motion region is less than zero. The results imply the chaotic threshold curve is correct.Figure 4
Chaotic threshold curve and the location of each point.Figure 5
The largest Lyapunov exponents (F = 1 N).Figure 6
The largest Lyapunov exponents (F = 1.5 N).Figure 7
The largest Lyapunov exponents (F = 1.75 N).Some typical points are shown in Figure4 with solid dots. The points from A to E are in the region of the periodic motion, and points F, G, and H are in the region of the chaotic motion. Among them, points A and G would occur when the excitation amplitude F = 1, points H, E, D, and B would occur when the excitation amplitude F = 1.5, while points F and C would occur when F = 1.75. The horizontal abscissa of the intersection points of the three straight lines and chaotic edge curves are 1.45, 1.86, and 2, respectively, and the Lyapunov exponents are equal to zero at these frequencies. In order to further explain Figure 4, the largest Lyapunov exponents and the system responses at each point are analyzed, respectively, in the following text.


### 3.1. The Largest Lyapunov Exponents

The Lyapunov exponents (LEs) measure the exponential growth, or decay, of infinitesimal phase space perturbations of a chaotic dynamical system. We use the Benettin–Wolf algorithm to determine the Lyapunov exponent. The main steps to determine numerically the Lyapunov exponent are as follows: using the numerical method of power series expansion (PSE) to numerically integrate the fractional-order differential equations to obtain the time series ofx(t) [41], the Gram–Schmidt procedure and picking up the exponents during the renormalization procedure, and the Lyapunov exponents being determined as the average of the logarithm of the stretching factor of each perturbation. The numerical simulation method for the largest Lyapunov exponents about the Duffing oscillator with the fractional-order derivative can be found in [44]. The largest Lyapunov exponents at F = 1 N, F = 1.5 N, and F = 1.75 N are shown in Figures 5–7.Figure5 shows that the maximum Lyapunov exponents of point G is greater than zero, the maximum Lyapunov exponents of point A is less than zero, and the critical excitation frequency for chaos is 1.51. Accordingly, the system performs periodic motions at point A and chaotic motions at point G.Figure6 shows that the maximum Lyapunov exponents of point H is greater than zero, the maximum Lyapunov exponents of points B, D, and E are less than zero, and the critical excitation frequency for chaos is 1.92. Accordingly, the system performs periodic motion at points B, D, and E and chaotic motion at point H.Figure7 shows that the maximum Lyapunov exponents of point F is greater than zero, the maximum Lyapunov exponents of point C is less than zero, and the critical excitation frequency for chaos is 2.06. Accordingly, the system performs periodic motions at point C and chaotic motions at point F.Compared with Figures4–7, it could be found that the critical excitation frequency based on the Melnikov method is little less than the numerical simulation method for the largest Lyapunov exponents (the difference value of critical excitation frequencies between the two methods is 0.06, 0.06, and 0.06, respectively). Although there is difference between the two results quantitatively, this difference is reasonable because the Melnikov method is a first-order approximate way to find chaotic motions [32]. This could also be useful to design or revise the similar dynamical system. From the observation of Figures 5–7, it could be found that increase of the critical excitation frequency would be against the chaos generation. In addition, it could be found that the maximum Lyapunov exponents in the chaotic region may appear as the negative value (when the system is in the state of periodic motion) from the region surrounded by the black circles 1, 2, and 3 in Figures 5–7, which show that the chaotic threshold curve determined based on Melnikov theory is the necessary condition for the chaotic motion of the system.


### 3.2. Analysis of the System Responses at Each Point

In order to further verify the conclusion, the typical dynamical responses, including the time history, frequency spectrograms, phase portraits, and Poincare maps, at different typical points are presented and analyzed in the following text. During the numerical simulation process [19, 41], the adopted numerical scheme is(24)Dpxtl≈h−p∑j=0lCjpxtl−j,where tl = lh is the time sample points, h is the sample step, and Cjp is the fractional binomial coefficient with the iterative relationship as(25)c0p=1,cjp=1−1+pjcj−1p.Here, we selecth = 0.01, the total computation time is 500 excitation periods, and 500 points were sampled in each period. The frontal transient response is omitted, and the posterior 100 periods are selected for analysis. In the following section, the points (in Figure 4) in the region of the periodic motion and chaotic motion are simulated and analyzed. Some typical samples in different ranges are taken to analyze the dynamical phenomena. Through the analysis of these dynamical characteristics at different typical points, the conclusion is that the curve obtained in Figure 4 is the boundary between chaos and periodic motion.


#### 3.2.1. Region of Periodic Motions

The graphs in Figure8 illustrate the periodic motion when ω = 2 rad/s and F = 1 N (point A).Figure 8
System response atω = 2 rad/s and F = 1 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure9 illustrate the periodic motion when ω = 3 rad/s and F = 1.5 N (point B).Figure 9
System response atω = 3 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure10 illustrate the periodic motion when ω = 2.5 rad/s and F = 1.75 N (point C).Figure 10
System response atω = 2.5 rad/s and F = 1.75 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figures8–10 show that, for different external excitation frequencies and excitation amplitudes, the phase portraits have a stable phase trajectory, and the Poincare maps have only one fixed point. The dominant frequency on the frequency spectrograms indicates that the system has forced vibration when the frequency is equal to the excitation frequency. Therefore, the system performs the steady-state periodic motion.The graphs in Figure11 illustrate the period-2 motion when ω = 2.2 rad/s and F = 1.5 N (point D).Figure 11
System response atω = 2.2 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figure11 shows that the Poincare maps have two fixed points, and there are dominant frequencies, a half frequency division, and weak frequency multiplication in the frequency spectrogram, showing that there was period-2 motion in this case.The graphs in Figure12 illustrate the period-4 motion when ω = 1.95 rad/s and F = 1.5 N (point E).Figure 12
System response atω = 1.95 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figure12 shows that the Poincare maps have four fixed points, and there are dominant frequencies, quarter frequency division, half frequency division, and weak frequency multiplication in the frequency spectrogram. The frequency multiplication will reduce the smoothness of the time-domain curve of the system, showing that there is period-4 motion in this case.


#### 3.2.2. Region of Chaotic Motions

The graphs in Figure13 illustrate the chaotic motion when ω = 1.85 rad/s and F = 1.75 N (point F).Figure 13
System response atω = 1.85 rad/s and F = 1.75 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure14 illustrate the chaotic motion when ω = 1 rad/s and F = 1 N (point G).Figure 14
System response atω = 1 rad/s and F = 1 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure15 illustrate the chaotic motion when ω = 1.5 rad/s and F = 1.5 N (point H).Figure 15
System response atω = 1.5 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)In Figures13 and 14, the time histories show that the motion of the system is irregular. In the frequency spectrograms, except for a dominant frequency that is equal to the excitation frequency, there are continuous frequency components. Phase trajectories show disorder. The Poincare maps show a collection of discrete points. These all show that the system has the chaotic motion in these cases. In Figure 15, the time-history response curve of the system shows the paroxysmal motion state with a random alternation between the regular and irregular motion. This state is one of the paths through which the system moves from the steady-state periodic motion to chaotic motion. With the change of system parameters, the time of the irregular motion becomes longer and longer. Finally, the system completely enters the chaotic motion state.The system responses for different typical points are shown in Tables1–3 when F = 1 N, F = 1.5 N, and F = 1.75 N.Table 1
System responses for different typical points (F = 1 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
A
2
−0.0022
Period-1
G
1
0.0080
ChaosTable 2
System responses for different typical points (F = 1.5 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
B
3
−0.0287
Period-1
D
2.2
−0.0033
Period-2
E
1.95
−0.0004
Period-4
H
1.5
0.0277
ChaosTable 3
System responses for different typical points (F = 1.75 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
C
2.5
−0.0023
Period-1
F
1.85
0.0067
ChaosBased on the above analytical results, we concluded that when the parameters of the system are above the chaotic threshold curve, the system shows the stable periodic motion, and when they are under the chaotic threshold curve, the system may have chaotic motions. From the dynamic response diagram of the system at the pointsB, D, E, and H, it can be seen that as the points gradually approach the chaotic threshold curve, the system changes from single periodic motion to period-2 motion and period-4 motion to chaotic motion. When it enters the chaotic motion region, it has reciprocating aperiodic motion, so there are abundant nonlinear dynamic phenomena in the region near the chaotic threshold curve. Accordingly, the numerical simulation results verify the correctness of the chaotic threshold curve determined by the previous calculations, and the results further illustrate the validity of the necessary condition obtained by the Melnikov method.


## 3.1. The Largest Lyapunov Exponents

The Lyapunov exponents (LEs) measure the exponential growth, or decay, of infinitesimal phase space perturbations of a chaotic dynamical system. We use the Benettin–Wolf algorithm to determine the Lyapunov exponent. The main steps to determine numerically the Lyapunov exponent are as follows: using the numerical method of power series expansion (PSE) to numerically integrate the fractional-order differential equations to obtain the time series ofx(t) [41], the Gram–Schmidt procedure and picking up the exponents during the renormalization procedure, and the Lyapunov exponents being determined as the average of the logarithm of the stretching factor of each perturbation. The numerical simulation method for the largest Lyapunov exponents about the Duffing oscillator with the fractional-order derivative can be found in [44]. The largest Lyapunov exponents at F = 1 N, F = 1.5 N, and F = 1.75 N are shown in Figures 5–7.Figure5 shows that the maximum Lyapunov exponents of point G is greater than zero, the maximum Lyapunov exponents of point A is less than zero, and the critical excitation frequency for chaos is 1.51. Accordingly, the system performs periodic motions at point A and chaotic motions at point G.Figure6 shows that the maximum Lyapunov exponents of point H is greater than zero, the maximum Lyapunov exponents of points B, D, and E are less than zero, and the critical excitation frequency for chaos is 1.92. Accordingly, the system performs periodic motion at points B, D, and E and chaotic motion at point H.Figure7 shows that the maximum Lyapunov exponents of point F is greater than zero, the maximum Lyapunov exponents of point C is less than zero, and the critical excitation frequency for chaos is 2.06. Accordingly, the system performs periodic motions at point C and chaotic motions at point F.Compared with Figures4–7, it could be found that the critical excitation frequency based on the Melnikov method is little less than the numerical simulation method for the largest Lyapunov exponents (the difference value of critical excitation frequencies between the two methods is 0.06, 0.06, and 0.06, respectively). Although there is difference between the two results quantitatively, this difference is reasonable because the Melnikov method is a first-order approximate way to find chaotic motions [32]. This could also be useful to design or revise the similar dynamical system. From the observation of Figures 5–7, it could be found that increase of the critical excitation frequency would be against the chaos generation. In addition, it could be found that the maximum Lyapunov exponents in the chaotic region may appear as the negative value (when the system is in the state of periodic motion) from the region surrounded by the black circles 1, 2, and 3 in Figures 5–7, which show that the chaotic threshold curve determined based on Melnikov theory is the necessary condition for the chaotic motion of the system.


## 3.2. Analysis of the System Responses at Each Point

In order to further verify the conclusion, the typical dynamical responses, including the time history, frequency spectrograms, phase portraits, and Poincare maps, at different typical points are presented and analyzed in the following text. During the numerical simulation process [19, 41], the adopted numerical scheme is(24)Dpxtl≈h−p∑j=0lCjpxtl−j,where tl = lh is the time sample points, h is the sample step, and Cjp is the fractional binomial coefficient with the iterative relationship as(25)c0p=1,cjp=1−1+pjcj−1p.Here, we selecth = 0.01, the total computation time is 500 excitation periods, and 500 points were sampled in each period. The frontal transient response is omitted, and the posterior 100 periods are selected for analysis. In the following section, the points (in Figure 4) in the region of the periodic motion and chaotic motion are simulated and analyzed. Some typical samples in different ranges are taken to analyze the dynamical phenomena. Through the analysis of these dynamical characteristics at different typical points, the conclusion is that the curve obtained in Figure 4 is the boundary between chaos and periodic motion.


### 3.2.1. Region of Periodic Motions

The graphs in Figure8 illustrate the periodic motion when ω = 2 rad/s and F = 1 N (point A).Figure 8
System response atω = 2 rad/s and F = 1 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure9 illustrate the periodic motion when ω = 3 rad/s and F = 1.5 N (point B).Figure 9
System response atω = 3 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure10 illustrate the periodic motion when ω = 2.5 rad/s and F = 1.75 N (point C).Figure 10
System response atω = 2.5 rad/s and F = 1.75 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figures8–10 show that, for different external excitation frequencies and excitation amplitudes, the phase portraits have a stable phase trajectory, and the Poincare maps have only one fixed point. The dominant frequency on the frequency spectrograms indicates that the system has forced vibration when the frequency is equal to the excitation frequency. Therefore, the system performs the steady-state periodic motion.The graphs in Figure11 illustrate the period-2 motion when ω = 2.2 rad/s and F = 1.5 N (point D).Figure 11
System response atω = 2.2 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figure11 shows that the Poincare maps have two fixed points, and there are dominant frequencies, a half frequency division, and weak frequency multiplication in the frequency spectrogram, showing that there was period-2 motion in this case.The graphs in Figure12 illustrate the period-4 motion when ω = 1.95 rad/s and F = 1.5 N (point E).Figure 12
System response atω = 1.95 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figure12 shows that the Poincare maps have four fixed points, and there are dominant frequencies, quarter frequency division, half frequency division, and weak frequency multiplication in the frequency spectrogram. The frequency multiplication will reduce the smoothness of the time-domain curve of the system, showing that there is period-4 motion in this case.


### 3.2.2. Region of Chaotic Motions

The graphs in Figure13 illustrate the chaotic motion when ω = 1.85 rad/s and F = 1.75 N (point F).Figure 13
System response atω = 1.85 rad/s and F = 1.75 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure14 illustrate the chaotic motion when ω = 1 rad/s and F = 1 N (point G).Figure 14
System response atω = 1 rad/s and F = 1 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure15 illustrate the chaotic motion when ω = 1.5 rad/s and F = 1.5 N (point H).Figure 15
System response atω = 1.5 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)In Figures13 and 14, the time histories show that the motion of the system is irregular. In the frequency spectrograms, except for a dominant frequency that is equal to the excitation frequency, there are continuous frequency components. Phase trajectories show disorder. The Poincare maps show a collection of discrete points. These all show that the system has the chaotic motion in these cases. In Figure 15, the time-history response curve of the system shows the paroxysmal motion state with a random alternation between the regular and irregular motion. This state is one of the paths through which the system moves from the steady-state periodic motion to chaotic motion. With the change of system parameters, the time of the irregular motion becomes longer and longer. Finally, the system completely enters the chaotic motion state.The system responses for different typical points are shown in Tables1–3 when F = 1 N, F = 1.5 N, and F = 1.75 N.Table 1
System responses for different typical points (F = 1 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
A
2
−0.0022
Period-1
G
1
0.0080
ChaosTable 2
System responses for different typical points (F = 1.5 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
B
3
−0.0287
Period-1
D
2.2
−0.0033
Period-2
E
1.95
−0.0004
Period-4
H
1.5
0.0277
ChaosTable 3
System responses for different typical points (F = 1.75 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
C
2.5
−0.0023
Period-1
F
1.85
0.0067
ChaosBased on the above analytical results, we concluded that when the parameters of the system are above the chaotic threshold curve, the system shows the stable periodic motion, and when they are under the chaotic threshold curve, the system may have chaotic motions. From the dynamic response diagram of the system at the pointsB, D, E, and H, it can be seen that as the points gradually approach the chaotic threshold curve, the system changes from single periodic motion to period-2 motion and period-4 motion to chaotic motion. When it enters the chaotic motion region, it has reciprocating aperiodic motion, so there are abundant nonlinear dynamic phenomena in the region near the chaotic threshold curve. Accordingly, the numerical simulation results verify the correctness of the chaotic threshold curve determined by the previous calculations, and the results further illustrate the validity of the necessary condition obtained by the Melnikov method.


## 3.2.1. Region of Periodic Motions

The graphs in Figure8 illustrate the periodic motion when ω = 2 rad/s and F = 1 N (point A).Figure 8
System response atω = 2 rad/s and F = 1 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure9 illustrate the periodic motion when ω = 3 rad/s and F = 1.5 N (point B).Figure 9
System response atω = 3 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure10 illustrate the periodic motion when ω = 2.5 rad/s and F = 1.75 N (point C).Figure 10
System response atω = 2.5 rad/s and F = 1.75 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figures8–10 show that, for different external excitation frequencies and excitation amplitudes, the phase portraits have a stable phase trajectory, and the Poincare maps have only one fixed point. The dominant frequency on the frequency spectrograms indicates that the system has forced vibration when the frequency is equal to the excitation frequency. Therefore, the system performs the steady-state periodic motion.The graphs in Figure11 illustrate the period-2 motion when ω = 2.2 rad/s and F = 1.5 N (point D).Figure 11
System response atω = 2.2 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figure11 shows that the Poincare maps have two fixed points, and there are dominant frequencies, a half frequency division, and weak frequency multiplication in the frequency spectrogram, showing that there was period-2 motion in this case.The graphs in Figure12 illustrate the period-4 motion when ω = 1.95 rad/s and F = 1.5 N (point E).Figure 12
System response atω = 1.95 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)Figure12 shows that the Poincare maps have four fixed points, and there are dominant frequencies, quarter frequency division, half frequency division, and weak frequency multiplication in the frequency spectrogram. The frequency multiplication will reduce the smoothness of the time-domain curve of the system, showing that there is period-4 motion in this case.


## 3.2.2. Region of Chaotic Motions

The graphs in Figure13 illustrate the chaotic motion when ω = 1.85 rad/s and F = 1.75 N (point F).Figure 13
System response atω = 1.85 rad/s and F = 1.75 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure14 illustrate the chaotic motion when ω = 1 rad/s and F = 1 N (point G).Figure 14
System response atω = 1 rad/s and F = 1 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)The graphs in Figure15 illustrate the chaotic motion when ω = 1.5 rad/s and F = 1.5 N (point H).Figure 15
System response atω = 1.5 rad/s and F = 1.5 N: (a) time history; (b) frequency spectrogram; (c) phase portrait; (d) Poincare maps.
(a)
(b)
(c)
(d)In Figures13 and 14, the time histories show that the motion of the system is irregular. In the frequency spectrograms, except for a dominant frequency that is equal to the excitation frequency, there are continuous frequency components. Phase trajectories show disorder. The Poincare maps show a collection of discrete points. These all show that the system has the chaotic motion in these cases. In Figure 15, the time-history response curve of the system shows the paroxysmal motion state with a random alternation between the regular and irregular motion. This state is one of the paths through which the system moves from the steady-state periodic motion to chaotic motion. With the change of system parameters, the time of the irregular motion becomes longer and longer. Finally, the system completely enters the chaotic motion state.The system responses for different typical points are shown in Tables1–3 when F = 1 N, F = 1.5 N, and F = 1.75 N.Table 1
System responses for different typical points (F = 1 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
A
2
−0.0022
Period-1
G
1
0.0080
ChaosTable 2
System responses for different typical points (F = 1.5 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
B
3
−0.0287
Period-1
D
2.2
−0.0033
Period-2
E
1.95
−0.0004
Period-4
H
1.5
0.0277
ChaosTable 3
System responses for different typical points (F = 1.75 N).
Points
Excitation frequency (rad/s)
The largest Lyapunov exponents
Motion state
C
2.5
−0.0023
Period-1
F
1.85
0.0067
ChaosBased on the above analytical results, we concluded that when the parameters of the system are above the chaotic threshold curve, the system shows the stable periodic motion, and when they are under the chaotic threshold curve, the system may have chaotic motions. From the dynamic response diagram of the system at the pointsB, D, E, and H, it can be seen that as the points gradually approach the chaotic threshold curve, the system changes from single periodic motion to period-2 motion and period-4 motion to chaotic motion. When it enters the chaotic motion region, it has reciprocating aperiodic motion, so there are abundant nonlinear dynamic phenomena in the region near the chaotic threshold curve. Accordingly, the numerical simulation results verify the correctness of the chaotic threshold curve determined by the previous calculations, and the results further illustrate the validity of the necessary condition obtained by the Melnikov method.


## 4. Effect of System Parameters on the Chaotic Threshold Curve

The main purpose of this section is to describe the effects of four parametersk1, c1, h, and p of the fractional-order Duffing system on the chaotic threshold curve. Based on the necessary condition determined by equation (21), the chaotic threshold curve can be obtained. By changing the value of each parameter, the changes for the chaotic threshold curves can be analyzed. The results are presented individually in Figures 16–19. The natural frequency of the system discussed below is ω0=k/m=1.Figure 16
Chaotic threshold curve with differentk1 values (m = 1, c1 = 0.5, k3 = 1, h = 0.5, and p=0.5).Figure 17
Chaotic threshold curve with differentp values (m = 1, k1 = 1, c1 = 0.5, k3 = 1, and h = 0.5).Figure 18
Chaotic threshold curve with differentc1 values (m = 1, k1 = 1, k3 = 1, h = 0.5, and p=0.5).Figure 19
Chaotic threshold curve with differenth values (m = 1, k1 = 1, c1 = 0.5, k3 = 1, and p=0.5).Figure16 shows that increasing linear stiffness coefficient k1 leads to smaller threshold values for generating chaos, so the generation of chaos becomes more difficult. In addition, when the external excitation frequency is greater than a certain value, the attenuation speed of the chaotic threshold curve decreases with the increasing system stiffness, and it eventually approaches to zero. When the external excitation frequency is relatively large, the chaotic motion of the system is not likely to occur.Figure18 shows that when the linear damping coefficient c1 increases, the peak value of the chaotic threshold curve decreases, the area where chaos occurs decreases, and the possibility of chaos in the system decreases.Figures16 and 18 show that, with increasing stiffness and damping of the system, the peak value of the chaotic threshold curve decreases, the area where chaos may occur decreases, and the possibility of chaos in the system decreases. Compared with the damping of the system, the stiffness of the system has a greater impact on the chaotic region of the system. Therefore, in practical work, the chaotic phenomena can be avoided by adjusting the stiffness of the system.Figure19 shows that increasing coefficient h of the fractional-order derivative leads to decreasing the peak value of the chaotic threshold curve. The area where chaos occurs decreases, and the possibility of chaos in the system decreases.Figure17 shows that increasing the order p of the fractional-order derivative leads to increasing threshold values for generating chaos, increasing areas where chaos occurs, and increasing likelihood of chaos generation.An analysis of Figures16 and 17 shows that the parameters of the fractional-order derivative have effects on the threshold value of chaos, especially the coefficient of the fractional-order derivative term. We also found that the system would be less chaotic with the smaller order and larger coefficient of the fractional-order derivative.The fractional-order derivative was introduced to more accurately describe the nonlinear dynamic behaviors of physical objects with fractional-order characteristics, such as viscoelastic materials, damping, dynamic characteristics of gear clearance, friction, and impact in mechanical systems [45–48]. The analysis shows that the parameters of the fractional-order derivative have an important effect on the chaotic threshold of the system.The above analysis shows that if the system is under stationary excitation amplitude with a small amplitude, chaotic motion will not occur, but under large excitation amplitude, such as impact and chaotic motion, may occur. Since chaotic motion may increase the uncertainty of the system, appropriate parameters should be selected to avoid the occurrence of chaotic phenomena.


## 5. Conclusions

In this paper, threshold for chaos of a Duffing oscillator with the fractional-order derivative is investigated. Furthermore, the effects of various parameters on the threshold value of chaotic motions are analyzed.(1)
The necessary condition for generating chaos in the sense of Smale horseshoes is obtained based on the Melnikov method. The relationship between each parameter is established, and the effect of each parameter on the chaotic threshold is studied, respectively.(2)
The typical dynamical responses, including the largest Lyapunov exponents, time histories, frequency spectrograms, phase portraits, and Poincare maps, are obtained by a numerical simulation method to analyze the dynamic response of the system. It could be found that the system has chaotic motions under certain conditions, and the transition from the steady periodic motion to chaotic motion is a gradual process, which verified the correctness of the necessary condition indirectly.(3)
The effects of the system parameters on the chaotic threshold curve are analyzed in detail. It concluded that the stiffness coefficient and damping coefficient of the system are not only effects on the chaotic threshold curve, and the fractional-order parameters also have a significant effect on the chaotic threshold curve.The corresponding results could provide a reference for analyzing dynamic systems with a fractional-order derivative, which is important to design or control this kind of system.



---

*Source: 1230194-2019-08-27.xml*