Spectral Collocation Method for Fractional Differential/Integral Equations with Generalized Fractional Operator

Generalized fractional operators are generalization of the Riemann-Liouville and Caputo fractional derivatives, which include Erdélyi-Kober and Hadamard operators as their special cases. Due to the complicated form of the kernel and weight function in the convolution, it is even harder to design high order numerical methods for differential equations with generalized fractional operators. In this paper, we first derive analytical formulas for α − th (α > 0) order fractional derivative of Jacobi polynomials. Spectral approximation method is proposed for generalized fractional operators through a variable transform technique. Then, operational matrices for generalized fractional operators are derived and spectral collocation methods are proposed for differential and integral equations with different fractional operators. At last, the method is applied to generalized fractional ordinary differential equation and Hadamard-type integral equations, and exponential convergence of the method is confirmed. Further, based on the proposed method, a kind of generalized grey Brownian motion is simulated and properties of the model are analyzed.


Introduction
During the last decade, fractional calculus has emerged as a model for a broad range of nonclassical phenomena in the applied sciences and engineering [1][2][3][4][5]. Along with the expansion of numerous and even unexpected recent applications of the operators of the classical fractional calculus, the generalized fractional calculus is another powerful tool stimulating the development of this field [6][7][8]. The notion "generalized operator of fractional integration" appeared first in the papers of the jubilarian professor S. L. Kalla in the years 1969-1979 [9,10]. A notable reference is the book by professor Kiryakova [11]. Later, generalized fractional operator was used successfully in papers of Mura and Mainardi [12] to model a class of self-similar stochastic processes with stationary increments, which provided models for both slowand fast-anomalous diffusion. A brief review of generalized fractional calculus is surveyed in [8] and further generalizations of fractional integrals and derivatives are presented by Agrawal in [13]. In the definition of generalized fractional operator, more peculiar kernels are used, which include the classical power kernel function as a special case. Due to the special formulation, all known fractional integrals and derivatives and other generalized integration and differential operators, such as Hadamard operator and the Erdélyi-Kober operator, in various areas of analysis happened to fall in the framework of this generalized fractional calculus. It is shown in [13] that many integral equations can be written and solved in an elegant way using the generalized fractional operator.
Differential and integral equations with generalized fractional operators have been investigated by many authors from theoretical and application aspects [6,11,12,[14][15][16]. In [12,17,18], analytical form solutions of some of these differential and integral equations are given based on transmutation method, Mittag-Leffler function, Mainardi function, and H Fox functions. However, the analytical form solutions are very complicated with infinity serials or integrations, which makes them not suitable for fast computing, while numerical methods are more practical for these equations in applications. In [19], Xu proposes a finite difference scheme for timefractional advection-diffusion equations with generalized 2 International Journal of Differential Equations fractional derivative [19]. Later, a finite difference scheme and an analytical solution are studied for generalized timefractional diffusion equation by Xu and Argrawal [20,21]. These schemes are based on finite difference discretization for first order derivative, which converge with order no more than 2.
Numerical methods with high order convergence have also been developed for fractional differential equations, e.g., spectral method [22], discontinuous Galerkin method [23], and wavelet method [24]. However, they have not yet been applied to fractional differential equations with generalized fractional operator. Among these, spectral method has been confirmed to be efficient and of high accuracy for some fractional differential equations. To name a few, in [25][26][27], Liu and his cooperators proposed a spectral approximation method to fractional derivatives and studied spectral methods for time-fractional Fokker-Planck equations, Riesz space fractional nonlinear reaction-diffusion equations. Li and Xu proposed a spectral method for time-space fractional diffusion equation [22]. Doha studied a spectral collocation method for multiterm fractional differential equations [28]. Zayernouri and Karniadakis proposed a spectral method for fractional ODEs based on polyfractonomials [29]. Zhao and Zhang studied the super-convergence property of spectral method for fractional differential equations [30]. Xu, Hesthaven, and Chen proposed multidomain spectral methods for space and time-fractional differential equations separately [23,31]. Recently, efficient spectral-Galerkin algorithms are developed by Mao and Shen to solve multidimensional fractional elliptic equations with variable coefficients in conserved form as well as nonconserved form [32]. The classical fractional derivative of any power function can be expressed analytically. This indicates that fractional derivative of any function in polynomial spaces can be evaluated exactly. However, the definition of generalized fractional derivative is more complicated than classical fractional derivative. It is far more difficult to construct a direct spectral approximation to the generalized fractional derivative. In this paper, we firstly study fractional derivative of Jacobi polynomials and derive a general formula for fractional derivative of Jacobi polynomial of any order. Through a suitable variable transform technique, spectral approximation formulas are proposed for generalized fractional operators based on Jacobi polynomials. Then, operational matrices are constructed and efficient spectral collocation methods are proposed for the generalized fractional differential and integral equations appearing in the references and applications.
The rest of the current paper is organized as follows: In Section 2, we introduce the definitions of different generalized fractional operators, and some of their properties are also given. In Section 3, a general formula for fractional derivative of Jacobi polynomial of any order is derived. A spectral approximation method for generalized fractional operators is proposed and operational matrices are constructed. In Section 4, collocation methods are proposed for several differential and integral equations. Numerical experiments are carried out to verify the accuracy and efficiency of the methods. Finally, we draw our conclusions in Section 5.

Notations and Definitions
We first introduce some definitions and properties of fractional operators.
Definition 1 (fractional integral [33]). The left fractional integral of order of a given function ( ) in [ , ] is defined as The right fractional integral of order of ( ) in [ , ] is defined as Here Γ(⋅) denotes the Gamma function.
Based on the fractional integral, Riemann-Liouville and Caputo derivatives can be defined in the following way.
Definition 2 (Riemann-Liouville derivative [33]). The left Riemann-Liouville derivative of order of function ( ) in [ , ] is defined as The right Riemann-Liouville derivative of order of ( ) in [ , ] is defined as where = ⌈ ⌉ is an integer.
Definition 3 (Caputo derivative [33]). The left Caputo derivative of order of function ( ) in [ , ] is given by The right Caputo derivative of order function ( ) in [ , ] is given by where = ⌈ ⌉ is an integer.
International Journal of Differential Equations 3 In investigations of dual integral equations in some applications, the modifications of Riemann-Liouville fractional integrals and derivatives are widely used. The important cases include Hadamard fractional operators and Erdélyi-Kober fractional operators.
Riemann-Liouville fractional integro-differentiation is formally a fractional power ( / ) of the differentiation operator / and is invariant relative to translation if considered on the whole axis. Hadamard suggested a construction of fractional integro-differentiation which is a fractional power of the type ( ( / )) . This construction is well suited to the case of the half-axis and invariant relative to dilation [7, §18.3]. Thus Hadamard introduced fractional integrals of the following form.
In 1993, Kilbas studied a weighted Hadamard fractional integral, also called Hadamard-type fractional integral, which extended the application of Hadamard operators.
In investigation of Hankel transform, Erdélyi and Kober proposed the Erdélyi-Kober (E-K) operators. They are generalizations of the classical Riemann-Liouville fractional operators. The left-sided E-K fractional integral of order is defined by the following formula.
In order to introduce the definition of E-K fractional derivative and its properties, we define a special space of functions that was first introduced in [36].
The E-K fractional derivative of order is defined in the following form.
For the functions from the space , ≥ − ( + 1), the left-sided E-K fractional derivative is a left-inverse operator to the left-sided E-K fractional integral (9) [38]; then the relation ( , , holds true for every ∈ . In order to unify these definitions, Agrawal [13] proposed a new definition which includes most of them as special cases. Definition 9 (generalized fractional integral [13]). The left/forward weighted/scaled fractional integral of order > 0 of a function ( ) with respect to another function ( ) and weight ( ) is defined as In this definition, if we set ( ) = , ( ) = 1, it reduces to the classical Riemann-Liouville fractional integral. Similarly, setting ( ) = log( ), ( ) = will lead to Hadamard integral and ( ) = , ( ) = will lead to E-K fractional integral with a factor .
Definition 10 (see [13]). The left/forward weighted/scaled derivative of integer order ⩾ 1 of a function ( ) with respect to another function ( ) and weight ( ) is defined as Definition 11 (generalized Riemann-Liouville derivative [13]). The left/forward weighted generalized Riemann-Liouville fractional derivative of order > 0 of a function ( ) with respect to another function ( ) and weight ( ) is defined as Definition 12 (generalized Caputo derivative [13]). The left/forward weighted generalized Caputo fractional derivative of order > 0 of a function ( ) with respect to another function ( ) and weight ( ) is defined as Remark 13. In the definitions of generalized fractional operators, more general kernels and weight functions are used. It generalized nearly all the existing fractional operators in one space dimension, such as the Riemann-Liouville derivative, the Grünwald-Letnikov derivative, the Caputo derivative, the Erdélyi-Kober-type fractional operator, and the Hadamardtype fractional operator.

Spectral Approximation of Generalized Fractional Operator
In this section, we will first study fractional derivative/integral of Jacobi polynomials and then derive a spectral 4 International Journal of Differential Equations approximation for generalized fractional operators based on Agrawal's definitions. Fractional derivatives/integrals of others type can be obtained as special cases.

Lemma 14.
For any ∈ N + , , ∈ R, > 0, the following relation holds: Proof. According to [39, 4.21], Jacobi polynomials with parameters , ∈ R are defined by Considering property of Jacobi polynomials, Jacobi polynomials can be rewritten in the following form: Define the following symbols: Through a series of calculation, we have The equality (19) is proved. Equality (18) can be proved similarly.
Through the relationship between Caputo derivative and Riemann-Liouville derivative, we immediately obtain the fractional derivative for Caputo derivative.
Then, fractional derivative of ( ) can be approximated as And we have the following lemma.
Combining Lemmas 18 and 19 and (37), approximation method can be obtained immediately for Caputo derivative and we have the following corollary.
Remark 23. These two lemmas establish important relation between classical and generalized fractional operators. With these lemmas, generalized fractional differential equations can be solved via classical fractional differential equations, and vice versa. Which kind of transform should be taken depends on characters of the problem to be solved.
In order to design a high order numerical approximation of the generalized fractional operator, we define a scaled space P [ , ] , such that where P is a polynomial space of up to order .
Define a inner product and norm for space P [ , ] , such that  (51) Let → ∞; then Φ ( ), = 0, 1, . . . , , . . ., form a 2 (Ω) space, and for any ( ) ∈ 2 (Ω), the projection ( ) can be written as Herê( = 0, . . . , ) are expansion coefficients such that The weight function ( ) plays an important role in the computational process and analysis of the method. Here, we choose a proper weight function to use properties of orthogonal polynomials and make the computation more efficient. Let ( ( )) = ( )Φ ( ), and note that The function ( ) can be expressed using both Jacobi polynomials and Lagrange polynomials. The following relation is derived: Considering the equivalence between Legendre basis and Lagrange basis, the following equality holds:

(65)
Example 26. Now we give an example to show the effectiveness and accuracy of the method. Assuming ( ) = √ 3 , ( ) = 1, we consider generalized fractional derivative of ( ) on the interval [0, 1] with the following form: The exact generalized fractional derivative of ( ) is Numerical approximation to generalized fractional derivative of ( ) can be obtained using (58). We consider the maximum absolute error  Figure 1. For this example, it is easy to check that ∈ P 10 , . From theory of spectral approximation, the error would decrease exponentially when < 10 and the numerical fractional derivative would be exact when ≥ 10. Our numerical result coincides with the theory exactly.
Example 27. In this example, we test the spectral approximation of Hadamard integral. Considering Hadamard-type fractional integral of ( ) = sin( ) on the interval [1,2], when ( ) = 1 and = 1, the Hadamard integral of ( ) is Sine integral function, Si( ) − Si( ); for more general ( ) and , the exact Hadamard-type integral of ( ) is unknown. Here we consider several pairs of ( ) and . Numerical Hadamard-type integral of ( ) would be computed using (57) with ( ) = log( ). To evaluate the approximation accuracy, for the case ( ) = 1, = 1, the exact Hadamard fractional integral is computed using MATLAB built-in function ; for other cases, "exact" Hadamard fractional integral is computed using (57) with large (e.g., = 50), which is treated as reference solution. The results for numerical Hadamard integral and approximation error are shown in Figures 2 and 3.
For this example is not in the space P , for any . Maximum error of approximated fractional integral converges exponentially until reaching machine accuracy.
In order to compute the fractional matrix M more efficiently, we define a few more matrices. L and L are ( + 1) × ( + 1) matrices such that L = L ( ( )), L = Evaluating the matrices multiplication for each element, the following relation is derived: The case < 0 can be proved similarly.

Remark 29. Collocation points and the interpolation points
are not necessary the same. To obtain a good approximation, interpolation points of Gauss-type are usually used. At the same time, collocation points should be chosen properly to guarantee stability properties of the method. In the following numerical examples, for computation and stability aim, both interpolation and collocation points are chosen based on Gauss-type points with respect to ( ).
Example 30. Consider the following example: When ( ) = 2/3 , ( ) = , the scaled polynomial space P , becomes For = 7 the error converges exponentially and reaches machine accuracy at = 12. It is faster than any finite difference method, while for = 6 solution converges algebraically as increases; however, it still reaches machine accuracy at = 23. The major reason for this is that ( ) = 7 ∈ P 12 , and ( ) = 6 ∉ P , for any ∈ N.  Remark 31. Since space P , is transformed from classical polynomial space with respect to ( ), the convergence of spectral collocation method for ODEs with generalized fractional operators depends not only on the smoothness of the solution itself, but also on the scale function ( ).

Hadamard-Type Integral Equations.
We consider the following Hadamard-type boundary value problem: International Journal of Differential Equations 11 In [42], Kilbas discussed the existence of the solution of (92). Explicit formulas for the solution ( ) were established in the following theorem.

Erdélyi-Kober Fractional Diffusion Equation.
In this subsection, we consider the following Erdélyi-Kober fractional diffusion equation [12]: Erdélyi-Kober fractional diffusion equation, which is also called stretched time-fractional diffusion equation, is the master equation of a kind of generalized grey Brownian motion (ggBm). The ggBm is a parametric class of stochastic processes that provides models for both fast and slow anomalous diffusion. This class is made up of self-similar processes , ( ) with stationary increments and it depends on two real parameters: 0 < ≤ 2 and 0 < ≤ 1. It includes the fractional Brownian motion when 0 < ≤ 2 and = 1, the time-fractional diffusion stochastic processes when 0 < = < 1, and the standard Brownian motion when = = 1. About the relationship between stochastic process , ( ) and stretched time-fractional diffusion equation, the following proposition is presented in [12].  Recalling the definition of generalized fractional integral and setting ( ) = / , ( ) = 1, the equation can be rewritten as ( , ) = 0 ( ) + 0+, [ , ] ( , ) .
Define space collocation matrix M 2 , such that M 2 = ( 2 / 2 )L ( ), and generalized fractional integral matrix, M . Matrix M is computed through Theorem 28 and the space-time collocation matrices are obtained using Kronecker product ⊗ . Suppose A and B are space-time collocation matrices with dimension ( + 1)( + 1) × ( + 1)( + 1) for the second order derivative and fractional integral of order separately. Then where, in the definition of 0 , each 0 ( ) is repeated + 1 times.
The matrix form discretized equation of (98) is obtained as In the discretized equation, initial condition is explicitly involved. After boundary condition added properly, numerical solution can be obtained by solving the matrix equation. Erdélyi-Kober diffusion equation characterizes the marginal density function of the process , ( ), ≥ 0. When = = 1, we recover the standard diffusion equation. When 0 < = < 1, we get the time-fractional diffusion equation of order . When = 1 and 0 < < 2, we have the equation of the fractional Brownian motion marginal density.
As shown in Figures 7 and 8, when 1 < < 2, the diffusion is fast, and the increments exhibit long-range dependence; when 0 < < 1, the diffusion is slow, and the increments form a stationary process which does not exhibit long-range dependence. The results coincide with theoretical analysis in [12,14].

Conclusion
In this paper, we propose a spectral collocation method for differential and integral equations with generalized fractional operators. To deal with the difficulty in designing spectral approximation scheme due to complexity of integral kernel and weight, a variable transform technique is applied to the generalized fractional operator, and a spectral approximation method is proposed for the generalized fractional operator. Operational matrices for generalized fractional operators are derived. Spectral collocation methods are designed for fractional ordinary differential equations, Hadamard-type integral equations, and Erdélyi-Kober diffusion equations separately. Numerical experiments are carried out to verify the accuracy and efficiency of the method and characteristics of the Erdélyi-Kober diffusion equation are analyzed based on numerical results.

Data Availability
(i) The programs used to support the findings of this study have been deposited in the GitHub repository (https://github .com/qinwuxu/SpectralGFPDE ). (ii) No data were used to support this study.

Disclosure
An earlier version of this work was presented at the "8th International Congress on Industrial and Applied Mathematics (ICIAM 2015)".

Conflicts of Interest
The authors declare that they have no conflicts of interest.