STABILITY ANALYSIS OF DELAYED FRACTIONAL INTEGRO-DIFFERENTIAL EQUATIONS WITH APPLICATIONS OF RLC CIRCUITS

This article presents the stability analysis of integro-differential equations with a delay and a fractional order derivative via some approximation techniques for the derived nonlinear terms of characteristic exponents. Based on these techniques, the existence of some analytical solutions at the neighborhood of their equilibrium points is proved. Stability charts are constructed and so both of the critical time delay and the critical frequency formulae are obtained. The impact of this work into general RLC circuit applications containing delays and fractional order derivatives is discussed. Integro-differential equations; Fractional calculus; Stability analysis; Periodic solutions:


INTRODUCTION
Delayed fractional dynamical systems are systems characterized by the existence of fractional order derivatives and time delay parameters depending on their past states. During the last decades, the combination of the past effect and the fractional order rate has a significant interest in the field of fractional order differential equations(FODEs), since it could be a closer step to describe real processes or a finer approach to describe real states of systems, cf. Bellen and Zennaro [6], Bhalekar et al. [7] and Chen and Moore [10]. Obviously, it can be noticed that the descriptive modeling using only the integer order derivative for such applied processes may differ significantly from their actual or experimental behaviors, cf. Deng et al. [11], Miller and Ross [29] and Podlubny et al. [31,32]. However, many physical models can display the two coupled effects such as anomalous diffusion systems, feedback amplifiers and capacitor systems, electrode-electrolyte interface models, fractional multipoles, viscoelasticity models, electrical circuital systems, electroanalytical chemistry, economical models, polymer models, and control systems, cf. Baleanu et al. [3], Kempfle et al. [18], Kilbas et al. [19], Machado et al. [24], Mainardi [25], Manabe et al. [26] and Rossikhin and Shitikova [36]. Even for fractional order models of happiness and love which have been developed in recent years, it is claimed that they give a better representation than the classical approach, cf. Song et al. [38]. Similar to classical differential systems, the study of stability is also considered the central task for fractional differential systems. Many works e.g. Chen and Moore [10], Qian et al. [34] and the references therein are studied the stability of FODEs (with/without retarded time) under a variety of methods with various directions. In this work, it is devoted to discuss the stability of integro-differential equations that have the combination of the delay and fractional derivative. In particularly, it is focused on disclosing some new qualitative aspects for the given physical systems via new approximation techniques. The model which has to be studied is where q is the fractional order derivative, k is a constant, y τ = y(t − τ ) and f : (−∞, ∞) × C → R, where C is the space of continuous functions, is a nonlinear function. To specify the solution completely, it is required an initial time t o ∈ R and the history φ ∈ C(t o ) where φ : [−τ, t o ] → R. The linearized form at the equilibrium points(y * ) is to D q t y(t) = a(y − y * ) + b(y τ − y * ) + F (y − y * , y τ − y * ) + k t t−τ y(z)dz.
Then, the linearized form is generalized with x(·) = y(·) − y * as follows: where a, b and k are constants in general, taking in your account the function of nonlinear terms F (x, x τ ) satisfies lim (x,xτ )→(0,0) The typical example to describe the linearized form (Eq. (3)) is the shown voltagecurrent circuit in FIGURE 1. This figure explains RLC circuit with a source (E s ), resistance (R a and R b ), capacitance (C), inductance (L) and a linear amplifier as shown with phase shifting network producing a constant time delay between the input and the output of the amplifier. The resulting delayed fractional integrodifferential equation takes the form where q is the fractional order derivatives and i(t) is the current. To regard the detailed of the classical circuit(q = 1), cf. Burton [9]. The common definitions that concerning with the fractional derivative being Grünwald-Letnikov fractional derivative, Riemann-Liouville fractional derivative, and Caputo's fractional derivative, are the mostly used to express the fractional dynamical systems, cf. Babakhani et al. [2], Jarad et al. [17], Miller and Roos [29], Oldham and Spanier [30] and Podlubny [32]. The definitions of Caputo and Riemann-Liouville fractional derivatives and their properties are given in the following definition.
Definition 1.1. For a function x(t) defined on an interval [t o , T ], the Riemann−Liouville fractional integral of x(t) of order q > 0 is defined by and Riemann-Liouville fractional derivative of x(t) of order q > 0 defined by where n − 1 < q < n and Γ(q) is the gamma function. The definition of Caputo fractional derivative of order q > 0 is The relation between Caputo and Riemann-Liouville fractional derivatives is illustrated by the following expression: Moreover, Caputo derivative of x(t) = c (constant) is always 0, whereas in the case of zero value of the lower terminal t o , Riemann-Liouville fractional derivative of However, by setting t o −→ −∞ in both definitions and requiring reasonable behavior of x(t) and its derivatives, we end up with the same formula for both of them as the following, From Parra-Hinojosa et al. [33] and Tavazoei [39], the fractional derivatives of sine or cosine functions in the Riemann-Liouville sense have no periodic behavior and they are, where E q1,q2 (t) is the two-parameter Mittag-Leffler function, , q 1 > 0, q 2 > 0.
When the time tends to go to infinity which is corresponding to Weyl's definition of fractional derivative, then Eq. (12) converges to the following periodic functions: By setting t o −→ −∞ in Eq. (7) which is corresponding to Liouville's definition of fractional derivative, it also gives a periodic result, i.e., −∞ D q t cos(ωt) = ω q cos(ωt + When q = 1 and k = 0 in Eq. follows: This linear equation (Eq.(16)) is widely taken as a prototype model of analytical stability of zero solution for several studies involving constant past history, cf. Diekmann et al. [12], Michiels and Niculescu [28] and Smith [37]. In Hayes [14], the stability of the zero solution is elucidated and the stability chart of the parameters space (a, b) is drawn. In particular, it is showed that in Barwell [4] and Koto [21], if a and b satisfies where {a} is the real part of a, then the zero solution of Eq.(16) is asymptotically stable for any τ ≥ 0. However, the stability analysis of Eq.(16) w.r.t. the parameters a and b is largely discussed, cf. Bellen and Maset [5], Breda [8], Diekmann et al. [13] and Insperger et al. [15], and the resulting stability chart or the transition curves exhibiting that (a, b)-plane is divided into stable and unstable regions as shown in FIGURE 2. The significant modification into Eq.(16) is to include an integral term with time delay turns into delay integro-differential equation introduced in the following formula where τ ≥ 0 and (a, b, k) ∈ R 3 . Obviously, by setting k = 0 turns Eq.(18) into Eq. (16). As mentioned in Koto [20,21], the zero solution of Eq.(18) is asymptotically stable for any τ ≥ 0 if and only if the space (a, b, k) ∈ R 3 satisfies the following: where λ is the characteristic exponent of the system. When the space (a, b, k) ∈ R 3 and k = 0, these conditions are reduced to the simple condition a < 0, k < 0 and b 2 < a 2 + 2k. In Insperger et al. [16], the fractional version of Hayes equation is tackled, having different treatments employed to obtain the transition curves on the first Riemann sheet for different values of fractional order derivative. However, delayed fractional integro-differential equations would provide valuable tools for modeling physical phenomena by giving explanations behaved alike real processes. Computational difficulties that face researchers in solving many of delayed fractional integro-differential equations have leaded them to use approximations and numerical methods. Therefore, different studies have been done to build up approximate methods for giving approximate solutions, cf. Li et al [22,23], Ren et al. [35] and Zarebnia [41]. Motivated by the previous arguments, in this work we consider the stability analysis of the extended fractional order linear integro-differential equations with time delay τ , where 0 < q ≤ 2 for τ ≥ 0 and the space (a, b, k) ∈ R 3 with one equilibrium point x * = 0. Eq.(21) can be easily transformed to the following delayed fractional differential system A ∈ R n×n and A τ ∈ R n×n . The initial value in Riemann-Liouville sense is given by and in Caputo's sense, it is given by where the history To consider the stability problem in this work, it is necessary to introduce the stability definition related to the system (Eq.(22)) as follows, cf. Qien et al. [34] and Zhang et al. [42].
The zero solution of the system (Eq. (22)) with orders q * = (q 1 , q 2 , ..., q n ) and q j ∈ (0, 1], j = 1, 2, ... is said to be i) stable if, for any initial value x o , there exists an > 0 such that the assumptions hold locally on D). iii) (globally) asymptotically stable if, in addition to being stable, ||x(t)|| → 0 as t → ∞ if x ∈ R n (i.e. the assumptions hold globally on R n ).
In the following, the main work is explicated in the section of main results which is organized as follows: the first two parts present general conditions of zero stability by Liapunov function for the classical problem and the properties of the characteristic function for the fractional version. The third part presents necessary conditions for the existence of the special analytical solutions. The fourth part presents the stability analysis of the fractional order version with different approximation techniques. The fifth part presents the construction of stability charts and transition curves. In the final part, the impact of the work in the shown RLC circuit is discussed. In the last, some concluded remarks are given.

