Construction of Fractional Pseudospectral Differentiation Matrices with Applications

: Differentiation matrices are an important tool in the implementation of the spectral collocation method to solve various types of problems involving differential operators. Fractional differentiation of Jacobi orthogonal polynomials can be expressed explicitly through Jacobi–Jacobi transformations between two indexes. In the current paper, an algorithm is presented to construct a fractional differentiation matrix with a matrix representation for Riemann–Liouville, Caputo and Riesz derivatives, which makes the computation stable and efficient. Applications of the fractional differentiation matrix with the spectral collocation method to various problems, including fractional eigenvalue problems and fractional ordinary and partial differential equations, are presented to show the effectiveness of the presented method.


Introduction
A differentiation matrix plays a key role in implementing the spectral collocation method (also known as the pseudospectral method) to obtain numerical solutions to problems such as partial differential equations.The efficient computation of the differentiation matrix of integer order is discussed in [1][2][3][4][5] and a Matlab suite is introduced in [6].
During the past several decades, the theory and applications of fractional calculus have developed rapidly [7][8][9][10][11].Researchers note that fractional calculus has broad and significant application prospects [8,[12][13][14].Numerous differential equation models that involve fractional calculus operators, referring to fractional differential equations, have emerged in various areas, including physics, chemistry, engineering and even finance and the social sciences.
Publications devoted to numerical solutions to fractional differential equations are numerous.Among them, the spectral collocation method is highlighted as one of the most important numerical methods.Zayernouri et al. [15] studied fractional spectral collocation methods for linear and nonlinear variable-order fractional partial differential equations.Jiao et al. [16] suggested fractional collocation methods using a fractional Birkhoff interpolation basis for fractional initial and boundary value problems.Dabiri et al. [17] presented a modified Chebyshev differentiation matrix and utilized it to solve fractional differential equations.Al-Mdallal et al. [18] presented a fractional-order Legendre collocation method to solve fractional initial value problems.Gholami et al. [19] derived a new pseudospectral integration matrix to solve fractional differential equations.Wu et al. [20] applied fractional differentiation matrices in solving Caputo fractional differential equations.Moreover, the spectral collocation method has been applied to solve tempered fractional differential equations [21][22][23], variable-order Fokker-Planck equations [24] and Caputo-Hadamard fractional differential equations [25].
The spectral collocation method approximates the unknown solution to fractional differential equations with classical orthogonal polynomials and equates the numerical solution at collocation points (usually Gauss-type quadrature nodes).Fractional differentiation in physical space is performed by the differentiation matrix as a discrete version of its continuous fractional derivative operator.Therefore, a crucial step to implement the spectral collocation method is to form the differentiation matrix.One way to construct a fractional differentiation matrix for the spectral collocation method is based on the following formula of Jacobi polynomials (see Equation (3.96) of p. 71, [26]) x − 1 2 k , which suffers from the roundoff error [27,28].Another way is based on the well-known three-term recurrence relationship of Jacobi polynomials (Equation ( 7) in next section), which gives a fast and stable evaluation of a fractional differentiation matrix [29].Similar recurrence relationships are derived to evaluate the differentiation matrix [21,[23][24][25] for various problems.In addition, there is little research on efficiently evaluating the fractional collocation differentiation matrix.
The key point in implementing the spectral collocation method is to stably and efficiently evaluate the collocation differentiation matrix and to solve it.The aim of this work is to present an algorithm to compute the collocation differentiation matrix for fractional integrals and derivatives.In effect, we present a representation of the fractional differentiation matrix that is a product of some special matrices.This representation gives a direct way to form differentiation matrices with relatively fewer operations, which is easy to program.Another benefit of this representation is that the inverses of the differentiation matrices can be obtained in the meantime, which makes the discrete system easy to solve.To show the effectiveness of the algorithm, we apply the fractional differentiation matrix with the Jacobi spectral collocation method to a fractional eigenvalue problem, fractional ordinary differential equations and fractional partial differential equations.In addition, our results provide an alternative option to compute fractional collocation differentiation matrices.We expect that our findings will contribute to further applications of the spectral collocation method to fractional-order problems.
The main contribution of this work is to represent some fractional differentiation matrices as a product of several special matrices.This representation gives not only a direct, fast and stable algorithm of fractional differentiation matrices, but also more information which can be used for inverse or preconditioning.
This paper is organized as follows: We recall several definitions of fractional calculus and Jacobi polynomials with some basic properties in Section 2. We describe the pseudospectral differentiation matrix and its properties in Section 3. We present a representation of the fractional differentiation matrix in Section 4. We also discuss the distribution of the spectrum of the fractional differentiation matrix in this section.We apply the fractional differentiation matrix to fractional eigenvalue problems in Section 5. We apply the fractional differentiation matrix with the Jacobi spectral collocation method to fractional differential equations in Section 6.We finish with the conclusions in Section 7.

Definitions of Fractional Calculus
In this subsection, we recall some definitions of fractional calculi that include the widely used Riemann-Liouville integral, Riemann-Liouville derivative, Caputo derivative and Riesz derivative.We also present some basic properties of these notions.Definition 1.For z > 0, Euler's Gamma function Γ(z) is defined as Definition 2 ([7,8]).For a function f (x) on (a, b) ⊆ R, the µth order left-and right-sided Riemann-Liouville integrals are defined, respectively, as Lemma 1 ([8]).For µ, ν > 0, the Riemann-Liouville integral operator has the following properties: Lemma 2 ([11]).For n > −1, we have Definition 3 ([7,8]).For a function f (x) on (a, b) ⊆ R, the µth order left-and right-sided Riemann-Liouville derivatives are defined as respectively.Here and in the subsequent sections, m is a positive integer satisfying m − 1 < µ < m, m ∈ N.
In Definition 2, the given function f (x) can be absolutely continuous or of bounded variation, or, even more generally, in L p (a, b) with p ≥ 1.As for f (x) in Definitions 3 and 4, it can be in the absolute continuous function space AC m (a, b).In the following discussions, we always assume that f (x) satisfies the aforementioned conditions such that its Riemann-Liouville integrals and Riemann-Liouville and Caputo derivatives exist.

Jacobi Polynomials
Let I := [−1, 1] and P N be the collection of all algebraic polynomials defined on I with degree at most N.The Pochhammer symbol, for c ∈ R and j ∈ N, is defined by Jacobi polynomials P α,β n (s), s ∈ I with parameters α, β ∈ R ( [33]) are defined as and The well-known three-term recurrence relationship of Jacobi polynomials where For α, β > −1, the Jacobi polynomials are orthogonal with respect to the weight where , and δ mn is the Kronecker symbol, i.e., Gauss-type quadratures hold for α, β > −1.Let {x α,β j , ω α,β j } N j=0 be Jacobi-Gauss-Lobatto nodes and weights (the definition and computation of these nodes and weights can be found in [26], pp.83-86).Then, Some additional properties of Jacobi polynomials can be found in [26,33].
then there exists a unique transform matrix T T and the (i, j)-th entry, denoted as t J 1 ,J 2 i,j , of T Proof.The existence and uniqueness of the transform T are clear.From the orthogonality (8) of the Jacobi polynomials, we have Hence, we have It is clear that t J 1 ,J 2 i,j = 0 if i > j.For i ≤ j(i > 0, j > 1), we have Then, we have For i = 0 and j > 1, we have Thus, It is easy to determine that t J 1 ,J 2 0,0 = 1, and Thus, this completes the proof.
Remark 1.The transform matrix T is an upper triangle, and It is clear that = I n is the identity matrix of (n + 1) × (n + 1).Since only the orthogonality of P J 1 n (x) with respect to ω J 1 is used in the derivation, the condition α 1 , β 1 > −1 ensures the transform is valid.
The following results are very useful.

Pseudospectral Differentiation/Integration Matrix
To set up a pseudospectral differentiation/integration matrix, let {x j } N j=0 ⊆ [a, b], and we interpolate the unknown function f (x) at nodes {x j } N j=0 with Lagrange basis functions as where The pseudospectral differentiation/integration matrix is derived by performing an integral or derivative operator on L j (x) at the nodes {x j } N j=0 .In order to impose an initial or boundary value condition, endpoints a and b could be included in {x j } N j=0 .For the purpose of efficient implementation of the pseudospectral method, Gauss-type quadrature nodes, e.g., Gauss-Lobatto or Gauss-Radau points (the definition and computation of these points can be found in [26], pp.80-86), are usually employed.Considering this idea, we take a = −1, b = 1 and x j = x α,β j as the Jacobi-Gauss-Lobatto nodes in increasing order Definition 7. Let F µ be a fractional operator (here, we mean one of the Riemann-Liouville integrals and derivatives, Caputo derivative and Riesz derivatives) with order µ > 0. The discretized matrix F µ corresponding to F µ is an (N + 1) × (N + 1) matrix whose (i, j)-th entry is given by and and Proof. Then, and with some b k .Hence, Let x → −1 + , and one has By considering the second expression of ( 15), the desired equalities on the right-handsided operators can be derived in a similar way.We collect the coefficients of (15) in two matrices H l and H r as Let us introduce three (N + 1)-element vectors as and two (N + 1) × (N + 1) matrices Theorem 3. The following representations of the pseudospectral integration matrix are valid, where diag(•) is a diagonal matrix with the vector in brackets as its diagonal line.Similarly, by replacing −µ with µ, the following representations of the pseudospectral differentiation matrix are found to be valid Proof.From the first expression of (15) and Lemma 2, we have Hence, we have the first equality of (19).The second equality of ( 19) can be obtained in a similar way.Because the two equalities (20) are clear.
In fact, since the Lagrange basis function L j (x) satisfies L j (x i ) = δ ij , we also have the inverse relation as where I is the identity matrix of (N + 1) × (N + 1).Noting that (diag(v , from Theorem 3 and the relation (21), the inverse of the pseudospectral integration/differentiation matrix may be formally written as and Remark 2. Some remarks are listed as follows.

•
From Lemma 1, we expect Similarly, according to Lemma 3, we can expect However, these equalities cannot be verified from Theorem 3 and the above-mentioned equalities.
In fact, the first two equalities are not true according to the numerical test.

•
It is worth noting that

•
We point out that the two matrices B l and B r are Vandermonde's type.

•
From Theorem 2, the first entry of v µ l and the last entry of v µ r are illogical for µ < 0, which indicates the endpoint singularity of the fractional differential operator.In order to avoid this issue, the nodes {x j } N j=0 are altered to the Jacobi-Gauss type.

•
As stated above, the matrices diag(v µ l ) and diag(v µ r ) are singular (or some of their entries are not well defined).Thus, the inverse of the singular matrix in ( 22) and ( 23) should be considered as a generalized inverse or a pseudo-inverse.Theorem 4. From Definition 7, the following relation holds: As stated in Remark 2, since the decompositions in Theorem 3 involve Vandermonde's matrix B l and B r (as well as their inverse L l and L r ), this causes problems in implementing the collocation method.The condition number of B l and B r grows rapidly as N increases.

Representation of Fractional Integration/Differentiation Matrices
In this section, we present a novel representation of the above-mentioned fractional integration/differentiation matrices, which also gives a stable method to compute these matrices and their inverses.Let us consider expanding L j (x) as where From Theorem 1, we have We also have, from [26], , where {x i , ω i } N i=0 denotes the nodes and weights of the Jacobi-Gauss-Lobatto quadrature.Additionally, it is easy to see that

Riemann-Liouville Fractional Integral and Derivative
Here, we have the following representations.
Theorem 5.The following representations of the pseudospectral integration matrices are valid for α, β > −1 Moreover, the differentiation matrices are valid for α, β > −1 Proof.From (24) and Lemma 5, one has Then, (26), the first equality of ( 28) is obtained.The second equality of ( 28) is derived in a similar way.From (24) and Lemma 5, one has Then, The first equality of ( 29) is derived from (26).The second equality of ( 29) is derived in a similar way.
In Theorem 5, the fractional integration and differentiation matrices are expressed as a product of five matrices with size (N + 1) × (N + 1).We emphasize that the two matrices diag(v where Proof.Differentiating ( 24) m times, one has Furthermore, we have .
By taking x = x i (the Jacobi-Gauss-Lobatto nodes as before) for i = 0, 1, • • • , N, we have Combining this with a = T 1.

2.
If µ = k is a positive integer, then m = µ, diag(v 0 l ), diag(v 0 r ) and diag(c 0 ) are all identity matrices (but of size (N − m + 1) × (N − m + 1)).Then, we can obtain the well-known relation of the differentiation matrix of integer order: We investigate the eigenvalues of the fractional differentiation matrices obtained as above, which are closely related to the stability of the numerical method.The interior part of the fractional differentiation matrices is employed with consideration of boundary conditions.Tests are performed for all kinds of the fractional derivatives defined in Section 2.1 with various µ, N, α, β.We report some results as follows.Special emphasis is on the Legendre α = β = 0 and Chebyshev α = β = −0.5 cases.
Since C D µ l is a real matrix, it has complex conjugate pair eigenvalues.Figure 1 shows the distribution of the spectrum of C D µ l for Chebyshev-Gauss-Lobatto points of N = 64 with different µ = 0.1, 0.3, 0.5, 0.7, 0.9, 1.It is observed that the real parts of the eigenvalues are positive for all cases.It is well known that the differentiation matrix C D 1 of first order (µ = 1) with Chebyshev collocation points is semi-positive definite for every N ∈ N (property 5.3 of [17]).Hence, it is reasonable to estimate that C D µ l is semi-positive definite for all 0 < µ ≤ 1, N ∈ N. Similar conclusions are obtained for Legendre-Gauss-Lobatto points, whose eigenvalues are plotted in Figure 2.For fixed µ = 0.5, Figures 3 and 4 demonstrate that, with increasing N, the eigenvalues are more scattered away from the axis in both Chebyshev-Gauss-Lobatto and Legendre-Gauss-Lobatto cases.
Figure 5 shows the distribution of the spectrum of RL D µ l with different µ =1.1:0.2:1.9 and 2 for Legendre-Gauss-Lobatto points of N = 64.It is observed that the real parts of the eigenvalues are negative for all cases.Hence, we estimate that RL D µ l is semi-negative definite for all 1 < µ ≤ 2. For fixed µ = 1.5, Figure 6 demonstrates that, with increasing N, the eigenvalues are also more scattered away from the axis.
It is verified that the Riesz fractional derivative is self-adjoint and positive definite [36].Figure 7 shows the distribution of the spectrum of RZ D µ with different µ = 1.1, 1.3, 1.5, 1.7, 1.9 for Jacobi-Gauss-Lobatto points of α = β = −0.8,N = 64.Figure 8 shows the distribution of the spectrum of RZ D µ for Jacobi-Gauss-Lobatto points of α = −0.7,β= 1.9, N = 64.It is observed that the real parts of the eigenvalues are negative for all cases.

Applications in Fractional Eigenvalue Problems
Eigenvalue problems are of importance in the theory and application of partial differential equations.Duan et al. [37] studied fractional eigenvalue problems.Reutskiy [38] presented a numerical method to solve fractional eigenvalue problems based on external excitation and the backward substitution method.He et al. [39] computed second-order fractional eigenvalue problems with the Jacobi-Davidson method.Gupta and Ranta [40] applied the Legendre wavelet method to solve fractional eigenvalue problems.
For µ = 1.1, 1.3, 1.5, the first six eigenvalues are listed in Tables 5, which are comparable with the results in [40].It is shown that at least 5-6 digits after the decimal point are consistent.The eigenvalues for other cases of µ are also comparable with the results in [40], which are similar, so we do not list them here.
From the results in Tables 1-5, it is shown that our method is reliable for the evaluation of the eigenvalues.
We solve the eigenvalue problem (34) using the Jacobi spectral collocation method.It is known that the Riesz fractional operator is self-adjoint and positive definite in a proper Sobolev space [36].This fact is confirmed by numerical tests.It is observed that the discrete Riesz fractional differentiation matrix is strictly diagonally dominant with a negative diagonal for every value of parameter The same problem has been studied using the Jacobi-Galerkin spectral method [36].The first five eigenvalues obtained with the proposed method are listed in Tables 6 and 7, together with the eigenvalues given in [36] for comparison.From the results in Tables 6 and 7, we conclude that the presented method is reliable in solving eigenvalue problems.The basic fractional initial value problem reads For the well posedness of (35), we refer to [11].For 0 < µ ≤ 1, the discrete system of (35) reads with the collocation point vector x = [x 1 , • • • , x N ] T , the unknown u that approximates the exact solution at the collocation vector, and the differentiation matrix D µ with respect to the fractional derivative C D µ a,x .We obtain the numerical solution u N (x) = u ≈ u(x) by solving the system.
While exact solutions are known in the following examples, the error E ∞ between the exact and numerical solutions is measured with where u N (t) is the numerical solution and t i denotes the mapped collocation points (Jacobi-Gauss-Lobatto points).The numerical convergence order is estimated with where N 1 and N 2 are two different degrees of freedom.We employ CO to measure numerically the convergence order, e.g., the convergence rate likes O(N −CO ).
Example 3. Let 0 < µ < 1.Consider the scalar linear fractional differential equation The right-hand-side term f (t) is chosen so that the exact solution satisfies one of the following cases: For case C11, the source function is exactly f (t) = ∑ 5 k=1 Γ(kσ+1)t kσ−µ kΓ(kσ−µ+1) .Since the solution u has low regularity if σ > 0 is a small non-integer, the convergence order is limited.The errors E ∞ and convergence orders CO are listed in Tables 8 and 9.It is shown that the convergence order for this case is clearly O(N −2σ ) if σ > 0 is a small non-integer.When σ is an integer, the solution u is smooth and exactly belongs to P 5σ .Hence, we only need N = 5σ to resolve this problem.Table 10 lists the errors E ∞ which show that machine precision is almost achieved for the case when N = 5σ.
For case C13, the solution u has low regularity since the term t µ is involved.Thus, the convergence order should be low.Tables 12 and 13 list the errors E ∞ and convergence orders CO, which show clearly that the convergence order is O(N −2µ ).
with the initial value u 1 (0) = u 2 (0) = 0.If µ = 1, the solution of ( 38) is u 1 (t) = 1 − cos(t), u 2 (t) = sin(t).We apply the Legendre spectral collocation method with N = 160 to produce maximum errors of E 1,∞ = 2.3412e − 12 and E 2,∞ = 2.2929e − 12.We plot the curves of the numerical solution u 1 (t) for different µ in Figure 9.We plot the portrait of the phase plane u 1 (t) and u 2 (t) in Figure 10.The results show good agreement with the figures in [17].

Fractional Boundary Value Problems
Let 1 < µ < 2. As a benchmark fractional boundary value problem, we consider the one-dimensional fractional Helmholtz equation as where the fractional derivative operator D µ x will be specified in the following examples.The discrete system of (39) in the spectral collocation method reads with the (N − 1) × (N − 1) identical matrix I, the fractional differentiation matrix D µ (whose first and last row and column are removed) corresponding to the fractional derivative operator, the undetermined unknowns u = u N (x) and the interior collocation points For case C21, since the solution u has low regularity if σ > 0 is a small non-integer, the convergence order is limited.The errors E ∞ are plotted in Figures 11 and 12 for different N, σ and µ. Figure 12 shows that the convergence order is of the first order.It is shown that the convergence order for this case depends on σ and µ.A possible estimate of the convergence order is O(N −2(σ−µ+1) ) if σ > 0 is a small non-integer.
For case C22, the source function is approximated by with a sufficiently large integer L (here, L = 50).Since the solution is smooth, we expect the convergence order to be exponential.Figures 13 and 14 plot the errors E ∞ , which show "spectral accuracy".
We solve the fractional Poisson problems with the Riesz derivative by employing the Legendre spectral collocation method (α = β = 0) with N = 64.The profiles of the numerical solutions are plotted in Figures 15 and 16.From a comparison with the curves in [41] (Figure 1 and Figure 4 (left)), the two groups of curves match very well.

Fractional Initial Boundary Value Problems
The fractional Burgers equation is a generalized model that describes weak nonlinearity propagation, such as in the acoustic phenomenon of a sound wave through a gas-filled tube.There are many publications that present analytical and numerical solutions to the Riemann-Liouville and Caputo fractional Burgers equations (see [27,42]).One of the fractional Burgers equations (FBEs) in one-dimensional form reads with the boundary condition u(a, t) = u(b, t) = 0 and the initial profile u(x, 0) = u 0 (x).
For the time discretization, we employ a semi-implicit time discretization scheme with step size τ, namely the two-step Crank-Nicolson/leapfrog scheme.Then, the full discretization scheme reads: where D is the first-order differentiation matrix.In the following examples, we always take α = β = 0, N = 360 and τ = 10 We first consider the initial profile with two peaks, i.e., case C41.The surface of the numerical solution is plotted in Figure 17 for µ = 1.5 and T = 1.The evolution of the numerical solution is observed.The numerical solutions of the FBE at time t = 1 are plotted in Figure 18 for different values of fractional order µ = 1.1, 1.2, 1.3, 1.5, 1.8.
For case C42, the surface of the numerical solution for µ = 1.5, T = 1 is plotted in Figure 19.The evolution of the numerical solution is observed.The numerical solutions of the FBE at time t = 1 are plotted in Figure 20 for different values of fractional order µ = 1.1, 1.2, 1.3, 1.5, 1.8.x .The numerical test is performed for the same two initial profiles as in (7).
We proceed in the same way as in the above example.We first consider the initial profile with two peaks, i.e., case C41.The surface of the numerical solution is plotted in Figure 21 for µ = 1.5, T = 1.The evolution of the numerical solution is observed.The numerical solutions of the FBE at time t = 1 are plotted in Figure 22 for different values of fractional order µ = 1.1, 1.2, 1.3, 1.5, 1.8.
For case C42, the surface of the numerical solution is plotted in Figure 23 for µ = 1.5, T = 1.The evolution of the numerical solution is observed.The numerical solutions of the FBE at time t = 1 are plotted in Figure 24 for different values of fractional order µ = 1.1, 1.2, 1.3, 1.5, 1.8.

Conclusions
A fractional differentiation matrix is constructed by employing the Jacobi-Jacobi transformation between two indexes (α 1 , β 1 ) and (α 2 , β 2 ).In effect, the fractional differentiation matrix is given as a product of some special matrices.With the aid of this representation, the fractional differentiation matrix can be evaluated in a stable, fast, and efficient manner.This representation gives a direct way to form differentiation matrices with relatively fewer operations, which is easy to program.Another benefit of this representation is that the inverses of the differentiation matrices can be obtained in the meantime, which makes the discrete system easy to solve.We develop applications of the fractional differentiation matrix with the Jacobi spectral collocation method to fractional problems and fractional initial and boundary value problems, as well as fractional partial differential equations.Our numerical experiments involve the Riemann-Liouville, Caputo and Riesz derivatives.All numerical experiments demonstrate that the algorithm is efficient.In addition, our results provide an alternative option to compute fractional collocation differentiation matrices.We expect that our findings can contribute to further applications of the spectral collocation method to fractional-order problems.
The effectiveness of the suggested method needs further exploration.We will investigate the complexity of the method and perform comparisons with other methods in the future.We also expect to apply the method to fractional differential equations in high-dimensional domains.

ThenTheorem 2 .
the left-and right-sided Riemann-Liouville fractional integration matrices by replacing F µ with RL D left-and right-sided fractional differentiation matrices in the sense of Riemann-Liouville and Caputo by replacing F µ with RL D µ −1,x , RL D µ x,1 , C D µ −1,x and C D µ x,1 ; and RZ D µ is the Riesz fractional differentiation matrix by replacing F µ with RZ D µ x , respectively.The first row of RL D −µ l , C D µ l and the last row of RL D −µ r , C D µ r are all zeros; e.g., for j The above theorem states that if 0 < µ < 1, the difference between RL D µ l and C D µ l lies only in the entries of the first column, and the difference between RL D µ r and C D µ r lies only in the entries of the last column.
µ l ) and diag(v µ r ) are not invertible because the first diagonal entry vanishes.In addition, the two matrices diag(v −µ l ) and diag(v −µ r ) are not well defined since the first

Example 5 .
Consider Equation(39) with the Caputo derivative:D µ x = C D µ a,x .The source term f (x) is chosen so that the exact solution satisfies one of the following cases: C21

Table 8 .
The error E ∞ and convergence order CO of C11 in Example 3, computed with the Chebyshev collocation method (α

Table 12 .
The error E ∞ and convergence order CO of C13 in Example 3, computed with the Chebyshev collocation method (α
Figure 10.Portrait of phase plane u 1 (t) and u 2 (t) of Example 4.