Fractional Ince equation with a RiemannâLiouville fractional derivative
Applied Mathematics and Computation 219 (2013) 10695–10705
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Fractional Ince equation with a Riemann–Liouville fractional derivative Alfredo Parra-Hinojosa, Julio C. Gutiérrez-Vega ⇑ Photonics and Mathematical Optics Group, Tecnológico de Monterrey, Mexico 64849, Mexico
a r t i c l e
i n f o
Keywords: Ince equation Fractional calculus Ince polynomials Eigenvalue curves Stability Hill equation
a b s t r a c t We extend the classical treatment of the Ince equation to include the effect of a fractional derivative term of order a > 0 and amplitude c. A Fourier expansion is used to determine the eigenvalue curves aðÞ in function of the parameter , the stability domains, and the periodic stable solutions of the fractional Ince equation. Two important observations are the detachment of the eigenvalue curves from the a-axis in the ð; aÞ-plane, as well as the appearance of degenerate eigenvalues for suitable selections of the parameters. The fractional solutions, valid for the steady state of the system, are not orthogonal and have no defined parity. We also introduce a discrete numerical method to evaluate the Riemann–Liouville fractional derivative with lower terminal at 1 for a class of functions. The case a ¼ 1 represents the Ince equation with an additional constant damping. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction In recent years, the study of fractional calculus [1–3] has given rise to physical models that describe real-world phenomena in a more precise way. Several examples can be found in a vast range of scientific fields, such as optics [4], nonlinear dynamics [5], signal processing [6], and electrochemistry [7]. Sabatier et al. [8] presented a comprehensive compilation of 37 papers where fractional calculus has been applied to the areas of particle physics, diffusive systems, electrical systems, viscoelastic media and control, among others. In the same way, recent studies have dealt with the properties of diverse linear fractional differential equations [1,2]. These include, for example, fractional versions of the linear oscillator equation [9], the Duffing equation [10], the van der Pol equation [11], the wave equation [12], and the Mathieu equation [13]. Although fractional calculus was born well over 300 years ago, still many scientists and engineers nowadays are unaware of the properties and potential applications of this rich mathematical field. The canonical form of the ordinary Ince equation with eigenvalue a and parameters and p as presented by Arscott [14] is given by
"
# 2 d d þ sin 2t þ a p cos 2t NðtÞ ¼ 0: dt dt 2
ð1Þ
This equation and its different variants have appeared in the description of numerous physical problems in the fields evolutionary game theory [15], laser optics [16,17], quantum mechanics [18] and vibrational mechanics [19,20], among others. ⇑ Corresponding author. E-mail addresses: [email protected] (A. Parra-Hinojosa), [email protected] (J.C. Gutiérrez-Vega). URL: http://optica.mty.itesm.mx/pmog/ (J.C. Gutiérrez-Vega). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.04.044
10696
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
In this study, we extend the treatment of the ordinary Ince equation to include the effect of an additional fractional derivative term of order a 2 ð0; 2. We determine the eigenvalue curves and discuss the properties of the periodic solutions of this fractional Ince equation. This work also extends the results obtained for the closely related Mathieu equation [13] to which the Ince equation can be reduced [21]. 2. Statement of the problem We are interested in studying the ordinary Ince equation (1) with an additional fractional derivative term
"
# 2 d d a þ sin 2t þ a p cos 2t þ c D 1 t NðtÞ ¼ 0; dt dt 2
ð2Þ
where c 2 R is an amplitude constant and 1 Dat denotes the Riemann–Liouville fractional derivative of order a P 0 with lower limit b ¼ 1 n
a b Dt f ðtÞ ¼
1 d Cðn aÞ dt n
Z b
t
f ðnÞ ðt nÞ1þan
dn;
a P 0;
ð3Þ
with n being the smallest integer exceeding a, i.e., ðn 1Þ 6 a < n. The fractional Ince equation (3) has three special cases: (a) When c ¼ 0, it reduces to the ordinary Ince equation (1). (b) The Ince equation (1) is often regarded as a generalized Mathieu equation with periodic damping, i.e. the first derivative term. Therefore when a ¼ 1 and c – 0 the fractional term becomes a first-derivative damping term with constant amplitude c. The strengths of the periodic and the constant damping terms are adjusted by the constants and c, respectively. In analogy with the damped Mathieu equation [22], we will refer this special case as damped Ince equation. (c) When a ¼ 2 and c – 0, the fractional term becomes a second-derivative oscillatory term of amplitude c that can be added to the already existent second-derivative term. Thus, this case reduces to a scaled version of the ordinary Ince equation (1). The ordinary Ince equation (1) can be reduced to the canonical form of the Hill equation
y00 þ ½k þ Q ðtÞy ¼ 0;
ð4Þ
through the transformation
N ¼ y exp
4
cos 2t ;
ð5Þ
where k is the eigenvalue parameter and Q ðtÞ is a real p-periodic function of t [23]. For this case, we have
QðtÞ ¼ ðp þ 1Þ cos 2t þ ð2 =8Þðcos 4t 1Þ:
k ¼ a;
ð6Þ
From Floquet theory, it is known that Hill’s equation admits solutions of period p or 2p. The same holds for the Ince equation, since the transformation (5) preserves the periodicity and boundedness of the solution, and so the following theorem due to Liapounoff and Haupt applies to the Ince equation [23]: Theorem 1. To every differential equation (4) there belong two monotonically increasing infinite sequences of real numbers
k0 ;
k1 ;
k2 ; ;
ð7aÞ
k01 ;
k02 ;
k03 ;
ð7bÞ
and
k04 ; ;
such that Eq. (4) has a solution of period p if and only if k ¼ kn ; n ¼ 0; 1; 2; . . ., and a solution of period 2p if and only if k ¼ k0n ; n ¼ 1; 2; 3; . . .. The kn and k0n satisfy the inequalities
k0 < k01 6 k02 < k1 6 k2 < k03 < k04 < ;
ð8Þ
and the relations
lim k1 n!1 n
¼ 0;
lim ðk0n Þ1 ¼ 0;
ð9Þ
n!1
The solutions of Eq. (4) are stable in the intervals
ðk0 ; k01 Þ;
ðk02 ; k1 Þ;
ðk2 ; k03 Þ;
ðk04 ; k3 Þ; :
ð10Þ
At the endpoints of these intervals the solutions of Eq. (4) are, in general, unstable. This is always true for k ¼ k0 . The solutions of (4) are stable for k ¼ k2nþ1 or k ¼ k2nþ2 if and only if k2nþ1 ¼ k2nþ2 , and they are stable for k ¼ k02nþ1 or k ¼ k02nþ2 if and only if k02nþ1 ¼ k02nþ2 . We will investigate the effect of the fractional derivative term on the shape of the regions of stability on the parameter– eigenvalue ð; aÞ plane. These regions are delimited by the eigenvalue curves, on which periodic solutions exist.
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
10697
3. Periodic solutions of the fractional Ince equation m Before we obtain the solutions of the fractional Ince equation (2), we first recall that the even C m p ðt; Þ and odd Sp ðt; Þ solutions of order p 2 R and degree m ¼ 0; 1; 2; . . .., of the ordinary Ince equation (1) (i.e. the fractional Ince equation with c ¼ 0) can be expanded in terms of Fourier series. The expansions fall into four classes, according to their symmetry or antisymmetry, about t ¼ 0 and t ¼ p=2, namely [14]
C 2k p ðt; Þ ¼
1 X Ar cos 2rt;
k ¼ 0; 1; . . . ;
r¼0 1 X
ðt; Þ ¼ C 2kþ1 p
Ar cosð2r þ 1Þt;
ð11aÞ
k ¼ 0; 1; . . . ;
ð11bÞ
r¼0
S2k p ðt; Þ ¼
1 X Br sin 2rt;
k ¼ 1; 2; . . . ;
r¼1 1 X
ðt; Þ ¼ S2kþ1 p
Br sinð2r þ 1Þt;
ð11cÞ
k ¼ 0; 1; . . .
ð11dÞ
r¼0 m
1
2
0 1 2 The respective eigenvalues am p ðÞ and bp ðÞ satisfy the relation ap < bp < ap < bp < ap for given values of p and m m that each function is associated with an eigenvalue ap or bp . They satisfy the normalization
Z p
p
2 ½C m p ðt; Þ dt ¼
Z p
p
2 ½Sm p ðt; Þ dt ¼ p
> 0, so ð12Þ
and the orthogonality relation [14]
m0 exp cos 2t C m p ðt; ÞC p ðt; Þdt ¼ 0; 2 p
Z p
m – m0 :
ð13Þ
If p is a non-negative integer the expressions in (11) become finite sums, and the Ince functions become the Ince polynomials, namely
C 2k 2n ðt; Þ ¼
n X Ar cos 2rt;
k ¼ 0; 1; . . . ; n;
ð14aÞ
r¼0
C 2kþ1 2nþ1 ðt; Þ ¼
n X Ar cosð2r þ 1Þt;
k ¼ 0; 1; . . . ; n;
ð14bÞ
r¼0
S2k 2n ðt; Þ ¼
n X Br sin 2rt;
k ¼ 1; 2; . . . ; n;
ð14cÞ
r¼1
S2kþ1 2nþ1 ðt; Þ ¼
n X Br sinð2r þ 1Þt;
k ¼ 0; 1; . . . ; n:
ð14dÞ
r¼0
Even order p ¼ 2n and odd order p ¼ 2n þ 1 Ince polynomials are p- and 2p-periodic, respectively. To each Ince polynomial 2k 2kþ1 2kþ1 we associate the corresponding finite sequences of real eigenvalues a2k 2n ; b2n ; a2nþ1 , and b2nþ1 [17]. We now look for periodic solutions of the fractional Ince equation (2). Consider a formal solution in the form of a Fourier series
NðtÞ ¼ A0 þ
1 X Ar cos rt þ Br sin rt:
ð15Þ
r¼1
The choice of setting b ¼ 1 in Eq. (3) allows the solution to be periodic. This choice is crucial and illustrates the fact that the fractional derivative is a nonlocal operator, since it depends on the values between the lower limit b and the point of interest [24]. The Ince equation often appears in problems where time is the independent variable. The lower limit must then be set to the initial value of t, namely, t ¼ 0. For this case, since we are interested in looking for periodic solutions of Eq. (2), we make use of the following results provided by Tavazoei [25] on the fractional derivative of periodic functions with lower limit b ¼ 0: Theorem 2. Suppose that f ðtÞ is ðn 1Þ-times continuously differentiable and f ðnÞ ðtÞ is bounded. If f ðtÞ is a non-constant periodic function with period T, the function 0 Dat f ðtÞ, where 0 < a R N and n is the first integer greater than a, cannot be a periodic function of period T. Corollary 1. A differential equation of fractional order in the form a
0 Dt
f ðtÞ þ W½f ðtÞ; f 1 ðtÞ; . . . ; f ðnÞ ðtÞ ¼ 0;
ð16Þ
10698
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
where 0 < a R N, cannot have any non-constant smooth periodic solution. For this problem, f ðtÞ is either a sine or a cosine function. Indeed, their fractional derivatives in the Riemann–Liouville sense are not periodic, and they are given by [1] a
cos rt ¼ ðrtÞa E2;1a ðr 2 t 2 Þ;
ð17aÞ
a
sin rt ¼ ðrtÞ1a E2;2a ðr 2 t 2 Þ;
ð17bÞ
0 Dt 0 Dt
where Ea;b ðzÞ is the two-parameter Mittag–Leffler function [1]
Ea;b ðzÞ ¼
1 X
zk ; C ð a k þ bÞ k¼0
a>0
ð18Þ
b>0
Nevertheless, for large values of t, Eq. (17) converge to the following periodic functions: a
cos rt ¼ ra cosðrt þ ap=2Þ;
t ! 1;
ð19aÞ
a
sin rt ¼ r a sinðrt þ ap=2Þ;
t ! 1:
ð19bÞ
0 Dt 0 Dt
Mathematically, setting b ¼ 1 in Eq. (3) gives the same result, i.e., a
1 Dt
cos rt ¼ r a cosðrt þ ap=2Þ;
ð20aÞ
a
a 1 Dt sin rt ¼ r sinðrt þ ap=2Þ:
ð20bÞ
In more physical terms, the solutions are valid for the steady state of the system in question. Also, for this choice, the Riemann–Liouville and Caputo definitions of the fractional derivative are equal [1]. Inserting (15) into (2) and collecting terms we obtain the following recurrence relations:
aA0 ð=2Þðp þ 2ÞA2 ¼ 0;
ð21aÞ
pA0 þ ½a 4 þ 2a c cosðap=2ÞA2 ð=2Þðp þ 4ÞA4 þ 2a c sinðap=2ÞB2 ¼ 0;
ð21bÞ
½a 4 þ 2a c cosðap=2ÞB2 ð=2Þðp þ 4ÞB4 2a c sinðap=2ÞA2 ¼ 0;
ð21cÞ
½a 1 ð=2Þðp þ 1Þ þ c cosðap=2ÞA1 ð=2Þðp þ 3ÞA3 þ c sinðap=2ÞB1 ¼ 0;
ð21dÞ
½a 1 þ ð=2Þðp þ 1Þ þ c cosðap=2ÞB1 ð=2Þðp þ 3ÞB3 c sinðap=2ÞA1 ¼ 0:
ð21eÞ
For r P 3
ð=2Þðr 2 pÞAr2 þ ½a r 2 þ r a c cosðap=2ÞAr ð=2Þðr þ 2 þ pÞArþ2 þ r a c sinðap=2ÞBr ¼ 0;
ð21fÞ
ð=2Þðr 2 pÞBr2 þ ½a r 2 þ r a c cosðap=2ÞBr ð=2Þðr þ 2 þ pÞBrþ2 r a c sinðap=2ÞAr ¼ 0;
ð21gÞ
a
where we used 1Dt A0 ¼ 0 for a > 0. Note that the Ar and Br coefficients are now coupled, and there are two independent systems, one involving only even subscripts (A0 ; B2 ; A2 ; . . .) and the other involving odd subscripts (B1 ; A1 ; B3 ; . . .). This means that there exist only two kinds of solutions, instead of four, as is the case of the ordinary Ince equation aÞ N2k;ð ðt; Þ ¼ A0 þ p
1 X A2r cos 2rt þ B2r sin 2rt;
k ¼ 0; 1; 2; . . . ;
ð22aÞ
r¼1
Np2kþ1;ðaÞ ðt; Þ ¼
1 X A2rþ1 cosð2r þ 1Þt þ B2rþ1 sinð2r þ 1Þt;
k ¼ 0; 1; 2; . . .
ð22bÞ
r¼0
The first solution (22a) is p-periodic, whereas (22b) is 2p-periodic, but none has a defined parity. We label the eigenvalues of 2k;ðaÞ
aÞ N 2k;ð ðt; Þ as ap p
2kþ1;ðaÞ
aÞ ðÞ, and those of N 2kþ1;ð ðt; Þ as ap p 0;ðaÞ ap
2;ðaÞ ap
4;ðaÞ ap
ðÞ, where the superscript ðaÞ refers to the order of the fractional 1;ðaÞ
3;ðaÞ
5;ðaÞ
< < and ap < ap < ap for > 0, and are normalized in the same derivative. They follow the order way as the transcendental Ince functions Eqs. (12). It is not difficult to show that, when p is a non-negative integer, both series (22) become finite, and thus the fractional Ince functions (22) become fractional Ince polynomials. In particular, for solution (22a) (to which we associate the recurrence relations (21) with even coefficients), let p be a positive integer (p ¼ 2n), and let us choose the parameter a such that A2nþ2 ¼ B2nþ2 ¼ 0. Then, for r ¼ 2n þ 2, Eqs. (21f) and (21g) give A2nþ4 ¼ A2nþ6 ¼ ¼ 0 and B2nþ4 ¼ B2nþ6 ¼ ¼ 0 respectively, which means that the solution is finite (ending in A2n and B2n ). This condition results in the following set of recurrence relations:
aA0 þ ðn þ 1ÞA2 ¼ 0;
ð23aÞ
2nA0 þ ½4 a 2a c cosðap=2ÞA2 þ ðn þ 2ÞA4 2a c sinðap=2ÞB2 ¼ 0;
ð23bÞ
½4 a 2a c cosðap=2ÞB2 þ ðn þ 2ÞB4 þ 2a c sinðap=2ÞA2 ¼ 0:
ð23cÞ
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
10699
For r ¼ 4; 6; . . . ; 2n 2,
ðn r=2 þ 1ÞAr2 þ ½r 2 a ra c cosðap=2ÞAr þ ðn þ r=2 þ 1ÞArþ2 r a c sinðap=2ÞBr ¼ 0;
ð23dÞ
ðn r=2 þ 1ÞBr2 þ ½r 2 a r a c cosðap=2ÞBr þ ðn þ r=2 þ 1ÞBrþ2 þ cra sinðap=2ÞAr ¼ 0:
ð23eÞ
Finally,
A2n2 þ ½4n2 a ð2nÞa c cosðap=2ÞA2n ð2nÞa c sinðap=2ÞB2n ¼ 0; B2n2 þ ½4n2 a ð2nÞa c cosðap=2ÞB2n þ ð2nÞa c sinðap=2ÞA2n ¼ 0:
ð23fÞ ð23gÞ
There are 2n þ 1 equations for the ð2n þ 1Þ coefficients A2r and B2r . The case c ¼ 0 (or a ¼ 0), B2r ¼ 0 reduces to the recurrence relations reported by Arscott in [14]. A similar analysis can be carried out for the 2p-periodic solution with p ¼ 2n þ 1, which gives the following relations:
½a 1 ðn þ 1Þ þ c cosðap=2ÞA1 ðn þ 2ÞA3 þ c sinðap=2ÞB1 ¼ 0;
ð24aÞ
½a 1 þ ðn þ 1Þ þ c cosðap=2ÞB1 ðn þ 2ÞB3 c sinðap=2ÞA1 ¼ 0:
ð24bÞ
For r ¼ 3; 5; . . . ; 2n 1,
ð2n r þ 3Þð=2ÞAr2 þ ½r 2 a cra cosðap=2ÞAr þ ð2n þ r þ 3Þð=2ÞArþ2 cra sinðap=2ÞBr ¼ 0;
ð24cÞ
ð2n r þ 3Þð=2ÞBr2 þ ½r 2 a cr a cosðap=2ÞBr þ ð2n þ r þ 3Þð=2ÞBrþ2 þ cra sinðap=2ÞAr ¼ 0:
ð24dÞ
Finally,
A2n1 þ ½ð2n þ 1Þ2 a cð2n þ 1Þa cosðap=2ÞA2nþ1 cð2n þ 1Þa sinðap=2ÞB2nþ1 ¼ 0; B2n1 þ ½ð2n þ 1Þ2 a cð2n þ 1Þa cosðap=2ÞB2nþ1 þ cð2n þ 1Þa sinðap=2ÞA2nþ1 ¼ 0:
ð24eÞ ð24fÞ
Finally, the expansions of the fractional Ince polynomials are given by 2k;ðaÞ
N2n
ðt; Þ ¼ A0 þ
n X A2r cos 2rt þ B2r sin 2rt;
k ¼ 0; 1; . . . ; n;
ð25aÞ
r¼1
2kþ1;ðaÞ
N2nþ1
ðt; Þ ¼
n X A2rþ1 cosð2r þ 1Þt þ B2rþ1 sinð2r þ 1Þt;
k ¼ 0; 1; . . . ; n;
ð25bÞ
r¼0 2k;ðaÞ
2kþ1;ðaÞ
with eigenvalues a2n ðÞ and a2nþ1 ðÞ, and where the Fourier coefficients are calculated with the system (24)?. The fractional Ince solutions are normalized in the same way as the ordinary Ince solutions. 4. Behavior of the eigenvalue curves An important problem associated with the fractional Ince equation is the determination of the eigenvalues aðÞ, for which the solutions are periodic in t. To make suitable comparisons, in Fig. 1 we first plot the eigenvalue curves aðÞ and bðÞ for the m ordinary Ince equation. Thick and thin lines denote the characteristic curves for even C m p ðt; Þ and odd Sp ðt; Þ Ince functions, respectively. Even-order and odd-order curves are symmetrical and antisymmetrical about the axis a, respectively. Except for the trivial case ¼ 0, the curves do not intersect and thus divide the ð; aÞ plane in regions where the solutions of the Ince equation are stable (white regions) or unstable (gray regions). For referring purposes, the unstable regions are numbered sequentially, as shown in Fig. 2. Note that the jth unstable region touches the axis q ¼ 0 at values a ¼ jjj2 . To get insight into the eigenvalue curves of the fractional Ince equation, let us analyze first the simple case when p is an even integer, say p ¼ 2n ¼ 2. The p-periodic solution (22a) [which fulfills the system (23)] is finite and has the form (25a), 2k;ðaÞ namely N 2 ¼ A0 þ B2 sin 2t þ A2 cos 2t, whereas the 2p-periodic solution (22b) has an infinite number of terms whose coefficients fulfill the odd system in Eq. (21). Recurrence relations (23) can be written in matrix form as follows:
2
32 3 2 3 0 0 2 A0 A0 6 76 7 6 7 a a 0 4 2 c cosð a p =2Þ 2 c sinð a p =2Þ B ¼ a 4 54 2 5 4 B2 5: a a 2 2 c sinðap=2Þ 4 2 c cosðap=2Þ A2 A2
ð26Þ
The matrix in Eq. (26) (denoted by M) is the matrix representation of the fractional differential operator ½d=dt 2 sin 2td=dt þ p cos 2t c1 Dat in the Fourier basis. The eigenvalues of the system are found by setting detðM aIÞ ¼ 0, which gives a relation between and a that defines the eigenvalue curves as a function of the parameter for fixed values of c and a. This consequence is drawn from Floquet theory, which states that periodic solutions (of period p or 2p) exist on the ð; aÞ curves that separate regions of stable and unstable solutions.1 For this case, the characteristic equation is given by the expression 1
A solution is stable if every independent solution is bounded.
10700
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
25
20
4
4
4 8
4 8
a8 :(C8 )
-4
b :(S )
15
10
a
3
3
3
5
4
4 8
4
b :(S )
3
b8 :(S8 ) 2
3
3
3
3
a8 :(C8 )
a8 :(C8 )
-3
4
4 8
a8 :(C8 )
Stable
b8 :(S8 ) 2
2
3
2
a8 :(C8 )
a8 :(C8 )
Unstable 2
b82 :(S82 )
0
1 8
a :(C )
1
b8 :(S8 )
a :(C )
-2
1
b8 :(S8 )
-5
a80 :(C80 )
-1 -10 -5
-4
-3
2
b8 :(S8 )
1 8
-2
0
-1
1 8
1 8
1
1
2
a80 :(C80 )
0
1
1
2
3
4
5
ε m Fig. 1. Eigenvalue curves of the ordinary Ince equation (c ¼ 0) with p ¼ 8. Thick and thin lines denote curves for even C m p ðt; Þ and odd Sp ðt; Þ Ince functions, respectively. The characteristic curves divide the plane into regions of stability (white) and instability (gray). The unstable regions are numbered sequentially for referring purposes.
25
20
-4
a88,(0.5):(N88,(0.5))
a88,(0.5):(N88,(0.5))
Stable
a86,(0.5):(N86,(0.5))
4
a86,(0.5):(N86,(0.5))
15 a85,(0.5):(N85,(0.5))
10
-3
a
a87,(0.5):(N87,(0.5))
5
a
4,(0.5) 8
a
0
-2
a87,(0.5):(N87,(0.5))
4,(0.5) 8
:(N
2,(0.5) 8
:(N
a
)
2,(0.5) 8
) a
0,(0.5) 8
:(N
4,(0.5) 8
4,(0.5) 8
2,(0.5) 8
2,(0.5) 8
) a
0,(0.5) 8
:(N
:(N
) )
Unstable
2
-5
a81,(0.5):(N81,(0.5))
-4
-3
a83,(0.5):(N83,(0.5)) a83,(0.5):(N83,(0.5))
a81,(0.5):(N81,(0.5))
-1 -10 -5
3
a85,(0.5):(N85,(0.5))
1
0
-2
-1
0
1
2
3
4
5
ε Fig. 2. Eigenvalue curves for p ¼ 8 with fractional coefficient c ¼ 2 and half-order a ¼ 0:5. The curves originally crossing the axis at a ¼ 0; 1; 4; 9; . . ., become detached. Also, some curves exhibit crossings, and the crossing points indicate double eigenvalues.
a½ð4 aÞ2 2ð4 aÞ2a c cosðap=2Þ þ 22a c2 þ 42 ½4 a 2a c cosðap=2Þ ¼ 0;
ð27Þ
which is cubic in a and whose solutions give three curves a ¼ aðÞ. Setting ¼ 0 in (27), the intersection points in the a axis can be obtained. The only real solution is a ¼ 0, and the other two solutions correspond to a ¼ 4 2a c expðiap=2Þ. Thus, there is only one intersection when c and a are not zero, and the curves originally intersecting the axis at a ¼ 4 (as shown in Fig. 1) become ‘‘detached’’ from the vertical axis. This detachment of the eigenvalue curves may be anticipated from the Mathieu equation when a fractional derivative term is present [13]. The next finite p-periodic solution (namely, n ¼ 2; p ¼ 4, 2k;ðaÞ so that N 4 ¼ A0 þ B2 sin 2t þ A2 cos 2t þ B4 sin 4t þ A4 cos 4t) gives a third intersection when c ¼ 0 (for a ¼ 16), which also disappears when c and a are not zero.
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
10701
40 6,(1.75)
a116,(1.75):(N11
7,(1.75)
a117,(1.75):(N11
)
30
)
4 a115,(1.75):(N115,(1.75))
20
4,(1.75)
a114,(1.75):(N11
a
)
10 2,(1.75)
a112,(1.75):(N11
)
3
0 3,(1.75)
a110,(1.75):
a113,(1.75):(N11
-10 (N110,(1.75))
0
)
a111,(1.75):(N111,(1.75))
0
2
1 5
10
15
ε Fig. 3. Eigenvalue curves for p ¼ 11 with fractional coefficient c ¼ 2 and fractional order a ¼ 1:75. The crossings that appear for a 2 ð0; 1Þ (c > 0) are no longer present.
In general the effect of the fractional term on the eigenvalue curves is illustrated in Fig. 2 for p ¼ 8; c ¼ 2, and a ¼ 0:5. Note that the eigenvalue curves that arise from the 2p-periodic solutions have also been included. Given that these solutions are not finite for p – 2n þ 1, the way to obtain the eigenvalue curves consists in truncating the series solution (22b) and using recurrence relations (21), which is the general procedure for any p 2 R. The larger the order of truncation, the more precise the curves are. Here, a 2p-periodic solution with 52 terms was used. The curves divide the ð; aÞ-plane into regions where the solutions of the fractional Ince equation are stable (white region) or unstable (gray region). The unstable regions are numbered sequentially following the same convention as in Fig. 1. Observation of Fig. 2 gives the following conclusions for c > 0 and a 2 ð0; 1Þ: 0;ð0:5Þ
Except for the lowest-order curve a8 ðÞ, each eigenvalue curve of the fractional Ince equation becomes ‘‘detached’’ from the a-axis. The upper and lower boundaries of the unstable regions correspond to the original curves aðÞ and bðÞ of the ordinary Ince equation and are depicted using different line-widths in Fig. 2. In the ‘‘detached’’ (i.e., stable) region, the eigenvalues become complex, and they appear in complex conjugate pairs (i.e., 2nþ1;ðaÞ 2nþ3;ðaÞ 2nþ2;ðaÞ 2nþ4;ðaÞ ap ¼ ap and ap ¼ ap ; n even). The upper and lower eigenvalue curves of an unstable region connect smoothly forming a single continuous boundarycurve. For > 0, the turning-point of the boundary-curve is the minimum value of for which the eigenvalue a is real. For each boundary-curve, the solutions above and below the turning-points are labeled with different superscripts. Periodic solutions corresponding to different points ð; aÞ lying on the same boundary-curve have the same periodicity, namely, p and 2p for boundary-curves of even and odd unstable regions, respectively. For any given value of , there is only a finite number of real eigenvalues a, which increases as increases. Consequently, unlike the ordinary Ince equation, the fractional Ince equation has a finite number of real periodic solutions for a given value of . There exist real values where the eigenvalue curves intersect, indicating the existence of degenerate eigenvalues (i.e., double points) a of the fractional Ince equation. In particular, for the choice a ¼ 0:5; c ¼ 2, the crossings occur for all pairs 2n;ðaÞ 2nþ1;ðaÞ 4n;ðaÞ 4nþ2;ðaÞ 4nþ1;ðaÞ 4nþ3;ðaÞ of curves fap ; ap g when > 0, and for the curves fap ; ap g and fap ; ap g when < 0 (n ¼ 0; 1; . . .). The adjacent curves become asymptotically parallel after the crossing, therefore they only cross once. Additionally, the effect of increasing p is to reduce the area of stability (the curves appear less detached), which means that the number of real eigenvalues a for a given value of increases. Nevertheless, for a 2 ½1; 2, we observed that eigenvalue curves do not intersect, i.e., the degenerate eigenvalues disappear, as shown in Fig. 3. As expected, when a ¼ 2, we recover the ordinary Ince equation (1) scaled by a factor 1 þ c and the curves become ‘‘reattached’’, as shown in Fig. 1. 5. Behavior of the solutions of fractional order Since the fractional Ince equation involves two additional parameters (c and a), a general description of the behavior of the eigenfunctions can be fairly complicated. In this section we shall illustrate the effect of the fractional term on the shape
10702
(a)
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
1
(b)
C13(t,4)
1.5
S33(t,3) 1
0.5 0.5
C13(t,0)=cos(t)
0
3 S11 (t,3)
0
-0.5
-0.5
-1 0
π/4
π/2
3π/4
0
π
π/4
t
(c)
(d)
1
0.5
0.5
0
0
-0.5
-0.5
π/2
π t
3π/2
2π
α=0 α=0.5 α=1
-1
α=0 α=0.5 α=1
0
π
1.5
1
-1.5
3π/4
t
1.5
-1
π/2
3π/2
2π
-1.5
0
π/2
π t
3;ðaÞ
Fig. 4. Plots of several Ince polynomials and fractional solutions. (a) C 13 ðt; Þ; 2 f0; 1; 2; 3; 4g; (b) S3p ðt; 3Þ; p 2 f3; 5; 7; 9; 11g; (c) N 3 5;ðaÞ reduces to C 13 ðt; 11Þ when a ¼ 0; (d) N 3 ðt; 11Þ with c ¼ 1. It reduces to S33 ðt; 11Þ when a ¼ 0.
ðt; 11Þ with c ¼ 1. It
2;ðaÞ
Fig. 5. N 6:5 ðt; 5Þ eigenfunction of the fractional Ince equation for c ¼ 2; a 2 ½0; 2. Solid curves correspond to a ¼ f0; 0:5; 1; 1:5; 2g.
of the solutions with a few examples. First, in Fig. 4(a) and (b) we show the ordinary Ince polynomials C 13 ðt; Þ and S3p ðt; 3Þ. 3;ðaÞ 3;ðaÞ Next, we consider the fractional Ince polynomial N 3 ðt; 11Þ with eigenvalue a3 . For c ¼ 1 and a ¼ f0:5; 1g, recurrence
10703
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
relations (23) can be used to obtain the Fourier coefficients and plot the functions, which reduce to C 13 ðt; 11Þ when a or c 5;ðaÞ equal to zero (Fig. 4(c)). The same procedure for N 3 ðt; 11Þ is illustrated in Fig. 4(d). 2;ðaÞ In Fig. 5 we show the behavior of the p-periodic normalized function N 6:5 ðt; 5Þ in the range a 2 ½0; 2. When a vanishes, the fractional Ince equation reduces to the ordinary Ince equation with a shifted eigenvalue a þ c, and the eigenfunction becomes the Ince function S26:5 ðt; 5Þ. Similarly, when a ¼ 2, the fractional Ince equation reduces again to the ordinary Ince equa2;ð2Þ tion but now scaled by 1 þ c, and the solution N 6:5 ðt; 5Þ becomes S26:5 ðt; 5=ð1 þ cÞÞ. In all the examples shown, when the order a varies, the periodicity of the solutions is preserved, while the amplitude is altered and zeros of the polynomials are slightly shifted, except for a ¼ 0; 2. Likewise, the solutions have no defined parity, as opposed to the non-fractional Ince functions. 6. Behavior of the degenerate solutions of fractional order Since the fractional Ince equation has a periodic solution for each point ð; aÞ lying on the eigenvalue curves in Fig. 2, the cross-points between adjacent curves define degenerate Ince functions with the same values of and a. Let us next illustrate the behavior of the degenerate solutions. Fig. 6(a) shows the first two eigenvalue curves aðÞ for a ¼ 0:6; p ¼ 5:5, and c ¼ 2. The crossing points at the values ¼ n indicate degenerate eigenvalues an , namely aÞ a2n2;ð ðn Þ ¼ ap2n1;ðaÞ ðn Þ an ; p
n ¼ 0; 1; 2; . . .
ð28Þ
These are marked with asterisks on the curves in Fig. 6(a) and their numerical values are given up to ten decimals of precision in Fig. 6(b).
(a)
1 0.5 0
(b)
-0.5
a
-1
(ε1*, a*) 1
n=1 n=2
-1.5
Crossing points ε*n 0.60919851809 2.8966572868
a*n −1.4509203434 −2.4053344921
-2
(ε2*, a*) 2
-2.5 0
(c)
0.5
1
1.5
ε
2
2.5
3
3.5
1.5
(d)
1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
-1.5
0
π/2
π
t
3π/2
2π
-1.5
0
π/2
π
3π/2
2π
t
Fig. 6. (a) Eigenvalue curves for a ¼ 0:6; c ¼ 2. Asterisks indicate the crossing points between adjacent curves. (b) Numerical values of the first two crossing 0;ð0:6Þ 1;ð0:6Þ points. (c) Eigenfunctions corresponding to the first degenerate eigenvalue a1 . N 5:5 ðt; 1 Þ (red-circle line) and N 5:5 ðt; 1 Þ (blue-point line); (d) 2;ð0:6Þ 3;ð0:6Þ Eigenfunctions corresponding to the second degenerate eigenvalue a2 . N 5:5 ðt; 2 Þ (red-circle line) and N 5:5 ðt; 2 Þ (blue-point line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
10704
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
The behavior of the first and second pairs of degenerate fractional Ince functions is depicted in Fig. 6(c) and (d). Consider, for 0;ð0:6Þ 1;ð0:6Þ instance, the pair ð1 ; a1 Þ. The degenerate functions N 5:5 ðt; 1 Þ and N 5:5 ðt; 1 Þ satisfy the same differential equation, namely
"
# 2 d d f0;1g;ð0:6Þ 0:6 þ 1 sin 2t þ a1 5:51 cos 2t þ 21 Dt N5:5 ¼ 0: dt dt 2
ð29Þ
Note that the even-order functions (red-circle lines) are p-periodic whereas odd-order functions (blue-point lines) are 2pperiodic. In general, in any pair of degenerate solutions one of them has period p and the other 2p. In particular, the p-periodic functions (red-circle lines) shown in Fig. 6(c) and (d) belong to the same boundary-curve of the second unstable region. 7. Numerical considerations In order to test our solutions of the fractional Ince equation we implemented accurate numerical routines based on suitable finite differences schemes. Firstly, solution of the Ince equation (1) and the fractional Ince equation (2) with integer a derivatives can be performed using standard numerical methods. However, the numerical solution of the fractional Ince equation (2) with arbitrary a is a more complicated task that deserve some discussion. For 0 < a < 1 the fractional Ince functions were verified by calculating the fractional derivative (3) with an extended scheme based on the following discrete form of (3) with lower limit b ¼ 0 due to Taukenova and Shkhanukov-Lafishev [26]:
i Xh 1 ðt j t s1 Þ1a ðt j t s Þ1a ft;s ; Cð2 aÞ s¼1 j
a
0 Dt j f ðt j Þ
¼
ð30Þ
where t j ¼ jh; j ¼ 0; 1; . . . ; j0 ; ft;s ¼ ðfs fs1 Þ=h, such that a
0 Dt j f
¼
1 Cð1 aÞ
Z 0
tj
f_ ðnÞðt j nÞa dn ¼ 0 Datj f þ OðhÞ:
ð31Þ
This scheme is obtained in the following way, assuming that f is Lipschitz continuous. Then
1 Cð1 aÞ
Z
tj 0
j Z ts j Z ts _ X X f ðts1=2 Þ þ €f ðgÞðn t s1=2 Þdn f_ ðnÞdn f_ ðnÞdn 1 1 a ¼ a ¼ C ð1 a Þ C ð1 a Þ ðt j nÞ ðt j nÞ ðt j nÞa s¼1 ts1 s¼1 t s1 j j Z ts € X X f ðgÞðn ts1=2 Þdn 1 1 2 ¼ ½ðt j t s1 Þ1a ðt j t s Þ1a ft;s þ þ Oðh Þ; ð32Þ Cð2 aÞ s¼1 Cð1 aÞ s¼1 ts1 ðt j nÞa
where g is an intermediate point between n and t s1=2 . Now, since j€f ðnÞj 6 M, we find that the second term in (32) is of order OðhÞ, and we obtain (31). To extend this scheme with the lower limit at 1, set t0 ! 1, and a numerical value for t0 is to be obtained. We have a
1 Dt j f ðt j Þ
¼
1 Cð1 aÞ
Z
tj
f_ ðnÞðt j nÞa dn ¼
1
1 Cð1 aÞ
Z
0
f_ ðnÞðt j nÞa dn þ
Z
tj
f_ ðnÞðt j nÞa dn :
ð33Þ
0
1
f ðtÞ being periodical, the infinite limit can be discretized as an integer number of periods nT, so that
Z
0
f_ ðnÞðt j nÞa dn ¼ lim
n!1
1
We call
Ik ¼
Z
kT
Z
0
f_ ðnÞðt j nÞa dn ¼ lim
n!1
nT
n1 Z X k¼0
kT
f_ ðnÞðt j nÞa dn:
ð34Þ
ðkþ1ÞT
f_ ðnÞðt j nÞa dn:
ð35Þ
ðkþ1ÞT
Since each Ik can be evaluated using the scheme (30) [which is precise to the order OðhÞ], we want to know the value of k for 0 which the contributions to (33) are negligible [of order OðhÞ]. We call this value k , such that 0
a
1 Dt j f ðt j Þ
¼
0
sk k k X X X 1 1 Ik þ 0 Datj f ðt j Þ þ OðhÞ ¼ ½ðt t s1 Þ1a ðt j t s Þ1a ft;s þ 0 Datj f ðt j Þ þ OðhÞ; Cð1 aÞ k¼0 Cð2 aÞ k¼0 s¼1 j
ð36Þ 0
where, for every sum in k, we have t s ¼ ðk þ 1ÞT þ sh, with sk h ¼ T, and the criterion for choosing the value of k is given by s
s
k k X X 1 1 0 0 ½ðt j t s1 Þ1a ðt j t s Þ1a ft;s ¼ ½ðt þ ðk þ 1ÞT ðs 1ÞhÞ1a ðt j þ ðk þ 1ÞT shÞ1a ft;s Cð2 aÞ s¼1 Cð2 aÞ s¼1 j 0
¼ OðhÞ:
0
ð37Þ
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
10705
The Lipschitz continuity of f gives a weaker but sufficient condition, s
k X 1 MT 1a 1 ½ðtj ts1 Þ1a ðtj ts Þ1a ft;s 6 ¼ OðhÞ: Cð2 aÞ s¼1 Cð1 aÞ k0a 0
ð38Þ 0
We can readily conclude that the scheme converges fast when a ! 1, while k becomes very large as a ! 0. The problem of convergence near a ¼ 0 also appears in [26] for a certain fractional differential equation. This scheme can be implemented with other finite difference methods to solve a variety of fractional differential equations. 8. Conclusions In this study, we extended the treatment of the ordinary Ince equation to include the effect of a fractional derivative term of order a > 0. The case a ¼ 1 represents the Ince equation with additional constant damping. The eigenvalue curves and the periodic stable solutions of the equation were determined using a suitable Fourier expansion. The regions of stability and instability on the ð; aÞ-plane were characterized. Two important observations are the detachment of the curves from the a-axis, as well as the appearance of degenerate eigenvalues for a given choice of parameters. The properties of the eigenfunctions of the fractional Ince equation were discussed, and a some plots were presented to illustrate their behavior. These fractional solutions, valid for the steady state of the system, are not orthogonal and have no defined parity. An extension of a discrete method to evaluate the fractional derivative with lower terminal at 1 for a class of functions was also presented. Acknowledgments This research was supported by Consejo Nacional de Ciencia y Tecnologia Mexico (Grant No. 82407) and by Tecnologico de Monterrey (Grant No. CAT141). 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]
I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1990. A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier, Amsterdam, 2006. K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, San Diego, 1974. J.C. Gutiérrez-Vega, Fractionalization of optical beams: I, Planar Analysis. Opt. Lett. 32 (2007) 1521–1523. J. Cao, C. Ma, Z. Jiang, S. Liu, Nonlinear dynamic analysis of fractional order rub-impact rotor system, Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 1443–1463. C.C. Tseng, Designs of fractional delay filter, Nyquist filter, lowpass filter and diamond-shaped filter, Signal Process 87 (2007) 584–601. K. Oldham, Fractional differential equations in electrochemistry, Adv. Eng. Softw. 41 (2010) 9–12. J. Sabatier, O.P. Agrawal, J.A. Tenreiro-Machado, Advances in Fractional Calculus, Springer, Dordrecht, 2007. K. Yonggang, Zhang Xiue, Some comparison of two fractional oscillators. Physica B. 405 (2010) 369–373. P. Wahl, A. Chatterjee, Averaging oscillations with small fractional damping and delayed terms, Nonlinear Dyn. 38 (2004) 3–22. F. Xie, X. Lin, Asymptotic solution of the van der pol oscillator with small fractional damping, Phys. Scripta T136 (2009) 014033. F. Mainardi, Fractional relaxation-oscillation and fractional diffusion-wave phenomena, Chaos Solitons Fract. 7 (1996) 1461–1477. R.H. Rand, S.M. Sah, M.K. Suchorsky, Fractional Mathieu equation, Commun. Nonlinear Sci. Numer. Simulat. 15 (2010) 3254–3262. F.M. Arscott, Periodic Differential Equations, Pergamon Press, London, 1964. R.H. Rand, M. Yazhbin, D.G. Rand, Evolutionary dynamics of a system with periodic coefficients, Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 3887–3895. M. Bandres, J.C. Gutiérrez-Vega, Ince–Gaussian beams, Opt. Lett. 29 (2004) 144–146. M. Bandres, J.C. Gutiérrez-Vega, Ince–Gaussian modes of the paraxial wave equation and stable resonators, J. Opt. Soc. Am. A 21 (2004) 873–880. R. Cordero-Soto, S. Suslov, The degenerate parametric oscillator and Ince’s equation, J. Phys. A: Math. Theor. 44 (2011) 015101. L. Ng, R.H. Rand, Nonlinear effects on coexistence phenomenon in parametric excitation, Nonlinear Dyn. 31 (2003) 73–89. Hartono, A.H.P. van der Burgh, A linear differential equation with a time-periodic damping coefficient: stability diagram and an application, J. Eng. Math. 49 (2004) 99–112. E.L. Ince, A linear differential equation with periodic coefficients, Proc. London Math. Soc. 23 (1923) 56–74. J.H. Taylor, K.S. Narendra, Stability regions for the damped Mathieu equation, SIAM J. Appl. Math. 17 (1969) 343–352. W. Magnus, S. Winkler, Hill’s Equation, Wiley & Sons, New York, 1966. M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations, J. Comput. Appl. Math. 172 (2004) 65–77. M.S. Tavazoei, A note on fractional-order derivatives of periodic functions, Automatica 46 (2010) 945–948. F.I. Taukenova, M.Kh. Shkhanukov-Lafishev, Difference methods for solving boundary value problems for fractional differential equations, Comput. Math. Math. Phys. 46 (2006) 1871–1888.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Fractional Ince equation with a Riemann–Liouville fractional derivative Alfredo Parra-Hinojosa, Julio C. Gutiérrez-Vega ⇑ Photonics and Mathematical Optics Group, Tecnológico de Monterrey, Mexico 64849, Mexico
a r t i c l e
i n f o
Keywords: Ince equation Fractional calculus Ince polynomials Eigenvalue curves Stability Hill equation
a b s t r a c t We extend the classical treatment of the Ince equation to include the effect of a fractional derivative term of order a > 0 and amplitude c. A Fourier expansion is used to determine the eigenvalue curves aðÞ in function of the parameter , the stability domains, and the periodic stable solutions of the fractional Ince equation. Two important observations are the detachment of the eigenvalue curves from the a-axis in the ð; aÞ-plane, as well as the appearance of degenerate eigenvalues for suitable selections of the parameters. The fractional solutions, valid for the steady state of the system, are not orthogonal and have no defined parity. We also introduce a discrete numerical method to evaluate the Riemann–Liouville fractional derivative with lower terminal at 1 for a class of functions. The case a ¼ 1 represents the Ince equation with an additional constant damping. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction In recent years, the study of fractional calculus [1–3] has given rise to physical models that describe real-world phenomena in a more precise way. Several examples can be found in a vast range of scientific fields, such as optics [4], nonlinear dynamics [5], signal processing [6], and electrochemistry [7]. Sabatier et al. [8] presented a comprehensive compilation of 37 papers where fractional calculus has been applied to the areas of particle physics, diffusive systems, electrical systems, viscoelastic media and control, among others. In the same way, recent studies have dealt with the properties of diverse linear fractional differential equations [1,2]. These include, for example, fractional versions of the linear oscillator equation [9], the Duffing equation [10], the van der Pol equation [11], the wave equation [12], and the Mathieu equation [13]. Although fractional calculus was born well over 300 years ago, still many scientists and engineers nowadays are unaware of the properties and potential applications of this rich mathematical field. The canonical form of the ordinary Ince equation with eigenvalue a and parameters and p as presented by Arscott [14] is given by
"
# 2 d d þ sin 2t þ a p cos 2t NðtÞ ¼ 0: dt dt 2
ð1Þ
This equation and its different variants have appeared in the description of numerous physical problems in the fields evolutionary game theory [15], laser optics [16,17], quantum mechanics [18] and vibrational mechanics [19,20], among others. ⇑ Corresponding author. E-mail addresses: [email protected] (A. Parra-Hinojosa), [email protected] (J.C. Gutiérrez-Vega). URL: http://optica.mty.itesm.mx/pmog/ (J.C. Gutiérrez-Vega). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.04.044
10696
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
In this study, we extend the treatment of the ordinary Ince equation to include the effect of an additional fractional derivative term of order a 2 ð0; 2. We determine the eigenvalue curves and discuss the properties of the periodic solutions of this fractional Ince equation. This work also extends the results obtained for the closely related Mathieu equation [13] to which the Ince equation can be reduced [21]. 2. Statement of the problem We are interested in studying the ordinary Ince equation (1) with an additional fractional derivative term
"
# 2 d d a þ sin 2t þ a p cos 2t þ c D 1 t NðtÞ ¼ 0; dt dt 2
ð2Þ
where c 2 R is an amplitude constant and 1 Dat denotes the Riemann–Liouville fractional derivative of order a P 0 with lower limit b ¼ 1 n
a b Dt f ðtÞ ¼
1 d Cðn aÞ dt n
Z b
t
f ðnÞ ðt nÞ1þan
dn;
a P 0;
ð3Þ
with n being the smallest integer exceeding a, i.e., ðn 1Þ 6 a < n. The fractional Ince equation (3) has three special cases: (a) When c ¼ 0, it reduces to the ordinary Ince equation (1). (b) The Ince equation (1) is often regarded as a generalized Mathieu equation with periodic damping, i.e. the first derivative term. Therefore when a ¼ 1 and c – 0 the fractional term becomes a first-derivative damping term with constant amplitude c. The strengths of the periodic and the constant damping terms are adjusted by the constants and c, respectively. In analogy with the damped Mathieu equation [22], we will refer this special case as damped Ince equation. (c) When a ¼ 2 and c – 0, the fractional term becomes a second-derivative oscillatory term of amplitude c that can be added to the already existent second-derivative term. Thus, this case reduces to a scaled version of the ordinary Ince equation (1). The ordinary Ince equation (1) can be reduced to the canonical form of the Hill equation
y00 þ ½k þ Q ðtÞy ¼ 0;
ð4Þ
through the transformation
N ¼ y exp
4
cos 2t ;
ð5Þ
where k is the eigenvalue parameter and Q ðtÞ is a real p-periodic function of t [23]. For this case, we have
QðtÞ ¼ ðp þ 1Þ cos 2t þ ð2 =8Þðcos 4t 1Þ:
k ¼ a;
ð6Þ
From Floquet theory, it is known that Hill’s equation admits solutions of period p or 2p. The same holds for the Ince equation, since the transformation (5) preserves the periodicity and boundedness of the solution, and so the following theorem due to Liapounoff and Haupt applies to the Ince equation [23]: Theorem 1. To every differential equation (4) there belong two monotonically increasing infinite sequences of real numbers
k0 ;
k1 ;
k2 ; ;
ð7aÞ
k01 ;
k02 ;
k03 ;
ð7bÞ
and
k04 ; ;
such that Eq. (4) has a solution of period p if and only if k ¼ kn ; n ¼ 0; 1; 2; . . ., and a solution of period 2p if and only if k ¼ k0n ; n ¼ 1; 2; 3; . . .. The kn and k0n satisfy the inequalities
k0 < k01 6 k02 < k1 6 k2 < k03 < k04 < ;
ð8Þ
and the relations
lim k1 n!1 n
¼ 0;
lim ðk0n Þ1 ¼ 0;
ð9Þ
n!1
The solutions of Eq. (4) are stable in the intervals
ðk0 ; k01 Þ;
ðk02 ; k1 Þ;
ðk2 ; k03 Þ;
ðk04 ; k3 Þ; :
ð10Þ
At the endpoints of these intervals the solutions of Eq. (4) are, in general, unstable. This is always true for k ¼ k0 . The solutions of (4) are stable for k ¼ k2nþ1 or k ¼ k2nþ2 if and only if k2nþ1 ¼ k2nþ2 , and they are stable for k ¼ k02nþ1 or k ¼ k02nþ2 if and only if k02nþ1 ¼ k02nþ2 . We will investigate the effect of the fractional derivative term on the shape of the regions of stability on the parameter– eigenvalue ð; aÞ plane. These regions are delimited by the eigenvalue curves, on which periodic solutions exist.
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
10697
3. Periodic solutions of the fractional Ince equation m Before we obtain the solutions of the fractional Ince equation (2), we first recall that the even C m p ðt; Þ and odd Sp ðt; Þ solutions of order p 2 R and degree m ¼ 0; 1; 2; . . .., of the ordinary Ince equation (1) (i.e. the fractional Ince equation with c ¼ 0) can be expanded in terms of Fourier series. The expansions fall into four classes, according to their symmetry or antisymmetry, about t ¼ 0 and t ¼ p=2, namely [14]
C 2k p ðt; Þ ¼
1 X Ar cos 2rt;
k ¼ 0; 1; . . . ;
r¼0 1 X
ðt; Þ ¼ C 2kþ1 p
Ar cosð2r þ 1Þt;
ð11aÞ
k ¼ 0; 1; . . . ;
ð11bÞ
r¼0
S2k p ðt; Þ ¼
1 X Br sin 2rt;
k ¼ 1; 2; . . . ;
r¼1 1 X
ðt; Þ ¼ S2kþ1 p
Br sinð2r þ 1Þt;
ð11cÞ
k ¼ 0; 1; . . .
ð11dÞ
r¼0 m
1
2
0 1 2 The respective eigenvalues am p ðÞ and bp ðÞ satisfy the relation ap < bp < ap < bp < ap for given values of p and m m that each function is associated with an eigenvalue ap or bp . They satisfy the normalization
Z p
p
2 ½C m p ðt; Þ dt ¼
Z p
p
2 ½Sm p ðt; Þ dt ¼ p
> 0, so ð12Þ
and the orthogonality relation [14]
m0 exp cos 2t C m p ðt; ÞC p ðt; Þdt ¼ 0; 2 p
Z p
m – m0 :
ð13Þ
If p is a non-negative integer the expressions in (11) become finite sums, and the Ince functions become the Ince polynomials, namely
C 2k 2n ðt; Þ ¼
n X Ar cos 2rt;
k ¼ 0; 1; . . . ; n;
ð14aÞ
r¼0
C 2kþ1 2nþ1 ðt; Þ ¼
n X Ar cosð2r þ 1Þt;
k ¼ 0; 1; . . . ; n;
ð14bÞ
r¼0
S2k 2n ðt; Þ ¼
n X Br sin 2rt;
k ¼ 1; 2; . . . ; n;
ð14cÞ
r¼1
S2kþ1 2nþ1 ðt; Þ ¼
n X Br sinð2r þ 1Þt;
k ¼ 0; 1; . . . ; n:
ð14dÞ
r¼0
Even order p ¼ 2n and odd order p ¼ 2n þ 1 Ince polynomials are p- and 2p-periodic, respectively. To each Ince polynomial 2k 2kþ1 2kþ1 we associate the corresponding finite sequences of real eigenvalues a2k 2n ; b2n ; a2nþ1 , and b2nþ1 [17]. We now look for periodic solutions of the fractional Ince equation (2). Consider a formal solution in the form of a Fourier series
NðtÞ ¼ A0 þ
1 X Ar cos rt þ Br sin rt:
ð15Þ
r¼1
The choice of setting b ¼ 1 in Eq. (3) allows the solution to be periodic. This choice is crucial and illustrates the fact that the fractional derivative is a nonlocal operator, since it depends on the values between the lower limit b and the point of interest [24]. The Ince equation often appears in problems where time is the independent variable. The lower limit must then be set to the initial value of t, namely, t ¼ 0. For this case, since we are interested in looking for periodic solutions of Eq. (2), we make use of the following results provided by Tavazoei [25] on the fractional derivative of periodic functions with lower limit b ¼ 0: Theorem 2. Suppose that f ðtÞ is ðn 1Þ-times continuously differentiable and f ðnÞ ðtÞ is bounded. If f ðtÞ is a non-constant periodic function with period T, the function 0 Dat f ðtÞ, where 0 < a R N and n is the first integer greater than a, cannot be a periodic function of period T. Corollary 1. A differential equation of fractional order in the form a
0 Dt
f ðtÞ þ W½f ðtÞ; f 1 ðtÞ; . . . ; f ðnÞ ðtÞ ¼ 0;
ð16Þ
10698
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
where 0 < a R N, cannot have any non-constant smooth periodic solution. For this problem, f ðtÞ is either a sine or a cosine function. Indeed, their fractional derivatives in the Riemann–Liouville sense are not periodic, and they are given by [1] a
cos rt ¼ ðrtÞa E2;1a ðr 2 t 2 Þ;
ð17aÞ
a
sin rt ¼ ðrtÞ1a E2;2a ðr 2 t 2 Þ;
ð17bÞ
0 Dt 0 Dt
where Ea;b ðzÞ is the two-parameter Mittag–Leffler function [1]
Ea;b ðzÞ ¼
1 X
zk ; C ð a k þ bÞ k¼0
a>0
ð18Þ
b>0
Nevertheless, for large values of t, Eq. (17) converge to the following periodic functions: a
cos rt ¼ ra cosðrt þ ap=2Þ;
t ! 1;
ð19aÞ
a
sin rt ¼ r a sinðrt þ ap=2Þ;
t ! 1:
ð19bÞ
0 Dt 0 Dt
Mathematically, setting b ¼ 1 in Eq. (3) gives the same result, i.e., a
1 Dt
cos rt ¼ r a cosðrt þ ap=2Þ;
ð20aÞ
a
a 1 Dt sin rt ¼ r sinðrt þ ap=2Þ:
ð20bÞ
In more physical terms, the solutions are valid for the steady state of the system in question. Also, for this choice, the Riemann–Liouville and Caputo definitions of the fractional derivative are equal [1]. Inserting (15) into (2) and collecting terms we obtain the following recurrence relations:
aA0 ð=2Þðp þ 2ÞA2 ¼ 0;
ð21aÞ
pA0 þ ½a 4 þ 2a c cosðap=2ÞA2 ð=2Þðp þ 4ÞA4 þ 2a c sinðap=2ÞB2 ¼ 0;
ð21bÞ
½a 4 þ 2a c cosðap=2ÞB2 ð=2Þðp þ 4ÞB4 2a c sinðap=2ÞA2 ¼ 0;
ð21cÞ
½a 1 ð=2Þðp þ 1Þ þ c cosðap=2ÞA1 ð=2Þðp þ 3ÞA3 þ c sinðap=2ÞB1 ¼ 0;
ð21dÞ
½a 1 þ ð=2Þðp þ 1Þ þ c cosðap=2ÞB1 ð=2Þðp þ 3ÞB3 c sinðap=2ÞA1 ¼ 0:
ð21eÞ
For r P 3
ð=2Þðr 2 pÞAr2 þ ½a r 2 þ r a c cosðap=2ÞAr ð=2Þðr þ 2 þ pÞArþ2 þ r a c sinðap=2ÞBr ¼ 0;
ð21fÞ
ð=2Þðr 2 pÞBr2 þ ½a r 2 þ r a c cosðap=2ÞBr ð=2Þðr þ 2 þ pÞBrþ2 r a c sinðap=2ÞAr ¼ 0;
ð21gÞ
a
where we used 1Dt A0 ¼ 0 for a > 0. Note that the Ar and Br coefficients are now coupled, and there are two independent systems, one involving only even subscripts (A0 ; B2 ; A2 ; . . .) and the other involving odd subscripts (B1 ; A1 ; B3 ; . . .). This means that there exist only two kinds of solutions, instead of four, as is the case of the ordinary Ince equation aÞ N2k;ð ðt; Þ ¼ A0 þ p
1 X A2r cos 2rt þ B2r sin 2rt;
k ¼ 0; 1; 2; . . . ;
ð22aÞ
r¼1
Np2kþ1;ðaÞ ðt; Þ ¼
1 X A2rþ1 cosð2r þ 1Þt þ B2rþ1 sinð2r þ 1Þt;
k ¼ 0; 1; 2; . . .
ð22bÞ
r¼0
The first solution (22a) is p-periodic, whereas (22b) is 2p-periodic, but none has a defined parity. We label the eigenvalues of 2k;ðaÞ
aÞ N 2k;ð ðt; Þ as ap p
2kþ1;ðaÞ
aÞ ðÞ, and those of N 2kþ1;ð ðt; Þ as ap p 0;ðaÞ ap
2;ðaÞ ap
4;ðaÞ ap
ðÞ, where the superscript ðaÞ refers to the order of the fractional 1;ðaÞ
3;ðaÞ
5;ðaÞ
< < and ap < ap < ap for > 0, and are normalized in the same derivative. They follow the order way as the transcendental Ince functions Eqs. (12). It is not difficult to show that, when p is a non-negative integer, both series (22) become finite, and thus the fractional Ince functions (22) become fractional Ince polynomials. In particular, for solution (22a) (to which we associate the recurrence relations (21) with even coefficients), let p be a positive integer (p ¼ 2n), and let us choose the parameter a such that A2nþ2 ¼ B2nþ2 ¼ 0. Then, for r ¼ 2n þ 2, Eqs. (21f) and (21g) give A2nþ4 ¼ A2nþ6 ¼ ¼ 0 and B2nþ4 ¼ B2nþ6 ¼ ¼ 0 respectively, which means that the solution is finite (ending in A2n and B2n ). This condition results in the following set of recurrence relations:
aA0 þ ðn þ 1ÞA2 ¼ 0;
ð23aÞ
2nA0 þ ½4 a 2a c cosðap=2ÞA2 þ ðn þ 2ÞA4 2a c sinðap=2ÞB2 ¼ 0;
ð23bÞ
½4 a 2a c cosðap=2ÞB2 þ ðn þ 2ÞB4 þ 2a c sinðap=2ÞA2 ¼ 0:
ð23cÞ
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
10699
For r ¼ 4; 6; . . . ; 2n 2,
ðn r=2 þ 1ÞAr2 þ ½r 2 a ra c cosðap=2ÞAr þ ðn þ r=2 þ 1ÞArþ2 r a c sinðap=2ÞBr ¼ 0;
ð23dÞ
ðn r=2 þ 1ÞBr2 þ ½r 2 a r a c cosðap=2ÞBr þ ðn þ r=2 þ 1ÞBrþ2 þ cra sinðap=2ÞAr ¼ 0:
ð23eÞ
Finally,
A2n2 þ ½4n2 a ð2nÞa c cosðap=2ÞA2n ð2nÞa c sinðap=2ÞB2n ¼ 0; B2n2 þ ½4n2 a ð2nÞa c cosðap=2ÞB2n þ ð2nÞa c sinðap=2ÞA2n ¼ 0:
ð23fÞ ð23gÞ
There are 2n þ 1 equations for the ð2n þ 1Þ coefficients A2r and B2r . The case c ¼ 0 (or a ¼ 0), B2r ¼ 0 reduces to the recurrence relations reported by Arscott in [14]. A similar analysis can be carried out for the 2p-periodic solution with p ¼ 2n þ 1, which gives the following relations:
½a 1 ðn þ 1Þ þ c cosðap=2ÞA1 ðn þ 2ÞA3 þ c sinðap=2ÞB1 ¼ 0;
ð24aÞ
½a 1 þ ðn þ 1Þ þ c cosðap=2ÞB1 ðn þ 2ÞB3 c sinðap=2ÞA1 ¼ 0:
ð24bÞ
For r ¼ 3; 5; . . . ; 2n 1,
ð2n r þ 3Þð=2ÞAr2 þ ½r 2 a cra cosðap=2ÞAr þ ð2n þ r þ 3Þð=2ÞArþ2 cra sinðap=2ÞBr ¼ 0;
ð24cÞ
ð2n r þ 3Þð=2ÞBr2 þ ½r 2 a cr a cosðap=2ÞBr þ ð2n þ r þ 3Þð=2ÞBrþ2 þ cra sinðap=2ÞAr ¼ 0:
ð24dÞ
Finally,
A2n1 þ ½ð2n þ 1Þ2 a cð2n þ 1Þa cosðap=2ÞA2nþ1 cð2n þ 1Þa sinðap=2ÞB2nþ1 ¼ 0; B2n1 þ ½ð2n þ 1Þ2 a cð2n þ 1Þa cosðap=2ÞB2nþ1 þ cð2n þ 1Þa sinðap=2ÞA2nþ1 ¼ 0:
ð24eÞ ð24fÞ
Finally, the expansions of the fractional Ince polynomials are given by 2k;ðaÞ
N2n
ðt; Þ ¼ A0 þ
n X A2r cos 2rt þ B2r sin 2rt;
k ¼ 0; 1; . . . ; n;
ð25aÞ
r¼1
2kþ1;ðaÞ
N2nþ1
ðt; Þ ¼
n X A2rþ1 cosð2r þ 1Þt þ B2rþ1 sinð2r þ 1Þt;
k ¼ 0; 1; . . . ; n;
ð25bÞ
r¼0 2k;ðaÞ
2kþ1;ðaÞ
with eigenvalues a2n ðÞ and a2nþ1 ðÞ, and where the Fourier coefficients are calculated with the system (24)?. The fractional Ince solutions are normalized in the same way as the ordinary Ince solutions. 4. Behavior of the eigenvalue curves An important problem associated with the fractional Ince equation is the determination of the eigenvalues aðÞ, for which the solutions are periodic in t. To make suitable comparisons, in Fig. 1 we first plot the eigenvalue curves aðÞ and bðÞ for the m ordinary Ince equation. Thick and thin lines denote the characteristic curves for even C m p ðt; Þ and odd Sp ðt; Þ Ince functions, respectively. Even-order and odd-order curves are symmetrical and antisymmetrical about the axis a, respectively. Except for the trivial case ¼ 0, the curves do not intersect and thus divide the ð; aÞ plane in regions where the solutions of the Ince equation are stable (white regions) or unstable (gray regions). For referring purposes, the unstable regions are numbered sequentially, as shown in Fig. 2. Note that the jth unstable region touches the axis q ¼ 0 at values a ¼ jjj2 . To get insight into the eigenvalue curves of the fractional Ince equation, let us analyze first the simple case when p is an even integer, say p ¼ 2n ¼ 2. The p-periodic solution (22a) [which fulfills the system (23)] is finite and has the form (25a), 2k;ðaÞ namely N 2 ¼ A0 þ B2 sin 2t þ A2 cos 2t, whereas the 2p-periodic solution (22b) has an infinite number of terms whose coefficients fulfill the odd system in Eq. (21). Recurrence relations (23) can be written in matrix form as follows:
2
32 3 2 3 0 0 2 A0 A0 6 76 7 6 7 a a 0 4 2 c cosð a p =2Þ 2 c sinð a p =2Þ B ¼ a 4 54 2 5 4 B2 5: a a 2 2 c sinðap=2Þ 4 2 c cosðap=2Þ A2 A2
ð26Þ
The matrix in Eq. (26) (denoted by M) is the matrix representation of the fractional differential operator ½d=dt 2 sin 2td=dt þ p cos 2t c1 Dat in the Fourier basis. The eigenvalues of the system are found by setting detðM aIÞ ¼ 0, which gives a relation between and a that defines the eigenvalue curves as a function of the parameter for fixed values of c and a. This consequence is drawn from Floquet theory, which states that periodic solutions (of period p or 2p) exist on the ð; aÞ curves that separate regions of stable and unstable solutions.1 For this case, the characteristic equation is given by the expression 1
A solution is stable if every independent solution is bounded.
10700
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
25
20
4
4
4 8
4 8
a8 :(C8 )
-4
b :(S )
15
10
a
3
3
3
5
4
4 8
4
b :(S )
3
b8 :(S8 ) 2
3
3
3
3
a8 :(C8 )
a8 :(C8 )
-3
4
4 8
a8 :(C8 )
Stable
b8 :(S8 ) 2
2
3
2
a8 :(C8 )
a8 :(C8 )
Unstable 2
b82 :(S82 )
0
1 8
a :(C )
1
b8 :(S8 )
a :(C )
-2
1
b8 :(S8 )
-5
a80 :(C80 )
-1 -10 -5
-4
-3
2
b8 :(S8 )
1 8
-2
0
-1
1 8
1 8
1
1
2
a80 :(C80 )
0
1
1
2
3
4
5
ε m Fig. 1. Eigenvalue curves of the ordinary Ince equation (c ¼ 0) with p ¼ 8. Thick and thin lines denote curves for even C m p ðt; Þ and odd Sp ðt; Þ Ince functions, respectively. The characteristic curves divide the plane into regions of stability (white) and instability (gray). The unstable regions are numbered sequentially for referring purposes.
25
20
-4
a88,(0.5):(N88,(0.5))
a88,(0.5):(N88,(0.5))
Stable
a86,(0.5):(N86,(0.5))
4
a86,(0.5):(N86,(0.5))
15 a85,(0.5):(N85,(0.5))
10
-3
a
a87,(0.5):(N87,(0.5))
5
a
4,(0.5) 8
a
0
-2
a87,(0.5):(N87,(0.5))
4,(0.5) 8
:(N
2,(0.5) 8
:(N
a
)
2,(0.5) 8
) a
0,(0.5) 8
:(N
4,(0.5) 8
4,(0.5) 8
2,(0.5) 8
2,(0.5) 8
) a
0,(0.5) 8
:(N
:(N
) )
Unstable
2
-5
a81,(0.5):(N81,(0.5))
-4
-3
a83,(0.5):(N83,(0.5)) a83,(0.5):(N83,(0.5))
a81,(0.5):(N81,(0.5))
-1 -10 -5
3
a85,(0.5):(N85,(0.5))
1
0
-2
-1
0
1
2
3
4
5
ε Fig. 2. Eigenvalue curves for p ¼ 8 with fractional coefficient c ¼ 2 and half-order a ¼ 0:5. The curves originally crossing the axis at a ¼ 0; 1; 4; 9; . . ., become detached. Also, some curves exhibit crossings, and the crossing points indicate double eigenvalues.
a½ð4 aÞ2 2ð4 aÞ2a c cosðap=2Þ þ 22a c2 þ 42 ½4 a 2a c cosðap=2Þ ¼ 0;
ð27Þ
which is cubic in a and whose solutions give three curves a ¼ aðÞ. Setting ¼ 0 in (27), the intersection points in the a axis can be obtained. The only real solution is a ¼ 0, and the other two solutions correspond to a ¼ 4 2a c expðiap=2Þ. Thus, there is only one intersection when c and a are not zero, and the curves originally intersecting the axis at a ¼ 4 (as shown in Fig. 1) become ‘‘detached’’ from the vertical axis. This detachment of the eigenvalue curves may be anticipated from the Mathieu equation when a fractional derivative term is present [13]. The next finite p-periodic solution (namely, n ¼ 2; p ¼ 4, 2k;ðaÞ so that N 4 ¼ A0 þ B2 sin 2t þ A2 cos 2t þ B4 sin 4t þ A4 cos 4t) gives a third intersection when c ¼ 0 (for a ¼ 16), which also disappears when c and a are not zero.
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
10701
40 6,(1.75)
a116,(1.75):(N11
7,(1.75)
a117,(1.75):(N11
)
30
)
4 a115,(1.75):(N115,(1.75))
20
4,(1.75)
a114,(1.75):(N11
a
)
10 2,(1.75)
a112,(1.75):(N11
)
3
0 3,(1.75)
a110,(1.75):
a113,(1.75):(N11
-10 (N110,(1.75))
0
)
a111,(1.75):(N111,(1.75))
0
2
1 5
10
15
ε Fig. 3. Eigenvalue curves for p ¼ 11 with fractional coefficient c ¼ 2 and fractional order a ¼ 1:75. The crossings that appear for a 2 ð0; 1Þ (c > 0) are no longer present.
In general the effect of the fractional term on the eigenvalue curves is illustrated in Fig. 2 for p ¼ 8; c ¼ 2, and a ¼ 0:5. Note that the eigenvalue curves that arise from the 2p-periodic solutions have also been included. Given that these solutions are not finite for p – 2n þ 1, the way to obtain the eigenvalue curves consists in truncating the series solution (22b) and using recurrence relations (21), which is the general procedure for any p 2 R. The larger the order of truncation, the more precise the curves are. Here, a 2p-periodic solution with 52 terms was used. The curves divide the ð; aÞ-plane into regions where the solutions of the fractional Ince equation are stable (white region) or unstable (gray region). The unstable regions are numbered sequentially following the same convention as in Fig. 1. Observation of Fig. 2 gives the following conclusions for c > 0 and a 2 ð0; 1Þ: 0;ð0:5Þ
Except for the lowest-order curve a8 ðÞ, each eigenvalue curve of the fractional Ince equation becomes ‘‘detached’’ from the a-axis. The upper and lower boundaries of the unstable regions correspond to the original curves aðÞ and bðÞ of the ordinary Ince equation and are depicted using different line-widths in Fig. 2. In the ‘‘detached’’ (i.e., stable) region, the eigenvalues become complex, and they appear in complex conjugate pairs (i.e., 2nþ1;ðaÞ 2nþ3;ðaÞ 2nþ2;ðaÞ 2nþ4;ðaÞ ap ¼ ap and ap ¼ ap ; n even). The upper and lower eigenvalue curves of an unstable region connect smoothly forming a single continuous boundarycurve. For > 0, the turning-point of the boundary-curve is the minimum value of for which the eigenvalue a is real. For each boundary-curve, the solutions above and below the turning-points are labeled with different superscripts. Periodic solutions corresponding to different points ð; aÞ lying on the same boundary-curve have the same periodicity, namely, p and 2p for boundary-curves of even and odd unstable regions, respectively. For any given value of , there is only a finite number of real eigenvalues a, which increases as increases. Consequently, unlike the ordinary Ince equation, the fractional Ince equation has a finite number of real periodic solutions for a given value of . There exist real values where the eigenvalue curves intersect, indicating the existence of degenerate eigenvalues (i.e., double points) a of the fractional Ince equation. In particular, for the choice a ¼ 0:5; c ¼ 2, the crossings occur for all pairs 2n;ðaÞ 2nþ1;ðaÞ 4n;ðaÞ 4nþ2;ðaÞ 4nþ1;ðaÞ 4nþ3;ðaÞ of curves fap ; ap g when > 0, and for the curves fap ; ap g and fap ; ap g when < 0 (n ¼ 0; 1; . . .). The adjacent curves become asymptotically parallel after the crossing, therefore they only cross once. Additionally, the effect of increasing p is to reduce the area of stability (the curves appear less detached), which means that the number of real eigenvalues a for a given value of increases. Nevertheless, for a 2 ½1; 2, we observed that eigenvalue curves do not intersect, i.e., the degenerate eigenvalues disappear, as shown in Fig. 3. As expected, when a ¼ 2, we recover the ordinary Ince equation (1) scaled by a factor 1 þ c and the curves become ‘‘reattached’’, as shown in Fig. 1. 5. Behavior of the solutions of fractional order Since the fractional Ince equation involves two additional parameters (c and a), a general description of the behavior of the eigenfunctions can be fairly complicated. In this section we shall illustrate the effect of the fractional term on the shape
10702
(a)
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
1
(b)
C13(t,4)
1.5
S33(t,3) 1
0.5 0.5
C13(t,0)=cos(t)
0
3 S11 (t,3)
0
-0.5
-0.5
-1 0
π/4
π/2
3π/4
0
π
π/4
t
(c)
(d)
1
0.5
0.5
0
0
-0.5
-0.5
π/2
π t
3π/2
2π
α=0 α=0.5 α=1
-1
α=0 α=0.5 α=1
0
π
1.5
1
-1.5
3π/4
t
1.5
-1
π/2
3π/2
2π
-1.5
0
π/2
π t
3;ðaÞ
Fig. 4. Plots of several Ince polynomials and fractional solutions. (a) C 13 ðt; Þ; 2 f0; 1; 2; 3; 4g; (b) S3p ðt; 3Þ; p 2 f3; 5; 7; 9; 11g; (c) N 3 5;ðaÞ reduces to C 13 ðt; 11Þ when a ¼ 0; (d) N 3 ðt; 11Þ with c ¼ 1. It reduces to S33 ðt; 11Þ when a ¼ 0.
ðt; 11Þ with c ¼ 1. It
2;ðaÞ
Fig. 5. N 6:5 ðt; 5Þ eigenfunction of the fractional Ince equation for c ¼ 2; a 2 ½0; 2. Solid curves correspond to a ¼ f0; 0:5; 1; 1:5; 2g.
of the solutions with a few examples. First, in Fig. 4(a) and (b) we show the ordinary Ince polynomials C 13 ðt; Þ and S3p ðt; 3Þ. 3;ðaÞ 3;ðaÞ Next, we consider the fractional Ince polynomial N 3 ðt; 11Þ with eigenvalue a3 . For c ¼ 1 and a ¼ f0:5; 1g, recurrence
10703
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
relations (23) can be used to obtain the Fourier coefficients and plot the functions, which reduce to C 13 ðt; 11Þ when a or c 5;ðaÞ equal to zero (Fig. 4(c)). The same procedure for N 3 ðt; 11Þ is illustrated in Fig. 4(d). 2;ðaÞ In Fig. 5 we show the behavior of the p-periodic normalized function N 6:5 ðt; 5Þ in the range a 2 ½0; 2. When a vanishes, the fractional Ince equation reduces to the ordinary Ince equation with a shifted eigenvalue a þ c, and the eigenfunction becomes the Ince function S26:5 ðt; 5Þ. Similarly, when a ¼ 2, the fractional Ince equation reduces again to the ordinary Ince equa2;ð2Þ tion but now scaled by 1 þ c, and the solution N 6:5 ðt; 5Þ becomes S26:5 ðt; 5=ð1 þ cÞÞ. In all the examples shown, when the order a varies, the periodicity of the solutions is preserved, while the amplitude is altered and zeros of the polynomials are slightly shifted, except for a ¼ 0; 2. Likewise, the solutions have no defined parity, as opposed to the non-fractional Ince functions. 6. Behavior of the degenerate solutions of fractional order Since the fractional Ince equation has a periodic solution for each point ð; aÞ lying on the eigenvalue curves in Fig. 2, the cross-points between adjacent curves define degenerate Ince functions with the same values of and a. Let us next illustrate the behavior of the degenerate solutions. Fig. 6(a) shows the first two eigenvalue curves aðÞ for a ¼ 0:6; p ¼ 5:5, and c ¼ 2. The crossing points at the values ¼ n indicate degenerate eigenvalues an , namely aÞ a2n2;ð ðn Þ ¼ ap2n1;ðaÞ ðn Þ an ; p
n ¼ 0; 1; 2; . . .
ð28Þ
These are marked with asterisks on the curves in Fig. 6(a) and their numerical values are given up to ten decimals of precision in Fig. 6(b).
(a)
1 0.5 0
(b)
-0.5
a
-1
(ε1*, a*) 1
n=1 n=2
-1.5
Crossing points ε*n 0.60919851809 2.8966572868
a*n −1.4509203434 −2.4053344921
-2
(ε2*, a*) 2
-2.5 0
(c)
0.5
1
1.5
ε
2
2.5
3
3.5
1.5
(d)
1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
-1.5
0
π/2
π
t
3π/2
2π
-1.5
0
π/2
π
3π/2
2π
t
Fig. 6. (a) Eigenvalue curves for a ¼ 0:6; c ¼ 2. Asterisks indicate the crossing points between adjacent curves. (b) Numerical values of the first two crossing 0;ð0:6Þ 1;ð0:6Þ points. (c) Eigenfunctions corresponding to the first degenerate eigenvalue a1 . N 5:5 ðt; 1 Þ (red-circle line) and N 5:5 ðt; 1 Þ (blue-point line); (d) 2;ð0:6Þ 3;ð0:6Þ Eigenfunctions corresponding to the second degenerate eigenvalue a2 . N 5:5 ðt; 2 Þ (red-circle line) and N 5:5 ðt; 2 Þ (blue-point line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
10704
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
The behavior of the first and second pairs of degenerate fractional Ince functions is depicted in Fig. 6(c) and (d). Consider, for 0;ð0:6Þ 1;ð0:6Þ instance, the pair ð1 ; a1 Þ. The degenerate functions N 5:5 ðt; 1 Þ and N 5:5 ðt; 1 Þ satisfy the same differential equation, namely
"
# 2 d d f0;1g;ð0:6Þ 0:6 þ 1 sin 2t þ a1 5:51 cos 2t þ 21 Dt N5:5 ¼ 0: dt dt 2
ð29Þ
Note that the even-order functions (red-circle lines) are p-periodic whereas odd-order functions (blue-point lines) are 2pperiodic. In general, in any pair of degenerate solutions one of them has period p and the other 2p. In particular, the p-periodic functions (red-circle lines) shown in Fig. 6(c) and (d) belong to the same boundary-curve of the second unstable region. 7. Numerical considerations In order to test our solutions of the fractional Ince equation we implemented accurate numerical routines based on suitable finite differences schemes. Firstly, solution of the Ince equation (1) and the fractional Ince equation (2) with integer a derivatives can be performed using standard numerical methods. However, the numerical solution of the fractional Ince equation (2) with arbitrary a is a more complicated task that deserve some discussion. For 0 < a < 1 the fractional Ince functions were verified by calculating the fractional derivative (3) with an extended scheme based on the following discrete form of (3) with lower limit b ¼ 0 due to Taukenova and Shkhanukov-Lafishev [26]:
i Xh 1 ðt j t s1 Þ1a ðt j t s Þ1a ft;s ; Cð2 aÞ s¼1 j
a
0 Dt j f ðt j Þ
¼
ð30Þ
where t j ¼ jh; j ¼ 0; 1; . . . ; j0 ; ft;s ¼ ðfs fs1 Þ=h, such that a
0 Dt j f
¼
1 Cð1 aÞ
Z 0
tj
f_ ðnÞðt j nÞa dn ¼ 0 Datj f þ OðhÞ:
ð31Þ
This scheme is obtained in the following way, assuming that f is Lipschitz continuous. Then
1 Cð1 aÞ
Z
tj 0
j Z ts j Z ts _ X X f ðts1=2 Þ þ €f ðgÞðn t s1=2 Þdn f_ ðnÞdn f_ ðnÞdn 1 1 a ¼ a ¼ C ð1 a Þ C ð1 a Þ ðt j nÞ ðt j nÞ ðt j nÞa s¼1 ts1 s¼1 t s1 j j Z ts € X X f ðgÞðn ts1=2 Þdn 1 1 2 ¼ ½ðt j t s1 Þ1a ðt j t s Þ1a ft;s þ þ Oðh Þ; ð32Þ Cð2 aÞ s¼1 Cð1 aÞ s¼1 ts1 ðt j nÞa
where g is an intermediate point between n and t s1=2 . Now, since j€f ðnÞj 6 M, we find that the second term in (32) is of order OðhÞ, and we obtain (31). To extend this scheme with the lower limit at 1, set t0 ! 1, and a numerical value for t0 is to be obtained. We have a
1 Dt j f ðt j Þ
¼
1 Cð1 aÞ
Z
tj
f_ ðnÞðt j nÞa dn ¼
1
1 Cð1 aÞ
Z
0
f_ ðnÞðt j nÞa dn þ
Z
tj
f_ ðnÞðt j nÞa dn :
ð33Þ
0
1
f ðtÞ being periodical, the infinite limit can be discretized as an integer number of periods nT, so that
Z
0
f_ ðnÞðt j nÞa dn ¼ lim
n!1
1
We call
Ik ¼
Z
kT
Z
0
f_ ðnÞðt j nÞa dn ¼ lim
n!1
nT
n1 Z X k¼0
kT
f_ ðnÞðt j nÞa dn:
ð34Þ
ðkþ1ÞT
f_ ðnÞðt j nÞa dn:
ð35Þ
ðkþ1ÞT
Since each Ik can be evaluated using the scheme (30) [which is precise to the order OðhÞ], we want to know the value of k for 0 which the contributions to (33) are negligible [of order OðhÞ]. We call this value k , such that 0
a
1 Dt j f ðt j Þ
¼
0
sk k k X X X 1 1 Ik þ 0 Datj f ðt j Þ þ OðhÞ ¼ ½ðt t s1 Þ1a ðt j t s Þ1a ft;s þ 0 Datj f ðt j Þ þ OðhÞ; Cð1 aÞ k¼0 Cð2 aÞ k¼0 s¼1 j
ð36Þ 0
where, for every sum in k, we have t s ¼ ðk þ 1ÞT þ sh, with sk h ¼ T, and the criterion for choosing the value of k is given by s
s
k k X X 1 1 0 0 ½ðt j t s1 Þ1a ðt j t s Þ1a ft;s ¼ ½ðt þ ðk þ 1ÞT ðs 1ÞhÞ1a ðt j þ ðk þ 1ÞT shÞ1a ft;s Cð2 aÞ s¼1 Cð2 aÞ s¼1 j 0
¼ OðhÞ:
0
ð37Þ
A. Parra-Hinojosa, J.C. Gutiérrez-Vega / Applied Mathematics and Computation 219 (2013) 10695–10705
10705
The Lipschitz continuity of f gives a weaker but sufficient condition, s
k X 1 MT 1a 1 ½ðtj ts1 Þ1a ðtj ts Þ1a ft;s 6 ¼ OðhÞ: Cð2 aÞ s¼1 Cð1 aÞ k0a 0
ð38Þ 0
We can readily conclude that the scheme converges fast when a ! 1, while k becomes very large as a ! 0. The problem of convergence near a ¼ 0 also appears in [26] for a certain fractional differential equation. This scheme can be implemented with other finite difference methods to solve a variety of fractional differential equations. 8. Conclusions In this study, we extended the treatment of the ordinary Ince equation to include the effect of a fractional derivative term of order a > 0. The case a ¼ 1 represents the Ince equation with additional constant damping. The eigenvalue curves and the periodic stable solutions of the equation were determined using a suitable Fourier expansion. The regions of stability and instability on the ð; aÞ-plane were characterized. Two important observations are the detachment of the curves from the a-axis, as well as the appearance of degenerate eigenvalues for a given choice of parameters. The properties of the eigenfunctions of the fractional Ince equation were discussed, and a some plots were presented to illustrate their behavior. These fractional solutions, valid for the steady state of the system, are not orthogonal and have no defined parity. An extension of a discrete method to evaluate the fractional derivative with lower terminal at 1 for a class of functions was also presented. Acknowledgments This research was supported by Consejo Nacional de Ciencia y Tecnologia Mexico (Grant No. 82407) and by Tecnologico de Monterrey (Grant No. CAT141). 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]
I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1990. A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier, Amsterdam, 2006. K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, San Diego, 1974. J.C. Gutiérrez-Vega, Fractionalization of optical beams: I, Planar Analysis. Opt. Lett. 32 (2007) 1521–1523. J. Cao, C. Ma, Z. Jiang, S. Liu, Nonlinear dynamic analysis of fractional order rub-impact rotor system, Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 1443–1463. C.C. Tseng, Designs of fractional delay filter, Nyquist filter, lowpass filter and diamond-shaped filter, Signal Process 87 (2007) 584–601. K. Oldham, Fractional differential equations in electrochemistry, Adv. Eng. Softw. 41 (2010) 9–12. J. Sabatier, O.P. Agrawal, J.A. Tenreiro-Machado, Advances in Fractional Calculus, Springer, Dordrecht, 2007. K. Yonggang, Zhang Xiue, Some comparison of two fractional oscillators. Physica B. 405 (2010) 369–373. P. Wahl, A. Chatterjee, Averaging oscillations with small fractional damping and delayed terms, Nonlinear Dyn. 38 (2004) 3–22. F. Xie, X. Lin, Asymptotic solution of the van der pol oscillator with small fractional damping, Phys. Scripta T136 (2009) 014033. F. Mainardi, Fractional relaxation-oscillation and fractional diffusion-wave phenomena, Chaos Solitons Fract. 7 (1996) 1461–1477. R.H. Rand, S.M. Sah, M.K. Suchorsky, Fractional Mathieu equation, Commun. Nonlinear Sci. Numer. Simulat. 15 (2010) 3254–3262. F.M. Arscott, Periodic Differential Equations, Pergamon Press, London, 1964. R.H. Rand, M. Yazhbin, D.G. Rand, Evolutionary dynamics of a system with periodic coefficients, Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 3887–3895. M. Bandres, J.C. Gutiérrez-Vega, Ince–Gaussian beams, Opt. Lett. 29 (2004) 144–146. M. Bandres, J.C. Gutiérrez-Vega, Ince–Gaussian modes of the paraxial wave equation and stable resonators, J. Opt. Soc. Am. A 21 (2004) 873–880. R. Cordero-Soto, S. Suslov, The degenerate parametric oscillator and Ince’s equation, J. Phys. A: Math. Theor. 44 (2011) 015101. L. Ng, R.H. Rand, Nonlinear effects on coexistence phenomenon in parametric excitation, Nonlinear Dyn. 31 (2003) 73–89. Hartono, A.H.P. van der Burgh, A linear differential equation with a time-periodic damping coefficient: stability diagram and an application, J. Eng. Math. 49 (2004) 99–112. E.L. Ince, A linear differential equation with periodic coefficients, Proc. London Math. Soc. 23 (1923) 56–74. J.H. Taylor, K.S. Narendra, Stability regions for the damped Mathieu equation, SIAM J. Appl. Math. 17 (1969) 343–352. W. Magnus, S. Winkler, Hill’s Equation, Wiley & Sons, New York, 1966. M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations, J. Comput. Appl. Math. 172 (2004) 65–77. M.S. Tavazoei, A note on fractional-order derivatives of periodic functions, Automatica 46 (2010) 945–948. F.I. Taukenova, M.Kh. Shkhanukov-Lafishev, Difference methods for solving boundary value problems for fractional differential equations, Comput. Math. Math. Phys. 46 (2006) 1871–1888.
