Fractional Differential Equations 2011

1 Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, QLD 4001, Australia 2 Department of Mechanical Engineering and Energy Processes, Southern Illinois University, Carbondale, IL 62901, USA 3 Department of Mathematics, The University of Jordan, Amman 11942, Jordan 4 School of Mathematics, Cardiff University, Cardiff CF24 4AG, UK 5 Department of Engineering Mechanics, Hohai University, Nanjing 210098, China

It is our pleasure to bring this special issue of the International Journal of Differential Equations dedicated to fractional differential equations FDEs .
In recent years, a growing number of work by many authors from various fields of science and engineering deal with dynamical systems described by fractional partial differential equations. Due to the extensive applications of FDEs in engineering and science, research in this area has grown significantly all around the world.
This special issue of Fractional Differential Equations consists of 13 original articles covering various aspects of FDEs and their applications by some of the prominent researchers in the field. These papers could be broadly grouped into three categories, namely, numerical and approximate schemes to solve fractional dynamical models 1st to 5th paper , existence and uniqueness of solutions of fractional differential equations and other theoretical results 6th to 11th paper , and application of fractional differential equations in various fields 12th and 13th paper . Other papers could also have been considered in this last category. However, because of their emphasis, they have been included in the first two categories.
The first paper introduces a new modified step variational iteration method for solving biochemical reaction model. The second and third papers use homotopy analysis method for solving a space-and time-fractional foam drainage equation and nonlinear coupled equations with parameters derivative, respectively. The fourth paper develops a new application of Mittag-Leffler function method and extends the application of the method to fractional linear differential equations. The fifth paper proposes an explicit numerical method for the

Introduction
The mathematical modelling of numerous phenomena in various areas of science and engineering using fractional derivatives naturally leads, in most cases, to what is called fractional differential equations FDEs . Although the fractional calculus has a long history and has been applied in various fields in real life, the interest in the study of FDEs and their applications has attracted the attention of many researchers and scientific societies beginning only in the last three decades 1, 2 . Since the exact solutions of most of the FDEs cannot be found easily, thus analytical and numerical methods must be used. For example, the ABMM is one of the most used methods to solve fractional differential equations [3][4][5] . Several of the other numerical analytical methods for solving fractional problems are the Adomian decomposition method ADM , the homotopy perturbation method HPM and the homotopy analysis method HAM . For example, Ray 6 and Abdulaziz et al. 7 used ADM to solve fractional diffusion equations and solve linear and nonlinear fractional differential equations, respectively. Hosseinnia et al. 8 presented an enhanced HPM to obtain an approximate solution of FDEs, and Abdulaziz et al. 9 extended the application of HPM to systems of FDEs. The HAM was applied to fractional KDV-Burgers-Kuromoto equations 10 , systems of nonlinear FDEs 11 , and fractional Lorenz system 12 .
Another powerful method which can also give explicit form for the solution is the variational iteration method VIM . It was proposed by He 13,14 , and other researchers have applied VIM to solve various problems 15-17 . For example, Song et al. 18 used VIM to obtain approximate solution of the fractional Sharma-Tasso-Olever equations. Yulita Molliq et al. 19, 20 solved fractional Zhakanov-Kuznetsov and fractional heat-and wavelike equations using VIM to obtain the approximate solution have shown the accuracy and efficiently of VIM. Nevertheless, VIM is only valid for short-time interval for solving the fractional system.
In this paper, we propose a modification of VIM to overcome this weakness of VIM. In particular, motivated by the work of 12 the procedure of dividing the time interval of solution in VIM to subintervals with the same step size Δt and the solution at each subinterval must necessary to satisfy the initial condition at each of the subinterval has been considered. Unfortunately, this idea does not give a good approximate solution when compared to the ABMM. Therefore, to obtain a good approximate solution which has a good agreement with ABMM, another idea is used: motivated by HAM, a nonzero auxiliary parameter is considered into the correction functional in VIM. This parameter was inserted to adjust and control the convergence region of the series solutions. In general, it is straightforward to choose a proper value of from the so-called -curve. We call this modification involving time step and auxiliary parameter the MoSVIM. Strictly speaking MoSVIM is a modification of our earlier proposed method-step variational iteration method-which is still under review 21 .
As an application, this paper investigates for the first time the applicability and effectiveness of MoSVIM to obtain the approximate solutions of the fractional version of the biochemical reaction model as studied in 22 for interval 0, T . The fractional biochemical reaction model shortly called FBRM is considered in the following form: subject to initial conditions u 0 1, v 0 0, 1.2 t − ξ q−θ−1 ∂ q u ξ ∂ξ q dξ, for q − 1 < θ < q, ∂ q u t ∂t q , for θ q ∈ .

2.4
For more information on the mathematical properties of fractional derivatives and integrals, one can consult 1, 2 .

Step Variational Iteration Method
The approximate solutions of fractional biochemical reaction model will be obtained in this paper. A simple way of ensuring validity of the approximations is solving under arbitrary initial conditions. In this case, 0, T is regarded as interval. From idea of Alomari et al. 12 , the 0, T interval is divided to subintervals with time step Δt, and the solution at each subinterval was obtained. So it is necessary to satisfy the initial condition at each of the subinterval. Thus the step technique can describe as the following formula: where λ i , for i 1, . . . , m, is a general Lagrange multiplier, L is linear operator, N is nonlinear operator, and g is inhomogeneous term. As knowledge, the optimal general Lagrange multiplier is obtained by constructing the correction functional as in VIM which is u i,n is considered as restricted variations, that is, δ u i,n 0. Accordingly, the initial values u 1,0 , u 2,0 , . . . , u m,0 will be changed for each subinterval, that is, u 1 t * c * In general, we do not have this information at our clearance except at the initial point t * t 0 0, but these values can be obtained by assuming that the new initial condition is the solution in previous interval i.e., if the solution in interval t j , t j 1 is necessary, then the initial conditions of this interval will be as follows: where c i , i 0, 1, . . . , m are the initial conditions in the interval t j , t j 1 .

Modified Step Variational Iteration Method
Furthermore, to implement the modification of SVIM, we consider / 0, a nonzero auxiliary parameter. Multiply by correction functional in 3.1 , yield where i 0, 1, 2, . . . , m, m ∈ and is the convergence-control parameter which ensures that this assumption can be satisfied. The subscript n denotes the nth iteration. Accordingly, the successive approximations u n t , n ≥ 0 of the solution u t will be readily obtained by selecting initial approximation u 0 that at least satisfies the initial conditions. The computations and plotting of figures for the algorithm, has been done using Maple package.

Application
In this section, we demonstrate the efficiency of MoSVIM od fractional biochemical reaction model in 1.1 . The correction functionals for the system 1.1 can be approximately constructed as used by VIM and 2.4 to find the general Lagrange multiplier in the following forms: where λ 1 and λ 2 are general Lagrange multipliers which can be identified optimally via variational theory. n denotes the nth iteration. u n , v n , and u n v n denote restricted variations, 6 International Journal of Differential Equations that is, δ u n 0, δ v n 0, and δu n v n 0. In this case, the general Lagrange multiplier can be easily determined by choosing the number of order q, that is, q 1. Thus, the following sets of stationary conditions was obtained as follows: 1 λ 1 t | ξ t 0, λ 1 ξ − λ 1 ξ 0, 1 λ 2 t | ξ t 0, βλ 2 ξ μ − λ 2 ξ 0.

5.3
Here, the general Lagrange multiplier in 5.3 is expanded by Taylor series and is chosen only one term in order to calculate, the general Lagrange multiplier can write as follows

5.4
Substituting the general Lagrange multipliers in 5.4 into the correction functional in 5.1 results in the following iteration formula:

5.5
Furthermore, we multiply the nonzero auxiliary parameter by 5.5 which yields:

5.6
Then, the interval 0, 2 is divided into subintervals with time step Δt, and we get the solution at each subinterval. In this case, the initial condition is regarded as initial approximation, Here, the iteration was chosen from previously research by Goh et al. 24 . Thus, the solution will be as follows: where t * start from t 0 0 until t J T 2. To carry out the solution on every subinterval of equal length Δt, the values of the following initial conditions is presented below: In general, we do not have this information at our clearance except at the initial point t * t 0 0, but we can obtain these values by assuming that the new initial condition is the solution in the previous interval i.e., If we need the solution in interval t j , t j 1 then the initial conditions of this interval will be as

5.10
where c 1 , c 2 are the initial conditions in the interval t j , t j 1 .

Result and Discussion
To investigate the influence of on convergence of the solution series, we plot the -curves of u 4 0.01 and v 4 0.01 using the fifth iteration of MoSVIM when θ 0.7, and θ 0.8 as shown in Figure 1. We found that the range of values for is between 0. 1  by Hashim et al. 25 . In this case, the computational algorithms for the system in 1.1 are written using the Maple software. A good solutions of fractional biochemical reaction model when 0.25 and θ 0.7 and θ 0.8 was presented in Tables 1 and 2, respectively. From the tables, MoSVIM is more accurate than SVIM in different value of θ, that is, θ 0.7 and θ 0.8. Figure 2

Conclusions
In this paper, an algorithm of fractional biochemical reaction model FBRM using step modified variational iteration method MoSVIM was developed. For computations and plots, the Maple package were used. We found that MoSVIM is a suitable technique to   solve the fractional problem. This modified method yields an analytical solution in iterations of a rapid convergent infinite power series with enlarged intervals. Comparison between MoSVIM, MoVIM and ABMM were made; the MoSVIM was found to be more accurate than the MoVIM. MoSVIM is easier in calculation yet powerful method and also is readily applicable to the more complex cases of fractional problems which arise in various fields of pure and applied sciences.

Introduction
Many phenomena in engineering, physics, chemistry, and other science can be described very successfully by models using the theory of derivatives and integrals of fractional order. Interest in the concept of differentiation and integration to noninteger order has existed since the development of the classical calculus 1-3 . By implication, mathematical modeling of many physical systems are governed by linear and nonlinear fractional differential equations in various applications in fluid mechanics, viscoelasticity, chemistry, physics, biology, and engineering. Since many fractional differential equations are nonlinear and do not have exact analytical solutions, various numerical and analytic methods have been used to solve these equations. The Adomian decomposition method ADM 4 , the homotopy perturbation method HPM 5 , the variational iteration method VIM 6 , and other methods have been used to provide analytical approximation to linear and nonlinear problems 7, 8 . However, the convergence region of the corresponding results is rather small.
In 1992, Liao 9-13 employed the basic ideas of the homotopy in topology to propose a general analytic method for nonlinear problems, namely, Homotopy Analysis Method 2 International Journal of Differential Equations HAM . This method has been successfully applied to solve many types of nonlinear problems in science and engineering, such as the viscous flows of non-Newtonian fluids 14 , the KdV-type equations 15 , higher-dimensional initial boundary value problems of variable coefficients 16 , and finance problems 17 . The HAM contains a certain auxiliary parameter h which provides us with a simple way to adjust and control the convergence region and rate of convergence of the series solution.
The HAM offers certain advantages over routine numerical methods. Numerical methods use discretization which gives rise to rounding off errors causing loss of accuracy and requires large computer memory and time. This computational method yields analytical solutions and has certain advantages over standard numerical methods. The HAM method is better since it does not involve discretization of the variables and hence is free from rounding off errors and does not require large computer memory or time.
The study of foam drainage equation is very significant for that the equation is a simple model of the flow of liquid through channels Plateau borders 18 and nodes intersection of four channels between the bubbles, driven by gravity and capillarity 19 . It has been studied by many authors 20-22 . The study for the foam drainage equation with time and space-fractional derivatives of this form has been investigated by the ADM and VIM method in 23, 24 . The fractional derivatives are considered in the Caputo sense. When α β 1, the fractional equation reduces to the foam drainage equation of the form In this paper, we extend the application of HAM to obtain analytic solutions to the spaceand time-fractional foam drainage equation. Two cases of special interest such as the timefractional foam drainage equation and the space-fractional foam drainage equation are discussed in details. Further, we give comparative remarks with the results obtained using ADM and VIM method see 23, 24 . The paper has been organized as follows. Notations and basic definitions are given in Section 2. In Section 3 the homotopy analysis method is described. In Section 4 we extend the method to solve the space-and time-fractional foam drainage equation. Discussion and conclusions are presented in Section 5.

Description on the Fractional Calculus
Definition 2.1. A real function f t , t > 0 is said to be in the space C μ , μ ∈ R if there exists a real number p > μ, such that f t t p f 1 t where f 1 ∈ 0, ∞ , and it is said to be in the space C μ n l if and only if h n ∈ C μ , n ∈ N. Clearly C μ ⊂ C ν if ν ≤ μ.

