Fractional calculus, the theory of differentiation and integration to non-integer order, is very useful for the description of various physical phenomena, such as damping laws, diffusion process, etc. Fractional differential equations extend prominent tools that are perfectly representing many engineering and physical problems. However, solving fractional differential equations are a challenging and stimulating area of research in mathematics and engineering since these fractional equations do not have exact and analytic solutions and that is the main reason made accurate numerical techniques preferable for solving fractional differential equations.
Many applications of shifted Legendre polynomials have been exemplified in research    . In this article, an application of Legendre polynomials to solve fractional differential equations is provided. The main aim is to generalize Legendre operational matrix to fractional calculus and based on using the operational matrix of an orthogonal function for solving differential equations. In the last decade or so, extensive research has been done in the numerical development methods that are numerically stable for both linear and nonlinear fractional differential equations  - . Orthogonal functions and polynomial series have received appreciable attention in dealing with various problems of dynamic systems. Legendre polynomials are well known family of orthogonal polynomials on the interval [−1, 1] that have many applications . They are widely used because of their good properties in the approximation of functions.
The proposed method is to obtain the numerical solution for FDE of the form  presented in Equation (I)-(III) based on the shifted Legendre method:
subject to the conditions
or the boundary conditions
where are constants.
and f(x) are the source terms.
The fractional derivatives are defined in the Caputo sense. The main idea in the current work is to apply the shifted Legendre polynomials and the operational matrix of fractional derivative together to discretize Equation (1) to get a satisfactory result. Figures 1-3 and results in next sections show the effectiveness of the proposed method in comparison with the analytical results.
The remainder of the article is organized as follows. In the next section, some mathematical preliminaries of the fractional calculus theory are introduced in addition to some relevant properties of the Legendre polynomials. Section 3 summarizes the application of the shifted Legendre method to the solution of problems (I)-(III). As a result, a system of algebraic equations is obtained and the solution of the considered problem is given. Section 4 illustrates applying the Legendre operational matrix of fractional derivative for solving multi-order fractional differential equation. In Section 5, the proposed method is applied to several examples. Finally, conclusion is given in Section 6. The numerical results are all computed using Mathematica Development Environment.
Figure 1. Comparison of y(x) Analytical and Approximate for m = 3.
Figure 2. Comparison of y(x) Analytical and Approximate for m = 5.
Figure 3. Comparison of y(x) Analytical and Approximate for m = 3.
2. LOM Preliminaries and Notations
In this section, some basic definitions and properties of fractional calculus theory are given that are further used in this article.
2.1. The Fractional Derivative in the Caputo Sense
The fractional calculus is a name for the theory of integrals and derivatives of arbitrary order, which unifies and generalizes the notions of integer-order differentiation and n-fold integration  . There are various definitions of fractional integration and differentiation, such as Grunwald_Letnikov’s definition and Riemann_Liouville’s definition. In the present work, the fractional derivatives are considered in the Caputo sense. The reason for adopting the Caputo definition, as pointed by  , is as follows: to solve differential equations (both classical and fractional), I need to specify additional conditions in order to produce a unique solution. For the case of the Caputo fractional differential equations, these additional conditions are just the traditional conditions, which are akin to those of classical differential equations, and are therefore familiar to us. In contrast, for the Riemann_Liouville fractional differential equations, these additional conditions constitute certain fractional derivatives (and/or integrals) of the unknown solution at the initial point x = 0, which are functions of x. For more details see .
Definition 2.1. The Caputo definition of the fractional-order derivative is defined as
where, is the order of the derivative and n is the smallest integer greater than For the Caputo derivative I have 
(C is a constant), (2)
I use the ceiling function, to denote the smallest integer greater than or equal to , and the floor function to denote the largest integer less than or equal to . Also and . Recall that for , the Caputo differential operator coincides with the usual differential operator of an integer order.
2.2. Properties of Shifted Legendre Polynomials
The well-known Legendre polynomials are defined on the interval and can be determined with the aid of the following recurrence formulae:
where and . In order to use these polynomials on the interval I define the so-called shifted Legendre polynomials by introducing the change of variable.
. Let the shifted Legendre polynomials be denoted by . Then can be obtained as follows:
where and . The analytic form of the shifted Legendre polynomial of degree i given by
Note that and The orthogonality condition is
A function , square integrable in , may be expressed in terms of shifted Legendre polynomials as
where the coefficients are given by
In practice, only the first - terms shifted Legendre polynomials are considered. Then I have
where the shifted Legendre coefficient vector C and the shifted Legendre vector are given by
The derivative of the vector can be expressed by
where is the operational matrix of derivative given by
for example for even m I had
3. Generalized Legendre Operational Matrix to Fractional Calculus
By using Equation (9), it is clear that
where and the superscript, in , denotes matrix powers. Thus
Lemma 1. Let be a shifted Legendre polynomial then
Proof. Using Equations (2), (3) in Equation (6) the lemma can be proved.
In the following theorem I generalize the operational matrix of derivative of shifted Legendre polynomials given in (9) for fractional derivative.
Theorem 1. Let be shifted Legendre vector defined in (8) and also suppose, then
where is the operational matrix of fractional derivative of order in the Caputo sense and is defined as follows:
where is given by
Note that in the first rows, are all zero.
Proof. Using Equations (3), (4) and (6) I have
Now, approximate by terms of shifted Legendre series, I have
Employing Equations (16)-(18) I get
where is given in Equation (15). Rewrite Equation (19) as a vector form I have
Also according to Lemma 1, I can write
A combination of Eqs. (20) and (21) leads to the desired result.
Remark. If Then Theorem 1 gives the same result as Equation (11).
4. Applications of the Operational Matrix of Fractional Derivative
Operational matrix of fractional derivative is applied to solve multi-order fractional differential equation. The existence and uniqueness and continuous dependence of the solution to this problem are discussed in .
4.1. Linear Multi-Order Fractional Differential Equation 
Consider the linear multi-order fractional differential equation
with initial conditions
where are real constant coefficients and also , and denotes the Caputo fractional derivative of order
To solve problem (22) and (23) I approximate and by the shifted Legendre polynomials as
where vector is known but is an unknown vector. By using Equations (13) and (24) I have
Employing Equations (24)-(27) the residual for Equation (22) can be written as
As in a typical tau method  I generate m − n linear equations by applying
Also, by substituting Equations (11) and (24) in Equation (23) I get
Equations (29) and (30) generate m − n and n +1 set of linear equations, respectively. These linear equations can be solved for unknown coefficients of the vector C. Consequently, y(x) given in Equation (24) can be calculated.
4.2. Treatment of Nonhomogeneous Boundary Conditions
To solve Equation (22) with respect to the following boundary conditions (for n is even),
The same technique described in Section 4.1 was applied, but the (n) set of linear equations resulting from (31) is changed to be obtained from
Equations (29) and (31) generate (N + 1) system of linear equations. This system can be solved to determine the unknown coefficients of the vector C.
4.3. Nonlinear Multi-Order Fractional Differential Equation
4.3.1. Consider the Nonlinear Multi-Order Fractional Differential Equation
with initial conditions
where and denotes the Caputo fractional derivative of order . It should be noted that F can be nonlinear in general.
In order to use shifted Legendre polynomials for this problem, I first approximate and for as Equations (24), (26) and (27) respectively. By substituting these equations in Equation (33) I get
Also, by substituting Equations (11) and (24) in Equation (34) I obtain
To find the solution y(x), I first collocate Equation (35) at m − n points. For suitable collocation points I use the first (m-n) shifted Legendre roots of . These equations together with Equation (36) generate a system of (.m + 1) nonlinear equations which can be solved using Newton’s iterative method. Consequently y(x) given in Equation (24) can be calculated.
4.3.2. Boundary Value Problem
Consider the nonlinear FDE (33) with boundary conditions (31). I apply the same technique described in Section 4.3.1, but Equation (36) shall be changed to be (32)., I have a system of (N + 1) nonlinear algebraic equations, which can be solved using Newton’s iterative method.
5. Numerical Results
The presented method in the previous two sections has been applied to solve some examples. In this section, the results for the examples are shown along with their figures and a comparison with the analytical results was also presented in order to show the effectiveness of the technique.
Example 1.  Consider the linear fractional-order IVP:
Subject to the initial conditions
where f(x) is chosen such that the exact solution of (37) is
By applying the technique described in Section 4.1 with m = 3, solution approximated as
Here, I have
Therefore using Equation (29) I obtain
Now, by applying Equation (30) I have
Finally by solving Equations (39)-(42) I got
; ; ;
Thus I can write
which is the exact solution.
Example 2. Consider the boundary value problem
subject to boundary conditions (43)
where the exact solution of this problem is .
This fractional boundary value problem is solved by applying the method described in Section 4.2 by using shifted Legendre expansion and its operational matrices of derivatives with m = 5. Using (43) four linear equations obtained, and by applying boundary condition I had two linear equations. By solving this linear system I got the unknown vector C. By substituting this vector in Equation (43), I had the exact solution.
Then, the 6 unknown coefficients will be in the form
; ; ; ; ;
Therefore, I can write
Numerical results will not be presented since the exact solution is obtained.
Example 3.  Consider the equation
Subject to initial conditions
where the exact solution of this problem is
By applying the technique described in Section 4.1 with m = 3, the approximate solution and the right hand side may be written in the form
Here, I have
Therefore, using Equation (4.8) I obtain
Now, by applying Equation (4.9), I have
Finally, by solving linear system of four Equations (45)-(48), I obtained
which is the exact solution.
It is clear that in Examples 1 - 3 the present method can be considered as an efficient method.
Shifted Legendre approximation method for solving higher order fractional differential equations has been presented. These equations are transformed to a system of algebraic equations to provide a matrix representation. The solution is expressed as a truncated Legendre series, and so it can be easily evaluated for arbitrary values by using computer program. From illustrative examples, it can be seen that this matrix approach can obtain very accurate and satisfactory results. The solution obtained is in very excellent agreement with the already existing ones and shows that this approach can solve the problems effectively. Comparisons between approximate solutions and analytical solutions illustrate the validity and the great potential of the technique.
 Saadatmandi, A. and Dehghan, M. (2006) A Tau Method for the One-Dimensional Parabolic Inverse Problem Subject to Temperature over Specification. Computers & Mathematics with Applications, 52, 933-940.
 Saadatmandi, A. and Dehghan, M. (2008) Numerical Solution of a Mathematical Model for Capillary Formation in Tumor Angiogenesis via the Tau Method. Communications in Numerical Methods in Engineering, 24, 1467-1474.
 Saadatmandi, A. and Dehghan, M. (2007) Numerical Solution of the One-Dimensional Wave Equation with an Integral Condition. Numerical Methods for Partial Differential Equations, 23, 282-292.
 Momani, S. and Shawagfeh, N.T. (2006) Decomposition Method for Solving Fractional Riccati Differential Equations. Applied Mathematics and Computation, 182, 1083-1092.
 Wang, Q. (2006) Numerical Solutions for Fractional KdV-Burgers Equation by Adomian Decomposition Method. Applied Mathematics and Computation, 182, 1048-1055.
 Inc, M. (2008) The Approximate and Exact Solutions of the Space and Time-Fractional Burgers Equations with Initial Conditions by Variational Iteration Method. Journal of Mathematical Analysis and Applications, 45, 476-484.
 Momani, S. and Odibat, Z. (2006) Analytical Approach to Linear Fractional Partial Differential Equations Arising in Fluid Mechanics. Physics Letters A, 355, 271-279.
 Odibat, Z. and Momani, S. (2006) Application of Variational Iteration Method to Nonlinear Differential Equations of Fractional Order. International Journal of Nonlinear Sciences and Numerical Simulation, 7, 271-279.
 Hashim, I., Abdulaziz, O. and Momani, S. (2009) Homotopy Analysis Method for Fractional IVPs. Communications in Nonlinear Science and Numerical Simulation, 14, 674-684.
 Liua, F., Anh, V. and Turner, I. (2004) Numerical Solution of the Space Fractional Fokker-Planck Equation. Journal of Computational and Applied Mathematics, 166, 209-219.
 Esmaeili, S. and Shamsi, M. (2011) A Pseudo-Spectral Scheme for the Approximate Solution of a Family of Fractional Differential Equations. Communications in Nonlinear Science and Numerical Simulation, 16, 3646-3654.
 Mohammadi, F. (2014) Numerical Solution of Bagley-Torvik Equation Using Chebyshev Wavelet Operational Matrix of Fractional Derivative. International Journal of Applied Mathematics and Mechanics, 2, 83-91.
 Momani, S. and Noor, M.A. (2006) Numerical Methods for Fourth-Order Fractional Integro-Differential Equations. Applied Mathematics and Computation, 182, 754-760.
 Diethelm, K., Ford, N.J., Freed, A.D. and Luchko, Yu. (2005) Algorithms for the Fractional Calculus: A Selection of Numerical Methods. Computer Methods in Applied Mechanics and Engineering, 194, 743-773.
 Diethelm, K. and Ford, N.J. (2004) Multi-Order Fractional Differential Equations and Their Numerical Solution. Applied Mathematics and Computation, 154, 621-640.
 Dohaa, E.H., Bhrawyb, A.H. and Ezz-Eldien, S.S. (2011) A Chebyshev Spectral Method Based on Operational Matrix for Initial and Boundary Value Problems of Fractional Order. Computers and Mathematics with Applications, 62, 2364-2373.
 Keshavarz, E., Ordokhani, Y. and Razzaghi, M. (2014) Bernoulli Wavelet Operational Matrix of Fractional Order Integration and Its Applications in Solving the Fractional Order Differential Equations. Applied Mathematical Modelling, 38, 6038-6051. https://doi.org/10.1016/j.apm.2014.04.064
 Ford, N.J. and Connolly, J.A. (2009) Systems-Based Decomposition Schemes for the Approximate Solution of Multi-Term Fractional Differential Equations. Journal of Computational and Applied Mathematics, 229, 382-391.