- Email: [email protected]
The residue harmonic balance for fractional order van der Pol like oscillators
Journal of Sound and Vibration 331 (2012) 1115–1126
Contents lists available at SciVerse ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
The residue harmonic balance for fractional order van der Pol like oscillators A.Y.T. Leung n, H.X. Yang, Z.J. Guo Department of Civil and Architectural Engineering, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong
a r t i c l e i n f o
abstract
Article history: Received 14 July 2011 Accepted 17 October 2011 Handling Editor: L.N. Virgin Available online 16 November 2011
When seeking a solution in series form, the number of terms needed to satisfy some preset requirements is unknown in the beginning. An iterative formulation is proposed so that when an approximation is available, the number of effective terms can be doubled in one iteration by solving a set of linear equations. This is a new extension of the Newton iteration in solving nonlinear algebraic equations to solving nonlinear differential equations by series. When Fourier series is employed, the method is called the residue harmonic balance. In this paper, the fractional order van der Pol oscillator with fractional restoring and damping forces is considered. The residue harmonic balance method is used for generating the higher-order approximations to the angular frequency and the period solutions of above mentioned fractional oscillator. The highly accurate solutions to angular frequency and limit cycle of the fractional order van der Pol equations are obtained analytically. The results that are obtained reveal that the proposed method is very effective for obtaining asymptotic solutions of autonomous nonlinear oscillation systems containing fractional derivatives. The influence of the fractional order on the geometry of the limit cycle is investigated for the first time. & 2011 Elsevier Ltd. All rights reserved.
1. Introduction Since the original work of Balthazar van der Pol [1], the van der Pol equation has always been a prototype for systems with self-excited, limited cyclic oscillations. Nowadays, the Van der Pol oscillator as well as its modified versions is still widely studied in various scientific fields, ranging from biology, chemistry and physics to engineering [2–10]. It is given by the following second-order differential equation: € _ þ uðtÞ ¼ 0, uðtÞ eð1u2 ðtÞÞuðtÞ
(1)
for which the elastic restoring force is a linear function of the dependence variable. Here the over-dot denotes differentiation with respect to time variable t, e is a control parameter that reflects the degree of nonlinearity of the system. For e ¼ 0 this is just a simple linear oscillator. For large values of e, which is for e b1, the system exhibits a relaxation oscillation [2]. The generalized van der Pol oscillators, in which the fractional power of dependent variable u or its derivative occurs, have been of great interest for many years [6,7,11–14]. In [10], with the help of the general harmonic balance method,
n
Corresponding author. Tel.: þ852 34437600. E-mail addresses: [email protected], [email protected] (A.Y.T. Leung).
0022-460X/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jsv.2011.10.023
1116
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
Mickens derived the expression for the limit cycle of following van der Pol type oscillator: € _ þ u1=3 ðtÞ ¼ 0: uðtÞ eð1u2 ðtÞÞuðtÞ
(2)
Mickens [11] also proposed an iteration technique to the oscillator (2) and obtained the same result, but the frequency decreased for about 15% with respect to the classical van der Pol oscillator. To avoid tedious computations in dealing with nonlinear algebras involving in harmonic balance equations for higher-order approximations, Lim and Lai [7] introduced a linearization technique prior to proceeding with harmonic balance procedure and studied (2) under the assumption that 0 r e 5 1. The transient solutions for the first-, second- and third-order approximations are given with the help of averaging method [13] upon taking into account the weak damping effect for (2). Recently, the applications of fractional calculus [15,16] have grown rapidly, which has attracted fairly broad research activity, including the viscoelasticity [17,18], chaotic dynamics [19,20], mechanics of non-Hamiltonian systems [21,22], quantum mechanics [23], physical kinetics [24] and so on [25]. Due to the occurrence of fractional oscillators in dynamical applications, various kinds of van der Pol equations containing fractional derivatives have attracted more and more attention, see for instance [19,20,26–37] and references therein. Two important physical parameters, i.e., the angular frequency and the amplitude of oscillation are involved in the investigations of nonlinear oscillators. Many valuable contributions were made to calculate these two physical parameters. However, the solution of the differential equation containing a fractional derivative is very involved. An effective and easyto-use method for solving such equations is needed. Combining the harmonic balance method with parameter expansion technique, an analytical technique, named residue harmonic balance method, is recently constructed [34–36]. The primary idea of this method is that a bookkeeping parameter p is introduced into the residue of the Fourier truncation series in harmonic equations. The unbalanced residues due to Fourier truncation are considered iteratively by solving linear algebraic equations to improve the accuracy of the solutions successively. Using of residue harmonic balance method, this paper is devoted to investigating the approximations of the fractional van der Pol oscillator with fractional restoring and damping forces. The highly accurate solutions to the angular frequency and limit cycle are obtained. It is shown that the proposed technique is very effective for obtaining asymptotic solutions for nonlinear fractional order systems. 2. Fractional van der Pol oscillator Recently, study on dynamical behaviors of the fractional order systems has attracted increasing attention. In these studies, it has been found that the fractional order systems like integer order systems can generate oscillations. Pereira et al. considered a fractional version of the van der Pol equation given by _ þ uðtÞ ¼ 0, Dlt uðtÞ þ e½uðtÞ2 1uðtÞ
1 o l o 2,
which is obtained by substituting the capacitance by a fractance in the nonlinear RLC circuit model [19,37]. Xie and Lin [28] introduced the following Van der Pol oscillator with a small fractional damping term and studied the asymptotic solution of Cauchy problem: € uðtÞ e½1uðtÞ2 Dlt uðtÞ þ uðtÞ ¼ 0, _ where uð0Þ ¼ a, uð0Þ ¼ b,0 o l o1. Barbosa et al. [26] also suggested the introduction of a fractional order time derivative for Eq. (1) of the standard van der Pol oscillator in the form: D1t þ l uðtÞe½1uðtÞ2 Dlt uðtÞ þ uðtÞ ¼ 0,
0 o l o 1:
(3)
In [35], the oscillatory region and asymptotic solutions of Eq. (3) are discussed by the residue harmonic balance method. In this study, motivated by Mickens [6], a modified version of fractional van der Pol oscillator is considered: € uðtÞ e½1uðtÞ2 Dlt uðtÞ þuðtÞ1=3 ¼ 0,
(4)
where u is a displacement, e is a control parameter that reflects the degree of nonlinearity of the system and e 40, 0 o l o1 is the order of the fractional derivative in the Caputo sense defined by Z t m 1 d Dat uðtÞ ¼ ðttÞma1 m uðtÞdt, (5) GðmaÞ 0 dt for m1 o a rm, m 2 N, t 40, u 2 C m 1 : When l ¼1, Eq. (4) is reduced to the modified van der Pol oscillator (2) with fractional restoring force [6,7]. 3. Residue harmonic balance method The principle of the residue harmonic balance method to differential equations analogous to the Newton iteration for nonlinear algebraic equations is illustrated in Appendix A using power series. To solve Eq. (4) using the residue harmonic
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
1117
balance method, we introduce a new time scale t ¼ ot, then (4) reads
o2 u00 ðtÞeol ð1u2 ðtÞÞDlt uðtÞ þu1=3 ðtÞ ¼ 0,
(6)
where the prime denotes differentiation of u with respect to t, and o is an unknown parameter to be determined. Due to the symmetric characteristic of Eq. (6), u(t) can be expressed by uðtÞ ¼ a0 cosðtÞ þ
1 X
ða2k þ 1 cos½ð2k þ 1Þt þ b2k þ 1 sin½ð2k þ1ÞtÞ:
(7)
k¼0
Therefore, instead of the function u(t), the coefficients a0 ,a2k þ 1 , b2k þ 1 ðk ¼ 0,1,2,. . .Þ become the main unknowns to be determined. Definition (5) of the fractional derivatives allows us to obtain the expression of Dbt uðtÞ, ð0 o b o1Þ, where Z t 1 d a cosðksÞ þbk sinðksÞ ds ðtsÞb Dbt ak cosðktÞ þbk sinðktÞ ¼ ds k Gð1bÞ 0 Z t k t b bk cosðktktÞak sinðktktÞ dt, (8) ¼ Gð1bÞ 0 According to the trigonometric identities cosðktktÞ ¼ cos kt cos kt þsin kt sin kt, sinðktktÞ ¼ sin kt cos ktcos kt sin kt, Eq. (8) becomes
" # Z kt Z kt b k b b cosðktÞ u ðbk cos u þ ak sin uÞdu þ sinðktÞ u ðbk sin uak cos uÞdu Eq: ð8Þ9 Gð1bÞ 0 0 Z Z kt ub cos u kt ub sin u b b du þ k ak cosðktÞ þ bk sinðktÞ du ¼ k bk cosðktÞak sinðktÞ Gð1bÞ Gð1bÞ 0 0 b
b
¼ k ðbk I1 þ ak I2 ÞcosðktÞ þ k ðbk I2 ak I1 ÞsinðktÞ,
(9)
where I1 9
Z kt b u cos u du, Gð1bÞ 0
I2 9
Z kt b u sin u du: Gð1bÞ 0
Although these integrals cannot be evaluated in closed form for general interval, we are able to evaluate them using the contour integration, when steady-state response is of interest, where pffiffiffiffi pG½ð1bÞ=2 bp I1 ¼ sin , I1 ¼ b 1 ¼ tlim -1 2 2 Gðb=2ÞGð1bÞ pffiffiffiffi pGð1ðb=2ÞÞ bp I1 ¼ cos : I2 ¼ b 2 ¼ tlim -1 2 2 G½ð1 þ bÞ=2Gð1bÞ 1 In the following calculation, we substitute I1 and I2 in Eq. (9) byI1 1 and I2 , respectively. Substituting Eq. (7) into Eq. (6) and using the Galerkin procedure result in the following harmonic balance equations: Z 2 p Rci ðo,a1 ,a3 ,. . .,a2n þ 1 ,b3 ,. . .,b2n þ 1 Þ ¼ ½ou€ þ aol ðu2 1ÞDlt u þ u1=3 cosð2i þ 1Þt dt ¼ 0,i ¼ 0,1,. . .,n, p 0 Z 2 p 2 Rsi ðo,a1 ,a3 ,. . .,a2n þ 1 ,b3 ,. . .,b2n þ 1 Þ ¼ ½o u€ þ aol ðu2 1ÞDlt u þu1=3 sinð2i þ 1Þt dt ¼ 0,i ¼ 0,1,. . .,n:
p
0
The direct integration leads to ½Kðq, mÞfqg ¼ fRg ¼ 0, where [K(q,m)] is the total stiffness matrix, Consider the initial approximation:
fRg ¼ ½Rc0 ,Rc1 ,. . .,Rcn ,Rs1 ,. . .,Rsn T
(10) T
and fqg ¼ ½o,a1 ,a3 ,. . .,a2n þ 1 ,b1 ,b3 ,. . .,b2n þ 1 .
o ¼ o0 , u0 ðtÞ ¼ a1,0 cosðtÞ:
(11)
We introduce a bookkeeping parameter pA[0,1], and rewrite u(t) and o as u(t,p) and o(p), respectively:
oðpÞ ¼ o0 þ po1 þ p2 o2 þ , uðt,pÞ ¼ u0 ðtÞ þ pu1 ðtÞ þ p2 u2 ðtÞ þ :
(12)
Expanding the unknown functions uk(t) into Fourier series, one has uk ðtÞ ¼ a1,k cosðtÞ þ
k X i¼1
fa2i þ 1,k cos½ð2i þ1Þt þb2i þ 1,k sin½ð2i þ 1Þtg, k ¼ 1,2,. . . :
(13)
1118
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
It can be shown that all these series are convergent at p¼1. Then the final solution will be ( o ¼ o0 þ o1 þ o2 þ ,
(14)
uðtÞ ¼ u0 ðtÞ þu1 ðtÞ þu2 ðtÞ þ : 3.1. The zeroth-order approximation
Just for the convenience of discussion below, we introduce a function F that represents the governing equation of fractional van der Pol oscillator:
Fðuðt,pÞ, oðpÞÞ ¼ : o2 ðpÞu00 ðt,pÞeol ðpÞð1u2 ðt,pÞÞDlt uðt,pÞ þ u1=3 ðt,pÞ,
(15)
where using the Taylor series, we have
1 2=3 1 2=3 1 5=3 u0 u1 p þ u0 u2 u0 u21 p2 þ , 3 3 3 l 1 ol ¼ ðo0 þ o1 p þ o2 p2 þ Þl ¼ ol0 þ lo0l1 o1 p þ l o0l1 o2 þ o0l2 o21 p2 þ : 2
u1=3 ¼ ðu0 þu1 p þu2 p2 þ Þ1=3 ¼ u0 1=3 þ
(16)
Substituting Eqs. (8), (12) and (13) into Eq. (15) and equating the coefficients of p0 lead to
F0 ðu0 , o0 Þ ¼: o20 u000 eol0 ð1u20 ÞDlt u0 þ u1=3 0 ¼ G1,0 cost þ U1,0 sin t þ
X
½G2i þ 1,0 cosð2i þ1Þt þ U2i þ 1,0 sinð2i þ 1Þt:
(17)
iZ1
The residue in Eq. (17) is separated into two parts: the first part conforms with the initial approximation in Eq. (11), the second part in the summation sign for the determination of the next higher approximations, the number of terms of which will depend on the nature of nonlinearity of u. Based on the harmonic balance, the right hand side of Eq. (17) should not contain the secure terms cosðtÞ and sinðtÞ. Letting their coefficients be zero, we have Z 2 p G1,0 ðo0 ,a10 Þ ¼ F cosðtÞdt ¼ 0, p 0 0 Z p 2 U1,0 ðo0 ,a10 Þ ¼ F0 sinðtÞdt ¼ 0: (18)
p
0
Then the unknown coefficients o0, a1,0 for the initial approximation result by simply solving the nonlinear algebraic equation (18), then we have uð0Þ ðtÞ ¼ a1,0 cosðtÞ, oð0Þ ¼ o0 :
(19)
However, the second part of residue under the summation sign of Eq. (17) cannot be identically zero and should be considered in the first approximation procedure. We rewrite it as X X i F0 ¼ G3,0 cosð3tÞ þ U3,0 sinð3tÞ þ ½G2i þ 1,0 cosð2i þ 1Þt þ U2i þ 1,0 sinð2i þ 1Þtpi1 ¼ R10 þ R0 pi1 : (20) iZ2
iZ2
3.2. The first-order approximation Substituting Eqs. (12) and (20) into Eq. (15) and equating the coefficients of p1 lead to h i F1 ðu1 , o1 Þ ¼ : o0 2 u001 e ol0 ð1u0 2 ÞDlt u1 2ol0 u0 u1 Dlt u0 þ lo0l1 o1 ð1u0 2 ÞDlt u0 þ 2o0 o1 u000 þ
1 X 1 2=3 u0 u1 þ R10 ¼ ½G2i þ 1,1 cosð2i þ 1Þt þ U2i þ 1,1 sinð2iþ 1Þt 3 i¼0 X ½G2i þ 1,1 cosð2i þ 1Þt þ U2i þ 1,1 sinð2iþ 1Þt: þ
(21)
iZ2
Obviously, u1(t) can be expanded according to Eq. (13): u1 ðtÞ ¼ a1,1 cosðtÞ þ a3,1 cosð3tÞ þ b3,1 sinð3tÞ:
(22)
Based on the harmonic balance, the right hand side of Eq. (21) has no secure terms cosðtÞ, sinðtÞ, cosð3tÞ and sinð3tÞ: That is Z 2 p G2i þ 1,1 ðo1 ,a1,1 ,a3,1 ,b3,1 Þ ¼ F cosð2i þ 1Þt dt ¼ 0, i ¼ 0,1, p 0 1 Z p 2 U2i þ 1,1 ðo1 ,a1,1 ,a3,1 ,b3,1 Þ ¼ F1 sinð2i þ1Þt dt ¼ 0, i ¼ 0,1: (23)
p
0
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
1119
Then, the four coefficients in (22) are simply given by solving four differential equations defined by (23). The first-order approximation to the angular frequency and limit cycle are therefore determined by uð1Þ ðtÞ ¼ ða1,0 þ a1,1 ÞcosðtÞ þa3,1 cosð3tÞ þb3,1 sinð3tÞ, oð1Þ ¼ o0 þ o1 :
(24)
It is easy to find that the residues in Eq. (21) are separated into two parts again. The second part is generally nonzero when we substitute Eq. (24) into Eq. (21). We denote them as X X i F1 ¼ G5,1 cosð5tÞ þ U5,1 sinð5tÞ þ ½G2i þ 1,0 cosð2i þ1Þt þ U2i þ 1,0 sinð2i þ 1Þtpi2 ¼ R21 þ R1 pi2 : (25) iZ3
iZ3
3.3. The second-order approximation Substituting Eqs. (12) and (25) into Eq. (15) and equating the coefficients of p2 yield F2 ðu2 , o2 Þ ¼ : e ol0 ð1u20 ÞDlt u2 2ol0 u0 u1 Dlt u1 ol0 ð2u0 u2 þ u21 ÞDlt u0 2lo0l1 o1 u0 u1 Dlt u0 l1 l2 2 þ lo0l1 o1 ð1u20 ÞDlt u1 þ l o0l1 o2 þ o0 o1 ð1u20 ÞDlt u0 2 1 2=3 1 5=3 2 00 00 2 00 þ o0 u2 þ 2o0 o1 u1 þ ð2o0 o2 þ o1 Þu0 þ u0 u2 u0 u21 þ R21 3 9 ! 2 X X ¼ þ ½G2i þ 1,2 cosð2i þ 1Þt þ U2i þ 1,2 sinð2iþ 1Þt: i¼0
(26)
iZ2
From Eq. (13), we have u2 ðtÞ ¼ a1,2 cosðtÞ þ a3,2 cosð3tÞ þ b3,2 sinð3tÞ þ a5,2 cosð5tÞ þ b5,2 sinð5tÞ:
(27)
Similarly, no secular terms cosðtÞ, cosð3tÞ, cosð5tÞ, sinðtÞ, sinð3tÞ and sinð5tÞ in Eq. (26) imply that their coefficients are all zeros, that is Z 2 p G2i þ 1,2 ðo2 ,a1,2 ,a3,2 ,b3,2 ,a5,2 ,b5,2 Þ ¼ F cosð2iþ 1Þt dt ¼ 0, i ¼ 0,1,2, p 0 2 Z p 2 U2i þ 1,2 ðo2 ,a1,2 ,a3,2 ,b3,2 ,a5,2 ,b5,2 Þ ¼ F2 sinð2iþ 1Þt dt ¼ 0, i ¼ 0,1,2: (28)
p
0
There are six linear algebraic equations for six unknown parameters o2, a1,2 a3,2, b3,2, a5,2 and b5,2, which can be easily solved with the help of symbolic software. Then the second-order improved harmonic balance approximations to the angular frequency and limit cycle can be obtained: uð2Þ ðtÞ ¼
2 X
a1,i cosðtÞ þ
i¼0
2 X
ða3,i cosð3tÞ þ b3,i sinð3tÞÞþ a5,2 cosð5tÞ þ b5,2 sinð5tÞ,
i¼1
oð2Þ ¼ o0 þ o1 þ o2 :
(29)
Substituting Eq. (29) into Eq. (26) yields the new residue defined by X X i F2 ¼ G7,2 cosð7tÞ þ U7,2 sinð7tÞ þ ½G2i þ 1,0 cosð2i þ1Þt þ U2i þ 1,0 sinð2i þ 1Þtpi3 ¼ R32 þ R2 pi3 : iZ4
(30)
iZ4
3.4. The third-order approximation Substituting Eqs. (12) and (30) into Eq. (15) and equating the coefficients of p3 yield
F3 ðu3 , o3 Þ ¼ : o20 u003 þ 2o0 o1 u002 þ ð2o0 o2 þw21 Þu001 þ 2ðo0 o3 þ o1 o2 Þu000 1 2=3 2 5=3 5 8=3 3 u0 u3 u0 u1 u2 þ u0 u1 þ eY þ R32 3 3 27 ! 3 X X þ ½G2i þ 1,3 cosð2iþ 1Þt þ U2i þ 1,3 sinð2i þ1Þt,
þ ¼
i¼0
with
(31)
iZ3
h
Y ¼ ol0 ðu20 1ÞDlt u3 þ2ðu1 u2 þ u0 u3 ÞDlt u0 þ2u0 u1 Dlt u2 þ ð2u0 u2 þ u21 ÞDlt u1 h i þ lo0l1 o1 ðu20 1ÞDlt u2 þ2u0 u1 Dlt u1 þ ð2u0 u2 þ u21 ÞDlt u0
i
1120
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
i l1 2 h 2 þ lo0l2 o0 o2 þ o1 ðu0 1ÞDlt u1 þ 2u0 u1 Dlt u0 2 ðl1Þðl2Þ 3 l3 2 o0 o3 þ ðl1Þo0 þ o3 ðu20 1ÞDlt u0 : þ lo0 6 Expanding u3(t) according to Eq. (13) yields u3 ðtÞ ¼ a1,3 cosðtÞ þa3,3 cosð3tÞ þb3,3 sinð3tÞ þ a5,3 cosð5tÞ þ b5,3 sinð5tÞ þ a7,3 cosð7tÞ þ b7,3 sinð7tÞ:
(32)
Similarly, no secular terms cosð2k þ 1Þt, sinð2k þ 1Þt, and ðk ¼ 0,1,2,3Þ in Eq. (31) imply that their coefficients are all zeros, that is Z 2 p G2i þ 1,3 ðo3 ,a1,3 ,a3,3 ,b3,3 ,a5,3 ,b5,3 ,a7,3 ,b7,3 Þ ¼ F cosð2iþ 1Þt dt ¼ 0, i ¼ 0,1,2,3, p 0 3 Z p 2 U2i þ 1,3 ðo3 ,a1,3 ,a3,3 ,b3,3 ,a5,3 ,b5,3 ,a7,3 ,b7,3 Þ ¼ F3 sinð2i þ 1Þt dt ¼ 0, i ¼ 0,1,2,3: (33)
p
0
There are eight linear algebraic equations for eight unknown parameters o3, a1,3 a3,3, b3,3, a5,3, b5,3, a7,3 and b7,3, which can be easily solved with the help of symbolic software. Then the third-order improved harmonic balance approximations to the angular frequency and limit cycle can be obtained: uð3Þ ðtÞ ¼
3 X
a1,i cosðtÞ þ
i¼0
3 X
ða3,i cosð3tÞ þ b3,i sinð3tÞÞ þ
i¼1
3 X
ða5,i cosð5tÞ þ b5,i sinð5tÞÞ þ a7,3 cosð7tÞ þ b7,3 sinð7tÞ,
i¼2
oð3Þ ¼ o0 þ o1 þ o2 þ o3 :
(34)
The procedure will be iteratively carried out to generate higher-order approximations by solving linear algebraic equations successively. 4. Results and discussion The classical van der Pol damped oscillator with fractional restoring force, governed by the nonlinear second differential Eq. (2), has been studied by Mickens using the first-order harmonic balance method [12] and averaging technique [14], and also by Lim and Lai [7] with the help of combined linearization and harmonic balance method. By Fourier expansion [12], we have ðcos tÞ1=3 ¼
1 X
a2n þ 1 cosð2n þ1Þt,
n¼0
with a2n þ 1 ¼
3Gð7=3Þ 24=3 Gðn þ 5=3ÞGð2=3nÞ
,
n ¼ 0,1,2. . .,
in which G is the Gamma function. Therefore, the first several terms are
cos 3t cos 5t 7 cos 7t cos 9t þ þ , ðcos tÞ1=3 ¼ w cos t 5 10 110 22
(35)
with
w¼
3Gð7=3Þ 24=3 Gð5=3ÞGð2=3Þ
¼
Gð1=3Þ 21=3 Gð2=3Þ2
1:159595266963929:
As is discussed in the third section, the residue for the zeroth-order approximation is separated into two parts. The first part conforms with the initial approximation and the second part comes from the nonlinearity of the function u1/3. In this case, infinite terms should be considered in the each stage of subsequent harmonic balance procedure. It is worthwhile to note that low accuracy of approximation is a direct consequence of discarding the second part. After substituting Eqs. (11) and (35) into Eq. (17), we find that the differential Eq. (18) is defined by
eol0 a10 3 3 lp 1=3 þ wa10 o0 2 a10 ¼ 0, eol0 a10 ¼ 0: (36) eol0 a10 a210 1 cos 4 2 4 The amplitude a10 of the zeroth-order approximation (19) is determined by the second equation in (36), which gives að0Þ ¼ a10 ¼ 2 while o0 depends on the solution of the first nonlinear equation in (36). Since higher-order approximations of the residue harmonic balance method are governed by linear equations only, solutions can easily be obtained to any desired accuracy. To show the effectiveness of the proposed method for fractional van der Pol oscillator (4), the frequency–amplitude response under the alternations of fractional order l and control parameter e will be illustratively validated by some figures, wherein the obtained results are compared with the numerical ones.
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
1121
4.1. Effect of fractional order on the frequency and amplitude of steady-state response The fractional order–frequency and fractional order–amplitude curves of limit cycles for e ¼1 are shown in Figs. 1 and 2, respectively. It is easy to find that the angular frequency of limit cycle in steady state increases at first, and then decreases as the fractional order varies in the interval (0,1]. However, the amplitude of limit cycle in steady state decreases at first, and then increases as the fractional order varies in the interval (0,1]. It is also clear that the obtained approximations are highly accurate and in good agreement with the numerically integrated results based on the phase portraits in Fig. 3 when e ¼ 1 for l ¼{0.5,0.8,1}. 4.2. Effect of control parameter on the frequency and amplitude of steady-state response When l ¼0.5, the fractional van der Pol oscillator (4) is reduced to a specific case: 1=3 € e½1uðtÞ2 D0:5 ¼ 0: uðtÞ t uðtÞ þuðtÞ
(37)
For this case, the comparisons between the asymptotic approximations to the angular frequency and amplitude of limit cycle and corresponding numerical results are illustrated by Figs. 4 and 5. Obviously, we can observe that the obtained approximations are highly accurate. Moreover, as the control parameter e varies in the interval (0,2], the angular frequency of the limit cycle are always increasing, while the amplitude decreases throughout the whole region of e. To further illustrate and verify the accuracy and effectiveness of the residue harmonic balance approximations, the comparison of the 1.7 1.6
Frequency of limit cycle
1.5 1.4 1.3 1.2 1.1 1 0.9 0.8 0.7 0
0.1
0.2
0.3
0.4 0.5 0.6 Fractinal order
0.7
0.8
0.9
1
Fig. 1. Comparisons of frequency between numerical solutions and results provided by the presented method, where solid lines and dashed lines denote the third-order and the zeroth-order approximations, respectively; dotted line denotes the numerical result. The control parameter is given as e ¼1.
2.25
Amplitude of limit cycle
2.2 2.15 2.1 2.05 2 1.95 1.9 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fractinal order Fig. 2. Comparisons of amplitude between numerical solutions and results provided by the presented method, where solid lines and dashed lines denote the third-order and the zeroth-order approximations, respectively; dotted line denotes the numerical result. The control parameter is given as e ¼1.
1122
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
3 2
Dλt u(t)
1 0 λ=1 λ=0.8 λ=0.5
-1 -2 -3 -2.5 -2
-1.5 -1
-0.5
0 u(t)
0.5
1
1.5
2
2.5
Fig. 3. Phase portraits of limit cycle, where solid lines and dotted lines denote the third-order analytical solution and numerical result, respectively. The control parameter and fractional derivative are given as e ¼ 1, l ¼{0.5,0.8,1}.
2.6
Frequency of limit cycle
2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0
0.2
0.4
0.6
1.2 1.4 0.8 1 Control parameter
1.6
1.8
2
Fig. 4. Comparisons of frequency between numerical solutions and results provided by the presented method, where solid lines and dashed lines denote the third-order and the zeroth-order approximations, respectively, dotted line denotes the numerical result. The fractional derivative is given as l ¼0.5.
2.04
Amplitude of limit cycle
2.02 2 1.98 1.96 1.94 1.92 1.9 0
0.2
0.4
0.6
0.8 1 1.2 1.4 Control parameter
1.6
1.8
2
Fig. 5. Comparisons of amplitude between numerical solutions and results provided by the presented method, where solid lines and dashed lines denote the third-order and the first-order approximations, respectively; dotted line denotes the numerical result. The fractional derivative is given as l ¼0.5.
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
1123
3 2 ε=1
D0.5 t u(t)
1
ε=2
ε=0.1
0 -1 -2 -3 -2.5 -2
-1.5 -1
-0.5
0 0.5 u(t)
1
1.5
2
2.5
Fig. 6. Phase portraits of limit cycle, where solid lines and dotted lines denote the third-order analytical solution and numerical result, respectively. The control parameter and fractional derivative are given as e ¼{0.1.1.2}, l ¼0.5.
Table 1 3rd-order analytical expressions of fractional vdP oscillator when l ¼0.5.
e
0.1
1
1.5
2
o
0.929607 2.02137 0.0217533 0.0205061 0.00658125 0.000356124 0.00218757 0.0000279856
1.50826 1.85721 0.0630316 0.0954636 0.000080662 0.00810594 0.00127273 0.000330784
1.77356 1.82287 0.0778078 0.110466 0.001301 0.0114029 0.00135085 0.000483207
2.01975 1.80172 0.086429 0.119604 0.00218938 0.0136131 0.00144587 0.000586111
cosðotÞ cosð3otÞ sinð3otÞ cosð5otÞ sinð5otÞ cosð7otÞ sinð7otÞ
phase portraits is presented in Fig. 6 for e ¼{0.1,1,2}. The analytical expressions of the limit cycle of oscillator (37) are tabulated in Table 1.
5. Conclusion The residue harmonic balance method to obtain the accurate approximations to the angular frequency and limit cycles of vibration systems is proposed. The fractional order van der Pol oscillator with fractional restoring force and damping is investigated successfully. It is shown that the approximations of frequencies, amplitudes and period solutions to the numerical ones are in excellent agreement. The effects of fractional order and the control parameter on the steady-state responses of angular frequency and amplitude of oscillator are illustrated by the several figures. It is worthwhile to point out that, for the fractional van der Pol oscillator with fractional restoring force, the residue, produced by the fractional power function, in the zeroth-order harmonic equation is separated into different stages and added to the subsequent harmonic procedures to improve the accuracy of the approximations. Therefore, the presented approach can be easily applied to other nonlinear oscillators with fractional power and derivative nonlinearities. Since the higher-order approximations of residue harmonic balance method are governed by linear equations, solutions can easily be obtained to any desired accuracy.
Appendix A. Newton iteration for differential equations The Newton method has been very successful to solve nonlinear algebraic equations. We extend the method to solve differential equations using the series here. When dealing with algebraic equations, one would expect quadratic convergence, that is, one iteration would double the accuracy of a root. When dealing with differential equations using series to express the solution function, one would expect to double the effective number of terms in the function expression. We shall illustrate the method in solving a simple ordinary differential equation of the form by power series: € u,u,tÞ _ Fðu, ¼ u€ þ u_ þ u2 t ¼ 0:
1124
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
Let the initial (zeroth) approximation be u0(t) ¼a0 þa1t þa2t2:
F0 ¼ Fðu€ 0 , u_ 0 ,u0 ,tÞ ¼ a0 2 þ a1 þ 2a2 þ ð2a2 þ 2a0 a1 1Þt þ ða21 þ 2a0 a2 Þt2 þ ð2a1 a2 Þt3 þ a22 t4 ¼ R0,0 þ R1,0 t þ R2,0 t 2 þ R3,0 t 3 þ R4,0 t 4 ¼
2 X
Ri,0 t i þ
i¼0
4 X
Ri,0 t i ¼
i¼3
2 X
Ri,0 t i þ R0 ðtÞ:
i¼0
P2 i _ R2,0 ¼ 0. The residue is separated in two parts: balanceable, e.g., Let R0,0 ¼ uð0Þ; R1,0 ¼ uð0Þ; i ¼ 0 Ri,0 t , and nonP4 i balanceable, e.g., i ¼ 3 Ri,0 t . There are 3 nonlinear equations in the balanceable part for the 3 unknowns a0, a1, a2, and u0(t) can be determined. However, the non-balanceable part, Ri,0a0; i ¼3,4. For R0,0 ¼0; R1,0 ¼0; R2,0 ¼0, one has 8 9 8 9 > < a0 > = > < 1:1040 > = 1:8367 : a ¼ a1 ¼ > > :a > ; > : ; 1:5278 2 Let the first approximation be u(1)(t)¼u0(t) þu1(t). u1 ðtÞ ¼ b0 þ b1 t þ b2 t 2 þ b3 t 3 þ b4 t 4 with coefficients bi to be determined. Substitute u(1) into
Fðu€ 0 þ u€ 1 , u_ 0 þ u_ 1 ,u0 þu1 ,tÞ ¼ ðu€ 0 þ u€ 1 Þ þ ðu_ 0 þ u_ 1 Þ þ ðu0 þ u1 Þ2 t ¼ 0: Retain only the linear terms in u1:
@ @ @ þ u_ 1 þ u1 F1 ¼ F0 þ u€ 1 þ u_ 1 þ 2u0 u1 ¼ R0 þ u€ 1 F90 ¼ R0 þ @u€ @u_ @u
2 X
uðkÞ 1
k¼0
! @ F90 ¼ R0 þD1 F90 ¼ 0 @uðkÞ
Here the number in the superscript brackets denotes the order of differentiation w.r.t. t. F90 indicates that the function is to be evaluated at u ¼u0 after partial differentiation:
F1 ¼ b1 þ 2b2 þ 2a0 b0 þ ð2b2 þ6b3 þ 2a0 b1 þ 2a1 b0 Þt þ ð3b3 þ12b4 þ 2a0 b2 þ2a1 b1 þ2a2 b0 Þt 2 þ ð4b4 þ 2a1 a2 þ 2a0 b3 þ2a1 b2 þ2a2 b1 Þt 3 þ ða22 þ 2b2 a2 þ 2a0 b4 þ 2a1 b3 Þt 4 þ ð2a1 b4 þ 2a2 b3 Þt 5 þ ð2a2 b4 Þt 6 ¼
4 X
Ri,1 t i þ R1 ðtÞ
i¼0
From the balanceable part, Ri,1 ¼0; i¼0, 1, 2, 3, 4, we have 5 linear equations for the 5 unknowns bi and u1(t) can be determined: 9 2 38 9 8 1 2 0 0 > b0 > > 0 > 2a0 > > > > > > > > >b > > 0 > 6 2a > 2 6 0 7 > > > > 1> > 6 1 2a0 7> = 6 7< = < 6 2a2 2a1 2a0 7 0 3 12 7 b2 þ ¼ 0: 6 > > > > >b > 6 0 > > > > 2a1 a2 > 4 7 2a2 2a1 2a0 > > > 4 5> 3> > > > > > > > > 0 0 2a2 2a1 2a0 : b4 ; : a22 ; Then bT ¼[0.5620, 0.7003, 0.9706, 0.2783, 0.1769]. Therefore, F1 ¼R1. However, the non-balanceable part Ri,1a0; i¼ 5,6. P Let the second approximation be u(2)(t)¼u0(t) þu1(t)þu2(t). u2 ðtÞ ¼ 6i ¼ 0 ci t i where ci are 7 constants to be determined. Substituted into the governing equation and retain only the linear terms in u2: 2 !2 3 2 2 6 X X @ 1 X ðkÞ @ 5F9 ¼ F2 ¼ R1 ðtÞ þ 4 uðkÞ þ u Ri,2 t i þ R2 ðtÞ ¼ 0, 2 1 0 ðkÞ ðkÞ 2! k ¼ 0 @u @u i¼0 k¼0
F2 ¼ ðb24 þ2a2 c6 Þt8 þ ð2a1 c6 þ 2a2 c5 þ 2b3 b4 Þt7 þ ðb23 þ2a2 b4 þ 2a0 c6 þ2a1 c5 þ 2a2 c4 þ 2b2 b4 Þt6 þ ð6c6 þ 2a1 b4 þ 2a2 b3 þ 2a0 c5 þ 2a1 c4 þ 2a2 c3 þ 2b1 b4 þ 2b2 b3 Þt 5 2
þðb2 þ 5c5 þ 30c6 þ 2a0 c4 þ 2a1 c3 þ2a2 c2 þ 2b0 b4 þ 2b1 b3 Þt 4 þ ð4c4 þ 20c5 þ 2a0 c3 þ 2a1 c2 þ 2a2 c1 þ2b0 b3 þ 2b1 b2 Þt 3 2
þ ðb1 þ3c3 þ 12c4 þ 2a0 c2 þ 2a1 c1 þ2a2 c0 þ 2b0 b2 Þt 2 þ ð2c2 þ 6c3 þ 2a0 c1 þ 2a1 c0 þ2b0 b1 Þt 2
þ b0 þc1 þ 2c2 þ 2a0 c0 ,
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
2
2a0
6 2a 6 1 6 6 2a2 6 6 6 6 6 6 6 4
1
2
2a0
2
2a1
2a0
6 3
2a2
2a1
2a0
4
20
2a2
2a1 2a2
2a0 2a1
5 2a0
2a2
2a1
12
1125
9 38 9 8 2 > c0 > > b0 > > > > > > > > > > > 7> > > > >c > > > 2b b 1 7> > > > 0 > > > > > > > 7> 2 > > > >c > > > 7> b1 þ 2b0 b2 > = = > < < 2> 7> 7 c 2b0 b3 þ 2b1 b2 ¼ 0: þ 3 7 > > 7> > > > > > 2 > > 30 7 b2 þ 2b0 b4 þ 2b1 b3 > > > > > c4 > > > > 7> > > > > > > >c > 6 7 > > > > 2a1 b4 þ 2a2 b3 þ 2b1 b4 þ2b2 b3 > > 5> 5> > > > > > > > ; > > : > 2 ; : 2a0 c6 b3 þ 2a2 b4 þ 2b2 b4
Let Ri,2 ¼ 0; i ¼ 0,1,. . .,6. We can determine ci and u2. By mathematical induction, one can prove that when the (r 1)th approximation is available, the rth approximation is given by uðrÞ ðtÞ ¼ u0 ðtÞ þ þur ðtÞ Let ur ðtÞ ¼
P2ðr þ 1Þ
i
ai,r t 2 2 X @ 1 Fr ¼ Rr1 ðtÞ þ 4 uðkÞ þ r ðkÞ 2! @u k¼0 i¼0
2 X k¼0
uðkÞ r1
@ @uðkÞ
!2
1 þ þ r!
2 X k¼0
uðkÞ 1
@ @uðkÞ
!r 3 2ðr þ 1Þ X 5F9 ¼ Ri,r t i þ Rr ðtÞ ¼ 0: 0 i¼0
Let the balanceable part be Ri,r ¼ 0; i ¼ 0,1,. . .,2ðr þ 1Þ. We can determine ai,r and ur. When we use sine and cosine functions instead of power series, we have the residue harmonic balance method. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]
B. van der Pol, A theory of the amplitude of free and forced triode vibrations, Radiology Review 1 (1920) 701–710; 754–762. B. van der Pol, On relaxation oscillations, Philosophical Magazine 7 (2) (1926) 978–992. J. Guckenheimer, K. Hoffman, W. Weckesser, Numerical computation of canards, International Journal of Bifurcation and Chaos 10 (2000) 2669–2687. Y.H. Kao, C.S. Wang, Analog study of bifurcation structures in a van der Pol oscillator with a nonlinear restoring force, Physical Review E 48 (4) (1993) 2514–2520. C. Letellier, V. Messager, R. Gilmore, From quasi-periodicity to toroidal chaos: analogy between the Curry–Yorke map and the van der Pol system, Physical Review E 77 (2008) 046203. R.E. Mickens, Fractional van der Pol equations, Journal of Sound and Vibration 259 (2003) 457–460. C.W. Lim, S.K. Lai, Accurate higher-order analytical approximate solutions non-conservative nonlinear oscillators and application to van der Pol damped oscillators, International Journal of Nonlinear Mechanics 48 (2006) 483–492. G.M. Mahmoud, A.A.M. Farghaly, Chaos control of chaotic limit cycles of real and complex van der Pol oscillators, Chaos Solitons and Fractals 21 (2004) 915–924. C.M.A. Pinto, J.A.T. Machado, Complex order van der Pol oscillator, Nonlinear Dynamics, doi: 10.1007/s11071-010-9886-0. Y.S. Chen, A.Y.T. Leung, Bifurcation and Chaos in Engineering, Springer-Verlag, London, 1998. R.E. Mickens, A generalized iteration procedure for calculating approximations to periodic solutions of ‘‘truly nonlinear oscillators’’, Journal of Sound and Vibration 287 (2005) 1045–1051. R.E. Mickens, Analysis of nonlinear oscillators having non-polynomial elastic terms, Journal of Sound and Vibration 255 (4) (2002) 789–792. R.E. Mickens, Iteration method solutions for conservative and limit-cycle x1/3 force oscillators, Journal of Sound and Vibration 292 (2006) 964–968. R.E. Mickens, A combined equivalent linearization and averaging perturbation method for non-linear oscillator equations, Journal of Sound and Vibration 264 (2003) 1195–1200. I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier, New York, 2006. N. Shimizu, W. Zhang, Fractional calculus approach to dynamic problems of viscoelastic materials, JSME International Journal Series C, Mechanical Systems, Machine Elements and Manufacturing 42 (4) (1999) 825–837. Y.A. Rossikhin, M.V. Shitikova, Application of fractional calculus for dynamic problems of solid mechanics: novel trends and recent results, Applied Mechanics Reviews 63 (2010) 010801–010852. I. Petras, Fractional Order Nonlinear Systems: Modeling, Analysis and Simulation, Springer-Verlag, London, 2011. G.M. Zaslavsky, Hamiltonian Chaos and Fractional Dynamics, Oxford University Press, Oxford, 2005. V.E. Tarasov, Fractional variations for dynamical systems: Hamilton and Lagrange approaches, Journal of Physics A: Mathematical and General 39 (2006) 8409–8425. V.E. Tarasov, Fractional systems and fractional Bogoliubov hierarchy equations, Physical Review E 71 (2005) 011102. N. Laskin, Fractional quantum mechanics, Physical Review E 62 (2000) 3135–3145. A.I. Saichev, G.M. Zaslavsky, Fractional kinetic equations: solutions and applications, Chaos 7 (1997) 753–764. R. Magina, M.D. Ortigueirab, Igor Podlubny, J. Trujillo, On the fractional signals and systems, Signal Processing 91 (3) (2011) 350–371. R.S. Barbosa, J.A.T. Machado, M.I. Ferreira, K.J. Tar, Dynamics of the fractional order Van der Pol oscillator, Proceedings of the Second IEEE International Conference on Computational Cybernetics (ICCC’04), Vienna University of Technology, Austria, 30 August–1 September 2004, pp. 373–378. T. Sardar, S.S. Ray, R.K. Bera, B.B. Biswas, The analytical approximate solution of the multi-term fractionally damped van der Pol equation, Physica Scripta 80 (2009) 025003–025006. F. Xie, X.Y. Lin, Asymptotic solution of the van der Pol oscillator with small fractional damping, Physica Scripta 136 (2009) 014033–014034. V. Gafiychuk, B. Datsko, V. Meleshko, Analysis of fractional order Bonhoeffer–van der Pol oscillator, Physica A 387 (2008) 418–424. M. Attari, M. Haeri, M.S. Tavazoei, Analysis of a fractional order van der Pol like oscillator via describing function method, Nonlinear Dynamics 61 (2010) 265–274. M.S. Tavazoei, M. Haeri, M. Attari, S. Bolouki, M. Siami, More details on analysis of fractional order Van der Pol oscillator, Journal of Vibration and Control 15 (6) (2009) 803–819. A. Konuralp, C. Konuralp, A. Yıldırım, Numerical solution to the van der Pol equation with fractional damping, Physica Scripta 136 (2009) 014034–014035. J.H. Chen, W.C. Chen, Chaotic dynamics of the fractionally damped van der Pol equation, Chaos Solitons and Fractals 35 (2008) 188–198. A.Y.T. Leung, Z.J. Guo, Forward residue harmonic balance for autonomous and non-autonomous systems with fractional derivative damping, Communications in Nonlinear Science and Numerical Simulations 16 (4) (2011) 2169–2218.
1126
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
[35] Z.J. Guo, A.Y.T. Leung, H.X. Yang, Oscillatory region and asymptotic solution of fractional van der Pol oscillator via residue harmonic balance technique, Applied Mathematical Modelling 35 (2011) 3918–3925. [36] A.Y.T. Leung, Z.J. Guo, Residue harmonic balance for discontinuous nonlinear oscillator with fractional power restoring force, International Journal of Nonlinear Science and Numerical Simulations 11 (9) (2010) 705–723. [37] E. Pereira, C.A. Monje, B.M. Vinagre, F. Gordillo, MATLAB toolbox for the analysis of fractional order systems with hard nonlinearities, Proceedings of the First IFAC Congress on Fractional Differentiation and its Applications, Bordeaux (France), July 2004, pp. 214–219.
Contents lists available at SciVerse ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
The residue harmonic balance for fractional order van der Pol like oscillators A.Y.T. Leung n, H.X. Yang, Z.J. Guo Department of Civil and Architectural Engineering, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong
a r t i c l e i n f o
abstract
Article history: Received 14 July 2011 Accepted 17 October 2011 Handling Editor: L.N. Virgin Available online 16 November 2011
When seeking a solution in series form, the number of terms needed to satisfy some preset requirements is unknown in the beginning. An iterative formulation is proposed so that when an approximation is available, the number of effective terms can be doubled in one iteration by solving a set of linear equations. This is a new extension of the Newton iteration in solving nonlinear algebraic equations to solving nonlinear differential equations by series. When Fourier series is employed, the method is called the residue harmonic balance. In this paper, the fractional order van der Pol oscillator with fractional restoring and damping forces is considered. The residue harmonic balance method is used for generating the higher-order approximations to the angular frequency and the period solutions of above mentioned fractional oscillator. The highly accurate solutions to angular frequency and limit cycle of the fractional order van der Pol equations are obtained analytically. The results that are obtained reveal that the proposed method is very effective for obtaining asymptotic solutions of autonomous nonlinear oscillation systems containing fractional derivatives. The influence of the fractional order on the geometry of the limit cycle is investigated for the first time. & 2011 Elsevier Ltd. All rights reserved.
1. Introduction Since the original work of Balthazar van der Pol [1], the van der Pol equation has always been a prototype for systems with self-excited, limited cyclic oscillations. Nowadays, the Van der Pol oscillator as well as its modified versions is still widely studied in various scientific fields, ranging from biology, chemistry and physics to engineering [2–10]. It is given by the following second-order differential equation: € _ þ uðtÞ ¼ 0, uðtÞ eð1u2 ðtÞÞuðtÞ
(1)
for which the elastic restoring force is a linear function of the dependence variable. Here the over-dot denotes differentiation with respect to time variable t, e is a control parameter that reflects the degree of nonlinearity of the system. For e ¼ 0 this is just a simple linear oscillator. For large values of e, which is for e b1, the system exhibits a relaxation oscillation [2]. The generalized van der Pol oscillators, in which the fractional power of dependent variable u or its derivative occurs, have been of great interest for many years [6,7,11–14]. In [10], with the help of the general harmonic balance method,
n
Corresponding author. Tel.: þ852 34437600. E-mail addresses: [email protected], [email protected] (A.Y.T. Leung).
0022-460X/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jsv.2011.10.023
1116
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
Mickens derived the expression for the limit cycle of following van der Pol type oscillator: € _ þ u1=3 ðtÞ ¼ 0: uðtÞ eð1u2 ðtÞÞuðtÞ
(2)
Mickens [11] also proposed an iteration technique to the oscillator (2) and obtained the same result, but the frequency decreased for about 15% with respect to the classical van der Pol oscillator. To avoid tedious computations in dealing with nonlinear algebras involving in harmonic balance equations for higher-order approximations, Lim and Lai [7] introduced a linearization technique prior to proceeding with harmonic balance procedure and studied (2) under the assumption that 0 r e 5 1. The transient solutions for the first-, second- and third-order approximations are given with the help of averaging method [13] upon taking into account the weak damping effect for (2). Recently, the applications of fractional calculus [15,16] have grown rapidly, which has attracted fairly broad research activity, including the viscoelasticity [17,18], chaotic dynamics [19,20], mechanics of non-Hamiltonian systems [21,22], quantum mechanics [23], physical kinetics [24] and so on [25]. Due to the occurrence of fractional oscillators in dynamical applications, various kinds of van der Pol equations containing fractional derivatives have attracted more and more attention, see for instance [19,20,26–37] and references therein. Two important physical parameters, i.e., the angular frequency and the amplitude of oscillation are involved in the investigations of nonlinear oscillators. Many valuable contributions were made to calculate these two physical parameters. However, the solution of the differential equation containing a fractional derivative is very involved. An effective and easyto-use method for solving such equations is needed. Combining the harmonic balance method with parameter expansion technique, an analytical technique, named residue harmonic balance method, is recently constructed [34–36]. The primary idea of this method is that a bookkeeping parameter p is introduced into the residue of the Fourier truncation series in harmonic equations. The unbalanced residues due to Fourier truncation are considered iteratively by solving linear algebraic equations to improve the accuracy of the solutions successively. Using of residue harmonic balance method, this paper is devoted to investigating the approximations of the fractional van der Pol oscillator with fractional restoring and damping forces. The highly accurate solutions to the angular frequency and limit cycle are obtained. It is shown that the proposed technique is very effective for obtaining asymptotic solutions for nonlinear fractional order systems. 2. Fractional van der Pol oscillator Recently, study on dynamical behaviors of the fractional order systems has attracted increasing attention. In these studies, it has been found that the fractional order systems like integer order systems can generate oscillations. Pereira et al. considered a fractional version of the van der Pol equation given by _ þ uðtÞ ¼ 0, Dlt uðtÞ þ e½uðtÞ2 1uðtÞ
1 o l o 2,
which is obtained by substituting the capacitance by a fractance in the nonlinear RLC circuit model [19,37]. Xie and Lin [28] introduced the following Van der Pol oscillator with a small fractional damping term and studied the asymptotic solution of Cauchy problem: € uðtÞ e½1uðtÞ2 Dlt uðtÞ þ uðtÞ ¼ 0, _ where uð0Þ ¼ a, uð0Þ ¼ b,0 o l o1. Barbosa et al. [26] also suggested the introduction of a fractional order time derivative for Eq. (1) of the standard van der Pol oscillator in the form: D1t þ l uðtÞe½1uðtÞ2 Dlt uðtÞ þ uðtÞ ¼ 0,
0 o l o 1:
(3)
In [35], the oscillatory region and asymptotic solutions of Eq. (3) are discussed by the residue harmonic balance method. In this study, motivated by Mickens [6], a modified version of fractional van der Pol oscillator is considered: € uðtÞ e½1uðtÞ2 Dlt uðtÞ þuðtÞ1=3 ¼ 0,
(4)
where u is a displacement, e is a control parameter that reflects the degree of nonlinearity of the system and e 40, 0 o l o1 is the order of the fractional derivative in the Caputo sense defined by Z t m 1 d Dat uðtÞ ¼ ðttÞma1 m uðtÞdt, (5) GðmaÞ 0 dt for m1 o a rm, m 2 N, t 40, u 2 C m 1 : When l ¼1, Eq. (4) is reduced to the modified van der Pol oscillator (2) with fractional restoring force [6,7]. 3. Residue harmonic balance method The principle of the residue harmonic balance method to differential equations analogous to the Newton iteration for nonlinear algebraic equations is illustrated in Appendix A using power series. To solve Eq. (4) using the residue harmonic
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
1117
balance method, we introduce a new time scale t ¼ ot, then (4) reads
o2 u00 ðtÞeol ð1u2 ðtÞÞDlt uðtÞ þu1=3 ðtÞ ¼ 0,
(6)
where the prime denotes differentiation of u with respect to t, and o is an unknown parameter to be determined. Due to the symmetric characteristic of Eq. (6), u(t) can be expressed by uðtÞ ¼ a0 cosðtÞ þ
1 X
ða2k þ 1 cos½ð2k þ 1Þt þ b2k þ 1 sin½ð2k þ1ÞtÞ:
(7)
k¼0
Therefore, instead of the function u(t), the coefficients a0 ,a2k þ 1 , b2k þ 1 ðk ¼ 0,1,2,. . .Þ become the main unknowns to be determined. Definition (5) of the fractional derivatives allows us to obtain the expression of Dbt uðtÞ, ð0 o b o1Þ, where Z t 1 d a cosðksÞ þbk sinðksÞ ds ðtsÞb Dbt ak cosðktÞ þbk sinðktÞ ¼ ds k Gð1bÞ 0 Z t k t b bk cosðktktÞak sinðktktÞ dt, (8) ¼ Gð1bÞ 0 According to the trigonometric identities cosðktktÞ ¼ cos kt cos kt þsin kt sin kt, sinðktktÞ ¼ sin kt cos ktcos kt sin kt, Eq. (8) becomes
" # Z kt Z kt b k b b cosðktÞ u ðbk cos u þ ak sin uÞdu þ sinðktÞ u ðbk sin uak cos uÞdu Eq: ð8Þ9 Gð1bÞ 0 0 Z Z kt ub cos u kt ub sin u b b du þ k ak cosðktÞ þ bk sinðktÞ du ¼ k bk cosðktÞak sinðktÞ Gð1bÞ Gð1bÞ 0 0 b
b
¼ k ðbk I1 þ ak I2 ÞcosðktÞ þ k ðbk I2 ak I1 ÞsinðktÞ,
(9)
where I1 9
Z kt b u cos u du, Gð1bÞ 0
I2 9
Z kt b u sin u du: Gð1bÞ 0
Although these integrals cannot be evaluated in closed form for general interval, we are able to evaluate them using the contour integration, when steady-state response is of interest, where pffiffiffiffi pG½ð1bÞ=2 bp I1 ¼ sin , I1 ¼ b 1 ¼ tlim -1 2 2 Gðb=2ÞGð1bÞ pffiffiffiffi pGð1ðb=2ÞÞ bp I1 ¼ cos : I2 ¼ b 2 ¼ tlim -1 2 2 G½ð1 þ bÞ=2Gð1bÞ 1 In the following calculation, we substitute I1 and I2 in Eq. (9) byI1 1 and I2 , respectively. Substituting Eq. (7) into Eq. (6) and using the Galerkin procedure result in the following harmonic balance equations: Z 2 p Rci ðo,a1 ,a3 ,. . .,a2n þ 1 ,b3 ,. . .,b2n þ 1 Þ ¼ ½ou€ þ aol ðu2 1ÞDlt u þ u1=3 cosð2i þ 1Þt dt ¼ 0,i ¼ 0,1,. . .,n, p 0 Z 2 p 2 Rsi ðo,a1 ,a3 ,. . .,a2n þ 1 ,b3 ,. . .,b2n þ 1 Þ ¼ ½o u€ þ aol ðu2 1ÞDlt u þu1=3 sinð2i þ 1Þt dt ¼ 0,i ¼ 0,1,. . .,n:
p
0
The direct integration leads to ½Kðq, mÞfqg ¼ fRg ¼ 0, where [K(q,m)] is the total stiffness matrix, Consider the initial approximation:
fRg ¼ ½Rc0 ,Rc1 ,. . .,Rcn ,Rs1 ,. . .,Rsn T
(10) T
and fqg ¼ ½o,a1 ,a3 ,. . .,a2n þ 1 ,b1 ,b3 ,. . .,b2n þ 1 .
o ¼ o0 , u0 ðtÞ ¼ a1,0 cosðtÞ:
(11)
We introduce a bookkeeping parameter pA[0,1], and rewrite u(t) and o as u(t,p) and o(p), respectively:
oðpÞ ¼ o0 þ po1 þ p2 o2 þ , uðt,pÞ ¼ u0 ðtÞ þ pu1 ðtÞ þ p2 u2 ðtÞ þ :
(12)
Expanding the unknown functions uk(t) into Fourier series, one has uk ðtÞ ¼ a1,k cosðtÞ þ
k X i¼1
fa2i þ 1,k cos½ð2i þ1Þt þb2i þ 1,k sin½ð2i þ 1Þtg, k ¼ 1,2,. . . :
(13)
1118
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
It can be shown that all these series are convergent at p¼1. Then the final solution will be ( o ¼ o0 þ o1 þ o2 þ ,
(14)
uðtÞ ¼ u0 ðtÞ þu1 ðtÞ þu2 ðtÞ þ : 3.1. The zeroth-order approximation
Just for the convenience of discussion below, we introduce a function F that represents the governing equation of fractional van der Pol oscillator:
Fðuðt,pÞ, oðpÞÞ ¼ : o2 ðpÞu00 ðt,pÞeol ðpÞð1u2 ðt,pÞÞDlt uðt,pÞ þ u1=3 ðt,pÞ,
(15)
where using the Taylor series, we have
1 2=3 1 2=3 1 5=3 u0 u1 p þ u0 u2 u0 u21 p2 þ , 3 3 3 l 1 ol ¼ ðo0 þ o1 p þ o2 p2 þ Þl ¼ ol0 þ lo0l1 o1 p þ l o0l1 o2 þ o0l2 o21 p2 þ : 2
u1=3 ¼ ðu0 þu1 p þu2 p2 þ Þ1=3 ¼ u0 1=3 þ
(16)
Substituting Eqs. (8), (12) and (13) into Eq. (15) and equating the coefficients of p0 lead to
F0 ðu0 , o0 Þ ¼: o20 u000 eol0 ð1u20 ÞDlt u0 þ u1=3 0 ¼ G1,0 cost þ U1,0 sin t þ
X
½G2i þ 1,0 cosð2i þ1Þt þ U2i þ 1,0 sinð2i þ 1Þt:
(17)
iZ1
The residue in Eq. (17) is separated into two parts: the first part conforms with the initial approximation in Eq. (11), the second part in the summation sign for the determination of the next higher approximations, the number of terms of which will depend on the nature of nonlinearity of u. Based on the harmonic balance, the right hand side of Eq. (17) should not contain the secure terms cosðtÞ and sinðtÞ. Letting their coefficients be zero, we have Z 2 p G1,0 ðo0 ,a10 Þ ¼ F cosðtÞdt ¼ 0, p 0 0 Z p 2 U1,0 ðo0 ,a10 Þ ¼ F0 sinðtÞdt ¼ 0: (18)
p
0
Then the unknown coefficients o0, a1,0 for the initial approximation result by simply solving the nonlinear algebraic equation (18), then we have uð0Þ ðtÞ ¼ a1,0 cosðtÞ, oð0Þ ¼ o0 :
(19)
However, the second part of residue under the summation sign of Eq. (17) cannot be identically zero and should be considered in the first approximation procedure. We rewrite it as X X i F0 ¼ G3,0 cosð3tÞ þ U3,0 sinð3tÞ þ ½G2i þ 1,0 cosð2i þ 1Þt þ U2i þ 1,0 sinð2i þ 1Þtpi1 ¼ R10 þ R0 pi1 : (20) iZ2
iZ2
3.2. The first-order approximation Substituting Eqs. (12) and (20) into Eq. (15) and equating the coefficients of p1 lead to h i F1 ðu1 , o1 Þ ¼ : o0 2 u001 e ol0 ð1u0 2 ÞDlt u1 2ol0 u0 u1 Dlt u0 þ lo0l1 o1 ð1u0 2 ÞDlt u0 þ 2o0 o1 u000 þ
1 X 1 2=3 u0 u1 þ R10 ¼ ½G2i þ 1,1 cosð2i þ 1Þt þ U2i þ 1,1 sinð2iþ 1Þt 3 i¼0 X ½G2i þ 1,1 cosð2i þ 1Þt þ U2i þ 1,1 sinð2iþ 1Þt: þ
(21)
iZ2
Obviously, u1(t) can be expanded according to Eq. (13): u1 ðtÞ ¼ a1,1 cosðtÞ þ a3,1 cosð3tÞ þ b3,1 sinð3tÞ:
(22)
Based on the harmonic balance, the right hand side of Eq. (21) has no secure terms cosðtÞ, sinðtÞ, cosð3tÞ and sinð3tÞ: That is Z 2 p G2i þ 1,1 ðo1 ,a1,1 ,a3,1 ,b3,1 Þ ¼ F cosð2i þ 1Þt dt ¼ 0, i ¼ 0,1, p 0 1 Z p 2 U2i þ 1,1 ðo1 ,a1,1 ,a3,1 ,b3,1 Þ ¼ F1 sinð2i þ1Þt dt ¼ 0, i ¼ 0,1: (23)
p
0
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
1119
Then, the four coefficients in (22) are simply given by solving four differential equations defined by (23). The first-order approximation to the angular frequency and limit cycle are therefore determined by uð1Þ ðtÞ ¼ ða1,0 þ a1,1 ÞcosðtÞ þa3,1 cosð3tÞ þb3,1 sinð3tÞ, oð1Þ ¼ o0 þ o1 :
(24)
It is easy to find that the residues in Eq. (21) are separated into two parts again. The second part is generally nonzero when we substitute Eq. (24) into Eq. (21). We denote them as X X i F1 ¼ G5,1 cosð5tÞ þ U5,1 sinð5tÞ þ ½G2i þ 1,0 cosð2i þ1Þt þ U2i þ 1,0 sinð2i þ 1Þtpi2 ¼ R21 þ R1 pi2 : (25) iZ3
iZ3
3.3. The second-order approximation Substituting Eqs. (12) and (25) into Eq. (15) and equating the coefficients of p2 yield F2 ðu2 , o2 Þ ¼ : e ol0 ð1u20 ÞDlt u2 2ol0 u0 u1 Dlt u1 ol0 ð2u0 u2 þ u21 ÞDlt u0 2lo0l1 o1 u0 u1 Dlt u0 l1 l2 2 þ lo0l1 o1 ð1u20 ÞDlt u1 þ l o0l1 o2 þ o0 o1 ð1u20 ÞDlt u0 2 1 2=3 1 5=3 2 00 00 2 00 þ o0 u2 þ 2o0 o1 u1 þ ð2o0 o2 þ o1 Þu0 þ u0 u2 u0 u21 þ R21 3 9 ! 2 X X ¼ þ ½G2i þ 1,2 cosð2i þ 1Þt þ U2i þ 1,2 sinð2iþ 1Þt: i¼0
(26)
iZ2
From Eq. (13), we have u2 ðtÞ ¼ a1,2 cosðtÞ þ a3,2 cosð3tÞ þ b3,2 sinð3tÞ þ a5,2 cosð5tÞ þ b5,2 sinð5tÞ:
(27)
Similarly, no secular terms cosðtÞ, cosð3tÞ, cosð5tÞ, sinðtÞ, sinð3tÞ and sinð5tÞ in Eq. (26) imply that their coefficients are all zeros, that is Z 2 p G2i þ 1,2 ðo2 ,a1,2 ,a3,2 ,b3,2 ,a5,2 ,b5,2 Þ ¼ F cosð2iþ 1Þt dt ¼ 0, i ¼ 0,1,2, p 0 2 Z p 2 U2i þ 1,2 ðo2 ,a1,2 ,a3,2 ,b3,2 ,a5,2 ,b5,2 Þ ¼ F2 sinð2iþ 1Þt dt ¼ 0, i ¼ 0,1,2: (28)
p
0
There are six linear algebraic equations for six unknown parameters o2, a1,2 a3,2, b3,2, a5,2 and b5,2, which can be easily solved with the help of symbolic software. Then the second-order improved harmonic balance approximations to the angular frequency and limit cycle can be obtained: uð2Þ ðtÞ ¼
2 X
a1,i cosðtÞ þ
i¼0
2 X
ða3,i cosð3tÞ þ b3,i sinð3tÞÞþ a5,2 cosð5tÞ þ b5,2 sinð5tÞ,
i¼1
oð2Þ ¼ o0 þ o1 þ o2 :
(29)
Substituting Eq. (29) into Eq. (26) yields the new residue defined by X X i F2 ¼ G7,2 cosð7tÞ þ U7,2 sinð7tÞ þ ½G2i þ 1,0 cosð2i þ1Þt þ U2i þ 1,0 sinð2i þ 1Þtpi3 ¼ R32 þ R2 pi3 : iZ4
(30)
iZ4
3.4. The third-order approximation Substituting Eqs. (12) and (30) into Eq. (15) and equating the coefficients of p3 yield
F3 ðu3 , o3 Þ ¼ : o20 u003 þ 2o0 o1 u002 þ ð2o0 o2 þw21 Þu001 þ 2ðo0 o3 þ o1 o2 Þu000 1 2=3 2 5=3 5 8=3 3 u0 u3 u0 u1 u2 þ u0 u1 þ eY þ R32 3 3 27 ! 3 X X þ ½G2i þ 1,3 cosð2iþ 1Þt þ U2i þ 1,3 sinð2i þ1Þt,
þ ¼
i¼0
with
(31)
iZ3
h
Y ¼ ol0 ðu20 1ÞDlt u3 þ2ðu1 u2 þ u0 u3 ÞDlt u0 þ2u0 u1 Dlt u2 þ ð2u0 u2 þ u21 ÞDlt u1 h i þ lo0l1 o1 ðu20 1ÞDlt u2 þ2u0 u1 Dlt u1 þ ð2u0 u2 þ u21 ÞDlt u0
i
1120
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
i l1 2 h 2 þ lo0l2 o0 o2 þ o1 ðu0 1ÞDlt u1 þ 2u0 u1 Dlt u0 2 ðl1Þðl2Þ 3 l3 2 o0 o3 þ ðl1Þo0 þ o3 ðu20 1ÞDlt u0 : þ lo0 6 Expanding u3(t) according to Eq. (13) yields u3 ðtÞ ¼ a1,3 cosðtÞ þa3,3 cosð3tÞ þb3,3 sinð3tÞ þ a5,3 cosð5tÞ þ b5,3 sinð5tÞ þ a7,3 cosð7tÞ þ b7,3 sinð7tÞ:
(32)
Similarly, no secular terms cosð2k þ 1Þt, sinð2k þ 1Þt, and ðk ¼ 0,1,2,3Þ in Eq. (31) imply that their coefficients are all zeros, that is Z 2 p G2i þ 1,3 ðo3 ,a1,3 ,a3,3 ,b3,3 ,a5,3 ,b5,3 ,a7,3 ,b7,3 Þ ¼ F cosð2iþ 1Þt dt ¼ 0, i ¼ 0,1,2,3, p 0 3 Z p 2 U2i þ 1,3 ðo3 ,a1,3 ,a3,3 ,b3,3 ,a5,3 ,b5,3 ,a7,3 ,b7,3 Þ ¼ F3 sinð2i þ 1Þt dt ¼ 0, i ¼ 0,1,2,3: (33)
p
0
There are eight linear algebraic equations for eight unknown parameters o3, a1,3 a3,3, b3,3, a5,3, b5,3, a7,3 and b7,3, which can be easily solved with the help of symbolic software. Then the third-order improved harmonic balance approximations to the angular frequency and limit cycle can be obtained: uð3Þ ðtÞ ¼
3 X
a1,i cosðtÞ þ
i¼0
3 X
ða3,i cosð3tÞ þ b3,i sinð3tÞÞ þ
i¼1
3 X
ða5,i cosð5tÞ þ b5,i sinð5tÞÞ þ a7,3 cosð7tÞ þ b7,3 sinð7tÞ,
i¼2
oð3Þ ¼ o0 þ o1 þ o2 þ o3 :
(34)
The procedure will be iteratively carried out to generate higher-order approximations by solving linear algebraic equations successively. 4. Results and discussion The classical van der Pol damped oscillator with fractional restoring force, governed by the nonlinear second differential Eq. (2), has been studied by Mickens using the first-order harmonic balance method [12] and averaging technique [14], and also by Lim and Lai [7] with the help of combined linearization and harmonic balance method. By Fourier expansion [12], we have ðcos tÞ1=3 ¼
1 X
a2n þ 1 cosð2n þ1Þt,
n¼0
with a2n þ 1 ¼
3Gð7=3Þ 24=3 Gðn þ 5=3ÞGð2=3nÞ
,
n ¼ 0,1,2. . .,
in which G is the Gamma function. Therefore, the first several terms are
cos 3t cos 5t 7 cos 7t cos 9t þ þ , ðcos tÞ1=3 ¼ w cos t 5 10 110 22
(35)
with
w¼
3Gð7=3Þ 24=3 Gð5=3ÞGð2=3Þ
¼
Gð1=3Þ 21=3 Gð2=3Þ2
1:159595266963929:
As is discussed in the third section, the residue for the zeroth-order approximation is separated into two parts. The first part conforms with the initial approximation and the second part comes from the nonlinearity of the function u1/3. In this case, infinite terms should be considered in the each stage of subsequent harmonic balance procedure. It is worthwhile to note that low accuracy of approximation is a direct consequence of discarding the second part. After substituting Eqs. (11) and (35) into Eq. (17), we find that the differential Eq. (18) is defined by
eol0 a10 3 3 lp 1=3 þ wa10 o0 2 a10 ¼ 0, eol0 a10 ¼ 0: (36) eol0 a10 a210 1 cos 4 2 4 The amplitude a10 of the zeroth-order approximation (19) is determined by the second equation in (36), which gives að0Þ ¼ a10 ¼ 2 while o0 depends on the solution of the first nonlinear equation in (36). Since higher-order approximations of the residue harmonic balance method are governed by linear equations only, solutions can easily be obtained to any desired accuracy. To show the effectiveness of the proposed method for fractional van der Pol oscillator (4), the frequency–amplitude response under the alternations of fractional order l and control parameter e will be illustratively validated by some figures, wherein the obtained results are compared with the numerical ones.
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
1121
4.1. Effect of fractional order on the frequency and amplitude of steady-state response The fractional order–frequency and fractional order–amplitude curves of limit cycles for e ¼1 are shown in Figs. 1 and 2, respectively. It is easy to find that the angular frequency of limit cycle in steady state increases at first, and then decreases as the fractional order varies in the interval (0,1]. However, the amplitude of limit cycle in steady state decreases at first, and then increases as the fractional order varies in the interval (0,1]. It is also clear that the obtained approximations are highly accurate and in good agreement with the numerically integrated results based on the phase portraits in Fig. 3 when e ¼ 1 for l ¼{0.5,0.8,1}. 4.2. Effect of control parameter on the frequency and amplitude of steady-state response When l ¼0.5, the fractional van der Pol oscillator (4) is reduced to a specific case: 1=3 € e½1uðtÞ2 D0:5 ¼ 0: uðtÞ t uðtÞ þuðtÞ
(37)
For this case, the comparisons between the asymptotic approximations to the angular frequency and amplitude of limit cycle and corresponding numerical results are illustrated by Figs. 4 and 5. Obviously, we can observe that the obtained approximations are highly accurate. Moreover, as the control parameter e varies in the interval (0,2], the angular frequency of the limit cycle are always increasing, while the amplitude decreases throughout the whole region of e. To further illustrate and verify the accuracy and effectiveness of the residue harmonic balance approximations, the comparison of the 1.7 1.6
Frequency of limit cycle
1.5 1.4 1.3 1.2 1.1 1 0.9 0.8 0.7 0
0.1
0.2
0.3
0.4 0.5 0.6 Fractinal order
0.7
0.8
0.9
1
Fig. 1. Comparisons of frequency between numerical solutions and results provided by the presented method, where solid lines and dashed lines denote the third-order and the zeroth-order approximations, respectively; dotted line denotes the numerical result. The control parameter is given as e ¼1.
2.25
Amplitude of limit cycle
2.2 2.15 2.1 2.05 2 1.95 1.9 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fractinal order Fig. 2. Comparisons of amplitude between numerical solutions and results provided by the presented method, where solid lines and dashed lines denote the third-order and the zeroth-order approximations, respectively; dotted line denotes the numerical result. The control parameter is given as e ¼1.
1122
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
3 2
Dλt u(t)
1 0 λ=1 λ=0.8 λ=0.5
-1 -2 -3 -2.5 -2
-1.5 -1
-0.5
0 u(t)
0.5
1
1.5
2
2.5
Fig. 3. Phase portraits of limit cycle, where solid lines and dotted lines denote the third-order analytical solution and numerical result, respectively. The control parameter and fractional derivative are given as e ¼ 1, l ¼{0.5,0.8,1}.
2.6
Frequency of limit cycle
2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0
0.2
0.4
0.6
1.2 1.4 0.8 1 Control parameter
1.6
1.8
2
Fig. 4. Comparisons of frequency between numerical solutions and results provided by the presented method, where solid lines and dashed lines denote the third-order and the zeroth-order approximations, respectively, dotted line denotes the numerical result. The fractional derivative is given as l ¼0.5.
2.04
Amplitude of limit cycle
2.02 2 1.98 1.96 1.94 1.92 1.9 0
0.2
0.4
0.6
0.8 1 1.2 1.4 Control parameter
1.6
1.8
2
Fig. 5. Comparisons of amplitude between numerical solutions and results provided by the presented method, where solid lines and dashed lines denote the third-order and the first-order approximations, respectively; dotted line denotes the numerical result. The fractional derivative is given as l ¼0.5.
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
1123
3 2 ε=1
D0.5 t u(t)
1
ε=2
ε=0.1
0 -1 -2 -3 -2.5 -2
-1.5 -1
-0.5
0 0.5 u(t)
1
1.5
2
2.5
Fig. 6. Phase portraits of limit cycle, where solid lines and dotted lines denote the third-order analytical solution and numerical result, respectively. The control parameter and fractional derivative are given as e ¼{0.1.1.2}, l ¼0.5.
Table 1 3rd-order analytical expressions of fractional vdP oscillator when l ¼0.5.
e
0.1
1
1.5
2
o
0.929607 2.02137 0.0217533 0.0205061 0.00658125 0.000356124 0.00218757 0.0000279856
1.50826 1.85721 0.0630316 0.0954636 0.000080662 0.00810594 0.00127273 0.000330784
1.77356 1.82287 0.0778078 0.110466 0.001301 0.0114029 0.00135085 0.000483207
2.01975 1.80172 0.086429 0.119604 0.00218938 0.0136131 0.00144587 0.000586111
cosðotÞ cosð3otÞ sinð3otÞ cosð5otÞ sinð5otÞ cosð7otÞ sinð7otÞ
phase portraits is presented in Fig. 6 for e ¼{0.1,1,2}. The analytical expressions of the limit cycle of oscillator (37) are tabulated in Table 1.
5. Conclusion The residue harmonic balance method to obtain the accurate approximations to the angular frequency and limit cycles of vibration systems is proposed. The fractional order van der Pol oscillator with fractional restoring force and damping is investigated successfully. It is shown that the approximations of frequencies, amplitudes and period solutions to the numerical ones are in excellent agreement. The effects of fractional order and the control parameter on the steady-state responses of angular frequency and amplitude of oscillator are illustrated by the several figures. It is worthwhile to point out that, for the fractional van der Pol oscillator with fractional restoring force, the residue, produced by the fractional power function, in the zeroth-order harmonic equation is separated into different stages and added to the subsequent harmonic procedures to improve the accuracy of the approximations. Therefore, the presented approach can be easily applied to other nonlinear oscillators with fractional power and derivative nonlinearities. Since the higher-order approximations of residue harmonic balance method are governed by linear equations, solutions can easily be obtained to any desired accuracy.
Appendix A. Newton iteration for differential equations The Newton method has been very successful to solve nonlinear algebraic equations. We extend the method to solve differential equations using the series here. When dealing with algebraic equations, one would expect quadratic convergence, that is, one iteration would double the accuracy of a root. When dealing with differential equations using series to express the solution function, one would expect to double the effective number of terms in the function expression. We shall illustrate the method in solving a simple ordinary differential equation of the form by power series: € u,u,tÞ _ Fðu, ¼ u€ þ u_ þ u2 t ¼ 0:
1124
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
Let the initial (zeroth) approximation be u0(t) ¼a0 þa1t þa2t2:
F0 ¼ Fðu€ 0 , u_ 0 ,u0 ,tÞ ¼ a0 2 þ a1 þ 2a2 þ ð2a2 þ 2a0 a1 1Þt þ ða21 þ 2a0 a2 Þt2 þ ð2a1 a2 Þt3 þ a22 t4 ¼ R0,0 þ R1,0 t þ R2,0 t 2 þ R3,0 t 3 þ R4,0 t 4 ¼
2 X
Ri,0 t i þ
i¼0
4 X
Ri,0 t i ¼
i¼3
2 X
Ri,0 t i þ R0 ðtÞ:
i¼0
P2 i _ R2,0 ¼ 0. The residue is separated in two parts: balanceable, e.g., Let R0,0 ¼ uð0Þ; R1,0 ¼ uð0Þ; i ¼ 0 Ri,0 t , and nonP4 i balanceable, e.g., i ¼ 3 Ri,0 t . There are 3 nonlinear equations in the balanceable part for the 3 unknowns a0, a1, a2, and u0(t) can be determined. However, the non-balanceable part, Ri,0a0; i ¼3,4. For R0,0 ¼0; R1,0 ¼0; R2,0 ¼0, one has 8 9 8 9 > < a0 > = > < 1:1040 > = 1:8367 : a ¼ a1 ¼ > > :a > ; > : ; 1:5278 2 Let the first approximation be u(1)(t)¼u0(t) þu1(t). u1 ðtÞ ¼ b0 þ b1 t þ b2 t 2 þ b3 t 3 þ b4 t 4 with coefficients bi to be determined. Substitute u(1) into
Fðu€ 0 þ u€ 1 , u_ 0 þ u_ 1 ,u0 þu1 ,tÞ ¼ ðu€ 0 þ u€ 1 Þ þ ðu_ 0 þ u_ 1 Þ þ ðu0 þ u1 Þ2 t ¼ 0: Retain only the linear terms in u1:
@ @ @ þ u_ 1 þ u1 F1 ¼ F0 þ u€ 1 þ u_ 1 þ 2u0 u1 ¼ R0 þ u€ 1 F90 ¼ R0 þ @u€ @u_ @u
2 X
uðkÞ 1
k¼0
! @ F90 ¼ R0 þD1 F90 ¼ 0 @uðkÞ
Here the number in the superscript brackets denotes the order of differentiation w.r.t. t. F90 indicates that the function is to be evaluated at u ¼u0 after partial differentiation:
F1 ¼ b1 þ 2b2 þ 2a0 b0 þ ð2b2 þ6b3 þ 2a0 b1 þ 2a1 b0 Þt þ ð3b3 þ12b4 þ 2a0 b2 þ2a1 b1 þ2a2 b0 Þt 2 þ ð4b4 þ 2a1 a2 þ 2a0 b3 þ2a1 b2 þ2a2 b1 Þt 3 þ ða22 þ 2b2 a2 þ 2a0 b4 þ 2a1 b3 Þt 4 þ ð2a1 b4 þ 2a2 b3 Þt 5 þ ð2a2 b4 Þt 6 ¼
4 X
Ri,1 t i þ R1 ðtÞ
i¼0
From the balanceable part, Ri,1 ¼0; i¼0, 1, 2, 3, 4, we have 5 linear equations for the 5 unknowns bi and u1(t) can be determined: 9 2 38 9 8 1 2 0 0 > b0 > > 0 > 2a0 > > > > > > > > >b > > 0 > 6 2a > 2 6 0 7 > > > > 1> > 6 1 2a0 7> = 6 7< = < 6 2a2 2a1 2a0 7 0 3 12 7 b2 þ ¼ 0: 6 > > > > >b > 6 0 > > > > 2a1 a2 > 4 7 2a2 2a1 2a0 > > > 4 5> 3> > > > > > > > > 0 0 2a2 2a1 2a0 : b4 ; : a22 ; Then bT ¼[0.5620, 0.7003, 0.9706, 0.2783, 0.1769]. Therefore, F1 ¼R1. However, the non-balanceable part Ri,1a0; i¼ 5,6. P Let the second approximation be u(2)(t)¼u0(t) þu1(t)þu2(t). u2 ðtÞ ¼ 6i ¼ 0 ci t i where ci are 7 constants to be determined. Substituted into the governing equation and retain only the linear terms in u2: 2 !2 3 2 2 6 X X @ 1 X ðkÞ @ 5F9 ¼ F2 ¼ R1 ðtÞ þ 4 uðkÞ þ u Ri,2 t i þ R2 ðtÞ ¼ 0, 2 1 0 ðkÞ ðkÞ 2! k ¼ 0 @u @u i¼0 k¼0
F2 ¼ ðb24 þ2a2 c6 Þt8 þ ð2a1 c6 þ 2a2 c5 þ 2b3 b4 Þt7 þ ðb23 þ2a2 b4 þ 2a0 c6 þ2a1 c5 þ 2a2 c4 þ 2b2 b4 Þt6 þ ð6c6 þ 2a1 b4 þ 2a2 b3 þ 2a0 c5 þ 2a1 c4 þ 2a2 c3 þ 2b1 b4 þ 2b2 b3 Þt 5 2
þðb2 þ 5c5 þ 30c6 þ 2a0 c4 þ 2a1 c3 þ2a2 c2 þ 2b0 b4 þ 2b1 b3 Þt 4 þ ð4c4 þ 20c5 þ 2a0 c3 þ 2a1 c2 þ 2a2 c1 þ2b0 b3 þ 2b1 b2 Þt 3 2
þ ðb1 þ3c3 þ 12c4 þ 2a0 c2 þ 2a1 c1 þ2a2 c0 þ 2b0 b2 Þt 2 þ ð2c2 þ 6c3 þ 2a0 c1 þ 2a1 c0 þ2b0 b1 Þt 2
þ b0 þc1 þ 2c2 þ 2a0 c0 ,
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
2
2a0
6 2a 6 1 6 6 2a2 6 6 6 6 6 6 6 4
1
2
2a0
2
2a1
2a0
6 3
2a2
2a1
2a0
4
20
2a2
2a1 2a2
2a0 2a1
5 2a0
2a2
2a1
12
1125
9 38 9 8 2 > c0 > > b0 > > > > > > > > > > > 7> > > > >c > > > 2b b 1 7> > > > 0 > > > > > > > 7> 2 > > > >c > > > 7> b1 þ 2b0 b2 > = = > < < 2> 7> 7 c 2b0 b3 þ 2b1 b2 ¼ 0: þ 3 7 > > 7> > > > > > 2 > > 30 7 b2 þ 2b0 b4 þ 2b1 b3 > > > > > c4 > > > > 7> > > > > > > >c > 6 7 > > > > 2a1 b4 þ 2a2 b3 þ 2b1 b4 þ2b2 b3 > > 5> 5> > > > > > > > ; > > : > 2 ; : 2a0 c6 b3 þ 2a2 b4 þ 2b2 b4
Let Ri,2 ¼ 0; i ¼ 0,1,. . .,6. We can determine ci and u2. By mathematical induction, one can prove that when the (r 1)th approximation is available, the rth approximation is given by uðrÞ ðtÞ ¼ u0 ðtÞ þ þur ðtÞ Let ur ðtÞ ¼
P2ðr þ 1Þ
i
ai,r t 2 2 X @ 1 Fr ¼ Rr1 ðtÞ þ 4 uðkÞ þ r ðkÞ 2! @u k¼0 i¼0
2 X k¼0
uðkÞ r1
@ @uðkÞ
!2
1 þ þ r!
2 X k¼0
uðkÞ 1
@ @uðkÞ
!r 3 2ðr þ 1Þ X 5F9 ¼ Ri,r t i þ Rr ðtÞ ¼ 0: 0 i¼0
Let the balanceable part be Ri,r ¼ 0; i ¼ 0,1,. . .,2ðr þ 1Þ. We can determine ai,r and ur. When we use sine and cosine functions instead of power series, we have the residue harmonic balance method. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]
B. van der Pol, A theory of the amplitude of free and forced triode vibrations, Radiology Review 1 (1920) 701–710; 754–762. B. van der Pol, On relaxation oscillations, Philosophical Magazine 7 (2) (1926) 978–992. J. Guckenheimer, K. Hoffman, W. Weckesser, Numerical computation of canards, International Journal of Bifurcation and Chaos 10 (2000) 2669–2687. Y.H. Kao, C.S. Wang, Analog study of bifurcation structures in a van der Pol oscillator with a nonlinear restoring force, Physical Review E 48 (4) (1993) 2514–2520. C. Letellier, V. Messager, R. Gilmore, From quasi-periodicity to toroidal chaos: analogy between the Curry–Yorke map and the van der Pol system, Physical Review E 77 (2008) 046203. R.E. Mickens, Fractional van der Pol equations, Journal of Sound and Vibration 259 (2003) 457–460. C.W. Lim, S.K. Lai, Accurate higher-order analytical approximate solutions non-conservative nonlinear oscillators and application to van der Pol damped oscillators, International Journal of Nonlinear Mechanics 48 (2006) 483–492. G.M. Mahmoud, A.A.M. Farghaly, Chaos control of chaotic limit cycles of real and complex van der Pol oscillators, Chaos Solitons and Fractals 21 (2004) 915–924. C.M.A. Pinto, J.A.T. Machado, Complex order van der Pol oscillator, Nonlinear Dynamics, doi: 10.1007/s11071-010-9886-0. Y.S. Chen, A.Y.T. Leung, Bifurcation and Chaos in Engineering, Springer-Verlag, London, 1998. R.E. Mickens, A generalized iteration procedure for calculating approximations to periodic solutions of ‘‘truly nonlinear oscillators’’, Journal of Sound and Vibration 287 (2005) 1045–1051. R.E. Mickens, Analysis of nonlinear oscillators having non-polynomial elastic terms, Journal of Sound and Vibration 255 (4) (2002) 789–792. R.E. Mickens, Iteration method solutions for conservative and limit-cycle x1/3 force oscillators, Journal of Sound and Vibration 292 (2006) 964–968. R.E. Mickens, A combined equivalent linearization and averaging perturbation method for non-linear oscillator equations, Journal of Sound and Vibration 264 (2003) 1195–1200. I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier, New York, 2006. N. Shimizu, W. Zhang, Fractional calculus approach to dynamic problems of viscoelastic materials, JSME International Journal Series C, Mechanical Systems, Machine Elements and Manufacturing 42 (4) (1999) 825–837. Y.A. Rossikhin, M.V. Shitikova, Application of fractional calculus for dynamic problems of solid mechanics: novel trends and recent results, Applied Mechanics Reviews 63 (2010) 010801–010852. I. Petras, Fractional Order Nonlinear Systems: Modeling, Analysis and Simulation, Springer-Verlag, London, 2011. G.M. Zaslavsky, Hamiltonian Chaos and Fractional Dynamics, Oxford University Press, Oxford, 2005. V.E. Tarasov, Fractional variations for dynamical systems: Hamilton and Lagrange approaches, Journal of Physics A: Mathematical and General 39 (2006) 8409–8425. V.E. Tarasov, Fractional systems and fractional Bogoliubov hierarchy equations, Physical Review E 71 (2005) 011102. N. Laskin, Fractional quantum mechanics, Physical Review E 62 (2000) 3135–3145. A.I. Saichev, G.M. Zaslavsky, Fractional kinetic equations: solutions and applications, Chaos 7 (1997) 753–764. R. Magina, M.D. Ortigueirab, Igor Podlubny, J. Trujillo, On the fractional signals and systems, Signal Processing 91 (3) (2011) 350–371. R.S. Barbosa, J.A.T. Machado, M.I. Ferreira, K.J. Tar, Dynamics of the fractional order Van der Pol oscillator, Proceedings of the Second IEEE International Conference on Computational Cybernetics (ICCC’04), Vienna University of Technology, Austria, 30 August–1 September 2004, pp. 373–378. T. Sardar, S.S. Ray, R.K. Bera, B.B. Biswas, The analytical approximate solution of the multi-term fractionally damped van der Pol equation, Physica Scripta 80 (2009) 025003–025006. F. Xie, X.Y. Lin, Asymptotic solution of the van der Pol oscillator with small fractional damping, Physica Scripta 136 (2009) 014033–014034. V. Gafiychuk, B. Datsko, V. Meleshko, Analysis of fractional order Bonhoeffer–van der Pol oscillator, Physica A 387 (2008) 418–424. M. Attari, M. Haeri, M.S. Tavazoei, Analysis of a fractional order van der Pol like oscillator via describing function method, Nonlinear Dynamics 61 (2010) 265–274. M.S. Tavazoei, M. Haeri, M. Attari, S. Bolouki, M. Siami, More details on analysis of fractional order Van der Pol oscillator, Journal of Vibration and Control 15 (6) (2009) 803–819. A. Konuralp, C. Konuralp, A. Yıldırım, Numerical solution to the van der Pol equation with fractional damping, Physica Scripta 136 (2009) 014034–014035. J.H. Chen, W.C. Chen, Chaotic dynamics of the fractionally damped van der Pol equation, Chaos Solitons and Fractals 35 (2008) 188–198. A.Y.T. Leung, Z.J. Guo, Forward residue harmonic balance for autonomous and non-autonomous systems with fractional derivative damping, Communications in Nonlinear Science and Numerical Simulations 16 (4) (2011) 2169–2218.
1126
A.Y.T. Leung et al. / Journal of Sound and Vibration 331 (2012) 1115–1126
[35] Z.J. Guo, A.Y.T. Leung, H.X. Yang, Oscillatory region and asymptotic solution of fractional van der Pol oscillator via residue harmonic balance technique, Applied Mathematical Modelling 35 (2011) 3918–3925. [36] A.Y.T. Leung, Z.J. Guo, Residue harmonic balance for discontinuous nonlinear oscillator with fractional power restoring force, International Journal of Nonlinear Science and Numerical Simulations 11 (9) (2010) 705–723. [37] E. Pereira, C.A. Monje, B.M. Vinagre, F. Gordillo, MATLAB toolbox for the analysis of fractional order systems with hard nonlinearities, Proceedings of the First IFAC Congress on Fractional Differentiation and its Applications, Bordeaux (France), July 2004, pp. 214–219.