Basic Idea of HAM
To describe the basic ideas of the HAM, we consider the following differential equation: where N is nonlinear operator, D α t stand for the fractional derivative and is defined as in 2.3 , x, t denotes independent variables, and u x, t is an unknown function, respectively.
By means of generalizing the traditional homotopy method, Liao 9 constructs the so-called zero-order deformation equation where q ∈ 0, 1 is the embedding parameter, h / 0 is a nonzero auxiliary parameter, H t / 0 is an auxiliary function, L is an auxiliary linear operator, u 0 x, t is initial guesse of u x, t , and φ x, t, q is unknown function. It is important that one has great freedom to choose auxiliary things in HAM. Obviously, when q 0 and q 1, it holds that respectively. Thus, as q increases from 0 to 1, the solution φ x, t, q varies from the initial guess u 0 x, t to the solution u x, t . Expanding φ x, t, q in Taylor series with respect to q, we have If the auxiliary linear operator, the initial guess, the auxiliary parameter h, and the auxiliary function are so properly chosen, the series 3.4 converges at q 1, then we have which must be one of solutions of original nonlinear equation, as proved by Liao 11 . As h −1 and H t 1, 3.2 becomes International Journal of Differential Equations 5 which is used mostly in the homotopy perturbation method 29 , whereas the solution obtained directly, without using Taylor series. According to definition 3.5 , the governing equation can be deduced from the zero-order deformation equation 3.2 . Define the vector Differentiating 3.2 m times with respect to the embedding parameter q and then setting q 0 and finally dividing them by m!, we have the so-called mth-order deformation equation

3.10
Applying the Riemann-Liouville integral operator J α on both side of 3.9 , we have It should be emphasized that u m x, t for m 1 is governed by the linear equation 3.9 under the linear boundary conditions that come from original problem, which can be easily solved by symbolic computation software such as Matlab. For the convergence of the above method we refer the reader to Liao's work. Liao 10 proved that, as long as a series solution given by the homotopy analysis method converges, it must be one of exact solutions. So, it is important to ensure that the solution series is convergent. Note that the solution series contain the auxiliary parameter h, which we can choose properly by plotting the so-called h-curves to ensure solution series converge.
Remark 3.1. The parameters α and β can be arbitrarily chosen as, integer or fraction, bigger or smaller than 1. When the parameter is bigger than 1, we will need more initial and boundary conditions such as u 0 x, 0 , u 0 x, 0 , . . . and the calculations will become more complicated correspondingly. In order to illustrate the problem and make it convenient for the readers, we only confine the parameter to 0, 1 to discuss.

Application
In this section we apply this method for solving foam drainage equation with time-and space-fractional derivatives.
where c is the velocity of wavefront 15 . The exact solution of 4.1 for the special case α β 1 is

4.3
For application of homotopy analysis method, in view of 4.1 and the initial condition given in 4.2 , it in convenient to choose as the initial approximate. We choose the linear operator with the property L c 0, where c is constant of integration. Furthermore, we define a nonlinear operator as

4.6
We construct the zeroth-order and the mth-order deformation equations where

4.7
International Journal of Differential Equations 7 We now successively obtain . . .

4.8
By taking α 1, h −1, we reproduce the solution of problem as follows:     For application of homotopy analysis method, in view of 4.10 and the initial condition given in 4.2 , it is inconvenient to choose Initial condition has been taken as the above polynomial to avoid heavy calculation of fractional differentiation. We choose the linear operator L φ x, t, q ∂φ x, t, q ∂t , 4.13 with the property L c 0, where c is constant of integration. Furthermore, we define a nonlinear operator as

4.14
We construct the zeroth-order and the mth-order deformation equations where

4.16
International Journal of Differential Equations

Conclusion
In this paper, we have successfully developed HAM for solving space-and time-fractional foam drainage equation. HAM provides us with a convenient way to control the convergence of approximation series by adapting h, which is a fundamental qualitative difference in analysis between HAM and other methods. The obtained results demonstrate the reliability of the HAM and its wider applicability to fractional differential equation. It, therefore, provides more realistic series solutions that generally converge very rapidly in real physical problems.
Matlab has been used for computations in this paper.

Introduction
Fractional differential equations have gained importance and popularity during the past three decades or so, mainly due to its demonstrated applications in numerous seemingly diverse fields of science and engineering. For example, the nonlinear oscillation of earthquake can be modeled with fractional derivatives, and the fluid-dynamic traffic model with fractional derivatives can eliminate the deficiency arising from the assumption of continuum traffic flow. The differential equations with fractional order have recently proved to be valuable tools to the modeling of many physical phenomena 1, 2 . This is because of the fact that the realistic modeling of a physical phenomenon does not depend only on the instant time, but also on the history of the previous time which can also be successfully achieved by using fractional calculus. Most nonlinear fractional equations do not have exact analytic solutions, so approximation and numerical techniques must be used. The Adomain decomposition method 3 , the homotopy perturbation method 4 , the variational iteration 2 International Journal of Differential Equations method 5 , and other methods have been used to provide analytical approximation to linear and nonlinear problems. However, the convergence region of the corresponding results is rather small. In 1992, Liao employed the basic ideas of the homotopy in topology to propose a general analytic method for nonlinear problems, namely, homotopy analysis method 6-10 . This method has been successfully applied to solve many types of nonlinear problems in science and engineering, such as the viscous flows of non-Newtonian fluids 11 , the KdV-type equations 12 , finance problems 13 , fractional Lorenz system 14 , and delay differential equation 15 . The HAM contains a certain auxiliary parameter h which provides us with a simple way to adjust and control the convergence region and rate of convergence of the series solution.
The HAM offers certain advantages over routine numerical methods. Numerical methods use discretization which gives rise to rounding off errors causing loss of accuracy and requires large computer memory and time. This computational method yields analytical solutions and has certain advantages over standard numerical methods. The HAM method is better since it does not involve discretization of the variables and hence is free from rounding off errors and does not require large computer memory or time.
In this paper, we extend the application of HAM to discuss the explicit numerical solutions of a type of nonlinear-coupled equations with parameters derivative in this form: where L i and N i i 1, 2 are the linear and nonlinear functions of u and v, respectively, α and β are the parameters that describe the order of the derivative. Different nonlinear coupled systems can be obtained when one of the parameters α or β varies. The study of 1.1 is very necessary and significant because when α and β are integers, it contains many important mathematical physics equations.
The paper has been organized as follows. Notations and basic definitions are given in Section 2. In Section 3 the homotopy analysis method is described. In Section 4 applying HAM for two famous coupled examples: Burgers and mKdV equations. Discussion and conclusions are presented in Section 5.

Description on the Fractional Calculus
Definition 2.2. The Riemann-Liouville fractional integral operator J α of order α ≥ 0, of a function f ∈ C μ , μ ≥ −1, is defined as

2.3
Here, we also need two basic properties about them:

2.4
Definition 2.4. The Mittag-Leffler function E α z with a > 0 is defined by the following series representation, valid in the whole complex plane:

Basic Idea of HAM
To describe the basic ideas of the HAM, we consider the operator form of 1.1 : where N is nonlinear operator, D α t and D β t stand for the fractional derivative and are defined as in 2.3 , t denotes an independent operator, and u x, t , v x, t are unknown functions.

International Journal of Differential Equations
By means of generalizing the traditional homotopy method, Liao 6 constructs the so-called zero-order deformation equations: where q ∈ 0, 1 is the embedding parameter, h / 0 is a non-zero auxiliary parameter, H t / 0 is an auxiliary function, L is an auxiliary linear operator, u 0 x, t , v 0 x, t are initial guesses of u x, t , v x, t and φ 1 x, t, q , φ 2 x, t, q are two unknown functions, respectively. It is important that one has great freedom to choose auxiliary things in HAM. Obviously, when q 0 and q 1, the following holds: 3.4 respectively. Thus, as q increases from 0 to 1, the solution φ 1 x, t, q , φ 2 x, t, q varies from the initial guess u 0 x, t , v 0 x, t to the solution u x, t , v x, t . Expanding φ 1 x, t, q , φ 2 x, t, q in Taylor series with respect to q, we have

3.6
If the auxiliary linear operator, the initial guess, the auxiliary parameter h, and the auxiliary function are so properly chosen, the series 3.5 converges at q 1, then we have

3.7
International Journal of Differential Equations 5 which must be one of solutions of original nonlinear equation, as proved by Liao 8 . As h −1 and H t 1, 3.2 and 3.3 become which is used mostly in the homotopy perturbation method 20 , where as the solution obtained directly, without using Taylor series. According to the definition 3.6 , the governing equation can be deduced from the zero-order deformation equation 3.2 . Define the vector Differentiating equations 3.2 and 3.3 m times with respect to the embedding parameter q and then setting q 0 and finally dividing them by m!, we have the so-called mth-order deformation equation:

3.11
Applying the Riemann-Liouville integral operator J α , J β on both side of 3.10 , we have

3.12
It should be emphasized that u m x, t , v m x, t for m 1 is governed by the linear equation 3.10 , under the linear boundary conditions that come from original problem, which can be easily solved by symbolic computation software such as MATLAB. For the convergence 6 International Journal of Differential Equations of the above method we refer the reader to Liao's work. Liao 7 proved that, as long as a series solution given by the homotopy analysis method converges, it must be one of exact solutions. So, it is important to ensure that the solution series is convergent. Note that the solution series contain the auxiliary parameter h, which we can choose properly by plotting the so-called h-curves to ensure solution series converge.
Remark 3.1. The parameters α and β can be arbitrarily chosen as, integer or fraction, bigger or smaller than 1. When the parameters are bigger than 1, we will need more initial and boundary conditions such as u 0 x, 0 , u 0 x, 0 , . . . and the calculations will become more complicated correspondingly. In order to illustrate the problem and make it convenient for the readers, we only confine the parameters to 0, 1 to discuss.

The Nonlinear Coupled Burgers Equations with Parameters Derivative
In order to illustrate the method discussed above, we consider the nonlinear coupled Burgers equations with parameters derivative in an operator form: where t > 0, L x ∂/∂x and the fractional operators D α t and D β t are defined as in 2.3 . Assuming the initial value as The exact solutions of 4.1 for the special case: α β 1 are For application of homotopy analysis method, in view of 4.1 and the initial condition given in 4.2 , it is convenient to choose as the initial approximate of 4.1 . We choose the linear operators International Journal of Differential Equations 7 with the property L c 0 where c is constant of integration. Furthermore, we define a system of nonlinear operators as

4.6
We construct the zeroth-order and the mth-order deformation equations where

4.7
We start with an initial approximation u x, 0 sin x , v x, 0 sin x , thus we can obtain directly the other components as The absolute error of the 6th-order HAM and exact solution with h −1 as shown in Figure 1. Also the absolute errors |u t − φ 6 t | have been calculated in Table 1. Figure 2 shows the numerical solutions of the nonlinear coupled Burgers equations with parameters derivative with h −1, α β 1. Figure 3 shows the explicit numerical solutions with h −1, α 1/4, and β 1/3 at t 0.02. As suggested by Liao 7 , the appropriate region for h is a horizontal line segment. We can investigate the influence of h on the convergence of the solution series gevin by the HAM, by plotting its curve versus h, as shown in Figure 4.

The Nonlinear Coupled mKdV Equations with Parameters Derivative
In order to illustrate the method discussed above, we consider the nonlinear coupled mKdV equations with parameters derivative in an operator form: with the initial conditions,  As we know, when α β 1 4.9 has the kink-type soliton solutions

4.11
constructed by Fan 22 , where ξ x 1/4 −4k 2 − 6λ 6kλ/b 3b 2 /k 2 t, k / 0, and b / 0. For application of homotopy analysis method, in view of 4.9 and the initial condition given in 4.10 , it in convenient to choose

4.13
with the property L c 0 where c is constant of integration. Furthermore, we define a system of nonlinear operators as

4.14
We construct the zeroth-order and the mth-order deformation equations where

4.15
International Journal of Differential Equations 11 We start with an initial approximation u x, 1, k 1/3, thus we can obtain directly the other components as follows: . . . The absolute error of the 6th-order HAM and exact solution with h −1 as shown in Figure 5. Also the absolute errors |u t − φ 6 t | have been calculated for in Table 2. Figure 6 shows the numerical solutions of the nonlinear coupled Burgers equations with parameters derivative with h −1, α β 1. Figure 7 shows the explicit numerical solutions with h −1, α 1/2, and β 2/3 at t 0.002.
As suggested by Liao 7 , the appropriate region for h is a horizontal line segment. We can investigate the influence of h on the convergence of the solution series gevin by the HAM, by plotting its curve versus h, as shown in Figure 8.

Conclusion
In this paper, based on the symbolic computation MATLAB, the HAM is directly extended to derive explicit and numerical solutions of the nonlinear coupled equations with parameters International Journal of Differential Equations derivative. HAM provides us with a convenient way to control the convergence of approximation series by adapting h, which is a fundamental qualitative difference in analysis between HAM and other methods. The obtained results demonstrate the reliability of the HAM and its wider applicability to fractional differential equation. It, therefore, provides more realistic series solutions that generally converge very rapidly in real physical problems. MATLAB has been used for computations in this paper.

