Fractional spectral differentiation matrices based on Legendre approximation

*Correspondence: aghorbani@um.ac.ir 1Department of Applied Mathematics, Faculty of Mathematical Sciences, Ferdowsi University of Mashhad, Mashhad, Iran Full list of author information is available at the end of the article Abstract A simple scheme is proposed for computing N×N spectral differentiation matrices of fractional order α for the case of Legendre approximation. The algorithm derived here is based upon a homogeneous three-term recurrence relation and is numerically stable. The matrices are then applied to numerically differentiate.


Introduction
In the last few decades, applied scientists have outstretched new models that include fractional derivatives, fractional differential equations (FDEs). The classical and modern dynamical systems modeled by FDEs in physics, engineering, signal processing, fluid mechanics, and bioengineering, manufacturing, systems engineering, and project management can be observed in the recently published book [1]. They indicated that non-integer order derivatives are very appropriate for the explanation of properties of various real materials. Also, it has been exhibited that new order-fractional models are more sufficient than previously applied order-integer ones. In [2][3][4][5][6], the existence and uniqueness of solutions to FDEs have been given. Besides, the authors of [7][8][9][10][11][12][13] have published some new results of fractional operators and their applications. In general, it is impossible to gain analytic solutions of FDEs. Therefore, we have to trust in linearization/discretization and numerical integration to attain approximate solutions, see, e.g., [6,[14][15][16][17]. Moreover, the applications of series solutions for solving some types of FDEs can be observed in [18][19][20].
The differentiation matrices have been extended in the last few decades. It has been demonstrated that they are a useful tool in numerical integrations [21,22]. Those could approximate derivatives by differentiating trial or cardinal basis functions through collocation points. It is noted that these matrices are both dense and very sensitive to rounding errors.
The explicit formulas for the order-integer derivatives of Legendre approximation already exist. Moreover, the authors of [17] proposed a fractional operational matrix based on the shifted Legendre polynomials on the interval [0, 1]. In this article, we present a numerical algorithm for computing Legendre spectral differentiation matrices of orderfractional α (LFDMs) on the interval [0, T]. This algorithm is based on a homogeneous three-term recurrence relation similar to the Legendre recurrence formula and is numerically stable. The proposed LFDMs are then utilized to numerically differentiate.

Preliminaries and definitions
First here we give some basic definitions and properties of the fractional calculus theory, the details of which can be found in [6,23].
The Riemann-Liouville fractional integral operator of order α ≥ 0, of a function f (x) is defined as follows: Properties of the fractional integral operator can be found in [6].
The fractional derivative of f (x) in the Caputo sense is defined as follows: for m -1 < α ≤ m. Also, we require here four of its main properties [6]: x β-α for β ∈ {0, 1, . . .} and β ≥ α , where α denotes the ceiling function, the smallest integer greater than or equal to α, and λ and μ are constants. The standard Legendre polynomials (LPs) are determined by the following three-term recurrence formula: where L 0 (x) = 1 and L 1 (x) = x. Here we define the shifted Legendre polynomials Until further notice x ∈ [0, T] and x = 2 T x -1 (for simplicity of statements). The general form of the shifted LP of order i is given by the sum: where i/2 denotes the floor function, the largest integer less than or equal to i/2. Let Λ = (0, T). The set of shifted LPs is the L 2 (Λ)-orthogonal system in Λ, i.e., where δ ij is the Kronecker symbol. Thus, for any v ∈ L 2 (Λ), we have that Here, the set of all algebraic polynomials of degree at most N , i.e., π N (Λ) = span{P 0 , . . . , P N }, is denoted by π N . We consider the orthogonal projection P N : L 2 (Λ) → π N . For any v ∈ L 2 (Λ), it is defined by [24] (v -P N v, φ) = 0 for any φ ∈ π N (11) or, equivalently, A complete proof on estimating the difference between v and P N v is given by Guo [24]. Now from (12) where the shifted Legendre coefficient vector V and the shifted LP vector Φ(x) are given by [ v 0 , . . . , v N-1 ] and [P 0 (x), . . . , P N-1 (x)], respectively. We can express by dΦ(x) dx = DΦ the derivative of the vector Φ(x), where D is the N × N sparse first order Legendre spectral differentiation matrix. The coefficients of the matrix D are given by

Legendre fractional differentiation matrices
First of all, we require a set of N collocation points x 0 , . . . , x N-1 ∈ [0, T] for computing LFDMs. In this work, we consider the Legendre-Gauss-Lobbato nodes (LGL) x k , k = 0, 1, . . . , N -1. In this case, x 0 = -1, x N-1 = 1, and x k , 1 ≤ k ≤ N -2 are the roots of the first derivative of the standard LP of (N -1) order. Since any interval [a, b] may be scaled into [-1, 1], therefore we can obtain the LGL nodes on the desired interval. Assume that we have the value v = (v(x 0 ), . . . , v(x N-1 )) and we want to approximate derivatives of order α at those points, i.e., v (α) = (v (α) (x 0 ), . . . , v (α) (x N-1 )). One scheme is to use the polynomial interpolation at the grid points x k and then analytically differentiate it and evaluate these derivatives at the same points. The other scheme is to replace polynomials by linear combinations of shifted LPs. In either case, then there is a matrix D α such that The concept of a collocation derivative operator based on trial basis functions is associated with a finite expansion such as where the functions φ j (x) form a complete set of the approximating space π N-1 . Here we consider v : [0, T] → R and φ j (x) = P j (x), j = 0, 1, . . . , N -1, i.e., the shifted LPs on the interval [0, T]. As a result, the function v(x) can be represented in the spectral form Now, using Eqs. (3)-(6) in Eq. (8), we have Also, applying the above-mentioned spectral differentiation matric D, we get First, let us consider α ∈ (0, 1). Substituting (19) into (2), according to (6), we get where In a similar way to [25], substituting N values of x k , k = 0, 1, . . . , N -1, into (20), we can obtain the fractional derivative formulation as follows: where D is as noted in (14) and P is the following matrix of the values of LPs P k , k = 0, 1, . . . , N -1, at the vector x ∈ [0, T]: which can be easily calculated based upon the recurrence formula of the standard LPs.
The remaining problem is the estimation of J k (x i ) for i, k = 0, 1, . . . , N -1. This can be done using a homogeneous linear recurrence relation, as we shall illustrate in the following theorem. Deriving this recurrence relation is based on the properties of LPs. k = 0, 1, . . . , N -1, for x ∈ [0, T] satisfies the following linear, homogeneous, three-term recurrence relation: with starting values where p = 2 T and q = -1.
Proof Let Using the identity we can write (27) in the form On the other hand, by means of the following properties of the Legendre polynomials and by integration by parts of equation (27), we gain Now, equating (29) and (31) yields the desired recurrence relation (24). The starting values can be readily calculated directly.
In view of (22), LFDM in a physical form is given by and also, for a spectral form, we have For the case α > 1, LFDMs will be where r = α and β = αα . Note that the Legendre fractional integration matrices can be derived in an analogous way. So, for α ≥ 0, the fractional integration matrices in physical and spectral forms will respectively be as follows: where the matrix J can be computed from (24)-(26) with α → 1α.

Application of LFDMs
Here the results of two numerical tests are given by using Matlab 7 in double precision. The fractional differentiation of the functions (x + 1) 3 and e x was selected to verify the correctness of LFDMs. That is because the fractional differentiation of those functions is given by [15] D α (x + 1) 3 and which is used to compare the outputs gained by the above-mentioned method. In the above expressions E n 1 ,n 2 and 2 F 1 are respectively the two-parameter Mittag-Leffler function and Gauss's hypergeometric function. When α = 0.5, 1.5 and T = 1, the absolute er- and where α ∈ (0, 1] and J ν (·) is a Bessel function of the first kind. As before, the numerical results of the fractional differentiation for different values of α are reported in Table 1.   15 3.00×10 -12 5.68×10 -12 6.79×10 -12 20 7.77×10 -16 2.55×10 -15 2.02×10 -15

Conclusion
Based on the shifted Legendre polynomials, we proposed an easy procedure for calculating spectral integration/differentiation matrices of the arbitrary order α. The developed algorithm is based on a three-term recurrence relation, which is numerically stable. In order to demonstrate the efficiency of the presented results, four examples of numerical differentiations were given, and the present method performed excellently. Besides, the promising results obtained in the present paper can further be used for solving fractional order problems.