Fractional differential equations are those equations in which an unknown function appears under a fractional (or non-integer) order differentiation operator. Interest in the applications of fractional order differential equations has been increasing steadily during the last few decades. Differential equations of fractional (non-integer) order and the related fractional integrals have found applications in diverse areas such as chemistry, biology, physics, finance, pharmacology, electrical circuits, heat transfer, diffusion processes, material science, visco- elasticity etc. One of the main advantages of a fractional (non-integer) order derivative is that it provides a mechanism for the incorporation of memory and hereditary properties of various materials and phenomena and this aspect is briefly considered in the next section. While an integer order derivative of a function considers the local behaviour of the function, a fractional (non-integer) order derivative depends on the previous history of the function values (see for instance the books and survey articles from  -  ).
Fractionalization of a dynamical system governed by integer order dif- ferential equations means a reformulation of such a system with a system containing one or more fractional order (non-integer order) derivatives. In the analysis of dynamical characteristics of a fractionalized system, we have to choose a suitable definition of a fractional derivative. In fractional calculus, there are several definitions (see the next section below) of a fractional deriva- tive usually named after some pioneers of fractional calculus; to name a few we have Grunwald-Leitnikov, Riemann-Liouville, Weyl, Riesz and Caputo deriva- tives. In the case of integer order calculus, the derivative of a constant is zero. We have a similar result for the Caputo derivative but not so in some other cases. Also in the case of initial value problems of integer order calculus, the initial values correspond to the respective initial values of the state variables and this is the case for the initial value problems with the Caputo fractional derivatives. However this is not the case with initial value problems with the Riemann-Liouville fractional derivatives. Riemann-Liouville derivatives are widely studied in the mathematical literature due mainly to the development of fractional calculus. We note that the fractional order Riemann-Liouville derivative of a constant is not zero. We remark that when the initial values of state variables are zero, the fractional Caputo derivative of a function is the same as that of the Riemann-Liouville dervative, under some smoothness of functions involved (see  ).
The purpose of this article is to formulate a fractionalized version of a simple dynamical system defined by a scalar semi-linear integer order differential equation. While there are many ways with which one can fractionalize an integer order differential equation, we propose one of the simplest methods by which an ordinary integer order derivative is replaced by a fractional (non-integer) order one (more details are presented in section 5 below) and then study the global existence of a unique solution along with the stability of its equilibrium solution when it exists. Our method is based on an application of a symbolic operational calculus to derive an analogue of the variation of constants formula well known in the theory of integer order differential equations. Such a variation of constants formula for the fractionalized initial value problem is shown to reduce to the well known formula of the integer order system when the fractional order approaches an integer; we use the Fourier integral theorem for this purpose. We apply a novel method for the existence and stability of solutions by exploiting the complete monotonicity of the relevant Mittag-Leffler functions. In the opinion of these authors this appears to be the first time that such a method has been proposed in the study of the stability of fractional (non-integer) order differential equations; we remark that a comparison principle and the positivity of certain Mittag-Leffler functions have been used by  in their stability considerations of fractional differential equations.
2. Some Preliminaries
In this section we recall and collect for the convenience of the reader a number of definitions and properties of the fractional order derivatives used in the analysis of our fractionalized differential equation of interest. For more details of fractional order differential equations and their application we refer to the works of       -  .
Definition 1 The fractional integral of order of a Lebesgue measurable function is defined by the operator given by
provided the integral exists pointwise on where denotes the gamma function.
We remark that a sentence similar to “fractional order integrals and fractional order derivatives provide excellent instruments for the description of memory and hereditary properties of various materials and processes” appears in numerous publications. (See for examples) The precise meaning of such a statement is not made clear satisfactorily and so we provide the following discussion of the memory aspects of the fractional integral operator.
The existence of memory in a dynamic system means that if at time s the state of a system is given by then the time evolution of the system whose magnitude at time is given by a convolutionary integral of the form
in which is known as the memory kernel (see for instance  ). We can write the fractional integral in (2.1) in the form
which expresses the integral as the weighted average of past states . The weight or memory kernel changes with the order of the integral operator. When and s is large with the weight is large implying that the distant past states enhance or increase the memory and this is not natural in realistic systems in the opinion of this author. On the other hand if , the distant past effects with s large and , are given smaller weights in comparison with the more recent ones. This illustrates that when , the memory is fading and the kernel demonstrates a fading memory”. When the memory kernel becomes a step function indicating that the integral has perfect memory which gives the same weight for all previous states in the integral. Now when the order of the integral tends to zero, the singularity of the term near zero is neutralized by the Gamma function with . Thus as , tends to a Dirac delta function leading to
which indicates that when , the system has no memory of past states implying that it has no memory of its past history except the current one. We have already seen that when , the system has partial or fading memory, in the sense that with each increment of the time the current state is given a maximum weight and the past states are given weights decreasing like that of the power law of the type . We conclude our discussion of the memory effects of a fractional integral with the comment that the order of a fractional integral is an index of memory” (see  ) of the respective fractional operator and the order interpolates between a perfect memory ( ) and absolutely no memory ( ) of the dynamic process.
Definition 2 The Caputo fractional derivative of order is defined by the operator where
We note that the Caputo derivative of a constant is zero as it is the case with integer order derivatives. One requires the integrability of the integer order derivative for the existence of the Caputo derivative. This is not the case for another fractional order derivative known as the Riemann-Liouville fractional derivative defined by
Riemann-Liouville derivative is more general and a function need not even be continuous for the existence of this derivative. We will not consider the Riemann-Liouville derivative in this article. It is found from the definitions of fractional derivatives, that the expressions of derivatives contain convolutional integrals with power law kernels; such convolutions account for the memory effects of the processes modelled by the fractional derivatives. The integer order derivatives have a local character whereas the fractional order derivatives have a non-local nature in the sense that the derivative at some value of the independent variable depends on the past history of the function values. This is one of the major differences of a fractional derivative from its integer order counterpart. As mentioned in the introduction, fractional order derivatives incorporate memory in the system it represents. We mean by memory that a current state of a system depends on some or all of the previous states of the system with certain weights attached to the past states. Fractional derivatives with respect to time are characterized by the power law decay of memory manifested by the kernels in the convolutionary integrals in their definitions. For example in the definition of
the Caputo derivative in (2) the memory kernel is given by and this
corresponds to a fading memory since the kernel decays with time. Thus a fractional derivative incorporates a fading memory although the remnants of the memory lingers on for ever with decreasing intensity and this has been interpreted as fractional derivatives having infinite memory. Memory is a characteristic of all living organisms including humans; however many (non- living) materials also have some kind of memory of their past behaviour and this has been examined by  in an article entitled Dead matter has memory. We will now show that the fractional order indicates the strength of memory in the sense that when , the fractional order derivative tends to that of the integer order derivative which has no memory.
We shall briefly examine the behavior of as where . We suppose that exists for in the following derivation.
Since the integrand in the integral of (4) is not singular, we have the following:
Thus if the fractional order corresponds to the memory in the fractional derivative, such memory is lost when as the fractional derivative becomes an integer order derivative. It is not known to this author how to relax the restrictive requirement of the existence of the second integer order derivative in the above derivation of the limiting consistency of the fractional order derivative.
The well known rules for integer order differentiation such as the product rule and chain rule are not applicable for the fractional order differentiation; although such rules exist for fractional order differentiation in the form of infinite series, they are not as useful as in integer order differentiation. In integer order differential and integral calculus, derivatives can represent a rate of change, direction of increase or decrease of a function; an integral can mean an area under the graph of a function. Such interpretations are not applicable for fractional order derivatives and integrals. For the convenience of readers uninitiated in fractional calculus, this article is formatted in a self-contained fashion at the cost of minimal repetition.
One of the properties of fractional integrals we will use in the proof of the equivalence of a fractional differential equations and an associated fractional integral equation is the following semi group property of fractional integrals.
This is easily seen as follows: we first recall a relation between Beta and Gamma functions given by
in which the familiar integral representation of the Gamma function
of elementary calculus is used. Using the definition of the fractional integral, we have for , and let,
The proof of (8) is complete. It follows from the definition of Caputo derivative that for arbitrary constants a and b,
indicating that is a linear operator.
Another important property of which we use later is that for ,
This can be verified as follows:
It is found from the above that is not a left inverse of if .
We shall now show that is the left inverse of and this fact will be useful later in the conversion of a fractional order differential equation to its equivalent integral equation. First we will verify the following:
In fact we have from the definition of ,
which leads to (12).
Lemma 1 The fractional integral and the Caputo-derivative satisfy the following:
and note that
Now operate on both sides of the above by so that
and the proof is complete.
Throughout the following, we consider only the Caputo-derivative and we use for convenience the notation
3. Mittag-Leffler Functions
The exponential function where is a con-
stant appears naturally in the study of integer order differential equations due to the easily verifiable property,
In the study of fractional order differential equations, Mittag-Leffler functions (see    ) play a prominent role like that of the exponential functions in the case of integer order differential equations; It also serve as a tool for discrete fractional modelling (see  ) the one parameter Mittag-Leffler function with real parameter is defined by means of a power series for complex z given by
with the Gamma function given by
The two parameter Mittag-Leffler function and the three parameter Mittag-Leffler function are respectively defined by
where denotes the Pochhammer symbol defined by
One can verify that the functions and are entire functions defined for all complex z by using the ratio test for convergence and the fact that the ratio of gamma functions satisfy an asymptotic estimation given by 
Our analysis below requires Laplace transforms of fractional derivatives and Mittag-Leffler functions. Let denote the Laplace transform of defined for and assumed to have a Laplace transform defined by
We derive a few known results on Laplace transforms which we will use below in the stability analysis of our fractionalized system.
We recall some elementary and known facts regarding Laplace transforms (see also  ). If is a differentiable function whose laplace transform exists then
We consider the Laplace transform of the Mittag-Leffler function by using the uniform convergence of the series in the definition of the function as follows:
One can derive in a similar way that
By analytic continuation the Laplace transform is usually extended for complex values of s over the entire complex plane except for a set of singular points and in some cases the function turns out to be multi-valued having branch point singularities as in the case of the function
If is known then one can get the function by means of a complex inversion formula given by (see  )
in which the integration is performed along a vertical line in the complex plane where ; the real number is chosen such that the line lies to the right of all the singularities including the branch points, poles, essential singularities; the line of integration can otherwise be arbitraty; the remarkable property of this integral is that the value of the integral is independent of the number . We choose a closed contour in the complex plane known as the Hankel-Bromwich contour shown in Figure 1 below so that the theory of residues can be used for the evaluation of the inversion integral above. We have from the Cauchy’s residue theorem that
Figure 1. Hankel-Bromwich contour.
since there are no singularities of the integrand in the interior of the contour as
discussed in the following. We note that the only poles of are on the
negative real axis and there are two branch points respectively at the origin and at . For instance let us set . Then
since , which implies that there are no poles on the upper half complex plane. Similarly for we have implying that and hence there are no poles in the lower half of the complex plane. For , we have which implies that
. Thus it follows that the only singularities of are
on the negative real axis of the complex plane. One can then determine the inverse Laplace transform by the method of residues and the Cauchy’s theorem
from complex analysis. Since the origin is a branch point of
due to the multivalued nature of the function we cannot use the usual Bromwich contour. We have to make a branch cut along the negative real axis beginning from just above the negative real axis going around the origin in a clockwise sense proceeding to just below the negative real axis and then follow a Bromwich type contour. Such a contour is known as a Hankel- Bromwich contour.
We proceed to evaluate the limits of certain line integrals involved in the computations of the inverse Laplace transform of by means of the complex contour integral. We begin with
On HK, we let and obtain
In a similar way one can show that
Next we consider the integral along the arc KA as follows:
In a similar way one can show that
Since the only singularities are poles on the negative real axis and this has been cut out of the complex plane and branch points of the integrand and are also excluded for the integration along the Hankel path in the contour integral above. Thus we have from the above equation that
One can show that the integrals along the arc BCD goes to zero as the radius of this arc goes to 0. Thus we are left with the contributions along the line segments AB and DE. Thus we have from the above that
Thus we have that
This formula is referred to as Titchmarsh’s theorem’ in the literature on fractional calculus (see   ). However in the literature on integral transforms there is a result known as Titchmarsh’s theorem’ (see  , p. 68) which is different from what appears in fractional calculus by that name.
It follows from the above that when
We shall briefly consider the limiting case of the formula (21) as . First we let in (21) and obtain
A direct attempt to evaluate the limit of the right side of (22) as leads to an indeterminacy. We will follow an indirect procedure to evaluate this limit. We use the formula
where and the Fourier integral formula (see Sneddon ) in the form
for a suitable function . One can have from (24) that
We note from (25) that the inner integral behaves like that of the Dirac delta function (see  ) and thus we can write
where denotes the Dirac delta function. For convenience in the following we define as follows
and obtain from (24) that
which together with (22) and (23) leads to
We are now ready to evaluate the limit of 21) as by using the representation in (29).
Thus we capture the otherwise known result from the definition of for that
4. Complete Monotonicity
It can be found from its definition that provides a generalization of the exponential function since we have from the definition of that
We will see below that preserves an interesting property of the exponential function for namely
for . The function is completely monotonic for and . In general a smooth function for is said to be completely monotonic if is positive and its integer order derivatives satisfy the condition
The following result characterises completely monotonic functions (see for instance      ).
Theorem 1 (Hausdorff-Bernstein-Widder) A smooth function for is completely monotonic if and only if
where is a Borel measure on ; that is there exists a non-negative distribution for satisfying the relation
in which the integral converges and is infinitely differentiable in the region of convergence of the Laplace transform.
Often one uses the above representability for proving the complete mono- tonicity of a function by finding the relevant function . In fact a sufficient condition for the complete monotonicity of a function is the existence of a positive locally integrable function such that a relation of the type in (35) holds (see  , pp. 61-72). Such a function is known as the spectral distribution of the function . In our analysis below we use the positivity of the functions and for . For this purpose we obtain the spectral represen- tation of these two functions from their respective inverse Laplace transforma- tions. For instance we have from the previous section that
It is not difficult to verify that for . It will follow that for and which by Theorem 4.1 implies that is completely monotonic.
Next we consider a special case of the two parameter Mittag-Leffler function and show that is completely monotonic. We begin from the observation that
which can be extended for complex s by analytic continuation. As before by using the Hankel-Bromwich contour, we can derive that for ,
where for we have the following:
The complete monotonicity of follows with its spectral distribution given by (43). In fact vanishes when is an integer and ; for all when . One of the consequences of the integral representation of is the derivation of asymptotic estimate of for large values of ; for instance we derive the following asymptotic estimation for ;
Thus we have that as for
on using the Euler’s reflection formula
which indicates that the Mittag-Leffler function decays for large values of t as a power law and the decay is not exponential as it will be in the case of integer order derivatives. We note that the asymptotic estimation in (46) can be rigorously justified by using Watson’s lemma for loop integrals as discussed at the end of this article (see also for instance pp. 82-83 of  ).
5. A Variation of Constants Formula
One of the well known and widely used methods of studying initial value problems for integer order differential equations is by converting the differential system into an equivalent Volterra integral system; subsequently one studies the integral system equivalent to the original differential system. This method has been followed by many authors in the analysis of non-integer order differential systems. However the equivalence of the fractional order differential system and the corresponding integral system is not satisfactorily established in many cases as it has been pointed out in the survey article by  . In this section we derive a variation of constants formula for a semi-linear fractionalized version of a class of dynamical systems governed by a class of semi-linear differential equations. It is elementary to show that the solution of the initial value problem for the integer order system
where are real constants and is a real valued continuous function of the real state variable is given by
The Equation (48) is known as a variation of constants formula for (47). The term appearing under the integral in(48) is a solution of the linear homo- geneous initial value problem
and the solution of (49) has the semi-group property
The solution of the fractional order initial value problem
when is given by the one parameter Mittag-Leffler function . For instance one can verify this as follows by using the definition of
and the uniform convergence of the series defining the Mittag-Leffler
We remark that we have used the following relation between Beta and Gamma functions
in the evaluation of the fractional derivative and subsequent simplification. When , the function does not have the semi-group property (50) and it is known that (see  )
As noted by Diethelm  that the absence of a formula of the type for Mittag-Leffler functions has been one of the stumbling obstacles for the development of a comprehensive theory for linear fractional differential equtions.
It is the purpose of this section to obtain a variation of constants formula of the type in (48) for the fractional dynamic system governed by the fractional order semi-linear initial value problem,
We will use a method of symbolic calculus for this purpose similar to that of Babenko’s method (see Podlubny ). We apply the fractional integral operator to both sides of (53) so that on using (11)
and hence we have an operational equation of the form
in which I denotes an identity operator. We use the following operational identity (see for instance the Equations (2.7) and (2.8) in  )
and observe that
We can verify (58) as follows by using (9)
The other identity
can be similarly established. We have from (55) that
where the order of summation and integration has been interchanged which is justified due to the fact that the infinite series is uniformly convergent for finite values of the terms in the infinite series when is assumed to be continuous. One can also show that if is a solution of the integral Equation (58), then is also a solution of the fractional order differential equation of the initial value problem (53) by reversing the above sequence of steps. A solution of the fractional order integral equation satisfies the initial value of the fractional order initial value problem and for this it is sufficient to verify that
One can verify by direct computation that
By the continuity of , the complete monotonicity of and the mean value theorem of integral calculus we have
where is some number in the interval (0,1). Since , any solution of the fractional integral equation satisfies the initial condition of the initial value problem (53)
It will follow that if is a solution of the more general initial value problem
then it is a solution of the the singular Volterra integral equation
and vice versa. The Formula (61) is the variation of constants formula for fractional order systems of the type (60). It is well known that if the forcing term is independent of the state variable, one uses the method of Laplace transforms to derive a formula similar to that of (61). We note that one can use the above procedure to derive a variation of constants formula for more general semi- linear systems.
It is of some interest to examine whether the Formula (61) leads to the usual variation of constants formula of the integer order differential equation corresponding to the case . We have already seen from (30) that
We have to examine the limiting value of as in the integral in (61). It has already been shown at the end of Section 2 that the
fractional derivative in (60) approaches the integer order derivative
when . It will follow from (40)-(43)
Proceeding as we examined the case of , in Section 3 (see Equation (3.17)) we can obtain that
so that the Formula (61) coincides with
of the integer order differential equation
Thus the variation of constants formula of the fractionalized semi-linear differential equation reduces to the usual variation of constants formula of the integer order differential equation when the order of the fractional derivative approaches an integer.
6. Fractionalization, Global Existence and Stability
In this section we consider the dynamics of a fractionalized version of the initial value problem
It is elementary that the above integer order differential equation is equivalent to the integral equation
One of the methods of fractionalization of the above initial value problem is to replace the above integral equation by a fractional integral equation of the form
which is the same as
and this is equivalent to by virtue of Lemma 2.3,
which leads to
It is found that this method reduces to the replacement of the integer order derivative by fractional order derivative. Such a method has been predominantly used in the literature on fractional order differential equations. We remark that the integral equation
can also be fractionalized by the fractional integral equation
where and and this fractional integral equation is the same as that of
and this is not equivalent to the replacement of integer order derivative by a fractional order derivative. Thus there is no unique way of fractionalization of initial value problems with integer order differential equations. In this article we study the dynamics of the fractionalized version where and the more general case where will be considered in a future article.
One of the methods used in the theory of integer order differential equations to establish the existence of globally defined solutions is based on an application of technique proposed by  . In this method one usually modifies a norm in a suitable Banach space and use the induced metric in the application of the well known contraction mapping principle. The modified norm is equivalent to the original norm and a contraction in the new metric is also a contraction in the original metric. This method has been used by several authors (see     ).
Let denote the linear space of real valued continuous functions defined on where T denotes an arbitrary positive number. For , we define a norm given by
for a suitable positive number which will be determined below. We remark that the norm is equivalent to the usual supremum norm defined by where
and hence with the metric induced by the space becomes a complete metric space. We recall that by the variation of constants formula, every solution of the fractional order initial value problem
is a solution of the fractional order integral equation
and every solution of (70) is a solution of (69). We can now proceed to formulate the following global existence of solutions of (69).
Theorem 2 Let a be a positive number and let . Let b be any real number. Suppose that is continuous and satisfies a global Lipschitz con- dition such that
for some positive number K. Then the integral equation (70) has a unique solution in B.
We define an operator F acting on the functions of B by the following:
From the properties of and the continuity of it will follow that F maps B into B. We have for all
By using the positivity of and the Lipschitz continuity of f, we have from (73) that
which leads to the following:
in which we have used the boundedness of the function on the positive real line so that
We choose the real number such that
from which it will follow that is a contraction and by the well known contraction mapping principle, has a unique fixed point which is a solution of (70) and hence that of (69). Since T is a positive, otherwise unrestricted number, the solution of (70) and hence that of (69) is defined on . This completes the proof.
The asymptotic behaviour as of solutions of differential equations and the stability of equilibria of dynamic systems have been well studied in the case of integer order differential equations. There are several articles devoted to the general analysis of stability of equilibria of non-integer order differen- tial equations (see for instance the articles by     -  and applications of delay fractional differentail equations to epidemic analysis(see   ) The asymptotic behaviour of solutions of fractional order systems can be quite different from the corresponding integer order systems indicating a change in the asymptotic behavior of a fractionalized dynamic system from its integer order counterpart. For instance an unstable integer order system can become a stable one when the system is fractionalized; also new asymptotic behaviour not found in the integer order system can emerge in fractionalized systems (see   We shall now proceed to consider the asymptotic be- haviour of solutions of the semi-linear fractional order system
We suppose that there exists a unique real constant say such that
Since the Caputo derivative of a constant is zero, the constant is an equilibrium of (78). Although several authors have investigated the stability of equilibrium solution of fractional order differential equations, there are not many investigations dealing with global asymptotic stability of nonlinear systems. In the following we use once again the complete monotonicity of the one and two parameter Mittag-Leffler functions and in our analysis of the stability of the equilibrium .
Definition 3 The equilibrium of (78) with is said to be globally stable if for arbitrary , every solution of (78) satisfies
The equilibrium is said to be globally asymptotically stable if it is stable and furthermore satisfies
In the case of fractional order differential equations there is another concept of stability which different from the exponential asymptotic stability common to the theory of integer-order differential equations.
Definition 4 The equilibrium of (78) is said to be Mittag-Leffler stable if every solution of (78) satisfies
where is locally Lipschitz with and p is some positive number.
Theorem 3 Suppose is a unique equilibrium of (78) with
Then the equilibrium of (78) is globally Mittag-Leffler stable.
Proof. We have from
and hence by the variation of constants formula
Now we use the non-negative and complete monotone character of the functions and and obtain from (84) that
Let be a real valued and non-negative function such that
By using the convolution theorem of Laplace transforms, we derive from (86) that
where denotes the Laplace transform operator and
We have from the known Laplace transforms of the one and two parameter Mittag-Leffler functions that
which simplifies to
We can invert the Laplace transforms in (90) and obtain
On invoking the complete monotonicity of and the non- negativity of , we derive from (91) that
which is the same as
from which the result will follow. The proof is complete.
We note that when the above inequality leads to the well known exponential asymptotic stability corresponding to the ordinary differential equation.
7. Some Remarks
We shall briefly consider the power law decay nature of the function for by using Watson’s lemma for loop contours (see for instance  , pp. 82-83 or  , pp. 250-252) suitable for asymptotic expansions. We let for convenience and note that
It will follow from the inverse of Laplace transforms for loop contours that
The notation is known as a loop integral where the path of integration starts at with , encircles the origin once in the counter- clockwise direction and returns to with . We have from (93) and (95) that
which for reads
The non-exponential and power law nature of the convergence of solutions to the equilibrium displayed by the above relation is ascribed to the long memory effect of the fractional derivative.
We conclude with the following observation regarding the fractionalization of systems of the form
For instance the following are some of the fractionalized versions of (96):
where , are constants; a further generalized fractionalization of (96) is given by
where is a continuous function. General results related to the existence and stability of solutions of equations of the form (97) and (98) will be useful and interesting developments to the existing literature on fractional order dynamic systems.
The project was supported by Small Research Grant, 2015-2016, Hong Kong Institute of Education.
 Butkovskii, A.G. (2010) Fractional Integrodifferential Calculus and Its Control Theoretic Applications: I. Mathematical Fundamentals and the Problem of Interpretation. Automation and Remote Control, 74, 543-574. https://doi.org/10.1134/S0005117913040012
 Butkovskii, A.G., Postnov, S.S. and Postnova, E.A. (2013) Fractional Integrodifferential Calculus and Control Theoretic Applications. II. Fractional Dynamic Systems: Modelling and Hardware Implementation. Automation and Remote Control, 74, 725-749.
 Debnath, L. (2003) Recent Applications of Fractional Calculus to Science and Engineering. International Journal of Mathematics and Mathematical Sciences, 54, 3413-3442.
 Monge, C.A., Chen, Y.Q., Vinagre, B.M., Xue, D. and Feliu, V. (2010) Fractional Order Systems and Controls, Advances in Industrial Control. Springer, Heidelberg. https://doi.org/10.1007/978-1-84996-335-0
 Hilfer, R., Luchko, Y. and Tomovski, Z. (2009) Operational Method for the Solution of Fractional Differential Equations with Generalized Riemann-Liouville Fractional Derivatives. Fractional Calculus and Applied Analysis, 12, 299-318.
 Machado, J.T., Kiryakova, V. and Mainardi, F. (2010) Recent History of Fractional Calculus. Communications in Nonlinear Science and Numerical Simulation, 16, 1140-1153.
 Rihan, F.A., Lakshmanan, S., Hashish, A.H., Rakkiyappan, R. and Ahmed, E. (2015) Fractional-Order Delayed Predator-Prey Systems with Holling Type-II Functional Response. Nonlinear Dynamics, 80, 777-789. https://doi.org/10.1007/s11071-015-1905-8
 Li, Y., Chen, Y.Q. and Podlubny, I. (2010) Stability of Fractional Order Nonlinear Dynamic Systems: Lyapunov Direct Method and Generalized Mittag-Leffler Stability. Computers and Mathematics with Applications, 59, 1810-1821. https://doi.org/10.1016/j.camwa.2009.08.019
 Bobylev, A.V. and Cercignani, C. (2002) The Inverse Laplace Transform of Some Analytic Functions with an Application to the Eternal Solutions of the Boltzmann Equation. Applied Mathematics Letters, 15, 807-813. https://doi.org/10.1016/S0893-9659(02)00046-0
 Pollard, H. (1948) The Completely Monotonic Character of the Mittag-Leffler Function . Bulletin of the American Mathematical Society, 54, 1115-1116. https://doi.org/10.1090/S0002-9904-1948-09132-7
 Gorenflo, R. and Mainardi, F. (1997) Fractional Calculus: Integral and Differential Equations of Fractional Order. Fractals and Fractional Calculus in Continuum Mechanics, 197, 223-276. https://doi.org/10.1007/978-3-7091-2664-6_5
 Bielecki, A. (1956) Une Remarque sur la Methode de Banach-Cacciopoli-Tikhonov dans la Theorie des Equations differentielles Ordinaires. Bulletin of the Polish Academy of Sciences. Mathematics, 4, 261-264.
 Baleanu, D. and Mustafa, O.G. (2010) On the Global Existence of Solutions to a Class of Fractional Differential Equations. Computers and Mathematics with Applications, 59, 1835-1841. https://doi.org/10.1016/j.camwa.2009.08.028
 El-Raheem, Z.F.A. (2003) Modification of the Application of a Contraction Mapping Method on a Class of Fractional Differential Equations. Applied Mathematics and Computation, 137, 371-374. https://doi.org/10.1016/S0096-3003(02)00136-4
 Klimek, M. (2011) On Contraction Principle Applied to Nonlinear Fractional Differential Equations with Derivatives of Order . Banach Center Publications, 95, 325-338.
 Baleanu, D., Ranjbar, A., Sadati, S.J., Delavari, R.H., Abdeljawad, T. and Gejji, V. (2011) Lyapunov-Krasovskii Stability Theorem for Fractional Systems with Delay. Romanian Journal of Physics, 56, 636-643.
 Delavari, H., Baleanu, D. and Sadati, J. (2011) Stability Analysis of Caputo Fractional-Order Nonlinear Systems Revisited. Nonlinear Dynamics, 67, 2433-2439.
 Li, C.P. and Zhang, F.R. (2011) A Survey on the Stability of Fractional Differential Equations. The European Physical Journal Special Topics, 193, 27-47. https://doi.org/10.1140/epjst/e2011-01379-1
 Li, L. and Wang, Z.J. (2013) Global Stability of Periodic Solutions for a Discrete Predator-Prey System with Functional Response. Nonlinear Dynamics, 72, 507-516.
 Qian, D., Li, C., Agarwal, R.P. and Wong, P.J.Y. (2010) Stability Analysis of Fractional Differential System with Riemann-Liouville Derivative. Mathematical and Computer Modelling, 52, 862-874. https://doi.org/10.1016/j.mcm.2010.05.016
 Rivero, M., Rogosin, S.V., Machado, J.A.T. and Trujillo, J.J. (2013) Stability of Fractional Order Systems. Mathematical Problems in Engineering, 2013, Article ID: 356215.
 Sadati, S.J., Baleanu, D., Ranjbar, A., Ghaderi, R. and Abdeljawad, T. (2010) Mittag-Leffler Stability Theorem for Fractional Nonlinear Systems with Delay. Abstract and Applied Analysis, 2010, Article ID: 108651. https://doi.org/10.1155/2010/108651
 Zhang, F., Li, C. and Chen, Y.Q. (2011) Asymptotical Stability of Nonlinear Fractional Differential Systems with Caputo Derivative. International Journal of Differential Equations, 2011, Article ID: 635165. https://doi.org/10.1155/2011/635165
 Arenas, A.J., Gonzalez-Parra, G. and Chen-Charpentier, B.M. (2010) A Nonstandard Numerical Scheme of Predictor-Corrector Type for Epidemic Models. Abstract and Applied Analysis, 2010, Article ID: 124812.
 Rihan, F.A., Baleanu, D., Lakshmanan, S. and Rakkiyappan, R. (2014) On Fractional SIRC Model with Salmonella Bacterial Infection. Abstract and Applied Analysis, 2014, 9.
 Ahmed, E., El-Sayed, A.M.A. and El-Saka, H.A.A. (2007) Equilibrium Points, Stability and Numerical Solution of Fractional Order Predator-Prey and Rabies Models. Journal of Mathematical Analysis and Applications, 325, 542-553. https://doi.org/10.1016/j.jmaa.2006.01.087