Introduction
Fractional differential equations have excited, in recent years, a considerable interest both in mathematics and in applications. They were used in modeling of many physical and chemical processes and engineering see, e.g., 1-6 . In its turn, mathematical aspects of fractional differential equations and methods of their solutions were discussed by many authors: the iteration method in 7 , the series method in 8 , the Fourier transform technique in 9, 10 , special methods for fractional differential equations of rational order or for equations of special type in 11-16 , the Laplace transform technique in 3-6, 16, 17 , and the operational calculus method in 18-23 . Recently, several mathematical methods including the Adomian decomposition method 18-25 , variational iteration method 23-26 and homotopy perturbation method 27, 28 have been developed to obtain the exact and approximate analytic solutions. Some of these methods use transformation in order to reduce equations into simpler equations or systems of equations, and some other methods give the solution in a series form which converges to the exact solution.
The reason of using fractional order differential FOD equations is that FOD equations are naturally related to systems with memory which exists in most biological systems. Also they are closely related to fractals which are abundant in biological systems. The results 2 International Journal of Differential Equations derived from the fractional system are of a more general nature. Respectively, solutions to the fractional diffusion equation spread at a faster rate than the classical diffusion equation and may exhibit asymmetry. However, the fundamental solutions of these equations still exhibit useful scaling properties that make them attractive for applications.
The concept of fractional or noninteger order derivation and integration can be traced back to the genesis of integer order calculus itself 29 . Almost all of the mathematical theory applicable to the study of noninteger order calculus was developed through the end of the 19th century. However, it is in the past hundred years that the most intriguing leaps in engineering and scientific application have been found. The calculation techniques in some cases meet the requirement of physical reality. The use of fractional differentiation for the mathematical modeling of real-world physical problems has been widespread in recent years, for example, the modeling of earthquake, the fluid dynamic traffic model with fractional derivatives, and measurement of viscoelastic material properties. Applications of fractional derivatives in other fields and related mathematical tools and techniques could be found in 30-41 . In fact, real-world processes generally or most likely are fractional order systems.
The derivatives are understood in the Caputo sense. The general response expression contains a parameter describing the order of the fractional derivative that can be varied to obtain various responses.

Fractional Calculus
There are several approaches to the generalization of the notion of differentiation to fractional orders, for example, the Riemann-Liouville, Grünwald-Letnikov, Caputo, and generalized functions approach 42 . The Riemann-Liouville fractional derivative is mostly used by mathematicians but this approach is not suitable for real-world physical problems since it requires the definition of fractional order initial conditions, which have no physically meaningful explanation yet. Caputo introduced an alternative definition, which has the advantage of defining integer order initial conditions for fractional order differential equations 42 . Unlike the Riemann-Liouville approach, which derives its definition from repeated integration, the Grünwald-Letnikov formulation approaches the problem from the derivative side. This approach is mostly used in numerical algorithms.
Here, we mention the basic definitions of the Caputo fractional-order integration and differentiation, which are used in the upcoming paper and play the most important role in the theory of differential and integral equation of fractional order.
The main advantages of Caputo approach are the initial conditions for fractional differential equations with the Caputo derivatives taking on the same form as for integer order differential equations.
International Journal of Differential Equations 3 For the Caputo derivative we have D α C 0, C is constant, For m to be the smallest integer that exceeds α, the Caputo fractional derivative of order α > 0 is defined as

Analysis of the Method
The Mittag-Leffler 1902-1905 functions E α and E α,β 42 , defined by the power series have already proved their efficiency as solutions of fractional order differential and integral equations and thus have become important elements of the fractional calculus theory and applications.
In this paper, we will explain how to solve some of differential equations with fractional level through the imposition of the generalized Mittag-Leffler function E α z . The generalized Mittag-Leffler method suggests that the linear term y x is decomposed by an infinite series of components: We will use the following definitions of fractional calculus: This is based on the Caputo fractional is derivatives. The convergence of the Mittag Leffler function discussed in 42 .

Applications and Results
In this section, we consider a few examples that demonstrate the performance and efficiency of the generalized Mittag-Leffler function method for solving linear differential equations with fractional derivatives.
Example 4.1. Consider the following fractional differential equation 43 : Combining the alike terms and replacing n by n 1 in the first sum, we assume the form With the coefficient of x nα equal to zero and identifying the coefficients, we obtain recursive a n 1 − Aa n 0 ⇒ a n 1 Aa n ,

4.5
The general solution is International Journal of Differential Equations 5 We can write the general solution in the Mittag-Leffler function form as As α 1, we have the exact solution: which is the exact solution of the standard form.
Example 4.2. Consider the fractional differential equation 44 By using 3.2 and 3.4 into 4.9 we find ∞ n 2 a n Combining the alike terms and replacing n by n 2 in the first sum, we assume the form ∞ n 0 a n 2 x nα Γ nα 1 − ∞ n 0 a n x nα Γ nα 1 0, ∞ n 0 a n 2 − a n x nα Γ nα 1 0.

4.11
With the Coefficient of x nα equal to zero and identifying the coefficients, we obtain recursive a n 2 a n .

4.12
Substituting into 3.2 , we find that: If a 1, we can write the general solution in the Mittag-Leffler function form as 14 which is the exact solution of the linear fractional differential equation 4.9 .
Combining the alike terms and replacing n by n 2 in the first sum, we assume the form With the coefficient of x nα equal to zero and identifying the coefficients, we obtain recursive a n 2 2a n − a n 1 .

4.18
Substituting into 3.2 , we find that: If a 1, we can write the general solution in the Mittag-Leffler function form as which is the solution of the linear fractional differential equation 4.15 .

Conclusions
A new generalization of the Mittag-Leffler function method has been developed for linear differential equations with fractional derivatives. The new generalization is based on the Caputo fractional derivative. It may be concluded that this technique is very powerful and efficient in finding the analytical solutions for a large class of linear differential equations of fractional order.

Introduction
Fractional calculus is a key tool for solving some relevant scientific problems in physics, engineering, biology, chemistry, hydrology, and so on 1-6 . A field of research in which the fractional formalism has been particularly useful is that related to anomalous diffusion processes 1, 7-13 . This kind of process is singularly abundant and important in biological media 14-16 . In this context, the electrodiffusion of ions in neurons is an anomalous diffusion problem to which the fractional calculus has recently been applied. The precise origin of the anomalous character of this diffusion process is not clear see 17 and references therein , but in any case the consideration of anomalous diffusion in the modeling of electrodiffusion of ions in neurons seems pertinent. This problem has been addressed recently by Langlands et al. 17,18 . An equation that plays a key role in their analysis is the following fractional cable or telegrapher's or Cattaneo equation model II : ∂u ∂t with n − 1 < γ < n and n integer, is the Riemann-Liouville fractional derivative. Here u is the difference between the membrane potential and the resting membrane potential, γ 1 is the exponent characterizing the anomalous flux of ions along the nerve cell, and γ 2 is the exponent characterizing the anomalous flux across the membrane 17, 18 . Some earlier fractional cable equations were discussed in 19, 20 .
A variety of analytical and numerical methods to solve many classes of fractional equations have been proposed and studied over the last few years 10, 21-30 . Of the numerical methods, finite difference methods have been particularly fruitful 31-38 . These methods can be broadly classified as explicit or implicit 39 . An implicit method for dealing with 1.1 has recently been considered by Liu et al. 38 . Although implicit methods are more cumbersome than explicit methods, they usually remain stable over a larger range of parameters, especially for large timesteps, which makes them particularly suitable for fractional diffusion problems. Nevertheless, explicit methods have some features that make them widely appreciated 32, 39 : flexibility, simplicity, small computational demand, and easy generalization to spatial dimensions higher than one. Unfortunately, they can become unstable in some cases, so that it is of great importance to determine the conditions under which these methods are stable. In this paper we will discuss an explicit finite difference scheme for solving the fractional cable equation, which is close to the methods studied in 32, 33 . We shall address two main questions: i whether this kind of method can cope with fractional equations involving different fractional derivatives, such as the fractional cable equation; ii whether the von Neumann stability analysis put forward in 32, 34 is suitable for this kind of equation.

The Numerical Method
Henceforth, we will use the notation x j jΔx, t m mΔt, and u x j , t m u m j U m j , where U m j is the numerical estimate of the exact solution u x, t at x x j and t t m .
In order to get the numerical difference algorithm, we discretize the continuous differential and integro-differential operators as follows. For the discretization of the fractional Riemann-Liouville derivative we use the Grünwald-Letnikov formula International Journal of Differential Equations 3 and ω α 0 1. These coefficients come from the generating function 40 To discretize the integer derivatives we use standard formulas: for the second-order spatial derivative we employ the three-point centered formula and for the first-order time derivative we use the forward derivative Inserting 2.1 , 2.5 , and 2.7 into 1.1 , one gets where, as can easily be proved, the truncating error T x, t is Neglecting the truncating error we get the finite difference scheme we are seeking: To test this algorithm, we solved 1.1 in the interval −L/2 ≤ x ≤ L/2, with absorbing boundary conditions, u x −L/2, t u x L/2, t 0, and initial condition given by a Dirac's delta function centered at x 0: u x, 0 δ x . The exact solution of this problem for L → ∞ is 17 where H denotes the Fox H function 10, 41 . In our numerical procedure, the exact initial condition u x, 0 δ x is approximated by

2.15
The explicit difference scheme 2.12 is tested by comparing the analytical solution with the numerical solution for several cases of the problem described following 2.13 with different values of γ 1 and γ 2 . We have computed the analytical solution by means of 2.14 truncating the series at k 20. The corresponding Fox H function was evaluated by means of the series expansion described in 10, 41 truncating the infinite series after the first 50 terms. In Figures  1 and 2 we show the analytical and numerical solution for two values of γ 1 and γ 2 at x 0 and x 0.5. The differences between the exact and the numerical solution are shown in Figures 3 and 4. One sees that, except for very short times, the agreement is quite good. The large value of the error for small times is due in part to the approximation of the Dirac's delta function at x 0 by 2.15 . This is clearly appreciated when noticing the quite different scales of Figures 3 and 4: the error is much smaller for x 0.5 than for x 0. For the cases with γ 1 1/2 we used a smaller value of Δt and, simultaneously, a larger value of Δx than for the cases with γ 1 1 in order to keep the numerical scheme stable. This issue will be discussed in Section 3.

Stability
As usual for explicit methods, the present explicit difference scheme 2.12 is not unconditionally stable, that is, for any given set of values of γ 1 , γ 2 , μ, and K there are choices of Δx and Δt for which the method is unstable. Therefore, it is important to determine the conditions under which the method is stable. To this end, here we shall employ the fractional von Neumann stability analysis or Fourier analysis put forward in 32 see also 33-35 . The question we address is to what extent this procedure is valid for fractional diffusion equations that involve fractional derivatives of different order.
Proceeding as 32 , we start by recognizing that the solution of our problem can be written as the linear combination of subdiffusive modes, u sum is over all the wave numbers q supported by the lattice. Therefore, following the von Neumann ideas, we reduce the problem of analyzing the stability of the solution to the problem of analyzing the stability of a single generic subdiffusion mode, ζ m e iqjΔx . Inserting this expression into 2.12 one gets The stability of the mode is determined by the behavior of ζ m . Writing  and assuming that the amplification factor ξ of the subdiffusive mode is independent of time, we get If |ξ| > 1 for some q, the temporal factor of the solution grows to infinity c.f., 3.2 , and the mode is unstable. Considering the extreme value ξ −1, one gets from 3.3 that the numerical method is stable if this inequality holds: Therefore, because S ≤ S, we find that a sufficient condition for the present method to be stable is that S ≤ S × . In Figures 5 and 6 we show two representative examples of the problem considered in Figure 2 but for two values of S, respectively, larger and smaller than the stability bound provided by 3.7 . One sees that the value of S that one chooses is crucial: when S is smaller than S × one is inside the stable region and gets a sensible numerical solution Figure 5 ; otherwise one gets an evidently unstable and nonsensical solution Figure 6 .

Numerical Check of the Stability Analysis
In this section we describe a comprehensive check of the validity of our stability analysis by using many different values of the parameters γ 1 , γ 2 , Δt, and Δx and testing whether the stability of the numerical method is as predicted by 3.7 . Without loss of generality, we assume μ K 1 in all cases. We proceed in the following way. First, we choose a set of values of γ 1 , γ 2 , Δx, and S and integrate the corresponding fractional cable equation. If for λ 10 within the first 1000 integrations, then we say the method is unstable; otherwise, we label the method as stable. We generated Figure 7 by starting the integration for values of S well below the theoretical stability limit given by 3.7 and kept increasing its value by 0.001 until condition 4.1 was first reached. The last value for which the method was stable is recorded and plotted in Figure 7. The limit value λ 10 is arbitrary, but choosing any other reasonable value does not significantly change these plots.