Stability Analysis Of The Classical
Version. General conditions of stability can be obtained in its classical case of Eq.(21) using the second method of Liapunov, cf. Merkin [27]. The necessary conditions to obtain the asymptotic stability of the zero solution and the existence of a periodic solution are stated in the following theorem.
Theorem 2.1. Necessary condition to obtain a stable or an asymptotically stable zero solution for is a + b + kτ ≤ 0. Moreover, a periodic solution might be existed for some values of time delay(τ ).
Proof. Now, consider the functional with u = t + v, Differentiate w.r.t time and then we obtain (27) Using the fact that Then, we obtain In order to take dV dt ≤ 0, choose that V (x) is bounded from below and from the Cauchy criteria for convergence of it follows that as t → ∞. Since dV dt ≤ 0 and this depends on (a, b, k) and τ such that Thus, the zero solution might be stable or asymptotically stable. In case of k = 0, the functional can be taken as Hence, the asymptotic stability can be obtained if If the solution is represented by x(t) = x(t, t o , ψ) and then integrate Eq.(36), so that one obtains and this yields Hence, in general, it is easily to conclude that In general, the basic rule to obtain asymptotically stable solution is that, the real parts of the eigenvalues(λ = α + iω) of the characteristic function( i.e. roots of Λ(λ) = 0) should be negative, cf. Wilson [40]. So that, for given values of τ and (a, b, k) a necessary condition for Hopf bifurcation to occur at the equilibrium solution of Eq.
Proof. If we substitute λ = iω into Eq. (40), we obtain where from Eq.(43) we can get the following Proposition 2.5. Suppose k = 0, then Eq.(21) and Eq.(40) turn into a given values of τ and (a, b) is a solution to Eq.(46b) if and only if the following pair of real equations hold: Proof. Substitute λ = iω into Eq.(46b), we obtain Now, we are going to find some special explicit solutions for Eq. (21) under certain conditions on the parameters: τ , q and the space (a, b, k).
(1) If ωτ = 2nπ and |b| > |a|, then there exists a real value ω such that x(t) = exp(iωt) is a periodic solution for Eq.
2.4. Stability Analysis. Some approximation techniques are employed to obtain some values for the characteristic exponents(λ s) of the characteristic equation Eq. (40), since there are infinite number of them due to the existence of time delay(τ ) and fractional order (q) under certain conditions. The strategy here is to find a suitable better approximation for λ q and exp(−λτ ) expansions. So that, we can obtain an approximate quadratic algebraic equation for Eq.(40) with two dominant eigenvalues and all the remaining eigenvalues are clustered around these two eigenvalues when time delay is going to be very small(0 < τ 1).
Theorem 2.7. Consider the delayed fractional Integro-differential equation Eq. (21) so that (0 < q ≤ 1). Suppose that τ 1, {λ} < 1 τ . If it is considered that, then the system represented by Eq. (21) is stable under the following condition If bτ ≈ 0 and q = 1 then, the classical case condition of stability holds Proof. Suppose that the eigenvalues could be obtained easily as follows where λ 1 and λ 2 are the dominant eigenvalues that all the remaining eigenvalues are clustered around both of them when the time delay is very small, after that we used both of Eq.(63) to have a quadratic equation. The system is stable when λ 2 < 0, and unstable when λ 2 > 0. From FIGURE 3, it is shown that when τ 1 and {λ} < 1 τ , and Eq.(59) is supposed to be a better approximation for λ q than the first approximation of λ q that was introduced in Atherton [1] as follows Now, we are going to do the same as discussed in the previous theorem but the Figure 3. The Comparison shows that Eq.(59) has a better approximation for λ q than Eq.(66).
case of k = 0 is illustrated in the following theorem.
Proof. First of all, we are going to expand approximations of exp(−λτ ) and ln(λτ ) to the second order in order to obtain a quadratic approximate algebraic characteristic equation, thus we have the following eigenvalues are obtained as in the following, where a better approximation for λ q than Eq.(66).

