- Email: [email protected]
Dynamical analysis of Mathieu equation with two kinds of van der Pol fractional-order terms
International Journal of Non-Linear Mechanics 84 (2016) 130–138
Contents lists available at ScienceDirect
International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm
Dynamical analysis of Mathieu equation with two kinds of van der Pol fractional-order terms Shaofang Wen a, Yongjun Shen b,n, Xianghong Li c, Shaopu Yang b a
Transportation Institute, Shijiazhuang Tiedao University, Shijiazhuang 050043, China Department of Mechanical Engineering, Shijiazhuang Tiedao University, Shijiazhuang 050043, China c Department of Mathematics and Physics, Shijiazhuang Tiedao University, Shijiazhuang 050043, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 19 June 2015 Received in revised form 4 May 2016 Accepted 4 May 2016 Available online 13 May 2016
In this paper the dynamics of Mathieu equation with two kinds of van der Pol (VDP) fractional-order terms is investigated. The approximately analytical solution is obtained by the averaging method. The steady-state solution, existence conditions and stability condition for the steady-state solution are presented, and it is found that the two kinds of VDP fractional coefficients and fractional orders remarkably affect the steady-state solution, which is characterized by the additional damping coefficient (ADC) and additional stiffness coefficient (ASC). The comparisons between the analytical and numerical solutions verify the correctness and satisfactory precision of the approximately analytical solution. The presented typical amplitude–frequency curves illustrate the important effects of two kinds of VDP fractional-order terms on system dynamics. The application of two VDP fractional-order terms in vibration control is discussed. At last, the detailed results are summarized and the conclusions are made. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Mathieu equation Fractional-order derivative Van der Pol oscillator Averaging method
1. Introduction Although fractional calculus had been proposed for more than 300 years, its applications to physics and engineering were just a recent research focus [1,2]. Comparing with the traditional integerorder counterpart, the fractional-order system is much closer to the real nature of the world, and has more advantages, such as strong ability of anti-noise, good robustness, high control precision and so on. The fractional derivatives are an adequate tool to model the frequency-dependent damping behavior of materials and physical systems, and which have played a very important role in various fields such as viscoelasticity, electrochemistry, bioengineering, mechanics, automatic control and signal processing. Accordingly, a lot of researchers in some relevant fields had applied the fractional-order models to solve the problems they met. For example, Gorenflo et al. [3], Jumarie et al. [4], Ishteva [5] and Agnieszka [6] et al. respectively studied the definitions and numerical methods of fractional-order calculus under different senses. Wang and Hu et al. [7,8] investigated a linear single degree-of-freedom oscillator with damping force of fractional-order derivative, and obtained the composition of the solution without external excitation. Shen et al. [9–12] investigated several linear and nonlinear fractional-order oscillators by averaging method, and found that the fractional-order derivatives had both n
Corresponding author. E-mail address: [email protected] (Y. Shen).
http://dx.doi.org/10.1016/j.ijnonlinmec.2016.05.001 0020-7462/& 2016 Elsevier Ltd. All rights reserved.
damping and stiffness effects on the dynamical response in those oscillators. Li et al. [13] discussed the properties of three kinds of fractional derivatives, and the sequential property of the Caputo derivative was also derived. Li et al. [14,15] had done a lot of researches in the mathematical theory of fractional-order calculus, and also established some efficient numerical algorithms. Wahi and Chatterjee [16] studied a special linear single degree-of-freedom oscillator with fractional-order derivative by average method, and analyzed the effects of the fractional-orders derivative. Xu and Li [17] combined Lindstedt–Poincare method with multi-scale method to study fractional-order Duffing oscillator subject to random excitation. Chen and Zhu [18–20] studied some nonlinear fractional-order system with different kinds of noise, and obtained some important statistic properties of the fractional-order system. Caputo et al. [21] presented a new definition of fractional-order derivative which took on two different representations for the temporal and spatial variables. Atangana [22–26] studied new fractional-order derivatives in typical nonlinear equations, and derived some new results about fractional-order derivatives. The well-known Mathieu equation is a linear differential equation with periodic coefficients, and had been applied in physics and engineering fields. Many scholars had studied the Mathieu equation and found that the fractional-order system could generate different dynamical properties from the integer-order counterpart [27–29]. Van der pol (VDP) oscillator could model the typical self-excited or self-sustained oscillation. Leung et al. [30], Sardar et al. [31], Xie and Lin [32] studied different type of
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
fractional-order VDP oscillators by different methods, and found some other important dynamical behaviors. Mathieu-VDP equation occurs in many physics and engineering fields, and it has complex dynamical properties. In recent years, many scholars had studied the Mathieu-VDP equation. For example, Momeni et al. [33] proposed a Mathieu-VDP nonlinear equation to govern the dust grain dynamics in the vicinity of resonance, and finally numerically solved by a fourth-order Runge–Kutta method. Pandey et al. [34] investigated the dynamics of a Mathieu-VDP equation and a Mathieu-VDP-Duffing equation, which was forced both parametrically and nonparametrically. It had been shown that the steady-state response could consist of either 1:1 frequency locking, or quasiperiodic motion. Kalas et al. [35] studied the generalized Mathieu-VDP equation with a small parameter, and the existence of periodic and quasiperiodic solutions was proved by the averaging method and phase space analysis of a derived autonomous equation. Belhaq and Fahsi [36] studied the frequencylocking area of 2:1 and 1:1 resonances in a fast harmonically excited VDP-Mathieu-Duffing oscillator. An averaging technique over the fast excitation was used to derive an equation governing the slow dynamic of the oscillator. Veerman and Verhulst [37] analyzed the VDP-Mathieu equation near and at 1:2 resonance by averaging method, and proved the existence of stable and unstable periodic solutions near the parametric resonance frequency. In this study, we shall consider the Mathieu equation with two kinds of fractional-order VDP terms by averaging method, where the fractional-order derivatives are classified based on their range and could cover all the cases. The paper is organized as follow. Section 2 presents the approximately analytical solution of the Mathieu equation with two kinds of fractional-order VDP terms, where the effects of the two kinds of fractional-order VDP terms are formulated as additional damping coefficient (ADC) and additional stiffness coefficient (ASC). In Section 3 the steady-state solution and the stability condition of the steady-state solution are analyzed. In Section 4, the comparisons between the approximately analytical solution and the numerical one verify the correctness and satisfactory precision of the analytical solution. Moreover, the effects of the parameters in the two kinds of fractional-order VDP terms are also given in this section. At last, the detail results are summarized and the conclusions are made.
2. Approximately analytical solution of Mathieu equation with two kinds of VDP fractional-order terms In this paper, we shall consider the Mathieu equation with two kinds of van der Pol (VDP) fractional-order terms as
x¨ + 2ζẋ + (δ + 2B cos γt ) x + K1 (x2 − 1) D p1 [x (t )] + K2 (x2 − 1) D p2 [x (t )] = 0,
K2 ( K2 > 0) are the fractional coefficients of two kinds of fractionalorder terms. The introduction of fractional-order terms in Mathieu equation lies in some dynamical devices could be modeled by fractional-order derivative, such as viscoelastic and fluid–solid-coupling device. There are several definitions for fractional-order derivative, such as Grünwald–Letnikov, Riemann–Liouville and Caputo definitions [1,2]. Under wide senses, they are equivalent for most mathematical functions. Accordingly, Caputo’s definition is adopted with the form as
∫0
t
x (n)( u)
( t − u) p − n + 1
du , (2)
where Γ ( y ) is Gamma function satisfying Γ ( y + 1) = yΓ ( y ), and the fractional order meets n − 1 < p < n while n is a natural number. Using the following transformations of coordinates
ζ = εμ , δ = ω02, B = εβ, γ = 2ω, K1 = εk1, K2 = εk2 Eq. (1) becomes
x¨ + 2εμẋ + (ω02 + 2εβ cos 2ωt ) x + εk1 (x2 − 1) D p1 [x (t )] + εk2 (x2 − 1) D p2 [x (t )] = 0.
(3)
We shall consider the parametric excitation frequency to be γ = 2ω0 + εσ , where ε < < 1 is a small real parameter, i.e. ω = ω0 + εσ /2. In order to obtain the first approximate solution, it could also be rewritten as ω2 = ω02 + εσω0 . So Eq. (3) becomes
x¨ + ω2x = ε {σω 0 x − 2μẋ − (2β cos 2ωt ) x − k1 (x2 − 1) D p1 [x (t )] − k2 (x2 − 1) D p2 [x (t )]}
(4)
The solution for Eq. (4) could be supposed as
x = a cos φ,
(5a)
ẋ = − aω sin φ,
(5b)
where φ = ωt + θ . Based on the averaging method [38,39], one could obtain the standard equations as
ȧ = −
1 Tω
∫0
1 Tω
aθ ̇ = −
T
∫0
⎡⎣ P1 ( a, θ ) + P2 ( a, θ ) + P3 ( a, θ ) ⎤⎦ sin φdφ,
T
⎡⎣ P1 ( a, θ ) + P2 ( a, θ ) + P3 ( a, θ ) ⎤⎦ cos φdφ,
(6a)
(6b)
where
P1 ( a, θ ) = ε [σω 0 a cos φ + 2μaω sin φ − (2β cos 2ωt ) a cos φ], P2 ( a, θ ) = − εk1 [(a cos φ)2 − 1] D p1 [ a cos φ], P3 ( a, θ ) = − εk2 [(a cos φ)2 − 1] D p2 [ a cos φ]. From the averaging method, one could select the time terminal T as T = 2π when Pi ( a, θ )( i = 1, 2, 3) is periodic function, or T = ∞ when Pi ( a, θ ) ( i = 1, 2, 3) is aperiodic one. Expanding the integral in Eq. (6), one could obtain
ȧ = a1̇ + a2̇ + a3̇ ,
(7a)
aθ ̇ = aθ1̇ + aθ 2̇ + aθ 3̇ ,
(7b)
(1)
where 2ζ , δ , and 2B cos γt are the system linear damping coefficient, constant stiffness coefficient and the periodic time-varying stiffness coefficient respectively. K1 (x2 − 1) D p1 ⎡⎣ x ( t ) ⎤⎦ and K2 (x2 − 1) D p2 ⎡⎣ x ( t ) ⎤⎦ are two kinds of VDP fractional-order terms, where the fractional orders are restricted as 0 < p1 < 1 and 1 < p2 < 2, and K1 ( K1 > 0) and
1 Γ ( n − p)
D p ⎡⎣ x ( t ) ⎤⎦ =
131
The first part in Eq. (7) could be simplified as
a1̇ = −
1 2π ω
∫0
1 2π ω
aθ1̇ = −
2π
∫0
P1 (a, θ ) sin φdφ =
2π
εaβ sin 2θ − εμa, 2ω
P1 (a, θ ) cos φdφ =
εaβ cos 2θ εaσω 0 − . 2ω 2ω
(8a)
(8b)
In order to calculate the second part and the third part for Eq. (7)
a2̇ = − lim
1
T →∞ Tω
= lim
εk1
T →∞ Tω
∫0
∫0 T
T
P2 ( a, θ ) sin φdφ
[(a cos φ)2 − 1] D p1 ⎡⎣ a cos φ⎤⎦ sin φdφ,
(9a)
132
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
T 1 P2 ( a, θ ) cos φdφ T →∞ Tω 0 T εk = lim 1 [(a cos φ)2 − 1] D p1 ⎡⎣ a cos φ⎤⎦ cos φdφ, T →∞ Tω 0
where
∫
aθ 2̇ = − lim
∫
T 1 a3̇ = − lim P3 ( a, θ ) sin φdφ T →∞ Tω 0 T εk = lim 2 [(a cos φ)2 − 1] D p2 ⎡⎣ a cos φ⎤⎦ sin φdφ, T →∞ Tω 0
K (p) = K1ω p1 − 1 cos (9b)
T 1 P3 ( a, θ ) cos φdφ T →∞ Tω 0 T εk = lim 2 [(a cos φ)2 − 1] D p2 ⎡⎣ a cos φ⎤⎦ cos φdφ, T →∞ Tω 0
ȧ = (10a)
∫
εβa a⎛ a2 ⎞ sin 2θ − εμa + ⎜ 1 − ⎟ C (p), 2ω 2⎝ 4⎠
aθ ̇ =
∫
aθ 3̇ = − lim
(10b)
one could expand the integral in Eqs. (9a) and (10a), and get
εk a2 T sin 3φD p1 ⎡⎣ a cos φ⎤⎦ dφ a2̇ = lim 1 T →∞ 4Tω 0 εk a2 T sin φD p1 ⎡⎣ a cos φ⎤⎦ dφ + lim 1 T →∞ 4Tω 0 T εk sin φD p1 ⎡⎣ a cos φ⎤⎦ dφ, − lim 1 T →∞ Tω 0
(11)
εk2 a2 T sin 3φD p2 ⎡⎣ a cos φ⎤⎦ dφ 4Tω 0 εk a2 T sin φD p2 ⎡⎣ a cos φ⎤⎦ dφ + lim 2 T →∞ 4Tω 0 T εk sin φD p2 ⎡⎣ a cos φ⎤⎦ dφ. − lim 2 0 T →∞ Tω
(12)
ȧ =
Ba a⎛ a2 ⎞ sin 2θ − ζa + ⎜ 1 − ⎟ C (p), 2ω 2⎝ 4⎠
aθ ̇ =
∫
∫
a3̇ = lim
T →∞
∫
∫
We introduce some important results about the definite integrals, which had been deduced in Ref. [12,39]
1 T →∞ Tω
∫0
1
∫0
lim
lim
T →∞ Tω
T
T
sin φD p [ a cos φ] dφ = −
cos φD p [ a cos φ] dφ =
⎛ pπ ⎞ aω p − 1 sin ⎜ ⎟, ⎝ 2⎠ 2
⎛ pπ ⎞ aω p − 1 cos ⎜ ⎟, ⎝ 2⎠ 2
(13a)
(13b)
for 0 < p < 1, and
1
T →∞ Tω
∫0
1 T →∞ Tω
∫0
lim
lim
T
T
sin φD p [ a cos φ] dφ = −
⎛p−1 ⎞ aω p − 1 π ⎟, cos ⎜ ⎝ 2 ⎠ 2
(14a)
cos φD p [ a cos φ] dφ = −
⎛p−1 ⎞ aω p − 1 π ⎟, sin ⎜ ⎝ 2 ⎠ 2
(14b)
for 1 < p < 2 respectively. As the similarly procedure, the first parts in Eqs. (11) and (12) would be zero when T → ∞, and the detail deduction procedure is presented in Appendix. One can get the second part and the third part for Eq. (7a)
a2̇ + a3̇ =
a⎛ a2 ⎞ ⎟ C (p), ⎜1 − 2⎝ 4⎠
(15a)
where
C (p) = K1ω p1 − 1 sin
⎛p −1 ⎞ p1π + K2 ω p2 − 1 cos ⎜ 2 π⎟ ⎠ ⎝ 2 2
(15b)
is defined as the additional damping coefficient (ADC). Through the similar procedure, the second part and the third part for Eq. (7b) could also be deduced as
aθ 2̇ + aθ 3̇ =
⎞ a ⎛ 3a2 − 1⎟ K (p), ⎜ ⎠ 2⎝ 4
(16a)
⎞ εβa εσaω 0 a ⎛ 3a2 + ⎜ − 1⎟ K (p). cos 2θ − ⎠ 2ω 2ω 2⎝ 4
(17a)
(17b)
Substituting the original parameters into Eq. (17), we can get
∫
∫
(16b)
is defined as the additional stiffness coefficient (ASC). Substituting Eqs. (8), (15) and (16) into Eq. (7), we obtain
∫
∫
⎛p −1 ⎞ p1π − K2 ω p2 − 1 sin ⎜ 2 π⎟ ⎠ ⎝ 2 2
⎞ εσaω 0 Ba a ⎛ 3a2 + ⎜ − 1⎟ K (p). cos 2θ − ⎠ 2ω 2ω 2⎝ 4
(18a)
(18b)
Eq. (18) shows that all two kinds of VDP fractional-order terms have crucial effects on the damping force and the spring force in Eq. (1), so that the parameters in the two kinds of VDP fractionalorder terms, i.e. the fractional coefficients K1 and K2, the fractional orders p1 and p2, are all important in system dynamics. From the analysis of Eqs. (15b) and (16b), it could be found that the fractional coefficients K1 and K2 are linear with ADC and ASC. ADC and ASC will increase monotonically with the increase of K1. Moreover, ADC will also increase monotonically with the increase of K2, but ASC will decrease monotonically if K2 increases. Accordingly the fractional coefficients K1 and K2 have important effect on the amplitude and the resonance frequency. It should be notable the fractional orders p1 and p2 are more important parameters influencing the system dynamics. It could be found that the fractional orders p1 and p2 affect ADC and ASC by an exponential function and a trigonometric function. From Eqs. (15b) and (16b), it could be found that ADC and ASC are all related with the parametric excitation frequency too. Due to ω0 ¼ 1 in the following parts in this paper, the effect of the parametric excitation frequency could not be considered herein. When K2 ¼0, i.e. there is only the first kind of VDP term in Eq. (1). If p1 varies from 0 to 1, it could be concluded that the increase of p1 would make ADC larger and ASC smaller. When p1 → 0, ADC and ASC will approach 0 and K1 respectively. At that moment the whole system damping coefficient becomes the minimum value, and the whole system stiffness coefficient will be the maximum value. When p1 → 1, ADC will be K1 and ASC will approach 0, so that the whole system damping coefficient will the maximum value and the whole system stiffness coefficient will become the minimum value. When K1 ¼0, i.e. there will exist only the second kind of VDP term in Eq. (1). If p2 varies from 1 to 2, it could be found that the increase of p2 would make ADC and ASC smaller. When p2 → 1, ADC will approach K2 and ASC will be 0. If p2 → 2, ADC will be 0 and ASC will approach −K2. Then the whole system damping coefficient becomes the minimum value and the maximum amplitude is relatively larger. From the analysis of Eq. (16b), it could be found that there is a special case, which is shown in Eq. (19). In this case the ASC is zero, and the whole system stiffness is almost the same as the original linear stiffness.
K2 =
p1π 2 p2 − 1 π 2
K1ω p1 − 1 cos ω p2 − 1 sin
(
)
(19)
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
3. The steady-state solutions and their stability conditions 3.1. The fixed point (0, 0) Eq. (18) shows that a¼ 0 is always one of the system solutions. In order to study its stability condition, one should determine the eigenvalues of the Jacobian matrix of Eq. (18) at a¼ 0, which can be written as
⎡ ∂f ⎢ J = ⎢ ∂a ⎢ ∂g ⎢⎣ ∂a
∂f ⎤ ⎥ ∂θ ⎥ ∂g ⎥ ∂θ ⎥⎦
a=0
(20)
where
⎞ Ba εσaω 0 a ⎛ 3a2 cos 2θ − + ⎜ − 1⎟ K (p). ⎠ 2ω 2ω 2⎝ 4
The eigenvalues of the above Jacobian matrix will be
λ1 = 0, λ2 =
C (p) − ζ. 2
(21)
If 2ζ > C (p), i.e. the linear damping force is greater than ADC, this solution will be stable. From the analysis of Eq. (15b), it could be found that C (p) increases linearly with the increase of K1 and K2. The increase of p1 would make C (p) larger, but the increase of p2 would make C (p) smaller. When p1 → 1 and p2 → 1, two kinds of VDP fractional-order terms can provide the effect of damping force simultaneously, so that C (p) will be the maximum value K1+ K2. In the vicinity of the fixed point (0, 0) a is a small real parameter, so a2
that − 4 C (p) will be about 0, and the damping term of VDP fractional-order terms is negative. That would result into 2ζ < C (p), so that the fixed point (0, 0) will be unstable. When p1 → 0 and p2 → 2, the damping term of two kinds of VDP fractional-order terms is about 0, which will lead to 2ζ > C (p) and the fixed point (0, 0) will be stable. The above analysis is identical with the result in the monograph [40]. 3.2. The existence conditions and stability conditions for the steadystate nonzero solution Now we study the steady-state nonzero solution, which is more important and meaningful in vibration engineering. Letting ȧ = 0 and θ ̇ = 0, and omitting the zero solution, Eq. (18) becomes
⎧ B 1⎛ a¯ 2 ⎞ ⎪ sin 2θ¯ − ζ + ⎜ 1 − ⎟ C (p) = 0 ⎪ 2ω 2⎝ 4⎠ ⎨ , ⎞ ⎪ B εσω 0 1 ⎛ 3a¯ 2 ¯ − + − ( ) = cos 2 K p θ 1 0 ⎜ ⎟ ⎪ ⎠ ω 2⎝ 4 ⎩ 2ω
(22)
where a¯ and θ¯ are the amplitude and phase of the steady-state response respectively. After eliminating θ¯ from Eq. (22), one could obtain the amplitude-frequency equation
⎡ εσω 0 ⎤2 B2 [(1 − ρ) C (p) − 2ζ ]2 + ⎢ (3ρ − 1) K (p) − ⎥⎦ = 2 , ⎣ ω ω where ρ =
a¯ 2 . 4
where
(23)
⎡ εσω0 ⎤2 B2 H3 = [2ζ − C (p)]2 + ⎢ K (p) + ⎥− . ⎣ ω ⎦ ω2 And the phase-frequency equation is
⎤ ⎡ (1 − ρ) C (p) − 2ζ ⎥ 1 θ¯ = arctan ⎢ εσω 0 . 2 ⎢⎣ (3ρ − 1) K (p) − ω ⎥⎦
(25)
Due to H1 > 0 , the existence conditions for Eq. (24) solution could be summarized into two cases. The first case is that Eq. (24) has only one real solution when H3 o 0. The other case means Eq. H22 4H1
.
From the analysis of Eq. (24), it could be found that H1, H2 and H3 are related to ADC and ASC. Moreover, H2 and H3 are related to the system linear damping coefficient, the periodic time-varying coefficient and parametric excitation frequency. Accordingly, the range of nonzero solutions has complex interaction with two kinds of VDP fractional-order terms and the system parameters. Now, we select the system parameters as ζ ¼0.1, ε ¼0.01, B ¼0.25, and study the effects of the two kinds of VDP fractionalorder parameters K1, K2, p1 and p2 on the existence range of nonzero solutions. The results are shown in Fig. 1(a)–(d), where the X-axis is for parametric excitation frequency ω , the Y-axis denote for K1, K2, p1, p2 respectively, and the Z-axis is for the existence boundary of the solution for Eq. (24). Obviously the parts of greater than zero are for the existence range of nonzero solutions. From the analysis of Fig. 1, it could be found that the increase of K1 and p2 would make the existence range of nonzero solutions larger, but the increase of p1 would make the existence range of nonzero solutions smaller. Furthermore, the increase of K2 would make the existence range of nonzero solutions larger at first, and then smaller if K2 is big enough. Next we study the stability of the steady-state nonzero solution. Letting a = a¯ + Δa , θ = θ¯ + Δθ , and substituting them into Eq. (18), one can obtain
¯ cos 2θ¯ dΔa ⎡ B sin 2θ¯ C (p) ⎛ aB 3a¯ 2 ⎞ ⎤ =⎢ −ζ+ Δθ , ⎜1 − ⎟ ⎥ Δa + ⎣ dt 2ω 2 ⎝ ω 4 ⎠⎦ (26a) ¯ sin 2θ¯ dΔθ 3K (p) a¯ aB = Δa + Δθ . dt 4 ω
(26b)
Based on Eq. (22), one could eliminate θ¯ from Eq. (26) and get the characteristic determinant as
⎡ ⎢ A1 − λ det ⎢ ⎢ 3a¯ K (p) ⎣ 4
A1 = −
εσω0 a¯ ω
⎤ − 1 K (p ) ⎥ ⎥ = 0,where ⎥ A2 − λ ⎦
− a¯
(
3a¯ 2 4
)
⎛ a¯ 2 a¯ 2 ⎞ C (p), A2 = − 2ζ + ⎜ 1 − ⎟ C (p). ⎝ 4 4⎠
Expanding the characteristic determinant and substituting
ρ=
a¯ 2 4
into it, we could get the characteristic equation as
⎡ εσω 0 ⎤ λ2 − (A1 + A2 ) λ + A1A2 − 3ρK (p) ⎢ − (3ρ − 1) K (p) ⎥ = 0. ⎣ ω ⎦
(27)
Based on Eq. (27), the stability conditions for the steady-state solution could be established as
Expanding Eq. (23), we obtain
H1ρ2 + H2 ρ + H3 = 0,
⎡ εσω0 ⎤ ⎥, H1 = C (p)2 + 9K (p)2 , H2 = 2C (p)[2ζ − C (p)] − 6K (p)⎢ K (p) + ⎢⎣ ω ⎥⎦
(24) will have two real solutions when H2 o0 and 0 o H3 o
Ba a⎛ a2 ⎞ f (a, θ ) = ȧ = sin 2θ − ζa + ⎜ 1 − ⎟ C (p) 2ω 2⎝ 4⎠ ⎡ ⎤ Ba 1 a2 = sin 2θ − a ⎢ ζ − C (p) + C (p) ⎥, 2ω 2 4 ⎣ ⎦ g (a, θ ) = aθ ̇ =
133
(24)
A1 + A2 < 0, and
(28a)
134
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
Fig. 1. The effects of the fractional parameters on the existence range of nonzero solutions: (a) the effects of K1 ( K2 ¼ 0, p1 ¼0.5), (b) the effects of p1 ( K1 ¼0.1, K2 ¼ 0), (c) the effects of K2 ( K1 ¼0, p2 ¼1.5) and (d) the effects of p2 ( K1 ¼0, K2 ¼ 0.2).
⎡ εσω 0 ⎤ − (3ρ − 1) K (p) ⎥ > 0. A1A2 − 3ρK (p) ⎢ ⎣ ω ⎦
(28b)
After expanding Eq. (28), the equivalent forms of the stability conditions for the steady-state solution are
T1 = 2ρC (p) + [2ζ − C (p)] > 0,
condition T1 40. Then one should analyze whether the condition T2 40 holds based on the existence conditions of Eq. (24). The results are summarized as follows: a. when H3 o0, Eq.(24) has one real solution
(29a)
ρ=
and
⎧ ⎡ εσω 0 ⎤ ⎫ T2 = [C (p)2 + 9K (p)2] ρ2 + ⎨ C (p)[2ζ − C (p)] − 3K (p) ⎢ K (p) + ⎥⎬ ⎣ ⎩ ω ⎦⎭ ρ > 0. The stability conditions are equal to
T2 = 2H1ρ + H2 > 0,
H22 − 4H1H3
−H2 +
2H1
.
(30)
Substituting Eq. (30) into Eq. (29b), one could obtain
T2 = H22 − 4H1H3 40. This will be the fourth cases, so that the solution is stable. H2 b. when H2 o0 and 0 o H3 o 2 , Eq. (24) has two real solutions 4H1
(29b)
where T1 may be defined as the equivalent nonlinear damping coefficient and T2 is the nonlinear parameter of the stability condition. T1 and T2 are generally used to judge the stability for the steadystate nonzero solution. It can be found that the stability conditions for the steady-state solution could be divided into four cases, i. e., 1) T1 o0 and T2 o0, 2) T1 o0 and T2 40, 3) T1 40 and T2 o0, 4) T1 40 and T2 4 0. When T1 o 0, the system response will be divergent and the nonzero steady-state solutions are unstable. When T1 4 0 and T2 o0, the system response will be unstable periodic vibration because T2 is affected by the stiffness and damping at the same time. If T1 40 and T2 40, the steady-state solution will be stable. This implies that only the solution corresponding to the forth case is stable. Now, the research on stability conditions could be focused on the two parameters. At first one should select appropriate system parameter ζ so as to satisfy the
ρ1 =
−H2 +
H22 − 4H1H3 2H1
and ρ2 =
−H2 −
H22 − 4H1H3 2H1
.
Substituting ρ1 into Eq. (29b), we could obtain T2 40, so that this solution is stable. On the contrary, substituting ρ2 into Eq. (29b) will result into T2 o0, so that this solution will unstable.
4. Numerical simulations 4.1. Comparison between the approximate analytical solutions and the numerical results An illustrative example system is studied herein as defined by the basic system parameter, p1 ¼0.5, p2 ¼ 1.5, ζ ¼ 0.1, ε ¼ 0.01 and B ¼0.25. Three group of fractional-order parameters, i.e., K1 ¼ 0.1 and K2 ¼0.1, K1 ¼ 0.2 and K2 ¼0.4, K1 ¼0.2 and K2 ¼0.5, are selected to study the different response modes. Based on the amplitude–
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
Fig. 4. The amplitude–frequency curves ( K1 ¼ 0.2, K2 ¼ 0.5).
Fig. 2. The amplitude–frequency curves ( K1 ¼ 0.1 , K2 ¼ 0.1).
Fig. 5. The effects of the fractional coefficient K1 on the amplitude–frequency curves ( p1 ¼ 0.5, K2 ¼ 0).
Fig. 3. The amplitude–frequency curves ( K1 ¼0.2, K2 ¼ 0.4).
frequency equation Eq. (24), one could analytically obtain the amplitude–frequency curves shown in Figs. 2–4, where the solid line is for the asymptotically stable solution and the dash line denotes the unstable one. In the three figures, the parametric excitation frequency parameter ω is for the horizontal axis, and the steady-state amplitude a is for the vertical axis. In order to verify the precision of the analytical solutions, we also present the numerical results. The numerical method is [1,2]
D p ⎡⎣ x ( tl ) ⎤⎦ ≈ h−p
l
∑ C jp x ( tl − j ), j=0
(31)
where tl = lh is the time sample points, h is the sample step and C jp is the fractional binomial coefficient with the iterative relationship as
⎛ 1 + p⎞ p C0p = 1, C jp = ⎜ 1 − ⎟C . ⎝ j ⎠ j−1
135
(32)
Here we select h ¼ 0.001, and the total computation time is 300 s. The amplitude–frequency curves by numerical integration are also shown in Figs. 2–4, where those circles are for the numerical results. From the observation of Figs. 2–4, one could conclude that the approximately analytical solutions agree very well
with the numerical results and could present satisfactory precision even if the amplitude–frequency curves are different topological structures. 4.2. Effects of VDP fractional-order terms on system dynamics When the coefficient K1 of the first kind of VDP fractional-order term is changed, the amplitude–frequency curves are shown in Fig. 5, where K2 ¼0. From the observation of Fig. 5, it could be found that the larger the fractional coefficient K1, the smaller the maximum amplitude. This is due to the fact that ADC becomes larger in this procedure. Furthermore, it could also be found that the increase of the fractional coefficient K1 would make ASC larger, which could result into the rightwards bending of the backbone of the amplitude–frequency curves, and make the existence range of the steady-state solution larger. When the order p1 of the first kind of VDP fractional-order term varies from 0 to 1, the amplitude–frequency curves are shown in Fig. 6, where K2 ¼0. From the observation of Fig. 6, it could be found that the larger the fractional order p1, the smaller the maximum amplitude. The reason is that ADC will also become larger in this procedure. Furthermore, it could also be found that the increase of the fractional coefficient p1 would make ASC smaller, which could result into the leftwards bending of the
136
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
Fig. 6. the effects of the fractional order p1 on the amplitude–frequency curves ( K1 ¼ 0.1, K2 ¼0).
Fig. 8. (a) the effects of the fractional order p2 on the amplitude–frequency curves ( K1 ¼0, K2 ¼ 0.2). (b) the effects of the fractional order p2 on the amplitude–frequency curves ( K1 ¼ 0.1, p1 ¼0.5, K2 ¼ 0.2).
Fig. 7. (a) The effects of the fractional coefficient K2 on the amplitude–frequency curves ( K1 ¼0, p2 ¼1.5). (b) The effects of the fractional coefficient K2 on the amplitude–frequency curves ( K1 ¼0.2, p1 ¼ 0.5, p2 ¼1.5).
backbone of the amplitude–frequency curves, and make the existence range of the steady-state nonzero solution smaller. Moreover, the changes of order p1 will also change the topology structures of the amplitude–frequency curves. The role of the second kind of VDP fractional-order term is similar to that of the first kind of VDP fractional-order term. The detail discussions about the effects of the second kind of VDP fractional-order term are analyzed in the following parts, and the results are shown in Figs. 7 to 8. In Fig. 7 the subfigure 7(a) and 7(b) are for K1 ¼0 and K1 ¼0.2 respectively. From the observation of those curves in Fig. 7, it could be found that the increase of the fractional coefficient K2 would make ADC larger, which will make the maximum amplitude smaller. Simultaneously, it is easy to find that the increase of the fractional coefficient K2 will make ASC smaller, which could result into the leftwards bending of the backbone of the amplitude–frequency curves, and make the existence range of the steady-state solution larger at first and then smaller. The changes of the topological structures could also be found in Fig. 7. Additionally, the difference between Figs. 7(a) and 7(b) is the existence of the first kind of the VDP fractional-order term, which decreases the maximum amplitude, and diminishes the leftwards bending of the backbone of the amplitude–frequency curves. It will affect the existence range of the steady-state solution and the topological
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
structure at the same time. When the order p2 of the second kind of VDP fractional-order term varies from 1 to 2, the results are shown in Fig. 8, where Fig. 8 (a) is for K1 ¼ 0 and Fig. 8(b) is for K1 ¼0.1 respectively. From the observation of Fig. 8, it could be found that the larger the fractional order p2, the larger the maximum amplitude. This is due to the fact that ADC becomes smaller in this procedure. The existence range of the steady-state solution becomes larger too. Moreover, the topological structures is changed, and the backbone of the amplitude–frequency curves will bend towards the left because ASC becomes smaller in this procedure according to Eq. (16b). When p2 → 2, one could conclude that the maximum amplitude is too large, and the backbone of the amplitude–frequency curve bends towards the left. When p2 → 1, it could be found that the maximum amplitude is small, and the bending degree of the backbone of the amplitude–frequency curve is little. Similar to Fig. 7, the difference between Fig. 8(a) and (b) is the existence of the first kind of the VDP fractional-order term. When p2 is the same value in Fig. 8(a) and (b), one could conclude that the existence of the first kind of the VDP fractional-order term could also decrease the maximum amplitude, and diminish the leftwards bending of the backbone of the amplitude–frequency curves. It will affect the existence range of the steady-state solution and the topological structure at the same time. From the observation of Fig. 8(b), it could be found that the maximum amplitude is smaller when two kinds VDP fractional-order terms are exist. The backbone of the amplitude–frequency curve bends from the right to the left, and the existence range of the steady-state solution is smaller at first and then larger. Moreover, it could be found that there is a special case. In this case the ASC is zero, and the ADC is almost the maximum, so that the maximum amplitude will be reduced significantly. This special case may be greatly important in vibration control engineering. Based on the above analysis, one could select appropriate fractional parameters to improve the dynamical characteristics of Mathieu equation with two kinds of VDP fractional-order terms. Comparing with the integral-order counterpart, the response modes of the Mathieu equation with VDP fractional-order terms are more complicated. Especially, two kinds of VDP fractional-order terms will bring better control performance than the single first kind of VDP fractional-order term.
5. Conclusions The Mathieu equation with two kinds of VDP fractional-order terms is studied based on the averaging method, and the approximately analytical solution and the amplitude–frequency equation are obtained. The steady-state solution and stability condition of the steady-state solution are researched. The effects of two kinds of VDP fractional coefficients and fractional orders on the steady-state solution are analyzed, and it is found that these two kinds of VDP fractional-order terms affect the linear damping and the linear stiffness, which are characterized by ADC and ASC. The effects of the two kinds of VDP fractional-order parameters on the existence range of nonzero solutions and the stability of the steady-state nonzero solution are studied. The comparisons of the analytical solution with the numerical one verify the correctness of the analytical solution. When two kinds of VDP fractional-order term exist, it would be found that there is a special case. In this case the ASC is almost zero, and the ADC is almost the maximum, so that the maximum amplitude will become almost the minimum. The special case may be greatly important in vibration control engineering. The results are meaningful to analyze and control the dynamical behavior of the similar fractional-order system, and could present useful reference.
137
Acknowledgment The authors are grateful to the support by National Natural Science Foundation of China (No. 11372198), the Cultivation plan for Innovation team and leading talent in Colleges and universities of Hebei Province (LJRC018), the Program for advanced talent in the universities of Hebei Province (GCC2014053), the Program for advanced talent in Hebei (A201401001), and the education department project of Hebei Province (QN2016258)
Appendix A Taking the first parts in Eq. (11) as an example, According to Caputo’s definition in Eq. (2), and ϕ = ωt + θ , the first parts in Eq. (11) becomes lim
T →∞
εk1a 2 4Tω
∫0
εk1a 2 T →∞ 4Tω
= lim
T
sin 3 (ωt + θ ) D p1 [ a cos (ωt + θ ) ] dφ
∫0
T
⎡ 1 sin 3 (ωt + θ ) ⎢ ⎢⎣ Γ ( 1 − p)
∫0
t
⎤ −aω sin (ωu + θ ) ⎥ du dt . p ( t − u) ⎦⎥
(A.1)
Letting s = t − u and du = − ds , Eq. (A.1) becomes lim
T →∞
εk1a 2 4Tω
∫0
T
⎤ −aω sin (ωu + θ ) ⎥ du dt p ⎥⎦ t − u ( )
⎡ 1 sin 3 (ωt + θ ) ⎢ ⎢⎣ Γ ( 1 − p)
∫0
T
⎡ sin 3 (ωt + θ ) ⎢ ⎣
∫0
T
⎡ ⎢ sin 3 (ωt + θ ) sin (ωt + θ ) ⎣
∫0
T
⎡ ⎢ sin 3 (ωt + θ ) cos (ωt + θ ) ⎣
∫0
=
1 −εk1a 3 lim 4Γ ( 1 − p) T →∞ T
∫0
=
1 −εk1a 3 lim 4Γ ( 1 − p) T →∞ T
∫0
+
1 εk1a 3 lim 4Γ ( 1 − p) T →∞ T
∫0
t
t
sin (ωt + θ − ωs) ⎤ ds⎥ dt ⎦ sp t
cos ωs ⎤ ds⎥ dt ⎦ sp
τ
sin ωs ⎤ ds⎥ dt . ⎦ sp
(A.2)
Defining the first part as A1 and the second part as A2 in Eq. (A.2), and integrating it by parts respectively, one could obtain −εk1a3 1 lim 4Γ ( 1 − p ) T →∞ T
A1 =
=
=
⎧ 1 ⎨ − [ cos 4 (ωt + θ ) − cos 2 (ωt + θ ) ] ⎩ 2
T
εk1a3 1 ⎡ sin 4 (ωt + θ ) sin 2 (ωt + θ ) ⎤ − lim ⎢ ⎥ ⎦ 4ω 2ω 8Γ ( 1 − p ) T →∞ T ⎣ −
A2 =
∫0
εk1a3 1 lim 8Γ ( 1 − p ) T →∞ T
1 εk1a 3 lim 4Γ ( 1 − p) T →∞ T
∫0
∫0
T
T
∫0
1 εk1a 3 lim 8Γ ( 1 − p) T →∞ T
cos ωs ⎫ ds⎬ dt ⎭ sp T
cos ωs ds sp 0
⎧1 ⎨ ⎡⎣ sin 4 (ωt + θ ) + sin 2 (ωt + θ ) ⎤⎦ ⎩2
∫0
t
⎡ sin 4 (ωt + θ ) τ sin 2 (ωt + θ ) ⎤ cos ωt − dt , ⎢ ⎥ ⎣ ⎦ tp 4ω 2ω
1 ⎡ cos 4 (ωt + θ ) cos 2 (ωt + θ ) ⎤ −εk1a 3 lim ⎢ − ⎥ ⎦ 4ω 2ω 8Γ ( 1 − p) T →∞ T ⎣ +
t
∫0
T
∫0
t
∫0
τ
(A.3)
sin ωs ⎫ ds⎬ dt ⎭ sp T
sin ωs ds sp 0
⎡ cos 4 (ωt + θ ) τ cos 2 (ωt + θ ) ⎤ sin ωt dt . − ⎥ ⎢ ⎦ tp ⎣ 4ω 2ω
(A.4)
For simplicity, one could define the first and second part in Eq. (A.3) as BA11 and BA12, the first and second part in Eq. (A.4) as BA21 and BA22 respectively. We introduce two important formulae, which had been deduced in Refs. [11,12]
∫ T →∞ 0 lim
lim
T →∞
∫0
T
T
⎛ pπ ⎞ sin (ωt ) dt = ω p − 1Γ (1 − p) cos ⎜ ⎟ ⎝ 2⎠ tp
(A.5a)
⎛ pπ ⎞ cos (ωt ) dt = ω p − 1Γ (1 − p) sin ⎜ ⎟ ⎝ 2⎠ tp
(A.5b)
Substituting Eq. (A.5b) into BA11 and integrating BA12, we obtain BA11 =
εk1a 3 1 ⎡ sin 4 (ωT + θ ) sin 2 (ωT + θ ) ⎤ − lim ⎢ ⎥ ⎦ 8Γ ( 1 − p) T →∞ T ⎣ 4ω 2ω pπ
=
εk1a 3ω p − 1 sin ( 2 ) 8
∫0
T
cos ωs ds sp
1 ⎡ sin 4 (ωT + θ ) sin 2 (ωT + θ ) ⎤ − lim ⎢ ⎥ = 0, ⎣ ⎦ 4ω 2ω
T →∞ T
(A.6a)
138
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
BA12 = −
εk1a3 sin 2 (ωt + θ ) ⎤ cos ωt 1 T ⎡ sin 4 (ωt + θ ) τ − lim ∫ ⎢ dt ⎥ ⎦ tp 4ω 2ω 8Γ ( 1 − p) T →∞ T 0 ⎣
= 0.
(A.6b)
Similarly, BA21 and BA22 could also be deduced BA21 ¼ 0, BA22 ¼0. So, the first parts in Eq. (11) would be zero when T → ∞. As the similar procedure, one could find that the first parts in Eq. (12) would be zero too when T → ∞.
References [1] I. Petras, Fractional-order Nonlinear Systems: Modeling, Analysis and Simulation, Higher Education Press, Beijing, 2011. [2] R. Caponetto, G. Dongola, L. Fortuna, I. Petras, Fractional Order Systems: Modeling and Control Applications, World Scientific, New Jersey, 2010. [3] R. Gorenflo, E.A. Abdel-Rehim, Convergence of the Grünwald–Letnikov scheme for time-fractional diffusion, Journal. Comput. Appl. Math. 205 (2) (2007) 871–881. [4] G. Jumarie, Modified Riemann–Liouville derivative and fractional Taylor series of nondifferentiable functions further results, Comput. Math. Appl. 51 (2006) 1367–1376. [5] M. Ishteva, R. Scherer, L. Boyadjiev, On the Caputo operator of fractional calculus and C-Laguerre functions, Math. Sci. Res. J. 9 (6) (2005) 161–170. [6] Agnieszka B. Malinowska, Delfim F.M. Torres, Fractional calculus of variations for a combined Caputo derivative, Fract. Calc. Appl. Anal. 14 (4) (2011) 523–537. [7] Z.H. Wang, H.Y. Hu, Stability of a linear oscillator with damping force of fractional-order derivative, Sci. China G: Phys. Mech. Astron. 39 (10) (2009) 1495–1502. [8] Z.H. Wang, M.L. Du, Asymptotical behavior of the solution of a SDOF linear fractionally damped vibration system, Shock. Vib. 18 (1–2) (2011) 257–268. [9] Y.J. Shen, P. Wei, S.P. Yang, Primary resonance of fractional-order van der Pol oscillator, Nonlinear Dyn. 77 (2014) 1629–1642. [10] Y.J. Shen, P. Wei, C.Y. Sui, S.P. Yang, Subharmonic resonance of Van Der Pol oscillator with fractional-order derivative, Math. Probl. Eng. (2014) 7380871–738087-17. [11] Y.J. Shen, S.P. Yang, H.J. Xing, et al., Primary resonance of Duffing oscillator with fractional-order derivative, Commun. Nonlinear Sci. Numer. Simul. 17 (7) (2012) 3092–3100. [12] Y.J. Shen, S.P. Yang, H.J. Xing, et al., Primary resonance of Duffing oscillator with two kinds of fractional-order derivatives, Int. J. Non-Linear Mech. 47 (9) (2012) 975–983. [13] C.P. Li, W.H. Deng, Remarks on fractional derivatives, Appl. Math. Comput. 187 (2) (2007) 777–784. [14] J.X. Cao, H.F. Ding, C.P. Li, Implicit difference schemes for fractional diffusion equations, Commun. Appl. Math. Comput. 27 (1) (2013) 61–74. [15] F.H. Zeng, C.P. Li, High-order finite difference methods for time-fractional subdiffusion equation, Chin. Journal. Comput. Phys. 30 (4) (2013) 491–500. [16] P. Wahi, A. Chatterjee, Averaging oscillations with small fractional damping and delayed terms, Nonlinear Dyn. 38 (1-4) (2004) 3–22. [17] Y. Xu, Y.G. Li, D. Liu, W.T. Jia, H. Huang, Responses of Duffing oscillator with fractional damping and random phase, Nonlinear Dyn. 74 (3) (2013) 745–753.
[18] L.C. Chen, W.Q. Zhu, Stochastic jump and bifurcation of Duffing oscillator with fractional derivative damping under combined harmonic and white noise excitations, Int. J. Non-linear Mech. 46 (12) (2011) 1324–1329. [19] L.C. Chen, H.F. Li, Z.S. Li, W.Q. Zhu, First passage failure of single-degree-offreedom nonlinear oscillators with fractional derivative, J. Vib. Control. 19 (14) (2013) 2154–2163. [20] L.C. Chen, W.H. Wang, Z.S. Li, W.Q. Zhu, Stationary response of Duffing oscillator with hardening stiffness and fractional derivative, Int. J. Non-linear Mech. 48 (1) (2013) 44–50. [21] M. Caputo, M. Fabrizio, A new definition of fractional derivative without singular kernel, Progr. Fract. Differ. Appl. 1 (2) (2015) 1–13. [22] A. Atangana, D. Baleanu, New fractional derivatives with nonlocal and nonsingular kernel: theory and application to heat transfer model, Therm. Sci. (2016), http://dx.doi.org/10.2298/TSCI160111018A. [23] A. Atanganaa, I. Kocab, On the new fractional derivative and application to nonlinear Baggs and Freedman model, J. Nonlinear Sci. Appl. 9 (2016) 2467–2480. [24] A. Atangana, I. Koca, Chaos in a simple nonlinear system with Atangana–Baleanu derivatives with fractional order, Chaos, Solitons Fractals (2016) 10, http://dx.doi.org/10.1016/j.chaos.2016.02.012. [25] A. Atangana, B.S.T. Alkahtani, Analysis of the Keller–Segel model with a fractional derivative without singular kernel, Entropy 17 (6) (2015) 4439–4453. [26] A. Atangana, On the new fractional derivative and application to nonlinear Fisher's reaction–diffusion equation, Appl. Math. Comput. 273 (2016) 948–956. [27] A. Ebaid, D.M.M. ElSayed, M.D. Aljoufi, Fractional calculus model for damped mathieu equation: approximate analytical solution, Appl. Math. Sci. 6 (82) (2012) 4075–4080. [28] R.H. Rand, S.M. Sah, M.K. Suchorsky, Fractional Mathieu equation, Commun. Nonlinear Sci. Numer. Simul. 15 (11) (2010) 3254–3262. [29] A.Y.T. Leung, Z.J. Guo, H.X. Yang, Transition curves and bifurcations of a class of fractional Mathieu-type equations, Int. Journal. Bifurc. Chaos 22 (2012) 1250275-1–1250275-13. [30] A.Y.T. Leung, H.X. Yang, Z.J. Guo, The residue harmonic balance for fractional order van der Pol like oscillators, Journal. Sound. Vib. 331 (2012) 1115–1126. [31] T. Sardar, R.S. Saha, R.K. Bera, B.B. Biswas, The analytical approximate solution of the multi-term fractionally damped van der Pol equation, Phys. Scr. 80 (2009) 1–6. [32] F. Xie, X.Y. Lin, Asymptotic solution of the van der Pol oscillator with small fractional damping, Phys. Scr. 136 (2009) 014033-1-3-4. [33] M. Momeni, I. Kourakis, M. Moslehi-Fard, P.K. Shukla, A Van der Pol-Mathieu equation for the dynamics of dust grain charge in dusty plasmas, J. Phys. A: Math. Theor. 40 (2007) 473–481. [34] M. Pandey, R.H. Rand, A.T. Zehnder, Frequency locking in a forced Mathieu-van der Pol-Duffing system, Nonlinear Dyn. 54 (1–2) (2008) 3–12. [35] J. Kalas, Z. Kadeřábek, Periodic solutions of a generalized Van der Pol–Mathieu differential equation, Appl. Math. Comput. 234 (2014) 192–202. [36] M. Belhaq, A. Fahsi, 2:1 and 1:1 frequency-locking in fast excited van der Pol– Mathieu–Duffing oscillator, Nonlinear Dyn. 53 (1–2) (2008) 139–152. [37] F. Veerman, F. Verhulst, Quasiperiodic phenomena in the van der Pol-Mathieu equation, J. Sound. Vib. 326 (1–2) (2009) 314–320. [38] J.A. Sanders, F. Verhulst, J. Murdock, Averaging Methods in Nonlinear Dynamical Systems, 2nd ed, Springer-Verlag, New York, 2007. [39] Y.J. Shen, S.P. Yang, C.Y. Sui, Analysis on limit cycle of fractional-order van der Pol oscillator, Chaos, Solitons Fractals 67 (2014) 94–102. [40] A.H. Nayfeh, Nonlinear Oscillations, Wiley, New York, 1973.
Contents lists available at ScienceDirect
International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm
Dynamical analysis of Mathieu equation with two kinds of van der Pol fractional-order terms Shaofang Wen a, Yongjun Shen b,n, Xianghong Li c, Shaopu Yang b a
Transportation Institute, Shijiazhuang Tiedao University, Shijiazhuang 050043, China Department of Mechanical Engineering, Shijiazhuang Tiedao University, Shijiazhuang 050043, China c Department of Mathematics and Physics, Shijiazhuang Tiedao University, Shijiazhuang 050043, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 19 June 2015 Received in revised form 4 May 2016 Accepted 4 May 2016 Available online 13 May 2016
In this paper the dynamics of Mathieu equation with two kinds of van der Pol (VDP) fractional-order terms is investigated. The approximately analytical solution is obtained by the averaging method. The steady-state solution, existence conditions and stability condition for the steady-state solution are presented, and it is found that the two kinds of VDP fractional coefficients and fractional orders remarkably affect the steady-state solution, which is characterized by the additional damping coefficient (ADC) and additional stiffness coefficient (ASC). The comparisons between the analytical and numerical solutions verify the correctness and satisfactory precision of the approximately analytical solution. The presented typical amplitude–frequency curves illustrate the important effects of two kinds of VDP fractional-order terms on system dynamics. The application of two VDP fractional-order terms in vibration control is discussed. At last, the detailed results are summarized and the conclusions are made. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Mathieu equation Fractional-order derivative Van der Pol oscillator Averaging method
1. Introduction Although fractional calculus had been proposed for more than 300 years, its applications to physics and engineering were just a recent research focus [1,2]. Comparing with the traditional integerorder counterpart, the fractional-order system is much closer to the real nature of the world, and has more advantages, such as strong ability of anti-noise, good robustness, high control precision and so on. The fractional derivatives are an adequate tool to model the frequency-dependent damping behavior of materials and physical systems, and which have played a very important role in various fields such as viscoelasticity, electrochemistry, bioengineering, mechanics, automatic control and signal processing. Accordingly, a lot of researchers in some relevant fields had applied the fractional-order models to solve the problems they met. For example, Gorenflo et al. [3], Jumarie et al. [4], Ishteva [5] and Agnieszka [6] et al. respectively studied the definitions and numerical methods of fractional-order calculus under different senses. Wang and Hu et al. [7,8] investigated a linear single degree-of-freedom oscillator with damping force of fractional-order derivative, and obtained the composition of the solution without external excitation. Shen et al. [9–12] investigated several linear and nonlinear fractional-order oscillators by averaging method, and found that the fractional-order derivatives had both n
Corresponding author. E-mail address: [email protected] (Y. Shen).
http://dx.doi.org/10.1016/j.ijnonlinmec.2016.05.001 0020-7462/& 2016 Elsevier Ltd. All rights reserved.
damping and stiffness effects on the dynamical response in those oscillators. Li et al. [13] discussed the properties of three kinds of fractional derivatives, and the sequential property of the Caputo derivative was also derived. Li et al. [14,15] had done a lot of researches in the mathematical theory of fractional-order calculus, and also established some efficient numerical algorithms. Wahi and Chatterjee [16] studied a special linear single degree-of-freedom oscillator with fractional-order derivative by average method, and analyzed the effects of the fractional-orders derivative. Xu and Li [17] combined Lindstedt–Poincare method with multi-scale method to study fractional-order Duffing oscillator subject to random excitation. Chen and Zhu [18–20] studied some nonlinear fractional-order system with different kinds of noise, and obtained some important statistic properties of the fractional-order system. Caputo et al. [21] presented a new definition of fractional-order derivative which took on two different representations for the temporal and spatial variables. Atangana [22–26] studied new fractional-order derivatives in typical nonlinear equations, and derived some new results about fractional-order derivatives. The well-known Mathieu equation is a linear differential equation with periodic coefficients, and had been applied in physics and engineering fields. Many scholars had studied the Mathieu equation and found that the fractional-order system could generate different dynamical properties from the integer-order counterpart [27–29]. Van der pol (VDP) oscillator could model the typical self-excited or self-sustained oscillation. Leung et al. [30], Sardar et al. [31], Xie and Lin [32] studied different type of
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
fractional-order VDP oscillators by different methods, and found some other important dynamical behaviors. Mathieu-VDP equation occurs in many physics and engineering fields, and it has complex dynamical properties. In recent years, many scholars had studied the Mathieu-VDP equation. For example, Momeni et al. [33] proposed a Mathieu-VDP nonlinear equation to govern the dust grain dynamics in the vicinity of resonance, and finally numerically solved by a fourth-order Runge–Kutta method. Pandey et al. [34] investigated the dynamics of a Mathieu-VDP equation and a Mathieu-VDP-Duffing equation, which was forced both parametrically and nonparametrically. It had been shown that the steady-state response could consist of either 1:1 frequency locking, or quasiperiodic motion. Kalas et al. [35] studied the generalized Mathieu-VDP equation with a small parameter, and the existence of periodic and quasiperiodic solutions was proved by the averaging method and phase space analysis of a derived autonomous equation. Belhaq and Fahsi [36] studied the frequencylocking area of 2:1 and 1:1 resonances in a fast harmonically excited VDP-Mathieu-Duffing oscillator. An averaging technique over the fast excitation was used to derive an equation governing the slow dynamic of the oscillator. Veerman and Verhulst [37] analyzed the VDP-Mathieu equation near and at 1:2 resonance by averaging method, and proved the existence of stable and unstable periodic solutions near the parametric resonance frequency. In this study, we shall consider the Mathieu equation with two kinds of fractional-order VDP terms by averaging method, where the fractional-order derivatives are classified based on their range and could cover all the cases. The paper is organized as follow. Section 2 presents the approximately analytical solution of the Mathieu equation with two kinds of fractional-order VDP terms, where the effects of the two kinds of fractional-order VDP terms are formulated as additional damping coefficient (ADC) and additional stiffness coefficient (ASC). In Section 3 the steady-state solution and the stability condition of the steady-state solution are analyzed. In Section 4, the comparisons between the approximately analytical solution and the numerical one verify the correctness and satisfactory precision of the analytical solution. Moreover, the effects of the parameters in the two kinds of fractional-order VDP terms are also given in this section. At last, the detail results are summarized and the conclusions are made.
2. Approximately analytical solution of Mathieu equation with two kinds of VDP fractional-order terms In this paper, we shall consider the Mathieu equation with two kinds of van der Pol (VDP) fractional-order terms as
x¨ + 2ζẋ + (δ + 2B cos γt ) x + K1 (x2 − 1) D p1 [x (t )] + K2 (x2 − 1) D p2 [x (t )] = 0,
K2 ( K2 > 0) are the fractional coefficients of two kinds of fractionalorder terms. The introduction of fractional-order terms in Mathieu equation lies in some dynamical devices could be modeled by fractional-order derivative, such as viscoelastic and fluid–solid-coupling device. There are several definitions for fractional-order derivative, such as Grünwald–Letnikov, Riemann–Liouville and Caputo definitions [1,2]. Under wide senses, they are equivalent for most mathematical functions. Accordingly, Caputo’s definition is adopted with the form as
∫0
t
x (n)( u)
( t − u) p − n + 1
du , (2)
where Γ ( y ) is Gamma function satisfying Γ ( y + 1) = yΓ ( y ), and the fractional order meets n − 1 < p < n while n is a natural number. Using the following transformations of coordinates
ζ = εμ , δ = ω02, B = εβ, γ = 2ω, K1 = εk1, K2 = εk2 Eq. (1) becomes
x¨ + 2εμẋ + (ω02 + 2εβ cos 2ωt ) x + εk1 (x2 − 1) D p1 [x (t )] + εk2 (x2 − 1) D p2 [x (t )] = 0.
(3)
We shall consider the parametric excitation frequency to be γ = 2ω0 + εσ , where ε < < 1 is a small real parameter, i.e. ω = ω0 + εσ /2. In order to obtain the first approximate solution, it could also be rewritten as ω2 = ω02 + εσω0 . So Eq. (3) becomes
x¨ + ω2x = ε {σω 0 x − 2μẋ − (2β cos 2ωt ) x − k1 (x2 − 1) D p1 [x (t )] − k2 (x2 − 1) D p2 [x (t )]}
(4)
The solution for Eq. (4) could be supposed as
x = a cos φ,
(5a)
ẋ = − aω sin φ,
(5b)
where φ = ωt + θ . Based on the averaging method [38,39], one could obtain the standard equations as
ȧ = −
1 Tω
∫0
1 Tω
aθ ̇ = −
T
∫0
⎡⎣ P1 ( a, θ ) + P2 ( a, θ ) + P3 ( a, θ ) ⎤⎦ sin φdφ,
T
⎡⎣ P1 ( a, θ ) + P2 ( a, θ ) + P3 ( a, θ ) ⎤⎦ cos φdφ,
(6a)
(6b)
where
P1 ( a, θ ) = ε [σω 0 a cos φ + 2μaω sin φ − (2β cos 2ωt ) a cos φ], P2 ( a, θ ) = − εk1 [(a cos φ)2 − 1] D p1 [ a cos φ], P3 ( a, θ ) = − εk2 [(a cos φ)2 − 1] D p2 [ a cos φ]. From the averaging method, one could select the time terminal T as T = 2π when Pi ( a, θ )( i = 1, 2, 3) is periodic function, or T = ∞ when Pi ( a, θ ) ( i = 1, 2, 3) is aperiodic one. Expanding the integral in Eq. (6), one could obtain
ȧ = a1̇ + a2̇ + a3̇ ,
(7a)
aθ ̇ = aθ1̇ + aθ 2̇ + aθ 3̇ ,
(7b)
(1)
where 2ζ , δ , and 2B cos γt are the system linear damping coefficient, constant stiffness coefficient and the periodic time-varying stiffness coefficient respectively. K1 (x2 − 1) D p1 ⎡⎣ x ( t ) ⎤⎦ and K2 (x2 − 1) D p2 ⎡⎣ x ( t ) ⎤⎦ are two kinds of VDP fractional-order terms, where the fractional orders are restricted as 0 < p1 < 1 and 1 < p2 < 2, and K1 ( K1 > 0) and
1 Γ ( n − p)
D p ⎡⎣ x ( t ) ⎤⎦ =
131
The first part in Eq. (7) could be simplified as
a1̇ = −
1 2π ω
∫0
1 2π ω
aθ1̇ = −
2π
∫0
P1 (a, θ ) sin φdφ =
2π
εaβ sin 2θ − εμa, 2ω
P1 (a, θ ) cos φdφ =
εaβ cos 2θ εaσω 0 − . 2ω 2ω
(8a)
(8b)
In order to calculate the second part and the third part for Eq. (7)
a2̇ = − lim
1
T →∞ Tω
= lim
εk1
T →∞ Tω
∫0
∫0 T
T
P2 ( a, θ ) sin φdφ
[(a cos φ)2 − 1] D p1 ⎡⎣ a cos φ⎤⎦ sin φdφ,
(9a)
132
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
T 1 P2 ( a, θ ) cos φdφ T →∞ Tω 0 T εk = lim 1 [(a cos φ)2 − 1] D p1 ⎡⎣ a cos φ⎤⎦ cos φdφ, T →∞ Tω 0
where
∫
aθ 2̇ = − lim
∫
T 1 a3̇ = − lim P3 ( a, θ ) sin φdφ T →∞ Tω 0 T εk = lim 2 [(a cos φ)2 − 1] D p2 ⎡⎣ a cos φ⎤⎦ sin φdφ, T →∞ Tω 0
K (p) = K1ω p1 − 1 cos (9b)
T 1 P3 ( a, θ ) cos φdφ T →∞ Tω 0 T εk = lim 2 [(a cos φ)2 − 1] D p2 ⎡⎣ a cos φ⎤⎦ cos φdφ, T →∞ Tω 0
ȧ = (10a)
∫
εβa a⎛ a2 ⎞ sin 2θ − εμa + ⎜ 1 − ⎟ C (p), 2ω 2⎝ 4⎠
aθ ̇ =
∫
aθ 3̇ = − lim
(10b)
one could expand the integral in Eqs. (9a) and (10a), and get
εk a2 T sin 3φD p1 ⎡⎣ a cos φ⎤⎦ dφ a2̇ = lim 1 T →∞ 4Tω 0 εk a2 T sin φD p1 ⎡⎣ a cos φ⎤⎦ dφ + lim 1 T →∞ 4Tω 0 T εk sin φD p1 ⎡⎣ a cos φ⎤⎦ dφ, − lim 1 T →∞ Tω 0
(11)
εk2 a2 T sin 3φD p2 ⎡⎣ a cos φ⎤⎦ dφ 4Tω 0 εk a2 T sin φD p2 ⎡⎣ a cos φ⎤⎦ dφ + lim 2 T →∞ 4Tω 0 T εk sin φD p2 ⎡⎣ a cos φ⎤⎦ dφ. − lim 2 0 T →∞ Tω
(12)
ȧ =
Ba a⎛ a2 ⎞ sin 2θ − ζa + ⎜ 1 − ⎟ C (p), 2ω 2⎝ 4⎠
aθ ̇ =
∫
∫
a3̇ = lim
T →∞
∫
∫
We introduce some important results about the definite integrals, which had been deduced in Ref. [12,39]
1 T →∞ Tω
∫0
1
∫0
lim
lim
T →∞ Tω
T
T
sin φD p [ a cos φ] dφ = −
cos φD p [ a cos φ] dφ =
⎛ pπ ⎞ aω p − 1 sin ⎜ ⎟, ⎝ 2⎠ 2
⎛ pπ ⎞ aω p − 1 cos ⎜ ⎟, ⎝ 2⎠ 2
(13a)
(13b)
for 0 < p < 1, and
1
T →∞ Tω
∫0
1 T →∞ Tω
∫0
lim
lim
T
T
sin φD p [ a cos φ] dφ = −
⎛p−1 ⎞ aω p − 1 π ⎟, cos ⎜ ⎝ 2 ⎠ 2
(14a)
cos φD p [ a cos φ] dφ = −
⎛p−1 ⎞ aω p − 1 π ⎟, sin ⎜ ⎝ 2 ⎠ 2
(14b)
for 1 < p < 2 respectively. As the similarly procedure, the first parts in Eqs. (11) and (12) would be zero when T → ∞, and the detail deduction procedure is presented in Appendix. One can get the second part and the third part for Eq. (7a)
a2̇ + a3̇ =
a⎛ a2 ⎞ ⎟ C (p), ⎜1 − 2⎝ 4⎠
(15a)
where
C (p) = K1ω p1 − 1 sin
⎛p −1 ⎞ p1π + K2 ω p2 − 1 cos ⎜ 2 π⎟ ⎠ ⎝ 2 2
(15b)
is defined as the additional damping coefficient (ADC). Through the similar procedure, the second part and the third part for Eq. (7b) could also be deduced as
aθ 2̇ + aθ 3̇ =
⎞ a ⎛ 3a2 − 1⎟ K (p), ⎜ ⎠ 2⎝ 4
(16a)
⎞ εβa εσaω 0 a ⎛ 3a2 + ⎜ − 1⎟ K (p). cos 2θ − ⎠ 2ω 2ω 2⎝ 4
(17a)
(17b)
Substituting the original parameters into Eq. (17), we can get
∫
∫
(16b)
is defined as the additional stiffness coefficient (ASC). Substituting Eqs. (8), (15) and (16) into Eq. (7), we obtain
∫
∫
⎛p −1 ⎞ p1π − K2 ω p2 − 1 sin ⎜ 2 π⎟ ⎠ ⎝ 2 2
⎞ εσaω 0 Ba a ⎛ 3a2 + ⎜ − 1⎟ K (p). cos 2θ − ⎠ 2ω 2ω 2⎝ 4
(18a)
(18b)
Eq. (18) shows that all two kinds of VDP fractional-order terms have crucial effects on the damping force and the spring force in Eq. (1), so that the parameters in the two kinds of VDP fractionalorder terms, i.e. the fractional coefficients K1 and K2, the fractional orders p1 and p2, are all important in system dynamics. From the analysis of Eqs. (15b) and (16b), it could be found that the fractional coefficients K1 and K2 are linear with ADC and ASC. ADC and ASC will increase monotonically with the increase of K1. Moreover, ADC will also increase monotonically with the increase of K2, but ASC will decrease monotonically if K2 increases. Accordingly the fractional coefficients K1 and K2 have important effect on the amplitude and the resonance frequency. It should be notable the fractional orders p1 and p2 are more important parameters influencing the system dynamics. It could be found that the fractional orders p1 and p2 affect ADC and ASC by an exponential function and a trigonometric function. From Eqs. (15b) and (16b), it could be found that ADC and ASC are all related with the parametric excitation frequency too. Due to ω0 ¼ 1 in the following parts in this paper, the effect of the parametric excitation frequency could not be considered herein. When K2 ¼0, i.e. there is only the first kind of VDP term in Eq. (1). If p1 varies from 0 to 1, it could be concluded that the increase of p1 would make ADC larger and ASC smaller. When p1 → 0, ADC and ASC will approach 0 and K1 respectively. At that moment the whole system damping coefficient becomes the minimum value, and the whole system stiffness coefficient will be the maximum value. When p1 → 1, ADC will be K1 and ASC will approach 0, so that the whole system damping coefficient will the maximum value and the whole system stiffness coefficient will become the minimum value. When K1 ¼0, i.e. there will exist only the second kind of VDP term in Eq. (1). If p2 varies from 1 to 2, it could be found that the increase of p2 would make ADC and ASC smaller. When p2 → 1, ADC will approach K2 and ASC will be 0. If p2 → 2, ADC will be 0 and ASC will approach −K2. Then the whole system damping coefficient becomes the minimum value and the maximum amplitude is relatively larger. From the analysis of Eq. (16b), it could be found that there is a special case, which is shown in Eq. (19). In this case the ASC is zero, and the whole system stiffness is almost the same as the original linear stiffness.
K2 =
p1π 2 p2 − 1 π 2
K1ω p1 − 1 cos ω p2 − 1 sin
(
)
(19)
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
3. The steady-state solutions and their stability conditions 3.1. The fixed point (0, 0) Eq. (18) shows that a¼ 0 is always one of the system solutions. In order to study its stability condition, one should determine the eigenvalues of the Jacobian matrix of Eq. (18) at a¼ 0, which can be written as
⎡ ∂f ⎢ J = ⎢ ∂a ⎢ ∂g ⎢⎣ ∂a
∂f ⎤ ⎥ ∂θ ⎥ ∂g ⎥ ∂θ ⎥⎦
a=0
(20)
where
⎞ Ba εσaω 0 a ⎛ 3a2 cos 2θ − + ⎜ − 1⎟ K (p). ⎠ 2ω 2ω 2⎝ 4
The eigenvalues of the above Jacobian matrix will be
λ1 = 0, λ2 =
C (p) − ζ. 2
(21)
If 2ζ > C (p), i.e. the linear damping force is greater than ADC, this solution will be stable. From the analysis of Eq. (15b), it could be found that C (p) increases linearly with the increase of K1 and K2. The increase of p1 would make C (p) larger, but the increase of p2 would make C (p) smaller. When p1 → 1 and p2 → 1, two kinds of VDP fractional-order terms can provide the effect of damping force simultaneously, so that C (p) will be the maximum value K1+ K2. In the vicinity of the fixed point (0, 0) a is a small real parameter, so a2
that − 4 C (p) will be about 0, and the damping term of VDP fractional-order terms is negative. That would result into 2ζ < C (p), so that the fixed point (0, 0) will be unstable. When p1 → 0 and p2 → 2, the damping term of two kinds of VDP fractional-order terms is about 0, which will lead to 2ζ > C (p) and the fixed point (0, 0) will be stable. The above analysis is identical with the result in the monograph [40]. 3.2. The existence conditions and stability conditions for the steadystate nonzero solution Now we study the steady-state nonzero solution, which is more important and meaningful in vibration engineering. Letting ȧ = 0 and θ ̇ = 0, and omitting the zero solution, Eq. (18) becomes
⎧ B 1⎛ a¯ 2 ⎞ ⎪ sin 2θ¯ − ζ + ⎜ 1 − ⎟ C (p) = 0 ⎪ 2ω 2⎝ 4⎠ ⎨ , ⎞ ⎪ B εσω 0 1 ⎛ 3a¯ 2 ¯ − + − ( ) = cos 2 K p θ 1 0 ⎜ ⎟ ⎪ ⎠ ω 2⎝ 4 ⎩ 2ω
(22)
where a¯ and θ¯ are the amplitude and phase of the steady-state response respectively. After eliminating θ¯ from Eq. (22), one could obtain the amplitude-frequency equation
⎡ εσω 0 ⎤2 B2 [(1 − ρ) C (p) − 2ζ ]2 + ⎢ (3ρ − 1) K (p) − ⎥⎦ = 2 , ⎣ ω ω where ρ =
a¯ 2 . 4
where
(23)
⎡ εσω0 ⎤2 B2 H3 = [2ζ − C (p)]2 + ⎢ K (p) + ⎥− . ⎣ ω ⎦ ω2 And the phase-frequency equation is
⎤ ⎡ (1 − ρ) C (p) − 2ζ ⎥ 1 θ¯ = arctan ⎢ εσω 0 . 2 ⎢⎣ (3ρ − 1) K (p) − ω ⎥⎦
(25)
Due to H1 > 0 , the existence conditions for Eq. (24) solution could be summarized into two cases. The first case is that Eq. (24) has only one real solution when H3 o 0. The other case means Eq. H22 4H1
.
From the analysis of Eq. (24), it could be found that H1, H2 and H3 are related to ADC and ASC. Moreover, H2 and H3 are related to the system linear damping coefficient, the periodic time-varying coefficient and parametric excitation frequency. Accordingly, the range of nonzero solutions has complex interaction with two kinds of VDP fractional-order terms and the system parameters. Now, we select the system parameters as ζ ¼0.1, ε ¼0.01, B ¼0.25, and study the effects of the two kinds of VDP fractionalorder parameters K1, K2, p1 and p2 on the existence range of nonzero solutions. The results are shown in Fig. 1(a)–(d), where the X-axis is for parametric excitation frequency ω , the Y-axis denote for K1, K2, p1, p2 respectively, and the Z-axis is for the existence boundary of the solution for Eq. (24). Obviously the parts of greater than zero are for the existence range of nonzero solutions. From the analysis of Fig. 1, it could be found that the increase of K1 and p2 would make the existence range of nonzero solutions larger, but the increase of p1 would make the existence range of nonzero solutions smaller. Furthermore, the increase of K2 would make the existence range of nonzero solutions larger at first, and then smaller if K2 is big enough. Next we study the stability of the steady-state nonzero solution. Letting a = a¯ + Δa , θ = θ¯ + Δθ , and substituting them into Eq. (18), one can obtain
¯ cos 2θ¯ dΔa ⎡ B sin 2θ¯ C (p) ⎛ aB 3a¯ 2 ⎞ ⎤ =⎢ −ζ+ Δθ , ⎜1 − ⎟ ⎥ Δa + ⎣ dt 2ω 2 ⎝ ω 4 ⎠⎦ (26a) ¯ sin 2θ¯ dΔθ 3K (p) a¯ aB = Δa + Δθ . dt 4 ω
(26b)
Based on Eq. (22), one could eliminate θ¯ from Eq. (26) and get the characteristic determinant as
⎡ ⎢ A1 − λ det ⎢ ⎢ 3a¯ K (p) ⎣ 4
A1 = −
εσω0 a¯ ω
⎤ − 1 K (p ) ⎥ ⎥ = 0,where ⎥ A2 − λ ⎦
− a¯
(
3a¯ 2 4
)
⎛ a¯ 2 a¯ 2 ⎞ C (p), A2 = − 2ζ + ⎜ 1 − ⎟ C (p). ⎝ 4 4⎠
Expanding the characteristic determinant and substituting
ρ=
a¯ 2 4
into it, we could get the characteristic equation as
⎡ εσω 0 ⎤ λ2 − (A1 + A2 ) λ + A1A2 − 3ρK (p) ⎢ − (3ρ − 1) K (p) ⎥ = 0. ⎣ ω ⎦
(27)
Based on Eq. (27), the stability conditions for the steady-state solution could be established as
Expanding Eq. (23), we obtain
H1ρ2 + H2 ρ + H3 = 0,
⎡ εσω0 ⎤ ⎥, H1 = C (p)2 + 9K (p)2 , H2 = 2C (p)[2ζ − C (p)] − 6K (p)⎢ K (p) + ⎢⎣ ω ⎥⎦
(24) will have two real solutions when H2 o0 and 0 o H3 o
Ba a⎛ a2 ⎞ f (a, θ ) = ȧ = sin 2θ − ζa + ⎜ 1 − ⎟ C (p) 2ω 2⎝ 4⎠ ⎡ ⎤ Ba 1 a2 = sin 2θ − a ⎢ ζ − C (p) + C (p) ⎥, 2ω 2 4 ⎣ ⎦ g (a, θ ) = aθ ̇ =
133
(24)
A1 + A2 < 0, and
(28a)
134
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
Fig. 1. The effects of the fractional parameters on the existence range of nonzero solutions: (a) the effects of K1 ( K2 ¼ 0, p1 ¼0.5), (b) the effects of p1 ( K1 ¼0.1, K2 ¼ 0), (c) the effects of K2 ( K1 ¼0, p2 ¼1.5) and (d) the effects of p2 ( K1 ¼0, K2 ¼ 0.2).
⎡ εσω 0 ⎤ − (3ρ − 1) K (p) ⎥ > 0. A1A2 − 3ρK (p) ⎢ ⎣ ω ⎦
(28b)
After expanding Eq. (28), the equivalent forms of the stability conditions for the steady-state solution are
T1 = 2ρC (p) + [2ζ − C (p)] > 0,
condition T1 40. Then one should analyze whether the condition T2 40 holds based on the existence conditions of Eq. (24). The results are summarized as follows: a. when H3 o0, Eq.(24) has one real solution
(29a)
ρ=
and
⎧ ⎡ εσω 0 ⎤ ⎫ T2 = [C (p)2 + 9K (p)2] ρ2 + ⎨ C (p)[2ζ − C (p)] − 3K (p) ⎢ K (p) + ⎥⎬ ⎣ ⎩ ω ⎦⎭ ρ > 0. The stability conditions are equal to
T2 = 2H1ρ + H2 > 0,
H22 − 4H1H3
−H2 +
2H1
.
(30)
Substituting Eq. (30) into Eq. (29b), one could obtain
T2 = H22 − 4H1H3 40. This will be the fourth cases, so that the solution is stable. H2 b. when H2 o0 and 0 o H3 o 2 , Eq. (24) has two real solutions 4H1
(29b)
where T1 may be defined as the equivalent nonlinear damping coefficient and T2 is the nonlinear parameter of the stability condition. T1 and T2 are generally used to judge the stability for the steadystate nonzero solution. It can be found that the stability conditions for the steady-state solution could be divided into four cases, i. e., 1) T1 o0 and T2 o0, 2) T1 o0 and T2 40, 3) T1 40 and T2 o0, 4) T1 40 and T2 4 0. When T1 o 0, the system response will be divergent and the nonzero steady-state solutions are unstable. When T1 4 0 and T2 o0, the system response will be unstable periodic vibration because T2 is affected by the stiffness and damping at the same time. If T1 40 and T2 40, the steady-state solution will be stable. This implies that only the solution corresponding to the forth case is stable. Now, the research on stability conditions could be focused on the two parameters. At first one should select appropriate system parameter ζ so as to satisfy the
ρ1 =
−H2 +
H22 − 4H1H3 2H1
and ρ2 =
−H2 −
H22 − 4H1H3 2H1
.
Substituting ρ1 into Eq. (29b), we could obtain T2 40, so that this solution is stable. On the contrary, substituting ρ2 into Eq. (29b) will result into T2 o0, so that this solution will unstable.
4. Numerical simulations 4.1. Comparison between the approximate analytical solutions and the numerical results An illustrative example system is studied herein as defined by the basic system parameter, p1 ¼0.5, p2 ¼ 1.5, ζ ¼ 0.1, ε ¼ 0.01 and B ¼0.25. Three group of fractional-order parameters, i.e., K1 ¼ 0.1 and K2 ¼0.1, K1 ¼ 0.2 and K2 ¼0.4, K1 ¼0.2 and K2 ¼0.5, are selected to study the different response modes. Based on the amplitude–
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
Fig. 4. The amplitude–frequency curves ( K1 ¼ 0.2, K2 ¼ 0.5).
Fig. 2. The amplitude–frequency curves ( K1 ¼ 0.1 , K2 ¼ 0.1).
Fig. 5. The effects of the fractional coefficient K1 on the amplitude–frequency curves ( p1 ¼ 0.5, K2 ¼ 0).
Fig. 3. The amplitude–frequency curves ( K1 ¼0.2, K2 ¼ 0.4).
frequency equation Eq. (24), one could analytically obtain the amplitude–frequency curves shown in Figs. 2–4, where the solid line is for the asymptotically stable solution and the dash line denotes the unstable one. In the three figures, the parametric excitation frequency parameter ω is for the horizontal axis, and the steady-state amplitude a is for the vertical axis. In order to verify the precision of the analytical solutions, we also present the numerical results. The numerical method is [1,2]
D p ⎡⎣ x ( tl ) ⎤⎦ ≈ h−p
l
∑ C jp x ( tl − j ), j=0
(31)
where tl = lh is the time sample points, h is the sample step and C jp is the fractional binomial coefficient with the iterative relationship as
⎛ 1 + p⎞ p C0p = 1, C jp = ⎜ 1 − ⎟C . ⎝ j ⎠ j−1
135
(32)
Here we select h ¼ 0.001, and the total computation time is 300 s. The amplitude–frequency curves by numerical integration are also shown in Figs. 2–4, where those circles are for the numerical results. From the observation of Figs. 2–4, one could conclude that the approximately analytical solutions agree very well
with the numerical results and could present satisfactory precision even if the amplitude–frequency curves are different topological structures. 4.2. Effects of VDP fractional-order terms on system dynamics When the coefficient K1 of the first kind of VDP fractional-order term is changed, the amplitude–frequency curves are shown in Fig. 5, where K2 ¼0. From the observation of Fig. 5, it could be found that the larger the fractional coefficient K1, the smaller the maximum amplitude. This is due to the fact that ADC becomes larger in this procedure. Furthermore, it could also be found that the increase of the fractional coefficient K1 would make ASC larger, which could result into the rightwards bending of the backbone of the amplitude–frequency curves, and make the existence range of the steady-state solution larger. When the order p1 of the first kind of VDP fractional-order term varies from 0 to 1, the amplitude–frequency curves are shown in Fig. 6, where K2 ¼0. From the observation of Fig. 6, it could be found that the larger the fractional order p1, the smaller the maximum amplitude. The reason is that ADC will also become larger in this procedure. Furthermore, it could also be found that the increase of the fractional coefficient p1 would make ASC smaller, which could result into the leftwards bending of the
136
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
Fig. 6. the effects of the fractional order p1 on the amplitude–frequency curves ( K1 ¼ 0.1, K2 ¼0).
Fig. 8. (a) the effects of the fractional order p2 on the amplitude–frequency curves ( K1 ¼0, K2 ¼ 0.2). (b) the effects of the fractional order p2 on the amplitude–frequency curves ( K1 ¼ 0.1, p1 ¼0.5, K2 ¼ 0.2).
Fig. 7. (a) The effects of the fractional coefficient K2 on the amplitude–frequency curves ( K1 ¼0, p2 ¼1.5). (b) The effects of the fractional coefficient K2 on the amplitude–frequency curves ( K1 ¼0.2, p1 ¼ 0.5, p2 ¼1.5).
backbone of the amplitude–frequency curves, and make the existence range of the steady-state nonzero solution smaller. Moreover, the changes of order p1 will also change the topology structures of the amplitude–frequency curves. The role of the second kind of VDP fractional-order term is similar to that of the first kind of VDP fractional-order term. The detail discussions about the effects of the second kind of VDP fractional-order term are analyzed in the following parts, and the results are shown in Figs. 7 to 8. In Fig. 7 the subfigure 7(a) and 7(b) are for K1 ¼0 and K1 ¼0.2 respectively. From the observation of those curves in Fig. 7, it could be found that the increase of the fractional coefficient K2 would make ADC larger, which will make the maximum amplitude smaller. Simultaneously, it is easy to find that the increase of the fractional coefficient K2 will make ASC smaller, which could result into the leftwards bending of the backbone of the amplitude–frequency curves, and make the existence range of the steady-state solution larger at first and then smaller. The changes of the topological structures could also be found in Fig. 7. Additionally, the difference between Figs. 7(a) and 7(b) is the existence of the first kind of the VDP fractional-order term, which decreases the maximum amplitude, and diminishes the leftwards bending of the backbone of the amplitude–frequency curves. It will affect the existence range of the steady-state solution and the topological
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
structure at the same time. When the order p2 of the second kind of VDP fractional-order term varies from 1 to 2, the results are shown in Fig. 8, where Fig. 8 (a) is for K1 ¼ 0 and Fig. 8(b) is for K1 ¼0.1 respectively. From the observation of Fig. 8, it could be found that the larger the fractional order p2, the larger the maximum amplitude. This is due to the fact that ADC becomes smaller in this procedure. The existence range of the steady-state solution becomes larger too. Moreover, the topological structures is changed, and the backbone of the amplitude–frequency curves will bend towards the left because ASC becomes smaller in this procedure according to Eq. (16b). When p2 → 2, one could conclude that the maximum amplitude is too large, and the backbone of the amplitude–frequency curve bends towards the left. When p2 → 1, it could be found that the maximum amplitude is small, and the bending degree of the backbone of the amplitude–frequency curve is little. Similar to Fig. 7, the difference between Fig. 8(a) and (b) is the existence of the first kind of the VDP fractional-order term. When p2 is the same value in Fig. 8(a) and (b), one could conclude that the existence of the first kind of the VDP fractional-order term could also decrease the maximum amplitude, and diminish the leftwards bending of the backbone of the amplitude–frequency curves. It will affect the existence range of the steady-state solution and the topological structure at the same time. From the observation of Fig. 8(b), it could be found that the maximum amplitude is smaller when two kinds VDP fractional-order terms are exist. The backbone of the amplitude–frequency curve bends from the right to the left, and the existence range of the steady-state solution is smaller at first and then larger. Moreover, it could be found that there is a special case. In this case the ASC is zero, and the ADC is almost the maximum, so that the maximum amplitude will be reduced significantly. This special case may be greatly important in vibration control engineering. Based on the above analysis, one could select appropriate fractional parameters to improve the dynamical characteristics of Mathieu equation with two kinds of VDP fractional-order terms. Comparing with the integral-order counterpart, the response modes of the Mathieu equation with VDP fractional-order terms are more complicated. Especially, two kinds of VDP fractional-order terms will bring better control performance than the single first kind of VDP fractional-order term.
5. Conclusions The Mathieu equation with two kinds of VDP fractional-order terms is studied based on the averaging method, and the approximately analytical solution and the amplitude–frequency equation are obtained. The steady-state solution and stability condition of the steady-state solution are researched. The effects of two kinds of VDP fractional coefficients and fractional orders on the steady-state solution are analyzed, and it is found that these two kinds of VDP fractional-order terms affect the linear damping and the linear stiffness, which are characterized by ADC and ASC. The effects of the two kinds of VDP fractional-order parameters on the existence range of nonzero solutions and the stability of the steady-state nonzero solution are studied. The comparisons of the analytical solution with the numerical one verify the correctness of the analytical solution. When two kinds of VDP fractional-order term exist, it would be found that there is a special case. In this case the ASC is almost zero, and the ADC is almost the maximum, so that the maximum amplitude will become almost the minimum. The special case may be greatly important in vibration control engineering. The results are meaningful to analyze and control the dynamical behavior of the similar fractional-order system, and could present useful reference.
137
Acknowledgment The authors are grateful to the support by National Natural Science Foundation of China (No. 11372198), the Cultivation plan for Innovation team and leading talent in Colleges and universities of Hebei Province (LJRC018), the Program for advanced talent in the universities of Hebei Province (GCC2014053), the Program for advanced talent in Hebei (A201401001), and the education department project of Hebei Province (QN2016258)
Appendix A Taking the first parts in Eq. (11) as an example, According to Caputo’s definition in Eq. (2), and ϕ = ωt + θ , the first parts in Eq. (11) becomes lim
T →∞
εk1a 2 4Tω
∫0
εk1a 2 T →∞ 4Tω
= lim
T
sin 3 (ωt + θ ) D p1 [ a cos (ωt + θ ) ] dφ
∫0
T
⎡ 1 sin 3 (ωt + θ ) ⎢ ⎢⎣ Γ ( 1 − p)
∫0
t
⎤ −aω sin (ωu + θ ) ⎥ du dt . p ( t − u) ⎦⎥
(A.1)
Letting s = t − u and du = − ds , Eq. (A.1) becomes lim
T →∞
εk1a 2 4Tω
∫0
T
⎤ −aω sin (ωu + θ ) ⎥ du dt p ⎥⎦ t − u ( )
⎡ 1 sin 3 (ωt + θ ) ⎢ ⎢⎣ Γ ( 1 − p)
∫0
T
⎡ sin 3 (ωt + θ ) ⎢ ⎣
∫0
T
⎡ ⎢ sin 3 (ωt + θ ) sin (ωt + θ ) ⎣
∫0
T
⎡ ⎢ sin 3 (ωt + θ ) cos (ωt + θ ) ⎣
∫0
=
1 −εk1a 3 lim 4Γ ( 1 − p) T →∞ T
∫0
=
1 −εk1a 3 lim 4Γ ( 1 − p) T →∞ T
∫0
+
1 εk1a 3 lim 4Γ ( 1 − p) T →∞ T
∫0
t
t
sin (ωt + θ − ωs) ⎤ ds⎥ dt ⎦ sp t
cos ωs ⎤ ds⎥ dt ⎦ sp
τ
sin ωs ⎤ ds⎥ dt . ⎦ sp
(A.2)
Defining the first part as A1 and the second part as A2 in Eq. (A.2), and integrating it by parts respectively, one could obtain −εk1a3 1 lim 4Γ ( 1 − p ) T →∞ T
A1 =
=
=
⎧ 1 ⎨ − [ cos 4 (ωt + θ ) − cos 2 (ωt + θ ) ] ⎩ 2
T
εk1a3 1 ⎡ sin 4 (ωt + θ ) sin 2 (ωt + θ ) ⎤ − lim ⎢ ⎥ ⎦ 4ω 2ω 8Γ ( 1 − p ) T →∞ T ⎣ −
A2 =
∫0
εk1a3 1 lim 8Γ ( 1 − p ) T →∞ T
1 εk1a 3 lim 4Γ ( 1 − p) T →∞ T
∫0
∫0
T
T
∫0
1 εk1a 3 lim 8Γ ( 1 − p) T →∞ T
cos ωs ⎫ ds⎬ dt ⎭ sp T
cos ωs ds sp 0
⎧1 ⎨ ⎡⎣ sin 4 (ωt + θ ) + sin 2 (ωt + θ ) ⎤⎦ ⎩2
∫0
t
⎡ sin 4 (ωt + θ ) τ sin 2 (ωt + θ ) ⎤ cos ωt − dt , ⎢ ⎥ ⎣ ⎦ tp 4ω 2ω
1 ⎡ cos 4 (ωt + θ ) cos 2 (ωt + θ ) ⎤ −εk1a 3 lim ⎢ − ⎥ ⎦ 4ω 2ω 8Γ ( 1 − p) T →∞ T ⎣ +
t
∫0
T
∫0
t
∫0
τ
(A.3)
sin ωs ⎫ ds⎬ dt ⎭ sp T
sin ωs ds sp 0
⎡ cos 4 (ωt + θ ) τ cos 2 (ωt + θ ) ⎤ sin ωt dt . − ⎥ ⎢ ⎦ tp ⎣ 4ω 2ω
(A.4)
For simplicity, one could define the first and second part in Eq. (A.3) as BA11 and BA12, the first and second part in Eq. (A.4) as BA21 and BA22 respectively. We introduce two important formulae, which had been deduced in Refs. [11,12]
∫ T →∞ 0 lim
lim
T →∞
∫0
T
T
⎛ pπ ⎞ sin (ωt ) dt = ω p − 1Γ (1 − p) cos ⎜ ⎟ ⎝ 2⎠ tp
(A.5a)
⎛ pπ ⎞ cos (ωt ) dt = ω p − 1Γ (1 − p) sin ⎜ ⎟ ⎝ 2⎠ tp
(A.5b)
Substituting Eq. (A.5b) into BA11 and integrating BA12, we obtain BA11 =
εk1a 3 1 ⎡ sin 4 (ωT + θ ) sin 2 (ωT + θ ) ⎤ − lim ⎢ ⎥ ⎦ 8Γ ( 1 − p) T →∞ T ⎣ 4ω 2ω pπ
=
εk1a 3ω p − 1 sin ( 2 ) 8
∫0
T
cos ωs ds sp
1 ⎡ sin 4 (ωT + θ ) sin 2 (ωT + θ ) ⎤ − lim ⎢ ⎥ = 0, ⎣ ⎦ 4ω 2ω
T →∞ T
(A.6a)
138
S. Wen et al. / International Journal of Non-Linear Mechanics 84 (2016) 130–138
BA12 = −
εk1a3 sin 2 (ωt + θ ) ⎤ cos ωt 1 T ⎡ sin 4 (ωt + θ ) τ − lim ∫ ⎢ dt ⎥ ⎦ tp 4ω 2ω 8Γ ( 1 − p) T →∞ T 0 ⎣
= 0.
(A.6b)
Similarly, BA21 and BA22 could also be deduced BA21 ¼ 0, BA22 ¼0. So, the first parts in Eq. (11) would be zero when T → ∞. As the similar procedure, one could find that the first parts in Eq. (12) would be zero too when T → ∞.
References [1] I. Petras, Fractional-order Nonlinear Systems: Modeling, Analysis and Simulation, Higher Education Press, Beijing, 2011. [2] R. Caponetto, G. Dongola, L. Fortuna, I. Petras, Fractional Order Systems: Modeling and Control Applications, World Scientific, New Jersey, 2010. [3] R. Gorenflo, E.A. Abdel-Rehim, Convergence of the Grünwald–Letnikov scheme for time-fractional diffusion, Journal. Comput. Appl. Math. 205 (2) (2007) 871–881. [4] G. Jumarie, Modified Riemann–Liouville derivative and fractional Taylor series of nondifferentiable functions further results, Comput. Math. Appl. 51 (2006) 1367–1376. [5] M. Ishteva, R. Scherer, L. Boyadjiev, On the Caputo operator of fractional calculus and C-Laguerre functions, Math. Sci. Res. J. 9 (6) (2005) 161–170. [6] Agnieszka B. Malinowska, Delfim F.M. Torres, Fractional calculus of variations for a combined Caputo derivative, Fract. Calc. Appl. Anal. 14 (4) (2011) 523–537. [7] Z.H. Wang, H.Y. Hu, Stability of a linear oscillator with damping force of fractional-order derivative, Sci. China G: Phys. Mech. Astron. 39 (10) (2009) 1495–1502. [8] Z.H. Wang, M.L. Du, Asymptotical behavior of the solution of a SDOF linear fractionally damped vibration system, Shock. Vib. 18 (1–2) (2011) 257–268. [9] Y.J. Shen, P. Wei, S.P. Yang, Primary resonance of fractional-order van der Pol oscillator, Nonlinear Dyn. 77 (2014) 1629–1642. [10] Y.J. Shen, P. Wei, C.Y. Sui, S.P. Yang, Subharmonic resonance of Van Der Pol oscillator with fractional-order derivative, Math. Probl. Eng. (2014) 7380871–738087-17. [11] Y.J. Shen, S.P. Yang, H.J. Xing, et al., Primary resonance of Duffing oscillator with fractional-order derivative, Commun. Nonlinear Sci. Numer. Simul. 17 (7) (2012) 3092–3100. [12] Y.J. Shen, S.P. Yang, H.J. Xing, et al., Primary resonance of Duffing oscillator with two kinds of fractional-order derivatives, Int. J. Non-Linear Mech. 47 (9) (2012) 975–983. [13] C.P. Li, W.H. Deng, Remarks on fractional derivatives, Appl. Math. Comput. 187 (2) (2007) 777–784. [14] J.X. Cao, H.F. Ding, C.P. Li, Implicit difference schemes for fractional diffusion equations, Commun. Appl. Math. Comput. 27 (1) (2013) 61–74. [15] F.H. Zeng, C.P. Li, High-order finite difference methods for time-fractional subdiffusion equation, Chin. Journal. Comput. Phys. 30 (4) (2013) 491–500. [16] P. Wahi, A. Chatterjee, Averaging oscillations with small fractional damping and delayed terms, Nonlinear Dyn. 38 (1-4) (2004) 3–22. [17] Y. Xu, Y.G. Li, D. Liu, W.T. Jia, H. Huang, Responses of Duffing oscillator with fractional damping and random phase, Nonlinear Dyn. 74 (3) (2013) 745–753.
[18] L.C. Chen, W.Q. Zhu, Stochastic jump and bifurcation of Duffing oscillator with fractional derivative damping under combined harmonic and white noise excitations, Int. J. Non-linear Mech. 46 (12) (2011) 1324–1329. [19] L.C. Chen, H.F. Li, Z.S. Li, W.Q. Zhu, First passage failure of single-degree-offreedom nonlinear oscillators with fractional derivative, J. Vib. Control. 19 (14) (2013) 2154–2163. [20] L.C. Chen, W.H. Wang, Z.S. Li, W.Q. Zhu, Stationary response of Duffing oscillator with hardening stiffness and fractional derivative, Int. J. Non-linear Mech. 48 (1) (2013) 44–50. [21] M. Caputo, M. Fabrizio, A new definition of fractional derivative without singular kernel, Progr. Fract. Differ. Appl. 1 (2) (2015) 1–13. [22] A. Atangana, D. Baleanu, New fractional derivatives with nonlocal and nonsingular kernel: theory and application to heat transfer model, Therm. Sci. (2016), http://dx.doi.org/10.2298/TSCI160111018A. [23] A. Atanganaa, I. Kocab, On the new fractional derivative and application to nonlinear Baggs and Freedman model, J. Nonlinear Sci. Appl. 9 (2016) 2467–2480. [24] A. Atangana, I. Koca, Chaos in a simple nonlinear system with Atangana–Baleanu derivatives with fractional order, Chaos, Solitons Fractals (2016) 10, http://dx.doi.org/10.1016/j.chaos.2016.02.012. [25] A. Atangana, B.S.T. Alkahtani, Analysis of the Keller–Segel model with a fractional derivative without singular kernel, Entropy 17 (6) (2015) 4439–4453. [26] A. Atangana, On the new fractional derivative and application to nonlinear Fisher's reaction–diffusion equation, Appl. Math. Comput. 273 (2016) 948–956. [27] A. Ebaid, D.M.M. ElSayed, M.D. Aljoufi, Fractional calculus model for damped mathieu equation: approximate analytical solution, Appl. Math. Sci. 6 (82) (2012) 4075–4080. [28] R.H. Rand, S.M. Sah, M.K. Suchorsky, Fractional Mathieu equation, Commun. Nonlinear Sci. Numer. Simul. 15 (11) (2010) 3254–3262. [29] A.Y.T. Leung, Z.J. Guo, H.X. Yang, Transition curves and bifurcations of a class of fractional Mathieu-type equations, Int. Journal. Bifurc. Chaos 22 (2012) 1250275-1–1250275-13. [30] A.Y.T. Leung, H.X. Yang, Z.J. Guo, The residue harmonic balance for fractional order van der Pol like oscillators, Journal. Sound. Vib. 331 (2012) 1115–1126. [31] T. Sardar, R.S. Saha, R.K. Bera, B.B. Biswas, The analytical approximate solution of the multi-term fractionally damped van der Pol equation, Phys. Scr. 80 (2009) 1–6. [32] F. Xie, X.Y. Lin, Asymptotic solution of the van der Pol oscillator with small fractional damping, Phys. Scr. 136 (2009) 014033-1-3-4. [33] M. Momeni, I. Kourakis, M. Moslehi-Fard, P.K. Shukla, A Van der Pol-Mathieu equation for the dynamics of dust grain charge in dusty plasmas, J. Phys. A: Math. Theor. 40 (2007) 473–481. [34] M. Pandey, R.H. Rand, A.T. Zehnder, Frequency locking in a forced Mathieu-van der Pol-Duffing system, Nonlinear Dyn. 54 (1–2) (2008) 3–12. [35] J. Kalas, Z. Kadeřábek, Periodic solutions of a generalized Van der Pol–Mathieu differential equation, Appl. Math. Comput. 234 (2014) 192–202. [36] M. Belhaq, A. Fahsi, 2:1 and 1:1 frequency-locking in fast excited van der Pol– Mathieu–Duffing oscillator, Nonlinear Dyn. 53 (1–2) (2008) 139–152. [37] F. Veerman, F. Verhulst, Quasiperiodic phenomena in the van der Pol-Mathieu equation, J. Sound. Vib. 326 (1–2) (2009) 314–320. [38] J.A. Sanders, F. Verhulst, J. Murdock, Averaging Methods in Nonlinear Dynamical Systems, 2nd ed, Springer-Verlag, New York, 2007. [39] Y.J. Shen, S.P. Yang, C.Y. Sui, Analysis on limit cycle of fractional-order van der Pol oscillator, Chaos, Solitons Fractals 67 (2014) 94–102. [40] A.H. Nayfeh, Nonlinear Oscillations, Wiley, New York, 1973.