Convergence Analysis
In this section we show that the present numerical method is convergent, that is, that the numerical solution converges towards the exact solution when the size of the spatiotemporal   Therefore, from 5.2 one gets ζ 1 χ 0 . But from 2.10 one knows that |T where |ζ {m} | is the maximum value of |ζ k | for k 0, . . . , m. Taking into account 2.4 , using the value z 1, and because ω α 0 1, it is easy to see that ∞ k 1 ω α k −1 or, equivalently, International Journal of Differential Equations fact, it is smaller than 2 . Using this result in 5.3 , together with |ζ k | ≤ C Δt Δx 2 and |χ k | ≤ C Δt Δx 2 , we find that Therefore the amplitude of the subdiffusive modes goes to zero when the spatiotemporal mesh goes to zero. Employing the Parseval relation, this means that the norm of the error goes to zero when Δt and Δx go to zero. This is what we aimed to prove.

Conclusions
An explicit method for solving a kind of fractional diffusion equation that involves several fractional Riemann-Liouville derivatives, which are approximated by means of the Grünwald-Letnikov formula, has been considered. The method was used to solve a class of equations of this type fractional cable equations with free boundary conditions, Dirac's delta initial condition, and different fractional exponents. The error of the numerical method is compatible with the truncating error, which is of order O Δt O Δx 2 . It was also proved that the method is convergent. Besides, it was also found that a fractional von-Neumann stability analysis, which provides very precise stability conditions for standard fractional diffusion equations, leads also to a very accurate estimate of the stability conditions for cable equations.

Introduction
Let X 1 0 , X 1 1 , . . . , X m 1 , X 2 0 , X 2 1 , . . . , X m 2 be 2m 2 vector fields on R d with bounded derivatives at each order. Let be an Hoermander's type operator on R 1 d . Let be a second Hoermander's operator on R 1 d . Bismut 1 considers the generator and the Markov semi-group exp tA . This semi-group has a probabilistic representation. We consider a Brownian motion t → z t independent of the others Brownian motions B i t . Bismut introduced the solution of the stochastic differential equation starting at x in Stratonovitch sense: 1.4 where t → B i t are m independent Brownian motions. Let us introduce the local time t → L t associated with t → z t and its right inverse t → A t see 2, 3 . Then, Such operator is classically related to the Dirichlet Problem 3 . Classically 4 , where x 1 t x is the solution of the Stratonovitch differential equation starting at x: The question is as following: is there an heat-kernel associated with the semi-group exp tL 1 ? This means that There are several approaches in analysis to solve this problem, either by using tools of microlocal analysis or tools of harmonic analysis. Malliavin 5 uses the probabilistic representation of the semi-group. Malliavin uses a heavy apparatus of functional analysis number operator on Fock space or equivalently Ornstein-Uhlenbeck operator on the Wiener space, Sobolev spaces on the Wiener space in order to solve this problem. Bismut 6 avoids using this machinery to solve this hypoellipticity problem. In particular, Bismut's approach can be adapted immediately to the case of the Poisson process 7 . The main difficulty to treat in the case of a Poisson process is the following: in general the solution of a stochastic differential equation with jumps is not a diffeomorphism when the starting point is moving see [8][9][10] .
The main remark of Bismut in 1 is that if we consider the jump process t → x 1 A t x , then it is a diffeomorphism almost surely in x. So, Bismut mixed the tools of the Malliavin International Journal of Differential Equations 3 Calculus for diffusion on the process t → x 1 t x and the tools of the Malliavin Calculus for Poisson process on the jump process t → A t in order to show that this is the problem if Developments on Bismut's idea was performed by Léandre in 9, 11 . Let us remark that this problem is related to study the regularity of the Dirichlet problem see 1, page 598 see 12-14 for related works .
Recently, we have translated into the language of semi-group theory the Malliavin Calculus of Bismut type for diffusion 15 . We have translated in semi-group theory a lot of tools on Poisson processes [16][17][18][19][20][21][22] . Especially, we have translated the Malliavin Calculus of Bismut type for Poisson process in semi-group theory in 17 . It should be tempting to translate in semi-group theory Bismut's Calculus on boundary process. It is the object of this work.
On the general problematic on this work, we refer to the review papers of Léandre 23-25 . It enters in the general program to introduce stochastic analysis tools in the theory of partial differential equation see 26-28 .

Statements of the Theorems
Let us recall some basis on the study of fractional powers of operators 29 . Let L be a generator of a Markovian semi-group P s . Then, The results of this paper could be extended to generators of the type where ∞ 0 g t ∧ 1dt < ∞ and g ≥ 0, but we have chosen the operator of the type 1.3 to be more closely related to the original intuition on Bismut's Calculus on boundary process. Let be On E d , we consider the vector fields:

4 International Journal of Differential Equations
We consider the Malliavin generator L 1 on E d : We consider the Malliavin semi-group P 1 t associated and − L 1 . We perform the same algebraic considerations on L 2 . We get L 2 , P 2 t , and − L 2 . Let us consider the total generator We get a theorem which enters in the framework of the Malliavin Calculus for heatkernel.

Theorem 2.1. Let one suppose that the Malliavin condition in x is checked:
holds for all p, then where q t s, y is the density of a probability measure on R 1 d .
is invertible in x, then the Malliavin condition holds in x.

Remark 2.3.
We give simple statements to simplify the exposition. It should be possible by the method of this paper to translate the results of 9, part III , got by using stochastic analysis as a tool.

Integration by Parts on the Underlying Diffusion
We consider the vector fields on R 1 d 1 , x is a convenient matrix on R m which depends smoothly on x and whose derivatives at each order are bounded. s, t → h j s,t does not depend on x, and h j s,t belong to R m . Let f be a smooth function on R 1 d 1 , D f denotes its gradient, and D 2 f denotes its Hessian.
We consider the generator L j,1 s,t acting on smooth functions on R 1 d 1 , We put It generates a semi-group P 1 t . Let us consider the Hoermander's type generator associated with the smooth Lipschitz vector fields on R 1 d d s, x, U on R 1 d d :

3.5
We consider the heat semi-group associated with L Let us consider the semi-group P 2 t associated with We get, with the same notations for s, x, u, U the following.

Theorem 3.2.
For f bounded continuous with compact support in s, x , one has the following relation: Proof. For the integrability conditions, we refer to the appendix. We remark that ∂/∂u commute with A 1 t , therefore with P 1 t . We deduce that By the method of variation of constants, In order to show that, we follow the lines of 2.17 and 2.18 in 15 . We apply A 1 t to 3.11 . By Lemma 3.1, Let us consider the vector fields on We consider the Hoermander's type operator associated with these vector fields:

3.15
International Journal of Differential Equations 7 We consider the generator It generates a semi-group P 3 t . By lemma 3.2 of 15 , we have

3.18
In 15, Equation 3.18 , we consider the semi-group P t instead of the semi-group P j,2 s,t and the test function Df instead as of the test function P 3 t DfU ·, I here. Y j i,v,t is considered as an element of R and not as a one-order differential operator: Therefore,

3.20
We write The Volterra expansion see 15, Equation 3.17 if it converges gives the following formula: International Journal of Differential Equations s,t fU s 0 , x 0 , U 0 is linear in u 0 . Therefore:

3.24
In this last formula, Y j i,s,t are considered as differential operators. Therefore, A 3 t fU s 0 , x 0 , U 0 does not depend on U 0 and is equal to But U 0 → P 3 s fU s 0 , x 0 , U 0 is linear. Therefore, s−v,t P 3 t fU s 0 , x 0 , I ds dv.

3.27
It remains to replace f by Df in this last equation and to compare 3.26 with 3.13 and 3.20 .
We consider the Malliavin generator A. We can perform the same algebraic construction as in Theorem 3.2. We get two semi-groups P 2 t and P 1 t . Y j i,s,t and Z j i,s,t are smooth with bounded derivatives in x x, U, U −1 , V . We get by the same procedure the following.

Theorem 3.3. If f is bounded with bounded derivatives and with compact support in s, then one gets
where one take does not derivative in the direction of s in D f.
We can perform the same improvements as in 15, page 512 . We define on R d × R d 1 × · · · × R d k some vectors fields: We can consider the generator A tot associated with these vector fields and perform the same algebraic computations as in Theorem 3.2. We get two semi-groups P 2,tot where D f tot does not include derivative in the direction of s.
We refer to the appendix for the proof and the subsequent estimates. According to the fact that f has compact support this means that we consider bounded values of A t , the transformed Brownian motion has an equivalent law through the Girsanov exponential to the original Brownian motions. The term in u in Theorem 3.2 come that from the fact we take the derivative in λ 0 of the Girsanov exponential. When we do this transformation, we get a random process x λ t x . Derivation of it in λ 0 is done classically according to the stochastic flow theorem, which leads to the study of generators of the type L j,2 s,t and of the type L j,3 .

Integration by Parts on the Subordinator
Let us consider diffusion type generator of the previous part: We have classically where the smooth vector fields are Lipschitz. Therefore, we can write We consider a diffeomorphsim f λ s of 0, ∞ with bounded derivative of first order in λ equal to s if s < and equals to s if s > 2 we suppose λ small . We can write We do this operation on the two operators on R 1 d giving A. We get a generator A λ . According the line of stochastic analysis, we consider a generator A λ,1 on R 1 d 1 . If L 1 is a generator on R 1 d with associated semi-group P s , then we consider A λ,1 the generator on R 1 d 1 , where J λ s is the Jacobian of the transformation s → f λ s . By doing this procedure in 1.3 , we deduce a global generator A λ,1 and a semi-group P λ,1 t associated with it. It is not clear that P λ,1 t is a Markovian semi-group. We decompose International Journal of Differential Equations

4.10
But by an elementary change of variable, The result holds by the unicity of the solution of the parabolic equation associated with A. To state the integrability of u, we refer to 16 .
where g λ are smooth and where the vector fields Y i are smooth Lipschitz on R d . We consider the semi-group P λ,· t associated with it. Let us introduce the vector fields on R d d : Proof. Let us introduce the vector fields on R d d : and the generator

4.21
But U 0 → P λ,2,· t df, U x, U 0 is linear in U 0 and therefore the quantity L λ,1 − L λ,2 P λ,2,· t−s n df, U x, U 0 does not depend on U 0 . Therefore the Volterra expansion reads Let us compute L λ,1 − L λ,2 . It is equal to Therefore,

4.26
We insert this formula in the right-hand side of 4. 23

4.30
It generates a semi-group P λ,2 t . In order to see that, we split the generator by keeping the values od s or s and we proceed as for A λ,1 see 4.10 .

International Journal of Differential Equations
We get the following.
Proof. We have

4.32
But by 15, Lemma 3.2. : Therefore DP λ,1 t satisfies the parabolic equation associated with P λ,2 t df, v, U s 0 , x 0 , u 0 , 1, I . Only the integrability of U puts any problem. It is solved by the appendix since f has compact support in s.

4.35
But U 0 → P λ,2,· t df, U x, U 0 is linear in U 0 and therefore the quantity L λ,1 − L λ,2 P λ,2,· t−s n df, U x, U 0 do not depend of U 0 . Therefore the Volterra expansion reads Namely, the generator of P j, √ f λ s t does not act on the u 0 component such that the two sides of 4.37 satisfy the same parabolic equality. Therefore, where J 0 s ∂/∂λ J 0 s , Therefore, Let A 3 be the generator on R 1 d 1 d : It generates a semi-group, P 3 t . We get the following.

4.42
Proof. If the Volterra expansion converges, then

4.43
But P 2 t−s n df, u, U s 0 , x 0 , u 0 , U 0 is linear in u 0 , U 0 . Let us explain the details of that. We have to compute By the technique of the beginning of the proof of Proposition 4.3, it does not depend on u 0 , U 0 . Therefore, the Volterra expansion reads: does not depend on u 0 , U 0 . Therefore, the right-hand side of formula 4.45 is equal to International Journal of Differential Equations 17 because u 0 , U 0 → P j t−v df, u, U s 0 , x 0 , u 0 , U 0 is linear in u 0 , U 0 and because the vector fields which give the generator of P j, √ s,2,· s are linear in u 0 , U 0 . Therefore,

4.49
But by an analog of Theorem 4.5, We can summarize the previous considerations in the next theorem.
By the appendix, if h is compact support in s. Let us consider the generator A 3, associated with f λ . If g df, u, U , then we have by Duhamel principle By the proof of Theorem 4.8, P 3 1−s g s 0 , x 0 , u 0 , U 0 is affine in u 0 , U 0 . Namely, in the proof of this theorem, we have removed the P 2 1−s g s 0 , x 0 , u 0 , U 0 which is equal to zero in u 0 0, U 0 0 because this expression is linear in u 0 , U 0 . Its component in u 0 , U 0 is smooth with bounded derivatives at each order. By Theorem 4.6, A 3, − A P 3 1−s g s 0 , x 0 , u 0 , U 0 is still affine in u 0 , U 0 and its components in u 0 , U 0 are smooth with bounded derivatives at each order. Moreover, if g 1 is affine in u 0 , U 0 with components in u 0 , U 0 smooth with bounded derivatives at each order, then we get that, for s ≤ 1, where C → 0 when → 0. This can be seen as an appliation of the Duhamel formula applied to the two semi-groups t → P j, √ s,1,· t g 1 and t → P ,j, √ s,1,· t g 1 . Then, the result arises from the Duhamel formula 4.55 .
We can consider vector fields at the manner of 3.30 and f λ s s λs 5 in a neighborhood of 0. We get a generator A tot and semi-groups P j, √ s,tot s and P 3,tot t . We have with the extension of Theorem 4.10 the following.