Stability
Charts. The construction of stability charts for Eq.(21) using the D-subdivision method for the space (a, b) under different values of k and q is illustrated. The construction is based on the numbers of stable/unstable characteristic exponents using the following analytic ways.
Proposition 2.9. If we set b = k = 0 in Eq. (21), then stability of the scalar fractional-order differential equation Eq.(71) alternates from stable to unstable and vice versa depending on the value of fractional order (q).
Proof. Assume that b = k = 0, then Eq.(21) reduces to the scalar fractional-order differential equation Eq.(71). With the characteristic function when the characteristic function Eq.(72) vanishes, we obtain its characteristic equation which is consequently, we have infinite characteristic exponents. In the following cases, we illustrate that the stability of the system introduced by Eq.(71) is affected by the value of the fractional order (q).

Case 1.:
If q = 1, then the system Eq.(71) reduces to a characteristic equation λ = a which is asymptotically stable when a < 0, and unstable when a > 0.

Case 2.:
If q = 1 2n , then Eq.(73) turns into λ = a 2n > 0, hence the system is unstable for all a, where n ∈ Z + . Case 3.: If q = 1 2n+1 , then Eq.(73) turns into λ = a 2n+1 , hence the system is asymptotically stable when a < 0, and unstable when a > 0, where n ∈ Z + . Case 4.: If q = m n , then Eq.(73) turns into λ = m √ a n . If m and n are even numbers, then a n > 0 and one of its m th roots is positive, one is negative and the rest (when m > 2) are complex. Hence the system is unstable. Case 5.: If q = 2m+1 2n , then Eq.(73) turns into λ = 2m+1 √ a 2n . As 2m + 1 is odd, then one of the (2m + 1) th roots is real and has the same sign as a 2n > 0, where the remaining 2m roots are not real. This implies that the system is unstable, where n and m ∈ Z + .
The presence of the two parameters b = 0 and k = 0 makes the stability properties of the system more complex. This is due to the nature of the delayed system, since it has an infinite-dimensional nature as well as the presence of the fractional order. Therefore, in the following theorem, we show that the effect of the changing of parameters b and k on the stability properties of Eq.(21). (1) If qω q−1 sin πq 2 + bτ cos(ωτ ) + k τ sin(ωτ ) ω + cos(ωτ ) ω sin(ωτ ) > 0, then there exists a number b cr such that the equilibrium x * is asymptotically stable for b < b cr and unstable for b > b cr .
ω > 0, then there exists a number k cr such that the equilibrium x * is asymptotically stable for k < k cr and unstable for k > k cr .
Because of the continuity of the function Λ(λ) with respect to changes in the system parameters (a, b, k, τ ), these curves divide the parameter plane (a, b) into stable and unstable regions separated by the transition curves. In some cases, it is noted that the numbers of unstable characteristic exponents are constant and typically the change of these numbers along the curves can be determined by the analysis of the exponent-crossing direction. This change is measured by the sign of the partial derivative of eigenvalue real part, α, with respect to one of the system parameters, cf. Inspeger et al. [15].
Normally, ω appears in the denominator of the transition curves in Eq.(77), thus ω can not vanish. Hence, we can not get a form for transition curve associated with a real critical characteristic exponent crossing the imaginary axis at 0. The transition curves given by Eq.(77) are associated with a complex conjugate pair of eigenvalues in the form λ = ±iω. Now, we are going to study the effect of b on stability by taking the partial derivatives of Eq.(75) and Eq.(76) with respect to b and considering that α = 0 along the transition curves to obtain − sin(ωτ ) = −qω q−1 cos πq 2 − bτ sin(ωτ ) + k τ cos(ωτ ) ω − sin(ωτ ) where α b and ω b are the partial derivatives of α and ω with respect to b. The solution of Eq.(78) and Eq.(79) for α b is where, Hence, sign qω q−1 sin πq 2 + bτ cos(ωτ ) + k τ sin(ωτ ) ω + cos(ωτ ) Obviously Eq.(82) tells us if α b > 0, all the roots crossing the imaginary axis at iω cross from left to right as b increases, i.e., α get closer to positivity and farther from negativity thus after a critical value of b stable characteristic exponent becomes unstable and this results in loss of stability. Now we are going to study the effect of k on stability by taking the partial derivative of Eq.(75) and Eq.(76) with respect to k and considering that α = 0 along the transition curves to obtain 1 ω (cos(ωτ ) − 1) = −qω q−1 cos πq 2 − bτ sin(ωτ ) + k τ cos(ωτ ) ω − sin(ωτ ) ω 2 α k + qω q−1 sin πq 2 + bτ cos(ωτ ) + k τ sin(ωτ ) ω + cos(ωτ ) where α k and ω k are the partial derivatives of α and ω with respect to k. The solution of Eq.(83) and Eq.(84) for α k is where β is defined in Eq.(81), hence sign [α k ] = sign sin(ωτ ) ω qω q−1 sin πq 2 + bτ cos(ωτ ) + k τ sin(ωτ ) ω + cos(ωτ ) Obviously Eq.(86) tells us if α k > 0, all the roots crossing the imaginary axis at iω cross from left to right as k increases, i.e., α get closer to positivity and farther from negativity thus after a critical value of k stable characteristic exponent becomes unstable and this results in loss of stability.
Remark 2.11. One can do the previous work and study the stability of Eq.(46a) w.r.t b by setting k = 0 to obtain the following; sign [α b ] = sign qω q−1 sin πq 2 + bτ cos(ωτ ) cos(ωτ ) if α b > 0, all the roots crossing the imaginary axis at iω cross from left to right as b increases, i.e., α gets closer to positivity and farther from negativity thus after a critical value of b stable characteristic exponent becomes unstable and this results in loss of stability.  13. In each figure, the angular frequency limits (ω) along the transition curves are shown. From these figures, one can notice that the angular frequency limits have a rather change due to the change of k and q, and it is correspondingly followed by the change of stable and unstable regions separated by transition curves. Particularly, it is resulted that, the stability state of a system can be influenced by the fractional order derivative terms and the transition curves as well.
2.6. The impact in fractional RLC circuits. By applying the deduced results to RLC shown in FIGURE 1 as an example to look for the impact of the work on the applied sciences. So then, it can be deduced the following states of the circuit under dissipative resistances R a and R b as follows, Case (i): The state of the circuit has stable behavior at q = 1 if Case (ii): The state of the circuit has stable behavior at q < 1, τ < 1 if Case (iii): The state of the circuit has periodic behavior at q = 1 if and τ < 1.

CONCLUDING REMARKS
In this work, the stability analysis of delayed linear integro-differential equations within the existence of fractional derivative operators and their impact in RLC circuit applications is studied. Such analysis is based on new approximations of fractional powered eigenvalue (λ q ) and the transcendental terms around the small values of time delay. It is resulted that, the deduced conditions to guarantee the stability generalize the classical cases conditions as well as the other corresponding results of fractional cases. The stability charts using D-division method are constructed and the influences of the fractional derivative on the transition curves are shown. Moreover, the necessary conditions for the existence of some explicit analytical solutions for the given system are concluded.