Theorem 4.11 Bismut . If f λ s s λs 5 in a neighborhood of 0 and is equal to 1 if s > 1, then one has the following integration by parts: let f tot be a function with compact support in s, bounded with bounded derivatives at each order. Then,
4.57

The Abstract Theorem
The proof of Theorem 2.1 follows the idea of Malliavin 5 . If there exist C l such that, for function f with compact support in 0, 1 × 0, l d , then the heat kernel q t s, y exists. There are two partial derivatives to treat: i the partial derivative in the time of the subordinator s, ii the partial derivatives in the space of the underlying diffusion x.
Let us begin by the most original part of Bismut's Calculus on boundary process, that is, the integration by parts in the time s.
We look at 4.42 . We remark see the next part that for all p. So, we take f tot s, x, u f s, x 1/u and we apply 4.42 for this convenient semigroup. We get R can be estimated by using the appendix by C l f ∞ for f with compact support in 0, l × 0, l d and by 5.2 .

20
International Journal of Differential Equations Proof. If f is a function with compact support depending only of s, x tot and V , we have We do the change of variable U → U and V → UV on the Malliavin generator A tot . By using Lemma 3.7 of 15 , it is transformed in A 2,tot f s, x tot , U, V where for A 2,tot we consider the same type of operator as A 2 but with the modified vector fields:

5.6
It remains to use the appendix to show the Lemma.
We consider Z j i t U −1 X j i . By the previous Lemma and Malliavin hypothesis, for all p if g s has compact support V is a matrix . After we consider a test function of the type of Bismut, we consider the component u i of U in 5.3 . We consider the Bismut function fV −1 u i /u . We integrate by parts as in Theorem 3.4. We deduce under Malliavin assumption that if f has compact support in 0, l × 0, l d . By the same way, we deduce that if f has compact support in 0, l × 0, l d then Therefore, the result is obtained .
Remark 5.2. We could do integration by parts to each order in order to show that the semigroup exp tA has a smooth heat-kernel under Malliavin assumption.
International Journal of Differential Equations 21

Inversion of the Malliavin Matrix
Proof of Theorem 2.2. Let s 1 < s 2 and let ξ be of modulus 1. Then, These two quantities are equal in t 0 when we consider the semi-group exp t A . Let us compute their derivative in time t. The derivative of the left-hand side is bigger than the derivative of the right-hand side because Proof. We consider a convex function decreasing from 0, ∞ into 0, 1 equal to 1 in 0 and tending to 0 at infinity. Let us introduce α s exp s A g |V ξ| I 0, 0, x, I, 0 . 6.5 In order to consider the derivative in s of α s , we study the expression We have only to consider by 6.3 the case where s is small enough, |x − x| is small enough, |U − I| is small enough, and the positive matrix V is small enough. For that we have to estimate γ u j P j,2 u g |V ξ| s , x , U , V − g |V ξ| 6.7

22
International Journal of Differential Equations for u between 0 and . The first derivative of γ u has an equivalent −C −1 when → 0, and its second derivative has a bound C −2 when → 0. Therefore, 0 ≥ γ u ≥ − Cu 6.8 on 0, and We deduce from that that Remark 6.2. We could improve 6.4 by showing that if s is small enough, |x − x| is small enough, |U − I| is small enough, and the positive matrix V is small enough. We consider a very small α. We slice the time interval 0, α in α−1 intervals of length . We have exp t A I 0,l I V ξ ≤ 0, x, I, 0 ≤ exp α A I 0,l I V ξ ≤ 0, x, I, 0 for a small C 0 . This last quantity is smaller than C p for all p by the previous lemma if α is small enough. The proof of Theorem 2. Proof. We remark that if we consider only functions of u, then where P 4 t is a Lévy semi-group with generator where g s 1 on a neighborhood of 0, is with compact support and is positive. The result follows from the adaptation in 17, 18 of the proof of 7 in semi-group theory. We remark that By using the adaptation in semi-group theory of the exponential martingales of Levy process of 17, 18 , we have The result holds from the Tauberian theorem of 7, 17, 18 .

Burkholder-Davies-Gundy Inequality
Theorem A.4. Let s 0 > 0 and p ∈ N. Then, Proof. Following the idea of 17, Appendix , we consider the auxiliary function International Journal of Differential Equations We get Let us consider an improvement of the Gronwall lemma: We remark that for K independent of C. Then, by the modified Gronwall lemma, where K does not depend on C. By Gronwall lemma, where K does not depend on C. The result arises by doing C → ∞.
By the same procedure, we get the following.
Theorem A.5. Let be s 0 > 0 and p ∈ N. Then and we get the following.
Theorem A.6. Let s 0 > 0 and p ∈ N: Remark A.7. We can show 6.3 by the same way.

Introduction
Fractional calculus is more than 300 years old, but it did not attract enough interest at the early stage of development. In the last three decades, fractional calculus has become popular among scientists in order to model various physical phenomena with anomalous decay, such as dielectric polarization, electrode-electrolyte polarization, electromagnetic waves, and viscoelastic systems 1 . Recent advances in fractional calculus have been reported in 2 .
Recently, stability of fractional differential systems has attracted increasing interest. In 1996, Matignon 3 firstly studied the stability of linear fractional differential systems with the Caputo derivative. Since then, many researchers have done further studies on the stability of linear fractional differential systems 4-11 . For the nonlinear fractional differential systems, the stability analysis is much more difficult and only a few are available.
Some authors 12, 13 studied the following nonlinear fractional differential system: , where m − 1 < q ≤ m. They discussed the continuous dependence of solution on initial conditions and the corresponding structural stability by applying Gronwall's inequality. In 14 the authors dealt with the following fractional differential system: where 0 < q ≤ 1, D q 0,t denotes either the Caputo, or the Riemann-Liouville fractional derivative operator. They proposed fractional Lyapunov's second method and firstly extended the exponential stability of integer order differential systems to the Mittag-Leffler stability of fractional differential systems. Moreover, the pioneering work on the generalized Mittag-Leffler stability and the generalized fractional Lyapunov direct method was proposed in 15 .
In this paper, we further study the stability of nonlinear fractional differential systems with Caputo derivative by utilizing a Lyapunov-like function. Taking into account the relation between asymptotical stability and generalized Mittag-Leffler stability, we are able to weaken the conditions assumed for the Lyapunov-like function. In addition, based on the comparison principle of fractional differential equations 16, 17 , we also study the stability of nonlinear fractional differential systems by utilizing the comparison method. Our contribution in this paper is that we have relaxed the condition of the Lyapunov-like function and that we have further studied the stability. The present paper is organized as follows. In Section 2, some definitions and lemmas are introduced. In Section 3, sufficient conditions on asymptotical stability and generalized Mittag-Leffler stability are given. The comparison method is applied to the analysis of the stability of fractional differential systems in Section 4. Conclusions are included in the last section.

Preliminaries and Notations
Let us denote by R the set of nonnegative real numbers, by R the set of real numbers, and by Z the set of positive integer numbers. Let 0 < q < 1 and set C q t 0 , T , x t ∈ C t 0 , T ×Ω, R }, where C t 0 , t , R denotes the space of continuous functions on the interval t 0 , t . Let us first introduce several definitions, results, and citations needed here with respect to fractional calculus which will be used later. As to fractional integrability and differentiability, the reader may refer to 18 .
Definition 2.1. The fractional integral with noninteger order q ≥ 0 of function x t is defined as follows: where Γ · is the Gamma function.
International Journal of Differential Equations 3 Definition 2.2. The Riemann-Liouville derivative with order q of function x t is defined as follows: where m − 1 ≤ q < m and m ∈ Z .
Definition 2.3. The Caputo derivative with noninteger order q of function x t is defined as follows: where m − 1 < q < m and m ∈ Z .

Definition 2.4. The Mittag-Leffler function is defined by
where α > 0, z ∈ R. The two-parameter Mittag-Leffler function is defined by where α > 0 and β ∈ R, z ∈ R.
Clearly E α z E α,1 z . The following definitions are associated with the stability problem in the paper.
Definition 2.5. The constant x eq is an equilibrium of fractional differential system D q t 0 ,t x t f t, x if and only if f t, x eq D q t 0 ,t x t | x t x eq for all t > t 0 , where D q t 0 ,t means either the Caputo or the Riemann-Liouville fractional derivative operator.
Throughout the paper, we always assume that x eq 0. Definition 2.6 see 15 . The zero solution of D q t 0 ,t x t f t, x t with order q ∈ 0, 1 is said to be stable if, for any initial value x 0 , there exists an ε > 0 such that x t ≤ ε for all t > t 0 . The zero solution is said to be asymptotically stable if, in addition to being stable, x t → 0 as t → ∞.
where t 0 is the initial time and x 0 is the corresponding initial value, q ∈ 0, 1 , λ ≥ 0, b > 0, m 0 0, m x ≥ 0, and m x is locally Lipschitz on x ∈ B ⊂ R n with the Lipschitz constant L 0 .

x t f t, x t is said to be generalized Mittag-Leffler stable if
where t 0 is the initial time and x 0 is the corresponding initial value, q ∈ 0, 1 , −q < γ ≤ 1 − q, λ ≥ 0, b > 0, m 0 0, m x ≥ 0, and m x is locally Lipschitz on x ∈ B ⊂ R n with the Lipschitz constant L 0 . Remark 2.9. Mittag-Leffler stability and generalized Mittag-Leffler stability both belong to algebraical stability, which also imply asymptotical stability see 15 .
Definition 2.10. A function α r is said to belong to class-K if α : R → R is continuous function such that α 0 0 and it is strictly increasing.
Definition 2.11 see 19 . The class-K functions α r and β r are said to be with local growth momentum at the same level if there exist r 1 > 0, k i > 0 i 1, 2 such that k 1 α r ≥ β r ≥ k 2 α r for all r ∈ 0, r 1 . The class-K functions α r and β r are said to be with global growth momentum at the same level if there exist k i > 0 i 1, 2 such that k 1 α r ≥ β r ≥ k 2 α r for all r ∈ R .
It is useful to recall the following lemmas for our developments in the sequel.

Lemma 2.12 see 20 .
Let v, w ∈ C 1−q t 0 , T , R be locally Hölder continuous for an exponent Remark 2.13. In Lemma 2.12, if we replace RL D q t 0 ,t by C D q t 0 ,t , but other conditions remain unchanged, then the same result holds.

Assume that both inequalities are nonstrict and h t, x is nondecreasing in x for each t. Further, suppose that h satisfies the standard Lipschitz condition
International Journal of Differential Equations 5 Remark 2. 15. In Lemmas 2.12 and 2.14, T can be ∞.

Stability of Nonlinear Fractional Differential Systems
Let us consider the following nonlinear fractional differential system 14, 15 : with the initial condition x 0 x t 0 , where f : t 0 , ∞ × Ω → R n is piecewise continuous in t and Ω ⊂ R n is a domain that contains the equilibrium point x eq 0, 0 < q < 1. Here and throughout the paper, we always assume there exists a unique solution x t ∈ C 1 t 0 , ∞ to system 3.1 with the initial condition x t 0 .
Recently, Li et al. 14, 15 investigated the Mittag-Leffler stability and the generalized Mittag-Leffler stability the asymptotic stability of system 3.1 by using the fractional Lyapunov's second method, where the following theorem has been presented.
In the following, we give a new proof for Theorem 3.1.

Proof of Theorem 3.1. From 3.2 and 3.3 , we can get
Obviously, for the initial value V 0, x 0 , the linear fractional differential equation has a unique solution V t, x t V 0, x 0 E p −α 3 /α 2 t p . Taking into account Remark 2.13 and the relationship between 3.4 and 3.5 , we obtain International Journal of Differential Equations where E p −α 3 /α 2 t p is a nonnegative function 21 . Substituting 3.6 in 3.2 yields where E p −α 3 /α 2 t p → 0 t → ∞ from the asymptotic expansion of Mittag-Leffler function 22 . Hence the proof is completed.
According to the above results, we have the following theorem.
where x ∈ D, p ∈ 0, 1 . Then x eq 0 is locally stable. If the assumptions hold globally on R n , then x eq 0 is globally stable.
Proof. Proceeding the same way as that in the proof of Theorem 3.1, it follows from 3.9 that V t, x t ≤ V t 0 , x t 0 . Again taking into account 3.8 , one can get where t ≥ t 0 . Therefore, the equilibrium point x eq 0 is stable. So the proof is finished.
In the above two theorems, the stronger requirements on function V have been assumed to ensure the existence of C D p t 0 ,t V t, x t . This undoubtedly increases the difficulty in choosing the function V t, x t . In fact, we can weaken the continuously differential function V t, x t as V t, x t ∈ C 1−p t 0 , ∞ × D, R . Here we give the corresponding results.

Theorem 3.3. Let x eq
0 be an equilibrium point of system 3.1 , and let D ⊂ R n be a domain containing the origin, V t, x t ∈ C 1−p t 0 , ∞ × D, R . Assume there exists a class-K function α such that where t > t 0 ≥ 0, x ∈ D, and p ∈ 0, 1 . Then x eq 0 is locally asymptotically stable. If the assumptions hold globally on R n , then x eq 0 is globally asymptotically stable.
Proof. Note that the linear fractional differential equation Taking into account Lemma 2.12 and the relationship between 3.12 and 3.13 , we obtain 3.14 Substituting 3.14 into 3.11 gives from the definition of class-K. This completes the proof.
where t > t 0 ≥ 0, x ∈ D, p ∈ 0, 1 , and a, b are arbitrary positive constants. Then x eq 0 is generalized Mittag-Leffler stable. If the assumptions hold globally on R n , then x eq 0 is globally generalized Mittag-Leffler stable.
Proof. In Theorem 3.3, by replacing α x by a x b , we can get so the conclusion holds.

3.18
ii there exists a > 0 such that α 1 r and r a have global growth momentum at the same level, where t > t 0 ≥ 0, x ∈ D, and p ∈ 0, 1 . Then x eq 0 is locally generalized Mittag-Leffler stable. If the assumptions hold globally on R n , then x eq 0 is globally generalized Mittag-Leffler stable.
Proof. It follows from condition i that there exists k 1 > 0 such that

3.19
On the other hand, the linear fractional differential equation has a unique solution for the initial value V 0 Γ p V t, x t t − t 0 1−p | t t 0 . Using 3.19 , 3.20 , and Lemma 2.12, we obtain where E p,p −k 1 t − t 0 p is a nonnegative function 23, 24 . In addition, using condition ii , one gets for all x ∈ D, where k 2 > 0. Substituting 3.23 into 3.22 , we finally obtain

3.24
International Journal of Differential Equations 9 Hence, the zero solution of system 3.1 is locally generalized Mittag-Leffler stable. If the assumptions hold globally on R n , then x eq 0 is globally generalized Mittag-Leffler stable. The proof is completed.
Remark 3.6. The nonnegative function E p,p −k 1 t − t 0 p tends to zero as t approaches infinity from the asymptotic expansion of two-parameter Mittag-Leffler function 22 , so the zero solution of system 3.1 satisfying the conditions of Theorem 3.5 is also asymptotically stable.

The Comparison Results on the Stability
It is well known that the comparison method is an effective way in judging the stability of ordinary differential systems. In this section, we will discuss similar results on the stability of fractional differential systems by using the comparison method.
In what follows, we consider system 3.1 with f t, 0 0 and the scalar fractional differential equation where the initial value u 0 ∈ R , u t ∈ C 1−q t 0 , ∞ , R , g ∈ C t 0 , ∞ × R, R is Lipschitz in u and g t, 0 0, 0 < q < 1. Also, we assume there exists a unique solution u t t ≥ t 0 for system 4.1 with the initial value u 0 . Theorem 4.1. For system 3.1 , let x eq 0 be an equilibrium point of system 3.1 , and let Ω ⊂ R n be a domain containing the origin. Assume that there exist a Lyapunov-like function V ∈ C 1−q t 0 , ∞ × Ω, R and a class-K function α such that V t, 0 0, V t, x ≥ α x , and V t, x satisfies the inequality Suppose further that g t, x is nondecreasing in x for each t.
i If the zero solution of 4.1 is stable, then the zero solution of system 3.1 is stable; ii if the zero solution of 4.1 is asymptotically stable, then the zero solution of system 3.1 is asymptotically stable, too.
Proof. Let x t x t, t 0 , x 0 denote the solution of system 3.1 with initial value x 0 ∈ Ω. Along the solution curve x t , V t, x t can be written as V t and Applying the fractional integral operator D Now, taking u 0 V 0 and applying Lemma 2.14 to inequalities 4.3 and 4.4 , one has V t ≤ u t , t > t 0 .
i If the zero solution of 4.1 is stable, then for any initial value u 0 ≥ 0, there exists > 0 such that |u t | < for all t > t 0 . Therefore, taking into account V t, x t ≥ α x , one gets that is, x t < α −1 , and the zero solution of system 3.1 is stable.
ii One can directly derive from the same argument in i . Then, taking the limit to both sides of 4.6 and combining with the definition of class-K function, one can obtain lim t → ∞ x t 0.
The proof is thus finished.

and V t, x is locally Lipschitz in x and satisfies the inequality
where m 0 0, m x ≥ 0 and m x is locally Lipschitz in x with Lipschitz constant L 0 .
International Journal of Differential Equations

11
Taking u 0 V 0 Γ q V t, x t − t 0 1−q | t t 0 and noting that V t, x ≤ u t holds from Theorem 4.1, then taking into account 4.8 and V t, x ≥ k x a , we obtain Furthermore,

4.10
Let x ≥ 0 and k > 0. In addition, M x is locally Lipschitz in x since m x and V t, x are locally Lipschitz in x. So, the zero solution of system 3.1 is generalized Mittag-Leffler stable. The proof is completed.

Conclusion
In this paper, we have studied the stability of the zero solution of nonlinear fractional differential systems with the Caputo derivative and the commensurate order 0 < q < 1 by using a Lyapunov-like function. Compared to 15 , we weaken the continuously differential function V t, x as V t, x ∈ C 1−p t 0 , ∞ ×D, R . Sufficient conditions on generalized Mittag-Leffler stability and asymptotical stability are derived. Meanwhile, comparison method is applied to the analysis of the stability of fractional differential systems by fractional differential inequalities.

Introduction
In the last century, notable contributions have been made to both the theory and applications of the fractional differential equations. For the theory part, Momani and Hadid have investigated the local and global existence theorem of both fractional differential equation and fractional integrodifferential equations; see 1-6 . Fractional-order differential equations have recently proved to be valuable tools in the modeling of many phenomena in various fields of science and engineering.
Integrodifferential equations with integral boundary conditions are often encountered in various applications; it is worthwhile mentioning the applications of those conditions in the study of population dynamics and cellular systems. For a detailed description of the integral boundary conditions, we refer the reader to a recent paper 7 . In 8 , Tidke studied the problem of existence of global solutions to nonlinear mixed Volttera-Fredholm integrodifferential equations with nonlocal condition.
Ahmad and Nieto 9 studied some existence results for boundary value problem involving a nonlinear integrodifferential equation of fractional order with integral equation.

International Journal of Differential Equations
Very recently N Guérékata 10 discussed the existence of solutions of fractional abstract differential equations with nonlocal initial condition. Anguraj et al. 11 studied the existence and uniqueness theorem for the nonlinear fractional mixed Volterra-Fredholm integrodifferential equation with nonlocal initial condition.
Motivated by these works, we study in this paper the existence of solution of boundary value problem for fractional integrodifferential equations in the case 1 < α ≤ 2 in Banach spaces by using Banach and Krasnosel'skii fixed-point theorems.

Preliminaries
First of all, we recall some basic definitions; see 12-15 .  where 1 < α ≤ 2, D α is the Caputo fractional derivative and the nonlinear functions f : 0, T × X × X × X → X, k, h 1 : 0, T × 0, T × X → X and g, h : X → X satisfy the following hypotheses: H3 there exists continuous functions p : 0, T → R 0, ∞ and p 1 : 0, T → R such that t 0 k t, s, x − k t, s, y ds ≤ p t x − y and t 0 k t, s, y ds ≤ p 1 t y , for every t, s ∈ 0, T and x, y ∈ X, H4 there exists continuous functions q : 0, T → R and q 1 : 0, T → R such that T 0 h 1 t, s, x − h 1 t, s, y ds ≤ q t x − y and T 0 h 1 t, s, y ds ≤ q 1 t y for every t, s ∈ 0, T andx, y ∈ X H5 there exists continuous function L : 0, T → R , and N 1 is positive constant such that f t, x 1 , y 1 , z 1 − f t, x 2 , y 2 , z 2 ≤ L t K x 1 − x 2 y 1 − y 2 z 1 − z 2 and N 1 sup t∈ 0,T f t, 0, 0, 0 , for every t ∈ 0, T and x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ∈ X, where K : R → 0, ∞ is continuous nondecreasing function satisfying K γ t x ≤ γ t K x , where γ is a continuous function γ : 0, T → R .

4 International Journal of Differential Equations
Proof. By Lemma 2.2, we reduce the problem 2.4 -2.5 to an equivalent integral equation

2.8
In view of the relations c D α I α y t y t and I α I β y t I α β y t , for α, β > 0, we obtain Applying the boundary condition 2.5 , we find that

Introduction
Fractional differential equations have been of great interest recently. It is caused both by the intensive development of the theory of fractional calculus itself and by the applications of such constructions in various sciences such as physics, mechanics, chemistry, and engineering. For details, see 1-6 and references therein. It should be noted that most of papers and books on fractional calculus are devoted to the solvability of linear initial fractional differential equations in terms of special functions 6-8 . Recently, there are some papers that deal with the existence and multiplicity of solution or positive solution of nonlinear initial fractional differential equation by the use of techniques of nonlinear analysis fixed-point theorems, Leray-Shauder theory, etc. ; see 9-17 . Recently, Bai and Lü 15 studied the existence of positive solutions of nonlinear fractional differential equation where 1 < α ≤ 2 is a real number, D α 0 is the standard Riemann-Liouville differentiation, and f : 0, 1 × 0, ∞ → 0, ∞ is continuous.
In this paper, we study the existence of positive solutions for fractional differential equation with nonlocal boundary condition We assume the following conditions hold throughout the paper:

The Preliminary Lemmas
For the convenience of the reader, we present here the necessary definitions from fractional calculus theory. These definitions can be found in the recent literature. where n α 1, provided the right side is pointwise defined on 0, ∞ .

Definition 2.3.
The map θ is said to be a nonnegative continuous concave functional on a cone P of a real Banach space E, provided that θ : P → 0, ∞ is continuous and for all x, y ∈ P and 0 ≤ t ≤ 1.
International Journal of Differential Equations 3 Remark 2.4. As a basic example, we quote for λ > −1, giving in particular D α 0 t α−m 0, m 1, 2, . . . , N, where N is the smallest integer greater than or equal to α.
From Definition 2.2 and Remark 2.4, we then obtain the following.

2.12
Proof. By applying Lemmas 2.6 and 2.7, we have

2.15
Lemma 2.9 see 15 . The function G t, s defined by 2.9 satisfies the following conditions: holds. Then, A has a fixed point in P ∩ Ω 2 \ Ω 1 .

Lemma 2.11 see 19 .
Let P be a cone in real Banach space E, P c {x ∈ P | x ≤ c}, θ a nonnegative continuous concave functional on P such that θ x ≤ x , for all x ∈ P c , and Suppose that A : P c → P c is completely continuous, and there exist constants 0 < a < b < d ≤ c such that Then, A has at least three fixed points x 1 , x 2 , x 3 with

2.17
Remark 2.12. If there holds d c, then condition C1 of Lemma 2.11 implies condition C3 of Lemma 2.11.

The Main Results
Let E C 0, 1 be endowed with the ordering u ≤ v if u t ≤ v t for all t ∈ 0, 1 , and the maximum norm, u max 0≤t≤1 |u t |. Define the cone P ⊂ E by P {u ∈ E | u t ≥ 0}. Let the nonnegative continuous concave functional θ on the cone P be defined by θ u min 1/4 ≤t≤ 3/4 |u t |.
then problem 1.2 has at least one positive solution u such that r 1 ≤ u ≤ r 2 .
Proof. By Lemmas 2.8 and 3.2, we know A : P → P is completely continuous, and problem 1.2 has a solution u u t if and only if u solves the operator equation u Au. In order to apply Lemma 2.10, we separate the proof into the following two steps.

3.4
International Journal of Differential Equations 7 So, Therefore, by ii of Lemma 2.10, we complete the proof.

3.7
With the use of Theorem 3.3, problem 3.6 has at least one positive solutions u such that 1/70 ≤ u ≤ 9/10 . Theorem 3.5. Assume (H1)-(H3) hold, and there exist constants 0 < a < b < c such that the following assumptions hold: Then, the boundary value problem 1.2 has at least three positive solutions u 1 , u 2 , u 3 with

8 International Journal of Differential Equations
Proof. We show that all the conditions of Lemma 2.9 are satisfied. If u ∈ P c , then u ≤ c. Assumption A3 implies f t, u t ≤ Mc for 0 ≤ t ≤ 1. Consequently, 3.9 Hence, A : P c → P c . In the same way, if u ∈ P a , then assumption A1 yields f t, u t < Ma, 0 ≤ t ≤ 1. Therefore, condition C2 of Lemma 2.11 is satisfied.
To check condition C1 of Lemma 2.11, we choose u t b c /2, 0 ≤ t ≤ 1. It is easy to see that u t b c /2 ∈ P θ, b, c , θ u θ b c /2 > b, and consequently, This shows that condition C1 of Lemma 2.11 is also satisfied. By Lemma 2.11 and Remark 2.12, the boundary value problem 1.2 has at least three positive solutions u 1 , u 2 , and u 3 with

3.11
The proof is complete.
International Journal of Differential Equations 9 Example 3.6. Consider the problem

3.14
With the use of Theorem 3.5, problem 3.12 has at least three positive solutions u 1 , u 2 , and u 3 with max 0≤t≤1 |u 1 t | < 3.15

Introduction
In this paper, we study the multipoint boundary value problem International Journal of Differential Equations f : 0, 1 × R 3 → R satisfying the Carathéodory conditions, e ∈ L 1 0, 1 . D α 0 and I α 0 are the standard Riemann-Liouville derivative and integral, respectively. We assume, in addition, that where Γ is the Gamma function. Due to condition 1.3 , the fractional differential operator in 1.1 , 1.2 is not invertible.
Fractional differential equation can describe many phenomena in various fields of science and engineering. Many methods have been introduced for solving fractional differential equations, such as the popular Laplace transform method, the iteration method. For details, see 1, 2 and the references therein.
Recently, there are some papers dealing with the solvability of nonlinear boundary value problems of fractional differential equation, by use of techniques of nonlinear analysis fixed-point theorems, Leray-Schauder theory, etc. , see, for example, 3-6 . But there are few papers that consider the fractional-order boundary problems at resonance. Very recently 7 , Y. H. Zhang and Z. B. Bai considered the existence of solutions for the fractional ordinary differential equation where n > 2 is a natural number, n − 1 < α ≤ n is a real number, f : 0, 1 × R n → R is continuous, and e ∈ L 1 0, 1 , σ ∈ 0, ∞ , and η ∈ 0, 1 are given constants such that ση α−1 1. D α 0 and I α 0 are the standard Riemann-Liouville derivative and integral, respectively. By the conditions, the kernel of the linear operator is one dimensional.
Motivated by the above work and recent studies on fractional differential equations 8-18 , in this paper, we consider the existence of solutions for multipoint boundary value problem 1.1 , 1.2 at resonance. Note that under condition 1.3 , the kernel of the linear operator in 1.1 , 1.2 is two dimensional. Our method is based upon the coincidence degree theory of Mawhin 18 . Now, we will briefly recall some notation and abstract existence result. Let Y, Z be real Banach spaces, let L : dom L ⊂ Y → Z be a Fredholm map of index zero, and let P : Y → Y, Q : Z → Z be continuous projectors such that Im P Ker P , Ker Q Im L , and Y Ker L ⊕ Ker P , Z Im L ⊕ Im Q . It follows that L| dom L ∩Ker P : dom L ∩ Ker P → Im L is invertible. We denote the inverse of the map by International Journal of Differential Equations 3 K P . If Ω is an open-bounded subset of Y such that dom L ∩Ω / ∅, the map N : Y → Z will be called L-compact on Ω if QN Ω is bounded and K P I − Q N : Ω → Y is compact.
The theorem that we used is Theorem 2.4 of 18 .
Theorem 1.1. Let L be a Fredholm operator of index zero and N be L-compact on Ω. Assume that the following conditions are satisfied: The rest of this paper is organized as follows. In Section 2, we give some notation and Lemmas. In Section 3, we establish an existence theorem of a solution for the problem 1.1 , 1.2 .

Background Materials and Preliminaries
For the convenience of the reader, we present here some necessary basic knowledge and definitions for fractional calculus theory, and these definitions can be found in the recent literature 1, 2 .
Definition 2.1. The fractional integral of order α > 0 of a function y : 0, ∞ → R is given by provided the right side is pointwise defined on 0, ∞ . And we let I 0 0 y t y t for every continuous y : 0, ∞ → R.
Definition 2.2. The fractional derivative of order α > 0 of a function y : 0, ∞ → R is given by where n α 1, provided the right side is pointwise defined on 0, ∞ .

International Journal of Differential Equations
We use the classical space C 0, 1 with the norm x ∞ max t∈ 0,1 |x t |. Given μ > 0 and N μ 1, one can define a linear space where x ∈ C 0, 1 and c i ∈ R, i 1, 2, . . . , N − 1. By means of the linear function analysis theory, one can prove that with the norm and equicontinuous means that for all ε > 0, ∃δ > 0 such that

2.6
Let Z L 1 0, 1 with the norm g 1 Definition 2.6. We say that the map f : 0, 1 × R → R satisfies the Carathéodory conditions with respect to L 1 0, 1 if the following conditions are satisfied: i for each z ∈ R, the mapping t → f t, z is Lebesgue measurable, ii for almost every t ∈ 0, 1 , the mapping z → f t, z is continuous on R, iii for each r > 0, there exists ρ r ∈ L 1 0, 1 , R such that, for a.e., t ∈ 0, 1 and every |z| ≤ r, we have |f t, z | ≤ ρ r t .
Define L to be the linear operator from dom L ∩ Y to Z with  hold. Then u t satisfies the boundary conditions 1.2 , that is, u ∈ dom L , and we have g ∈ Z | g satisfies 2.11 ⊆ Im L .

2.12
Let u ∈ dom L , then for D α 0 u ∈ Im L , we have which, due to the boundary value condition 1.2 , implies that D α 0 u satisfies 2.11 . In fact, from International Journal of Differential Equations Therefore, Im L g ∈ Z | g satisfies 2.11 .

2.17
Consider the continuous linear mapping Q 1 : Z → Z and Q 2 : Z → Z defined by 2.18 Using the above definitions, we construct the following auxiliary maps R 1 , R 2 : Z → Z:

2.19
Since the condition 1.4 holds, the mapping Q : Z → Z defined by Qy t R 1 g t t α−1 R 2 g t t α−2 2.20 is well defined. Recall 1.4 and note that

2.21
International Journal of Differential Equations 7 and similarly we can derive that

2.22
So, for g ∈ Z, it follows from the four relations above that that is, the map Q is idempotent. In fact, Q is a continuous linear projector. Note that g ∈ Im L implies Qg 0. Conversely, if Qg 0, then we must have R 1 g R 2 g 0; since the condition 1.4 holds, this can only be the case if Q 1 g Q 2 g 0, that is, g ∈ Im L . In fact, Im L Ker Q . Take g ∈ Z in the form g g − Qg Qg, so that g − Qg ∈ Im L Ker Q and Qg ∈ Im Q . Thus, Z Im L Im Q . Let g ∈ Im L ∩ Im Q and assume that g s as α−1 bs α−2 is not identically zero on 0, 1 , then, since g ∈ Im L , from 2.11 and the condition 1.4 , we derive a b 0, which is a contradiction. Hence, Im L ∩ Im Q {0}; thus, Z Im L ⊕ Im Q . Now, dim Ker L 2 co dim Im L , and so L is a Fredholm operator of index zero.
Let P : Y → Y be defined by

2.24
Note that P is a continuous linear projector and It is clear that Y Ker L ⊕ Ker P . Note that the projectors P and Q are exact. Define K P : Im L → dom L ∩ Ker P by t − s g s ds, 2.27 K P g ∞ ≤ g 1 , and thus In fact, if g ∈ Im L , then LK P g D α 0 I α 0 g g. Also, if u ∈ dom L ∩ Ker P , then

2.31
With arguments similar to those of 7 , we obtain the following Lemma.

The Main Results
Assume that the following conditions on the function f t, x, y, z are satisfied: H1 there exists a constant A > 0, such that for u ∈ dom L \ Ker L satisfying |D α−1 0 u t | |D α−2 0 u t | > A for all t ∈ 0, 1 , we have H2 there exist functions a, b, c, d, r ∈ L 1 0, 1 and a constant θ ∈ 0, 1 such that for all x, y, z ∈ R 3 and a.e., t ∈ 0, 1 , one of the following inequalities is satisfied: 3.2 H3 there exists a constant B > 0 such that for every a, b ∈ R satisfying a 2 b 2 > B, then either Remark 3.1. R 1 N at α−1 bt α−2 and R 2 N at α−1 bt α−2 from H3 stand for the images of u t at α−1 bt α−2 under the maps R 1 N and R 2 N, respectively. Proof. Set then for u ∈ Ω 1 , Lu λNu; thus, λ / 0, Nu ∈ Im L Ker Q , and hence QNu t 0 for all t ∈ 0, 1 . By the definition of Q, we have Q 1 Nu t Q 2 Nu t 0. It follows from H1 that there exists t 0 ∈ 0, 1 such that |D α−1 International Journal of Differential Equations

3.8
Now by 3.8 , we have

3.9
Note that I − P u ∈ Im K P dom L ∩ Ker P for u ∈ Ω 1 , then, by 2.28 and 2.30 ,

3.10
International Journal of Differential Equations 11 Using 3.9 and 3.10 , we obtain If the first condition of H2 is satisfied, then we have

3.13
and consequently, International Journal of Differential Equations Note that θ ∈ 0, 1 and a 1 b 1 c 1 < 1/τ, so there exists M 1 > 0 such that D α−1 0 u ∞ ≤ M 1 for all u ∈ Ω 1 . The inequalities 3.14 and 3.15 show that there exist M 2 , M 3 , that is, Ω 1 is bounded given the first condition of H2 . If the other conditions of H2 hold, by using an argument similar to the above, we can prove that Ω 1 is also bounded. Let We define the isomorphism J : Im Q → Ker L by If the first part of H3 is satisfied, let For every at α−1 bt α−2 ∈ Ω 3 , 3.20 If λ 1, then a b 0, and if a 2 b 2 > B, then by H3 , λ a 2 b 2 1 − λ aR 1 N at α−1 bt α−2 bR 2 N at α−1 bt α−2 < 0, 3.21 which, in either case, obtain a contradiction. If the other part of H3 is satisfied, then we take and, again, obtain a contradiction. Thus, in either case,

3.23
International Journal of Differential Equations 13 for all u ∈ Ω 3 , that is, Ω 3 is bounded.
In the following, we will prove that all the conditions of Theorem 1.1 are satisfied. Set Ω to be a bounded open set of Y such that U 3 i 1 Ω i ⊂ Ω. by Lemma 2.8, the operator K P I − Q N : Ω → Y is compact; thus, N is L-compact on Ω, then by the above argument, we have i Lu / λNx, for every u, λ ∈ dom L \ Ker L ∩ ∂Ω × 0, 1 , ii Nu / ∈ Im L , for every u ∈ Ker L ∩ ∂Ω.
Finally, we will prove that iii of Theorem 1.1 is satisfied. Let H u, λ ±Iu 1 − λ JQNu, where I is the identity operator in the Banach space Y . According to the above argument, we know that

Introduction
Mathematical model of the movement of gas in a channel surrounded by a porous environment was described by parabolic-hyperbolic equation. This was done in the fundamental work of Gel'fand 1 . Modeling of heat transfer processes in composite environment with finite and infinite velocities leads to boundary value problems BVPs for parabolic-hyperbolic equations 2 . Omitting the huge amount of works devoted to studying these kinds of equations, we refer the readers to 3, 4 . We would like to note works 5-10 , devoted to the studying of BVPs for parabolichyperbolic equations, involving fractional derivatives. In turn, applications of Fractionalorder differential equations can be found in the monographs 11-15 . We also note some recent papers 16-18 , related to the fractional diffusion and diffusion-wave equations.
BVP for parabolic-hyperbolic equations with integral gluing condition for the first time was investigated by Kapustin and Moiseev 19 and was generalized for this kind of equation, but with parameters, in the work 20 . Another motivation of the usage of integral gluing conditions comes from the appearance of them in heat exchange processes 21 .
The consideration of equations with parameters was interesting because of the possibility of studying some multidimensional analogues of the main BVP via reducing them by Fourier transformation to the BVP for equations with parameters. On the other hand, consideration of equations with parameters will give possibility to study some spectral properties of BVPs for this kind of equations such as the existence of nontrivial solutions for corresponding homogeneous problem at some values of parameters 22 .

Analog of the Tricomi Problem
Consider an equation is the αth Riemann-Liouville fractional-order derivative of a function f given on interval a, b , where n α 1 and α is the integer part of α, and Γ · is the Euler gamma function defined by For λ > 0 and 0 < α ≤ 1 given, we formulate the following problem called the analog of the Tricomi problem.

Problem AT
To find a solution of 2.1 , which belongs to the class of functions

2.7
Here ω x , ψ i y i 1, 2 are given functions such as lim y → 0 y 1−α ψ 1 y ω 0 . Solution of the Cauchy problem for 2.1 in Ω 2 defined as where J k · is the first-kind Bessel function of the order k, τ − y u 0 − , y , ν − y u x 0 − , y . We calculate u −y/2, y/2 in order to use condition 2.5 : u −y/2, y/2

2.9
Considering the condition 2.5 and the following integral operator 23 m, n 0, 1, 2.10 equality 2.9 can be written as follows 2.11 Now we use an integral operator Considering the following properties of operators 2.10 and 2.12 Taking gluing conditions 2.7 into account, we have Let us consider the following auxiliary problem: where Ψ y ψ * 1 y − P 0, y , K 2 y, η is the resolvent of the kernel K 1 y, η .

6
International Journal of Differential Equations Once we have obtained ν y , considering 2.18 or 2.25 we find function τ y . Then using gluing conditions 2.7 we find functions τ − y , ν − y . Finally, we can define solution of the considered problem by the formula 2.23 in the domain Ω 1 , by formula 2.8 in the domain Ω 2 .
Hence, we prove the following theorem.

Analog of the Gellerstedt Problem
We would like to note some related works. Regarding the consideration of Gellerstedt problem for parabolic-hyperbolic equations with constant coefficients we refer the readers to 3 and for loaded parabolic-hyperbolic equations work by Khubiev 27 , and also for Lavrent'ev-Bitsadze equation 28 .
Consider an equation where Φ 0 is a domain, bounded by segments AA 0 , BB 0 , A 0 B 0 of straight lines x 0, x 1, y 1, respectively; Φ 1 is a domain, bounded by the segment AE of the axe x and by characteristics of 3.1 AC 1 : x y 0, EC 1 : x − y r; Φ 2 is a domain, bounded by the segment EB of the axe x and by characteristics of 3.1 EC 2 : x − y r, BC 2 : x − y 1; I 0 is an interval 0 < x < 1, I 1 is an interval 0 < x < r, and I 2 is an interval r < x < 1.

Problem AG
To find a solution of 3.1 from the class of functions satisfying boundary conditions u 0, y ϕ 1 y , u 1, y ϕ 2 y , 0 ≤ y ≤ 1, 3.3 International Journal of Differential Equations together with gluing conditions lim y → 0 Here ϕ j · j 1, 4 are given functions such as lim y → 0 y 1−α ϕ 1 y ϕ 3 0 , lim y → 0 y 1−α u r, y ϕ 4 r .

Theorem 3.1. If the following conditions
are fulfilled, then the Problem AG has a unique solution.
Proof. Introduce the following designations: 3.9 Solution of the Cauchy problem for 3.1 in the domain Φ i i 1, 2 in case, when λ ≥ 0 has a form u x, y 1 2 3.10 8 International Journal of Differential Equations Using boundary conditions 3.4 , 3.5 , and gluing conditions 3.6 , 3.7 , from 3.10 we obtain 3.13 According to 10 , tending y to 0, from 3.1 we get In order to prove the uniqueness of the solution for the Problem AG, we need estimate the following integral: Considering homogeneous case of the condition 3.3 and taking designation 3.9 into account, after some evaluations we derive If λ ≥ 0, then I ≤ 0. On the other hand, if we consider homogeneous cases of 3.11 and 3.12 , one can easily be sure that I ≥ 0. Hence, we get that I ≡ 0. Based on 3.16 we can conclude that τ x 0 for all x ∈ I 0 . Due to the solution of the first boundary problem 24 we can conclude that u x, y ≡ 0 in Φ 0 . Further, according to the gluing conditions and the solution of Cauchy problem, we have u x, y ≡ 0 in Φ.

3.18
International Journal of Differential Equations 9 The problems 3.17 and 3.18 are model problems and can be solved directly. After the finding function τ x for all x ∈ I 0 , functions ν x and τ − x , ν − x can be defined by formulas 3.14 and 3.6 , 3.7 , respectively. Finally, solution of the Problem AG can be recovered by formulas 3.10 and 2.23 in the domains Φ i i 1, 2 and Φ 0 , respectively, but only with some changes in 2.23 , precisely, Green function G x, y, ξ, η should be replaced by which is the Green function of the first boundary problem for the 3.1 in Φ 0 24 . Theorem 3.1 is proved.

Introduction
There are many fluids in industry and technology whose behavior cannot be explained by the classical linearly viscous Newtonian model. The departure from the Newtonian behavior manifests itself in a variety of ways: non-Newtonian viscosity shear thinning or shear thickening , stress relaxation, nonlinear creeping, development of normal stress differences, and yield stress 1 . The Navier-Stokes equations are inadequate to predicted the behavior of such type of fluids; therefore, many constitutive relations of non-Newtonian fluids are proposed 2 . These constitutive relations give rise to the differential equations, which, in general, are more complicated and higher order than the Navier-Stokes equations. Therefore, it is difficult to obtain exact analytical solutions for non-Newtonian fluids 3 .
Modeling of the equation governing the behaviors of non-Newtonian fluids in different circumstance is important from many points of view. For examples, plastics and polymers are extensively handeled by the chemical industry, whereas biological and biomedical devices like hemodialyser make use of the rheological behavior 4 . In general, the analysis of the behavior of the fluid motion of non-Newtonian fluids tends to be much more complicated and subtle in comparison with that of the Newtonian fluids 5 .
The fractional calculus, almost as old as the standard differential and integral one, is increasingly seen as an efficient tool and subtle frame work within which useful generalization is quite long and arguments almost yearly. It includes fractal media, fractional wave diffusion, fractional Hamiltonian dynamics, and biopolymer dynamics as well as many other topics in physics. Fractional calculus is useful in the field of biorheology and bioengineering, in part, because many tissue-like materials polymers, gels, emulsions, composites, and suspensions exhibit power-law responses to an applied stress or strain 6, 7 . An example of such power-law behavior in elastic tissue was observed recently for viscoelastic measurements of the aorta, both in vivo and in vitro 8, 9 , and the analysis of these data was most conveniently performed using fractional order viscoelastic models. The starting point of the fractional derivative model of non-Newtonian model is usually a classical differential equation which is modified by replacing the time derivative of an integer order by the so-called Riemann-Liouville/Caputo fractional calculus operators. This generalization allows one to define precisely noninteger order integrals or derivatives. In general, fractional model of viscoelastic fluids is derived from well-known ordinary model by replacing the ordinary time derivatives, to fractional order time derivatives and this plays an important role to study the valuable tool of viscoelastic properties. We include here some investigation 10-18 in which the fractional calculus approach has been adopted for the flows of non-Newtonian fluids. Furthermore, the one-dimensional fractional derivative Maxwell model has been found very useful in modeling the linear viscoelastic response of some polymers in the glass transition and the glass state 19 . In other cases, it has been shown that the governing equations employing fractional derivatives are also linked to molecular theories 20 . The use of fractional derivatives within the context of viscoelasticity was firstly proposed by Germant  A general view of the literature shows that the slip effects on the flows of non-Newtonian fluids has been given not much attention. Especially, polymer melts exhibit a macroscopic wall slip. The fluids exhibiting a boundary slip are important in technological applications, for example, the polishing of artificial heart valves, rarefied fluid problems, and flow on multiple interfaces. In the study of fluid-solid surface interactions, the concept of slip of a fluid at a solid wall serves to describe macroscopic effects of certain molecular phenomena. When the molecular mean free path length of the fluid is comparable to the distance between the plates as in nanochannels or microchannels, the fluid exhibits noncontinuum effects such as slipflow, as demonstrated experimentally by Derek et al. 25 . Experimental observations show that 26-28 non-Newtonian fluids, such as polymer melts, often exhibit macroscopic wall slip, which, in general, is described by a nonlinear and nonmonotone relation between the wall slip velocity and the traction. A more realistic class of slip flows are those in which the magnitude of the shear stress reaches some critical value, here called the slip yield stress, before slip occurs. In fact, some experiments show that the onset slip and slip velocity may also depend on the normal stress at the boundary 26, 29 . Much of the research involving slip presumes that the slip velocity depends on the shear stress. The slip condition is an important factor in sharskin, spurt, and hysteresis effects, but the existing theory for non-Newtonian fluids with wall slippage is scant. We mention here some recent attempt regarding exact analytical solutions of non-Newtonian fluids with slip effects 30-36 .
The objective of this paper is twofold. Firstly, is to give few more exact analytical solutions for viscoelastic fluids with fractional derivative approach, which is more natural and appropriate tool to describe the complex behavior of such fluids. Secondly, is to study the slip effects on viscoelastic fluid flows, which is important due to their practical applications. More precisely, our aim is to find the velocity field and the shear stress corresponding to the motion of a Maxwell fluid due to a sudden moved plate, where no-slip assumption is no longer valid. However, for completeness, we will determine exact solutions for a larger class of such fluids. Consequently, motivated by the above remarks, we solve our problem for Maxwell fluids with fractional derivatives. The general solutions are obtained using the discrete Laplace transforms. They are presented in series form in terms of the Wright generalized hypergeometric functions p Ψ q and presented as sum of the slip contribution and the corresponding no-slip contributions. The similar solutions for ordinary Maxwell fluids can easily be obtained as limiting cases of general solutions. The Newtonian solutions are also obtained as special cases of fractional and ordinary Maxwell fluids. Furthermore, the solutions for fractional and ordinary Maxwell fluid for no-slip condition also obtained as a special cases, and they are similar with previously known results in the literature. Finally, the influence of the material, slip and fractional parameters on the motion of fractional and ordinary Maxwell fluids is underlined by graphical illustrations. The difference among fractional Maxwell, ordinary Maxwell, and Newtonian fluid models is also highlighted.

The Differential Equations Governing the Flow
The equations governing the flow of an incompressible fluid include the continuity equation and the momentum equation. In the absence of body forces, they are where ρ is the fluid density, V is the velocity field, t is the time, and ∇ represents the gradient operator. The Cauchy stress T in an incompressible Maxwell fluid is given by 10, 11, 14-17 where −pI denotes the indeterminate spherical stress due to the constraint of incompressibility, S is the extrastress tensor, L is the velocity gradient, A L L T is the first Rivlin Ericsen tensor, μ is the dynamic viscosity of the fluid, λ is relaxation time, the superscript T indicates the transpose operation, and the superposed dot indicates the material time derivative. The model characterized by the constitutive equations 2.2 contains as special

Numerical Results and Conclusions
In this paper, the unsteady flow of fractional Maxwell fluid over an infinite plate, where the no-slip assumption between the wall and the fluid is no longer valid, is studied by means of the discrete Laplace transforms. The motion of the fluid is due to the plate that at time t 0 is suddenly moved with a constant velocity U in its plane. Closed-form solutions are obtained for the velocity u y, t and the shear stress τ y, t in series form in terms of the Wright generalized hypergeometric functions. These solutions, presented as a sum of the slip contribution and the corresponding no-slip contributions, satisfy all imposed initial and boundary conditions. The corresponding solutions for ordinary Maxwell fluids are also obtained from general solutions for α → 1. In the special case when θ → 0, the general solution reduces to previously known results for Stokes' first problem. In order to reveal some relevant physical aspects of the obtained results, the diagrams of the velocity field u y, t and the shear stress τ y, t have been drawn against y and t for different values of t, material constants ν, λ, slip parameter θ, and fractional parameter α. From all figures, it is clear that increasing the slip parameter at the wall the velocity decreases at the wall. Figures 1 and 2   The influence of relaxation time and kinematic viscosity ν on fluid motion are presented in Figures 5-8. As expected, the two material parameter have opposite effects on fluid motion. For instance, the velocity and shear stress decreases with respect to λ. More important for us is to see the effects of fractional and slip parameter on fluid motion. It is observed that velocity and shear stress either slip effects present or not are increasing functions of fractional parameter, as shown in Figures 9 and 10. The effect of slip parameter is clear from Figure 11. Finally, for comparison, the velocity field and the shear stress corresponding to the three models fractional Maxwell, ordinary Maxwell, and Newtonian are together depicted in  Figure 13: Profiles of the velocity field u y, t and the shear stress τ y, t for fractional Maxwell, ordinary Maxwell and Newtonian fluids, for U 1, y 1, ν 0.379, μ 33, λ 1.5, α 0.2, t 2 s, and θ 0.2. moving plate. For higher values of slip parameter the fractional Maxwell fluid is swiftest and the ordinary Maxwell is slowest and shear stress corresponding to fractional Maxwell fluid is largest on the whole flow domain as it is clear from Figure 14. It is important to note the difference between fractional and ordinary Maxwell fluid that, when slip effect is not present, the ordinary Maxwell fluid have oscillating behavior near the moving plate as

Antisynchronization between Financial and Chen Systems of Fractional Order
Assuming that Chen system is antisynchronized with Financial system; define the drive system as

5.4
The linear functions V 4 , V 5 , V 6 are given by

5.5
With the values given in 5.4 and 5.5 , the error system 5. It can be observed that the coefficient matrix of the error system 5.6 has eigenvalues −1, −1, −1. So the system is stable and antisynchronization is achieved.

Simulations and Results
We take parameters for fractional-order Chen system as a 1 35, b 1 3, c 1 27. Parameters for the Financial system are same as given in Section 4.1. Experiments are done for fixed value of fractional-order α 0.95, which is same for drive and response system 5.1 and 5.2 . The initial conditions for the systems 5.1 and 5.2 are x 1 0 2, y 1 0 3, z 1 0 2 and x 2 0 10, y 2 0 25, z 2 0 36, respectively. For the error system 5.6 , the initial conditions turns out to be e 1 0 12, e 2 0 28, e 3 0 38. The simulation results are summarized in Figure 2. Antisynchronization between fractional Financial and Chen system is shown in Figure 2

Antisynchronization between Fractional Lü and Financial System
In this case, consider Lü system as the drive system 6.2