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

Authors: Qinwu Xu; Zhoushun Zheng
Journal: International Journal of Differential Equations
 (2019)

---

## Abstract

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.


---

## Body

## 1. 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–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–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 slow- and 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–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 time-fractional advection-diffusion equations with generalized fractional derivative [19]. Later, a finite difference scheme and an analytical solution are studied for generalized time-fractional 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–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 Section2, 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.


## 2. 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 u(t) in [a,b] is defined as(1)Iatαut≔1Γα∫0tt-sα-1usds.The right fractional integral of order α of u(t) in [a,b] is defined as(2)Itbαut≔1Γα∫tbs-tα-1usds.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 u(t) in [a,b] is defined as(3)Datαut≔dndtnIatn-αut=1Γn-αdndtn∫att-sn-1-αusds.The right Riemann-Liouville derivative of order α of u(t) in [a,b] is defined as(4)Dtbαut≔-ddtnItbn-αut=1Γn-α-ddtn∫tbs-tn-1-αusds,where n=⌈α⌉ is an integer.Definition 3 (Caputo derivative [33]).
The left Caputo derivative of orderα of function u(t) in [a,b] is given by(5)D0ctαut≔I0tn-αdnutdtn=1Γn-α∫0tt-sn-1-αdnusdsnds.The right Caputo derivative of order α function u(t) in [a,b] is given by(6)Dtcbαut≔-1nItbn-αdnutdtn=-1nΓn-α∫tbs-tn-1-αdnusdsnds,where n=⌈α⌉ is an integer.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 powerd/dxα of the differentiation operator d/dx 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 (xd/dx)α. 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.Definition 4 (Hadamard fractional integral [34]).
(7) H a + , 0 α u t ≔ 1 Γ α ∫ a t log ⁡ t s α - 1 u s d s s ,  t > a > 0 .
In 1993, Kilbas studied a weighted Hadamard fractional integral, also called Hadamard-type fractional integral, which extended the application of Hadamard operators.Definition 5 (Hadamard-type fractional integral [34]).
(8) H a + , μ α u t ≔ 1 Γ α ∫ a t s t μ log ⁡ t s α - 1 u s d s s ,  t > a > 0 .
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.Definition 6 (Erdélyi-Kober fractional integral [35]).
(9) I β γ , α u t ≔ β t - β γ + α Γ α ∫ 0 t t β - s β α - 1 s β γ + 1 - 1 u s d s ,  α > 0 ,   β = 1,2 ,   γ ∈ R .
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].Definition 7 (see [36]).
The function spaceCη,η∈R consists of all functions f(x),x>0, that can be represented in the form f(x)=xpf1(x) with p>η and f1∈C([0,∞]).The E-K fractional derivative of orderα is defined in the following form.Definition 8 (Erdélyi-Kober fractional derivative [11, 37]).
(10) D β γ , α u t ≔ ∏ j = 1 n γ + j + 1 β t d d t I β γ + α , n - α u t ,where n=⌈α⌉ is an integer.For the functions from the spaceCη,η≥-β(γ+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(11)Dβγ,αIβγ,αft=ftholds true for every f∈Cη.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 f(t) with respect to another function z(t) and weight w(t) is defined as(12)Ia+;z;wαfx=wx-1Γα∫axwtz′tftzx-zt1-αdt.In this definition, if we setz(s)=s,w(s)=1, it reduces to the classical Riemann-Liouville fractional integral. Similarly, setting z(s)=log⁡(s), w(s)=sμ will lead to Hadamard integral and z(s)=sβ,w(s)=sβγ will lead to E-K fractional integral with a factor tβα.Definition 10 (see [13]).
The left/forward weighted/scaled derivative of integer orderm⩾1 of a function f(t) with respect to another function z(t) and weight w(t) is defined as(13)Dz,w,Lmfx=wx-11z′xDxmwxfxx.Definition 11 (generalized Riemann-Liouville derivative [13]).
The left/forward weighted generalized Riemann-Liouville fractional derivative of orderα>0 of a function f(t) with respect to another function z(t) and weight w(t) is defined as(14)Da+;z,w,1αfx=Dz,w,LmIa+;z;wm-αfx.Definition 12 (generalized Caputo derivative [13]).
The left/forward weighted generalized Caputo fractional derivative of orderα>0 of a function f(t) with respect to another function z(t) and weight w(t) is defined as(15)Da+;z,w,2αfx=Ia+;z;wm-αDz,w,Lmfx.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 Hadamard-type fractional operator.


## 3. Spectral Approximation of Generalized Fractional Operator

In this section, we will first study fractional derivative/integral of Jacobi polynomials and then derive a spectral approximation for generalized fractional operators based on Agrawal’s definitions. Fractional derivatives/integrals of others type can be obtained as special cases.


### 3.1. Fractional Derivative of Orthogonal Polynomials

Denote byJjβ,γ(x) the j-th order Jacobi polynomial with index (β,γ) defined on [-1,1].As a set of orthogonal polynomials,{Jjβ,γ(x)}j=0N satisfies the following three-term-recurrence relation [39]:(16)J0β,γx=1,J1β,γx=β+γ+2x+β-γ,Jj+1β,γx=Ajβ,γx-Bjβ,γJjβ,γx-Cjβ,γJj-1β,γx,1⩽j⩽N-1,where the recursive coefficients are defined as (17)Ajβ,γ=2j+β+γ+12j+β+γ+22j+1j+β+γ+1,Bjβ,γ=γ2-β22j+β+γ+12j+1j+β+γ+12j+β+γ,Cjβ,γ=j+βj+γ2j+β+γ+2j+1j+β+γ+12j+β+γ.In order to derive fractional derivative of Jacobi polynomials, we introduce some useful lemmas first.Lemma 14.
For anyn∈N+,β,γ∈R,α>0, the following relation holds:(18)Jn-1-α,γ+α+1x=-n+γ+11-x2n-αJn-11-α,γ+α+1x-αJn-α,γ+αxn-α;(19)Jnβ+α+1,-1-αx=n+β+11+x2n-αJn-1β+α+1,1-αx-αJnβ+α,-αxn-α.Proof.
According to [39, 4.21], Jacobi polynomials with parameters β,γ∈R are defined by(20)Jnβ,γx=∑j=0n1j!n-j!Γn+j+β+γ+1Γn+β+γ+1Γn+β+1Γj+β+1x-12j.Considering property of Jacobi polynomials,(21)Jnβ,γx=-1nJnγ,β-x.Jacobi polynomials can be rewritten in the following form:(22)Jnβ,γx=∑j=0n-1j+nj!n-j!Γn+j+β+γ+1Γn+β+γ+1Γn+γ+1Γj+γ+1x+12j.
Define the following symbols:(23)PA=n+β+11+xJn-1β+α+1,1-αx;PB=2αJnβ+α,-αx;PC=2n-αJnβ+α+1,-1-αx.
Through a series of calculation, we have(24)PB+PC=∑j=0n2α-1j+nj!n-j!Γn+j+β+1Γn+β+1Γn-α+1Γj-α+1x+12j+∑j=0n2n-α-1j+nj!n-j!Γn+j+β+1Γn+β+1Γn-αΓj-αx+12j=∑j=0n-1j+nj!n-j!Γn+j+β+1Γn+β+1Γn-α+1Γj-α+1x+12j2j.PA=∑j=1n2-1j+nj-1!n-j!Γn+j+β+1Γn+β+1Γn-α+1Γj-α+1x+12j=PB+PC.The equality (19) is proved. Equality (18) can be proved similarly.Lemma 15 (see [39, P.96]).
Forα>0,-1<x<1,β>-1,γ∈R,(25)1-xβ+αJnβ+α,γ-αxJnβ+α,γ-α1=Γβ+α+1Γβ+1Γα∫x11-yβy-x1-αJnβ,γyJnβ,γ1dy.For α>0,-1<x<1,γ>-1,β∈R,(26)1+xγ+αJnβ-α,γ+αxJnβ-α,γ+α1=Γγ+α+1Γγ+1Γα∫-1x1+yγx-y1-αJnβ,γyJnβ,γ1dy.Lemma 16 (see [29, 40]).
Forα>0,-1<x<1,β,γ∈R,(27)Ix1αJn0,γx=1Γα∫x1Jn0,γyy-x1-αdy=n!1-xαΓn+α+1Jnα,γ-αx;(28)I-1xαJnβ,0x=1Γα∫-1xJnβ,0yx-y1-αdy=n!1+xαΓn+α+1Jnβ-α,αx.Lemma 17.
∀ α > 0 , α ∉ N , - 1 < x < 1 , β , γ ∈ R,(29)Dx1αJn0,γx=n!1-x-αΓn+1-αJn-α,γ+αx;(30)D-1xαJnβ,0x=n!1+x-αΓn+1-αJnβ+α,-αx.Proof.
First, for the case0<α<1, formulas (29) and (30) can be obtained by setting β=-α in formula (25), setting γ=-α in formula (26), and applying Riemann-Liouville fractional derivative operator to both sides of them. We refer to [22] for detailed discussion.
For the caseα>1,α∉N, the formulas cannot be obtained from Lemma 15 due to the constrains β>-1 and γ>-1 in Lemma 15. Here we prove the lemma inductively.
For left Riemann-Liouville derivative, from case 1, equality (30) holds for ⌊α⌋=0. Assume the equality holds for ⌊α⌋=k-1,k∈N+, when ⌊α⌋=k, and let α~=α-1; from the assumption, it holds that(31)D-1xα~Jnβ,0x=n!1+x-α~Γn+1-α~Jnβ+α~,-α~x.We apply d/dx to both sides of (31), and then (32)D-1xαJnβ,0x=D-1xα~+1Jnβ,0x=n!Γn+1-α~-α~1+x-α~-1Jnβ+α~,-α~x+1+x-α~ddxJnβ+α~,-α~x=n!1+x-α~-1Γn+1-α~-α~Jnβ+α~,-α~x+n+β+11+x2Jn-1β+α~+1,1-α~x.
Applying Lemma14, we have (33)D-1xαJnβ,0x=n!1+x-α~-1Γn+1-α~-α~n-αJnβ+α~+1,-α~-1x=n!1+x-α~-1Γn+1-α~n-α~Jnβ+α~+1,-α~-1x=n!1+x-αΓn+1-αJnβ+α,-αx.Equality (30) is proved. Equality (29) can be proved similarly.Through the relationship between Caputo derivative and Riemann-Liouville derivative, we immediately obtain the fractional derivative for Caputo derivative.Lemma 18.
Forα>0,-1<x<1, (34)Dxc1αJn0,γx=n!1-x-αΓn+1-αJn-α,γ+αx-∑j=0m-1∂xjJn0,γ11-xj-α-1jΓ1+j-α,D-1cxαJnβ,0x=n!1+x-αΓn+1-αJnβ+α,-αx-∑j=0m-1∂xjJnβ,0-11+xj-αΓ1+j-α,where m=⌈α⌉.Given a functionu(x)∈Hωm[-1,1] and polynomials space PN, the projection of u(x) in space PN, πNu(x), satisfies the following relation: (35)u-πNu,vω=0,∀v∈PN.According properties of space PN and Jacobi polynomials, uN can be expressed as (36)πNux=∑j=0Nu~jJjβ,γx,where u~j=ux,Jβ,γxω/Jβ,γxω2, ω=ω(x) is the weight function.Then, fractional derivative ofu(x) can be approximated as(37)D-1xαux≈∑j=0Nu~jD-1xαJjβ,γx.And we have the following lemma.Lemma 19 (see [41]).
∀ u ∈ H ω m + 1 [ - 1,1 ] ,   0 ≤ l ≤ m , l , m ∈ N , 0 < α < 1, uN is a spectral approximation to u in polynomial space such that (u-πNu,v)ω=0∀v∈PN, and then there exists a constant C such that(38)D-1xl+αux-πNuxωβ+α+l,α+l≤CNl-mD-1xm+αuωβ+α+m,α+m.where ω(β,γ)=(1-x)β(1+x)γ.Combining Lemmas18 and 19 and (37), approximation method can be obtained immediately for Caputo derivative and we have the following corollary.Corollary 20.
∀ u ∈ H ω m + 1 [ - 1,1 ] ,   0 ≤ l ≤ m , l , m ∈ N , 0 < α < 1, uN is a spectral approximation to u in polynomial space such that (u-πNu,v)ω=0∀v∈PN, and then there exists a constant C such that(39)D-1cxl+αux-πNuxω≤CNl-mD-1xm+αuω.Based on above lemmas, fractional integral and derivative of Jacobi polynomials in the standard interval[-1,1] can be expressed explicitly. Fractional integrals and derivatives in shifted interval [a,b] can be obtained through proper variable transforms.For anyy∈[a,b], we assume Jn,∗β,γ(y) is defined in [a,b]. Let h=b-a,x=-1+2y-a/h; then x∈[-1,1]. Substituting x=-1+2y-a/b-a into equation (27)-(28) and (29)-(30), we obtained (40)IybαJn,∗0,γy=h2αIx1αJn0,γx,IayαJn,∗β,0y=h2αI-1xαJnβ,0x,DycbαJn,∗0,γy=2hαDxc1αJn0,γx,DacyαJn,∗β,0y=2hαD-1cxαJnβ,0x.


### 3.2. Spectral Approximation to Generalized Fractional Operators

We assumez(x) and w(x) are positive monotone functions and z,w∈C[a,b]. Obviously, z(x) and w(x) are invertible. The following lemmas can be derived for generalized fractional operators.Lemma 21.
Letg(ζ)=w(z-1(ζ))f(z-1(ζ)), 0<α<1; then the generalized fractional integral operator is equivalent to the following classical fractional integral: (41)Ia+;z;wαfx=wx-1Iζ0ζαgzx.Here ζ0=z(a).Proof.
From the definition of generalized fractional integral, we have(42)Ia+;z;wαfx=wx-1Γα∫axwtz′tftzx-zt1-αdt.Since ζ=z(t) is positive monotone function, then ζ=z(t) is invertible, and we have t=z-1ζ, (43)Ia+;z;wαfx=wx-1Γα∫zazxwz-1ζz′z-1ζfz-1ζzx-ζ1-αdz-1ζ=wx-1Γα∫zazxwz-1ζfz-1ζzxk-ζ1-αdζ.Let g(ζ)=w(z-1(ζ))f(z-1(ζ)); then the generalized fractional integral of f(x) is converted to classical fractional integral of g(ζ) in the following form: (44)Ia+;z;wαfx=wx-1Iζ0ζαgzx.Lemma is proved.Lemma 22.
Letg(ζ)=w(z-1(ζ))f(z-1(ζ)),0<α<1; then the generalized fractional derivative of order α is equivalent to the following classical fractional derivative: (45)Da+;z;w,2αfx=wx-1Dζ0cζαgzx.Here ζ0=z(a).Proof.
Similar to the proof of Lemma21, we have (46)Da+;z;w,2αfx=wx-1Γ1-α∫zazx∂/∂ζwz-1ζfz-1ζ∂ζ/∂tzx-ζαdz-1ζ=wx-1Γ1-α∫zazx∂/∂ζwz-1ζfz-1ζzx-ζαdζ.
Letg(ζ)=w(z-1(ζ))f(z-1(ζ)); then the generalized fractional derivative of f(x) is expressed through the classical fractional derivative of g(ζ), (47)Da+;z;w,2αfx=wx-1Dζ0cζαgzx.Lemma is proved.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 spaceP[w,z]n, such that(48)Pw,zn=vx=wx-1gzx:gx∈Pn,x∈Ω=a,b,where Pn is a polynomial space of up to order n.Define a inner product and norm for spaceP[w,z]n, such that(49)ux,vxω=∫abuxvxωxdx,andvω=v,vω.Define projection QN into space Pw,zN, such that for any function u(x)(50)u-QNu,vω=0,∀v∈Pw,zN.SupposeΦj(x),j=0,…,N, are a set of orthogonal basis functions in space Pz,wN satisfying (51)Φix,Φjxω=1,i=j;0,i≠j.Let N→∞; then Φj(x),j=0,1,…,N,…, form a Lω2(Ω) space, and for any u(x)∈Lω2(Ω), the projection QNu(x) can be written as(52)QNux=∑j=0Nu^jΦjx.Here u^j(j=0,…,N) are expansion coefficients such that (53)u^j=ux,ΦjxωΦjxω2.The weight functionω(x) 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 pj(z(x))=w(x)Φj(x), and note that(54)QNux,Φixω=∫ab∑j=0Nu^jw-1xpjzxw-1xpizxωxdx.Take ω(x)=w2xz′(x); then(55)QNux,Φixω=∫zazb∑j=0Nu^jpjzpizdz.Since pj(z) is polynomial of order j, the computation can be carried out easily through properties of orthogonal polynomials.SupposePj,∗(x) is shifted Legendre polynomial defined on [z(a),z(b)]. Then, given any function u(x)∈Lω2 with weight ω(x)=w2xz′(x), we have(56)QNux=wx-1∑j=0Nu^jPj,∗zx.Recalling Lemmas21 and 22, the generalized fractional integral and derivative can be obtained in the form(57)Ia+;z;wαQNux=wx-1∑j=0Nu^jIζ0ζαPj,∗ζζ=zx;(58)Da+;z;w,2αQNux=wx-1∑j=0Nu^jDζ0cζαPj,∗ζζ=zx.From Lemma19, the following corollary can be obtained immediately.Corollary 24.
ζ ( x ) = z ( x ) is a monotone increasing function, w(x)>0. 0<α<1, QN is a projection into space Pw,zN. u(z-1(ζ))∈Cm(Ω),m∈N; then there exists a constant Cα such that(59)Da+,z,w,2αux-QNuxω2≤CαN1-mDζ0zm+αuζ-1zω1.Here the weight function ω1(ζ)=(ζ-a)-α(b-ζ)α,ω2(x)=ω1zxw2xz′(x).Remark 25.
Unlike the convergence theory in classical polynomial space, the convergence order in spacePw,zN depends on regularity of the function u(x) with respect to z(x). We will illustrate this through numerical examples.Most of the time, it is more convenient to consider the problems in nodal form. We assume the given interpolation points areζj=z(xj)(j=0,1,…,N); then the Lagrange basis functions Lj(x)(j=0,1,…,N) can be defined as follows: (60)Liζ=∏j=0,…,N;j≠iζ-ζjζi-ζj.The functionu(x) can be expressed using both Jacobi polynomials and Lagrange polynomials. The following relation is derived:(61)ux=wx-1∑j=0Ng^jPj,∗zx=wx-1∑j=0NwjujLjzx.Considering the equivalence between Legendre basis and Lagrange basis, the following equality holds:(62)Ljz=∑k=0NlkjPk,∗z,k=0,1,…,N.where lij=(Lj(z),Pi,∗(z))/Pi,∗zL22.Then the nodal form expansion ofu(x) is obtained(63)ux=wx-1∑j=0N∑k=0NwjujlkjPk,∗zx.From (57) and (58), the corresponding nodal form of generalized fractional integral and derivative are obtained(64)Ia+;z;wαux=wx-1∑j=0N∑k=0NwjujlkjIζ0ζαPk,∗ζζ=zx;(65)Da+;z;w,2αux=wx-1∑j=0N∑k=0NwjujlkjDζ0cζαPk,∗ζζ=zx.Example 26.
Now we give an example to show the effectiveness and accuracy of the method. Assumingz(x)=x3,w(x)=1, we consider generalized fractional derivative of y(x) on the interval [0,1] with the following form: (66)yx=∑j=15-1j+1x3j2j!.
The exact generalized fractional derivative ofy(x) is (67)D0+,z,w,2αyx=∑j=15-1j+1x3j-3α/2Γ2j+1-α.
Numerical approximation to generalized fractional derivative ofy(x) can be obtained using (58). We consider the maximum absolute error (68)e=maxx∈0,1⁡D0+,z,w,2αyx-D0+,z,w,2αQNyxof the numerical derivative. Results for α=0.2,0.5,0.8 are shown in Figure 1.
For this example, it is easy to check thaty∈Pw,z10. From theory of spectral approximation, the error would decrease exponentially when N<10 and the numerical fractional derivative would be exact when N≥10. Our numerical result coincides with the theory exactly.Figure 1
Error of numerical approximation to generalized fractional derivative ofy(x).Example 27.
In this example, we test the spectral approximation of Hadamard integral. Considering Hadamard-type fractional integral ofy(x)=sin⁡(πx) on the interval [1,2], when w(x)=1 and α=1, the Hadamard integral of y(x) is Sine integral function, Si(πx)-Si(π); for more general w(x) and α, the exact Hadamard-type integral of y(x) is unknown. Here we consider several pairs of w(x) and α. Numerical Hadamard-type integral of y(x) would be computed using (57) with z(x)=log⁡(x). To evaluate the approximation accuracy, for the case w(x)=1,α=1, the exact Hadamard fractional integral is computed using MATLAB built-in function sinint; for other cases, “exact” Hadamard fractional integral is computed using (57) with large N (e.g., N=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 exampleu is not in the space Pw,zn for any n. Maximum error of approximated fractional integral converges exponentially until reaching machine accuracy.Figure 2
Hadamard fractional integral for different weightw(x) and α.Figure 3
Approximation error of Hadamard fractional integral for different weightw and α.


### 3.3. Fractional Integral/Differential Matrices

Supposeu(x)∈Pz,wN, xi(i=0,1,…,N) are the interpolation points, and yj(j=0,1,…,N) are the collocation points; we define the following symbols:(69)U=ux0,ux1,…,uxNT;(70)Uα=Da+;z;wαuy0,Da+;z;wαuy1,…,Da+;z;wαuyNT,α>0;(71)Uα=Ia+;z;w-αuy0,Ia+;z;w-αuy1,…,Ia+;z;w-αuyNT,α<0.We define a generalized fractional differential/integral matrixMα, such that(72)Mijα=∑k=0NwxjlkjwyiDζ0cζαPk,∗ζζ=zyiforα>0;(73)Mijα=∑k=0NwxjlkjwyiIζ0ζαPk,∗ζζ=zyiforα<0.In order to compute the fractional matrixMα more efficiently, we define a few more matrices. L and Lα are (N+1)×(N+1) matrices such that Lij=Lj(z(yi)), Lijα=Dζ0ζαLj(z)z=z(yi),i,j=0,1,…,N. V is defined based on values of Legendre polynomials at interpolation points z(xi) such that(74)Vij=Pj,∗zxi,i,j=0,1,…,N.Dα is defined as fractional derivative/integral of Legendre polynomials at collocation points ζi=z(yi):(75)Dijα=Dζ0cζαPj,∗ζi,α>0;(76)Dijα=Iζ0ζαPj,∗ζi,α<0.W l and Wr are the weight matrices defined as follows:(77)Wijr=wxj,i=j,0,otherwise;(78)Wijl=1wyj,i=j>1,0,otherwise.Theorem 28.
Foru(x)∈Pz,wN, vectors defined in (69)–(71), and matrices defined in (72)–(78), the following relation holds:(79)Uα=MαU=WlDαV-1WrU.Proof.
From the definition (69)–(71), (72)–(78), and (63)-(64), it is easy to obtain (80)Uα=MαU.Next, we prove that Mα=WlDαV-1Wr.
From the definition of the interpolation functionLj(z), (81)Pj,∗z=∑i=0NPj,∗zxiLiz.Suppose P→(z)=(P0,∗(z),P1,∗(z),…,PN,∗(z)), L→(z)=(L0(z),L1(z),…,LN(z)); for α>0, we have (82)P→z=L→zV,Dζ0ζαP→z=Dζ0ζαL→zV.
Evaluating the matrices multiplication for each element, the following relation is derived:(83)WlDαV-1Wrij=WlLαWrij=wxjwyiLjαzyi(84)=∑k=0NwxjlkjwyiDζ0ζαPk,∗ζζ=zyi=Mijα.The case α<0 can be proved similarly.Remark 29.
Collocation pointsyj and the interpolation points xi 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 z(x).


## 3.1. Fractional Derivative of Orthogonal Polynomials

Denote byJjβ,γ(x) the j-th order Jacobi polynomial with index (β,γ) defined on [-1,1].As a set of orthogonal polynomials,{Jjβ,γ(x)}j=0N satisfies the following three-term-recurrence relation [39]:(16)J0β,γx=1,J1β,γx=β+γ+2x+β-γ,Jj+1β,γx=Ajβ,γx-Bjβ,γJjβ,γx-Cjβ,γJj-1β,γx,1⩽j⩽N-1,where the recursive coefficients are defined as (17)Ajβ,γ=2j+β+γ+12j+β+γ+22j+1j+β+γ+1,Bjβ,γ=γ2-β22j+β+γ+12j+1j+β+γ+12j+β+γ,Cjβ,γ=j+βj+γ2j+β+γ+2j+1j+β+γ+12j+β+γ.In order to derive fractional derivative of Jacobi polynomials, we introduce some useful lemmas first.Lemma 14.
For anyn∈N+,β,γ∈R,α>0, the following relation holds:(18)Jn-1-α,γ+α+1x=-n+γ+11-x2n-αJn-11-α,γ+α+1x-αJn-α,γ+αxn-α;(19)Jnβ+α+1,-1-αx=n+β+11+x2n-αJn-1β+α+1,1-αx-αJnβ+α,-αxn-α.Proof.
According to [39, 4.21], Jacobi polynomials with parameters β,γ∈R are defined by(20)Jnβ,γx=∑j=0n1j!n-j!Γn+j+β+γ+1Γn+β+γ+1Γn+β+1Γj+β+1x-12j.Considering property of Jacobi polynomials,(21)Jnβ,γx=-1nJnγ,β-x.Jacobi polynomials can be rewritten in the following form:(22)Jnβ,γx=∑j=0n-1j+nj!n-j!Γn+j+β+γ+1Γn+β+γ+1Γn+γ+1Γj+γ+1x+12j.
Define the following symbols:(23)PA=n+β+11+xJn-1β+α+1,1-αx;PB=2αJnβ+α,-αx;PC=2n-αJnβ+α+1,-1-αx.
Through a series of calculation, we have(24)PB+PC=∑j=0n2α-1j+nj!n-j!Γn+j+β+1Γn+β+1Γn-α+1Γj-α+1x+12j+∑j=0n2n-α-1j+nj!n-j!Γn+j+β+1Γn+β+1Γn-αΓj-αx+12j=∑j=0n-1j+nj!n-j!Γn+j+β+1Γn+β+1Γn-α+1Γj-α+1x+12j2j.PA=∑j=1n2-1j+nj-1!n-j!Γn+j+β+1Γn+β+1Γn-α+1Γj-α+1x+12j=PB+PC.The equality (19) is proved. Equality (18) can be proved similarly.Lemma 15 (see [39, P.96]).
Forα>0,-1<x<1,β>-1,γ∈R,(25)1-xβ+αJnβ+α,γ-αxJnβ+α,γ-α1=Γβ+α+1Γβ+1Γα∫x11-yβy-x1-αJnβ,γyJnβ,γ1dy.For α>0,-1<x<1,γ>-1,β∈R,(26)1+xγ+αJnβ-α,γ+αxJnβ-α,γ+α1=Γγ+α+1Γγ+1Γα∫-1x1+yγx-y1-αJnβ,γyJnβ,γ1dy.Lemma 16 (see [29, 40]).
Forα>0,-1<x<1,β,γ∈R,(27)Ix1αJn0,γx=1Γα∫x1Jn0,γyy-x1-αdy=n!1-xαΓn+α+1Jnα,γ-αx;(28)I-1xαJnβ,0x=1Γα∫-1xJnβ,0yx-y1-αdy=n!1+xαΓn+α+1Jnβ-α,αx.Lemma 17.
∀ α > 0 , α ∉ N , - 1 < x < 1 , β , γ ∈ R,(29)Dx1αJn0,γx=n!1-x-αΓn+1-αJn-α,γ+αx;(30)D-1xαJnβ,0x=n!1+x-αΓn+1-αJnβ+α,-αx.Proof.
First, for the case0<α<1, formulas (29) and (30) can be obtained by setting β=-α in formula (25), setting γ=-α in formula (26), and applying Riemann-Liouville fractional derivative operator to both sides of them. We refer to [22] for detailed discussion.
For the caseα>1,α∉N, the formulas cannot be obtained from Lemma 15 due to the constrains β>-1 and γ>-1 in Lemma 15. Here we prove the lemma inductively.
For left Riemann-Liouville derivative, from case 1, equality (30) holds for ⌊α⌋=0. Assume the equality holds for ⌊α⌋=k-1,k∈N+, when ⌊α⌋=k, and let α~=α-1; from the assumption, it holds that(31)D-1xα~Jnβ,0x=n!1+x-α~Γn+1-α~Jnβ+α~,-α~x.We apply d/dx to both sides of (31), and then (32)D-1xαJnβ,0x=D-1xα~+1Jnβ,0x=n!Γn+1-α~-α~1+x-α~-1Jnβ+α~,-α~x+1+x-α~ddxJnβ+α~,-α~x=n!1+x-α~-1Γn+1-α~-α~Jnβ+α~,-α~x+n+β+11+x2Jn-1β+α~+1,1-α~x.
Applying Lemma14, we have (33)D-1xαJnβ,0x=n!1+x-α~-1Γn+1-α~-α~n-αJnβ+α~+1,-α~-1x=n!1+x-α~-1Γn+1-α~n-α~Jnβ+α~+1,-α~-1x=n!1+x-αΓn+1-αJnβ+α,-αx.Equality (30) is proved. Equality (29) can be proved similarly.Through the relationship between Caputo derivative and Riemann-Liouville derivative, we immediately obtain the fractional derivative for Caputo derivative.Lemma 18.
Forα>0,-1<x<1, (34)Dxc1αJn0,γx=n!1-x-αΓn+1-αJn-α,γ+αx-∑j=0m-1∂xjJn0,γ11-xj-α-1jΓ1+j-α,D-1cxαJnβ,0x=n!1+x-αΓn+1-αJnβ+α,-αx-∑j=0m-1∂xjJnβ,0-11+xj-αΓ1+j-α,where m=⌈α⌉.Given a functionu(x)∈Hωm[-1,1] and polynomials space PN, the projection of u(x) in space PN, πNu(x), satisfies the following relation: (35)u-πNu,vω=0,∀v∈PN.According properties of space PN and Jacobi polynomials, uN can be expressed as (36)πNux=∑j=0Nu~jJjβ,γx,where u~j=ux,Jβ,γxω/Jβ,γxω2, ω=ω(x) is the weight function.Then, fractional derivative ofu(x) can be approximated as(37)D-1xαux≈∑j=0Nu~jD-1xαJjβ,γx.And we have the following lemma.Lemma 19 (see [41]).
∀ u ∈ H ω m + 1 [ - 1,1 ] ,   0 ≤ l ≤ m , l , m ∈ N , 0 < α < 1, uN is a spectral approximation to u in polynomial space such that (u-πNu,v)ω=0∀v∈PN, and then there exists a constant C such that(38)D-1xl+αux-πNuxωβ+α+l,α+l≤CNl-mD-1xm+αuωβ+α+m,α+m.where ω(β,γ)=(1-x)β(1+x)γ.Combining Lemmas18 and 19 and (37), approximation method can be obtained immediately for Caputo derivative and we have the following corollary.Corollary 20.
∀ u ∈ H ω m + 1 [ - 1,1 ] ,   0 ≤ l ≤ m , l , m ∈ N , 0 < α < 1, uN is a spectral approximation to u in polynomial space such that (u-πNu,v)ω=0∀v∈PN, and then there exists a constant C such that(39)D-1cxl+αux-πNuxω≤CNl-mD-1xm+αuω.Based on above lemmas, fractional integral and derivative of Jacobi polynomials in the standard interval[-1,1] can be expressed explicitly. Fractional integrals and derivatives in shifted interval [a,b] can be obtained through proper variable transforms.For anyy∈[a,b], we assume Jn,∗β,γ(y) is defined in [a,b]. Let h=b-a,x=-1+2y-a/h; then x∈[-1,1]. Substituting x=-1+2y-a/b-a into equation (27)-(28) and (29)-(30), we obtained (40)IybαJn,∗0,γy=h2αIx1αJn0,γx,IayαJn,∗β,0y=h2αI-1xαJnβ,0x,DycbαJn,∗0,γy=2hαDxc1αJn0,γx,DacyαJn,∗β,0y=2hαD-1cxαJnβ,0x.


## 3.2. Spectral Approximation to Generalized Fractional Operators

We assumez(x) and w(x) are positive monotone functions and z,w∈C[a,b]. Obviously, z(x) and w(x) are invertible. The following lemmas can be derived for generalized fractional operators.Lemma 21.
Letg(ζ)=w(z-1(ζ))f(z-1(ζ)), 0<α<1; then the generalized fractional integral operator is equivalent to the following classical fractional integral: (41)Ia+;z;wαfx=wx-1Iζ0ζαgzx.Here ζ0=z(a).Proof.
From the definition of generalized fractional integral, we have(42)Ia+;z;wαfx=wx-1Γα∫axwtz′tftzx-zt1-αdt.Since ζ=z(t) is positive monotone function, then ζ=z(t) is invertible, and we have t=z-1ζ, (43)Ia+;z;wαfx=wx-1Γα∫zazxwz-1ζz′z-1ζfz-1ζzx-ζ1-αdz-1ζ=wx-1Γα∫zazxwz-1ζfz-1ζzxk-ζ1-αdζ.Let g(ζ)=w(z-1(ζ))f(z-1(ζ)); then the generalized fractional integral of f(x) is converted to classical fractional integral of g(ζ) in the following form: (44)Ia+;z;wαfx=wx-1Iζ0ζαgzx.Lemma is proved.Lemma 22.
Letg(ζ)=w(z-1(ζ))f(z-1(ζ)),0<α<1; then the generalized fractional derivative of order α is equivalent to the following classical fractional derivative: (45)Da+;z;w,2αfx=wx-1Dζ0cζαgzx.Here ζ0=z(a).Proof.
Similar to the proof of Lemma21, we have (46)Da+;z;w,2αfx=wx-1Γ1-α∫zazx∂/∂ζwz-1ζfz-1ζ∂ζ/∂tzx-ζαdz-1ζ=wx-1Γ1-α∫zazx∂/∂ζwz-1ζfz-1ζzx-ζαdζ.
Letg(ζ)=w(z-1(ζ))f(z-1(ζ)); then the generalized fractional derivative of f(x) is expressed through the classical fractional derivative of g(ζ), (47)Da+;z;w,2αfx=wx-1Dζ0cζαgzx.Lemma is proved.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 spaceP[w,z]n, such that(48)Pw,zn=vx=wx-1gzx:gx∈Pn,x∈Ω=a,b,where Pn is a polynomial space of up to order n.Define a inner product and norm for spaceP[w,z]n, such that(49)ux,vxω=∫abuxvxωxdx,andvω=v,vω.Define projection QN into space Pw,zN, such that for any function u(x)(50)u-QNu,vω=0,∀v∈Pw,zN.SupposeΦj(x),j=0,…,N, are a set of orthogonal basis functions in space Pz,wN satisfying (51)Φix,Φjxω=1,i=j;0,i≠j.Let N→∞; then Φj(x),j=0,1,…,N,…, form a Lω2(Ω) space, and for any u(x)∈Lω2(Ω), the projection QNu(x) can be written as(52)QNux=∑j=0Nu^jΦjx.Here u^j(j=0,…,N) are expansion coefficients such that (53)u^j=ux,ΦjxωΦjxω2.The weight functionω(x) 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 pj(z(x))=w(x)Φj(x), and note that(54)QNux,Φixω=∫ab∑j=0Nu^jw-1xpjzxw-1xpizxωxdx.Take ω(x)=w2xz′(x); then(55)QNux,Φixω=∫zazb∑j=0Nu^jpjzpizdz.Since pj(z) is polynomial of order j, the computation can be carried out easily through properties of orthogonal polynomials.SupposePj,∗(x) is shifted Legendre polynomial defined on [z(a),z(b)]. Then, given any function u(x)∈Lω2 with weight ω(x)=w2xz′(x), we have(56)QNux=wx-1∑j=0Nu^jPj,∗zx.Recalling Lemmas21 and 22, the generalized fractional integral and derivative can be obtained in the form(57)Ia+;z;wαQNux=wx-1∑j=0Nu^jIζ0ζαPj,∗ζζ=zx;(58)Da+;z;w,2αQNux=wx-1∑j=0Nu^jDζ0cζαPj,∗ζζ=zx.From Lemma19, the following corollary can be obtained immediately.Corollary 24.
ζ ( x ) = z ( x ) is a monotone increasing function, w(x)>0. 0<α<1, QN is a projection into space Pw,zN. u(z-1(ζ))∈Cm(Ω),m∈N; then there exists a constant Cα such that(59)Da+,z,w,2αux-QNuxω2≤CαN1-mDζ0zm+αuζ-1zω1.Here the weight function ω1(ζ)=(ζ-a)-α(b-ζ)α,ω2(x)=ω1zxw2xz′(x).Remark 25.
Unlike the convergence theory in classical polynomial space, the convergence order in spacePw,zN depends on regularity of the function u(x) with respect to z(x). We will illustrate this through numerical examples.Most of the time, it is more convenient to consider the problems in nodal form. We assume the given interpolation points areζj=z(xj)(j=0,1,…,N); then the Lagrange basis functions Lj(x)(j=0,1,…,N) can be defined as follows: (60)Liζ=∏j=0,…,N;j≠iζ-ζjζi-ζj.The functionu(x) can be expressed using both Jacobi polynomials and Lagrange polynomials. The following relation is derived:(61)ux=wx-1∑j=0Ng^jPj,∗zx=wx-1∑j=0NwjujLjzx.Considering the equivalence between Legendre basis and Lagrange basis, the following equality holds:(62)Ljz=∑k=0NlkjPk,∗z,k=0,1,…,N.where lij=(Lj(z),Pi,∗(z))/Pi,∗zL22.Then the nodal form expansion ofu(x) is obtained(63)ux=wx-1∑j=0N∑k=0NwjujlkjPk,∗zx.From (57) and (58), the corresponding nodal form of generalized fractional integral and derivative are obtained(64)Ia+;z;wαux=wx-1∑j=0N∑k=0NwjujlkjIζ0ζαPk,∗ζζ=zx;(65)Da+;z;w,2αux=wx-1∑j=0N∑k=0NwjujlkjDζ0cζαPk,∗ζζ=zx.Example 26.
Now we give an example to show the effectiveness and accuracy of the method. Assumingz(x)=x3,w(x)=1, we consider generalized fractional derivative of y(x) on the interval [0,1] with the following form: (66)yx=∑j=15-1j+1x3j2j!.
The exact generalized fractional derivative ofy(x) is (67)D0+,z,w,2αyx=∑j=15-1j+1x3j-3α/2Γ2j+1-α.
Numerical approximation to generalized fractional derivative ofy(x) can be obtained using (58). We consider the maximum absolute error (68)e=maxx∈0,1⁡D0+,z,w,2αyx-D0+,z,w,2αQNyxof the numerical derivative. Results for α=0.2,0.5,0.8 are shown in Figure 1.
For this example, it is easy to check thaty∈Pw,z10. From theory of spectral approximation, the error would decrease exponentially when N<10 and the numerical fractional derivative would be exact when N≥10. Our numerical result coincides with the theory exactly.Figure 1
Error of numerical approximation to generalized fractional derivative ofy(x).Example 27.
In this example, we test the spectral approximation of Hadamard integral. Considering Hadamard-type fractional integral ofy(x)=sin⁡(πx) on the interval [1,2], when w(x)=1 and α=1, the Hadamard integral of y(x) is Sine integral function, Si(πx)-Si(π); for more general w(x) and α, the exact Hadamard-type integral of y(x) is unknown. Here we consider several pairs of w(x) and α. Numerical Hadamard-type integral of y(x) would be computed using (57) with z(x)=log⁡(x). To evaluate the approximation accuracy, for the case w(x)=1,α=1, the exact Hadamard fractional integral is computed using MATLAB built-in function sinint; for other cases, “exact” Hadamard fractional integral is computed using (57) with large N (e.g., N=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 exampleu is not in the space Pw,zn for any n. Maximum error of approximated fractional integral converges exponentially until reaching machine accuracy.Figure 2
Hadamard fractional integral for different weightw(x) and α.Figure 3
Approximation error of Hadamard fractional integral for different weightw and α.


## 3.3. Fractional Integral/Differential Matrices

Supposeu(x)∈Pz,wN, xi(i=0,1,…,N) are the interpolation points, and yj(j=0,1,…,N) are the collocation points; we define the following symbols:(69)U=ux0,ux1,…,uxNT;(70)Uα=Da+;z;wαuy0,Da+;z;wαuy1,…,Da+;z;wαuyNT,α>0;(71)Uα=Ia+;z;w-αuy0,Ia+;z;w-αuy1,…,Ia+;z;w-αuyNT,α<0.We define a generalized fractional differential/integral matrixMα, such that(72)Mijα=∑k=0NwxjlkjwyiDζ0cζαPk,∗ζζ=zyiforα>0;(73)Mijα=∑k=0NwxjlkjwyiIζ0ζαPk,∗ζζ=zyiforα<0.In order to compute the fractional matrixMα more efficiently, we define a few more matrices. L and Lα are (N+1)×(N+1) matrices such that Lij=Lj(z(yi)), Lijα=Dζ0ζαLj(z)z=z(yi),i,j=0,1,…,N. V is defined based on values of Legendre polynomials at interpolation points z(xi) such that(74)Vij=Pj,∗zxi,i,j=0,1,…,N.Dα is defined as fractional derivative/integral of Legendre polynomials at collocation points ζi=z(yi):(75)Dijα=Dζ0cζαPj,∗ζi,α>0;(76)Dijα=Iζ0ζαPj,∗ζi,α<0.W l and Wr are the weight matrices defined as follows:(77)Wijr=wxj,i=j,0,otherwise;(78)Wijl=1wyj,i=j>1,0,otherwise.Theorem 28.
Foru(x)∈Pz,wN, vectors defined in (69)–(71), and matrices defined in (72)–(78), the following relation holds:(79)Uα=MαU=WlDαV-1WrU.Proof.
From the definition (69)–(71), (72)–(78), and (63)-(64), it is easy to obtain (80)Uα=MαU.Next, we prove that Mα=WlDαV-1Wr.
From the definition of the interpolation functionLj(z), (81)Pj,∗z=∑i=0NPj,∗zxiLiz.Suppose P→(z)=(P0,∗(z),P1,∗(z),…,PN,∗(z)), L→(z)=(L0(z),L1(z),…,LN(z)); for α>0, we have (82)P→z=L→zV,Dζ0ζαP→z=Dζ0ζαL→zV.
Evaluating the matrices multiplication for each element, the following relation is derived:(83)WlDαV-1Wrij=WlLαWrij=wxjwyiLjαzyi(84)=∑k=0NwxjlkjwyiDζ0ζαPk,∗ζζ=zyi=Mijα.The case α<0 can be proved similarly.Remark 29.
Collocation pointsyj and the interpolation points xi 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 z(x).


## 4. Collocation Methods for Fractional Differential and Integral Equations

### 4.1. Fractional Ordinary Differential Equations

In this subsection, we consider collocation method for the generalized fractional ordinary differential equation of the form(85)D0+,z,w,2αux=λxux+fx,x∈Ω=0,b;(86)u0=u0.Here 0<α<1,w(x)>0,z(x)>0 and z(x) is a monotone function in Ω.We assumeuN(x)∈Pz,wN is a numerical solution of the equation, xi(i=0,1,…,N) are the chosen interpolation points, uNj=uN(xj). The following discretized equation is obtained:(87)wx-1∑j=0N∑k=0NwxjuNjlkjDζ0cζαPk,∗zx=λxuNx+fx.Let (87) hold on collocation points yj(j=1,…,N); the matrix form is obtained:(88)Mα-ΛU=F.Here Λ is a diagonal matrix with Λii=λ(yi), F=(u0,f(y1),…,f(yN))T.Considering the initial condition, we setS=Mα-Λ with the first row replaced by (1,0,…,0), F′=F with its first element replaced by u0. Then the solution U is obtained by solving the matrix equation SU=F′.Example 30.
Consider the following example:(89)D0+,z,w,2αux=1+xux+Γ1.5r+2.5xr-2/3αΓ1.5r+2.5-α-xr1+x,x∈0,1;(90)u0=0.Here z(x)=x2/3, w(x)=x, r is an arbitrary positive number.
The exact solution of the ordinary differential equation isu(x)=xr. Maximum absolute errors of numerical solutions for r=6,7 and α=0.3,0.6,0.9 are shown in Figures 4 and 5.
Whenz(x)=x2/3,w(x)=x, the scaled polynomial space Pz,wN becomes (91)Pz,wN=spanx2n/3-1,n=1,2,…,N.
Forr=7 the error converges exponentially and reaches machine accuracy at N=12. It is faster than any finite difference method, while for r=6 solution converges algebraically as N increases; however, it still reaches machine accuracy at N=23. The major reason for this is that u(x)=x7∈Pz,w12 and u(x)=x6∉Pz,wN for any N∈N.Figure 4
Log-log plot of the maximum error forr=6.Figure 5
Semilog plot of the maximum error forr=7.Remark 31.
Since spacePw,zN is transformed from classical polynomial space with respect to z(x), 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 z(t).


### 4.2. Hadamard-Type Integral Equations

We consider the following Hadamard-type boundary value problem:(92)1Γα∫axsxμlog⁡xsα-1fsdss=gx,x∈Ω=a,b.In [42], Kilbas discussed the existence of the solution of (92). Explicit formulas for the solution f(t) were established in the following theorem.Theorem 32 (see [42]).
Ifxμg(x)∈AC[a,b], then the Hadamard-type integral equation (92) with 0<α<1 is solvable in Xμ(a,b) and its solution may be represented in the form(93)fx=x-μΓ1-αaμgalog⁡xa-α+∫axlog⁡xs-αsμgs′ds.Here AC[a,b] is the set of absolutely continuous functions on [a,b] and Xμ(a,b) is space of those Lebesgue measurable functions f on [a,b] for which xμ-1f(x) is absolutely integrable.Solution of the Hadamard-type integral equation is exactly the same as generalized fractional derivative of Riemann-Liouville type withz(s)=log⁡(s),w(s)=sμ. The relationship between (93) and Caputo type generalized fractional derivative is(94)fx=aμgax-μΓ1-αlog⁡xa-α+Da+,z,w,2αgx.Supposexj,j=0,1,…,N, are interpolation points; the discretized solution of equation (92) is(95)fNx=aμgax-μΓ1-αlog⁡xa-α+x-μ∑j=0N∑k=0Nxjμgxjlkj(DζacζαPk,∗log⁡x.Rewriting (95) in matrix form, we have(96)F¯=Ga¯+MαG¯,where Mα is generalized fractional differential matrix; G¯=(g(x0),g(x1),…,g(xN))T, F=(fN(x0),fN(x1),…,fN(xN))T. G¯a is a vector about the initial condition of the integral equation, defined by (97)Ga¯=aμgaΓ1-αx0-μlog⁡x0a-α,x1-μlog⁡x1a-α,…,xN-μlog⁡xNa-αT.Example 33.
Assumeg(x)=sin⁡(x-1),μ=1/3,Ω=[1,10]. Solutions of (92) for α=0.3,0.6,0.9 are shown in Figure 6.Figure 6
Solution of (92).


### 4.3. Erdélyi-Kober Fractional Diffusion Equation

In this subsection, we consider the following Erdélyi-Kober fractional diffusion equation [12]:(98)ux,t=u0x+1Γβ∫0tαβsα/β-1tα/β-sα/ββ-1∂2∂x2ux,sds.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 processesBα,β(t) 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 Bα,β(t) and stretched time-fractional diffusion equation, the following proposition is presented in [12].Proposition 34.
The marginal probability density functionfα,β(x,t) of the process Bα,β(t), t≥0, is the fundamental solution of the stretched time-fractional diffusion equation:(99)ux,t=u0x+1Γβ∫0tαβsα/β-1tα/β-sα/ββ-1∂2∂x2ux,sds.Recalling the definition of generalized fractional integral and settingz(t)=tα/β,w(t)=1, the equation can be rewritten as(100)ux,t=u0x+I0+,z,wβuxxx,t.We use collocation method for both space and time discretization. We choose Legendre-Gauss-Lobatto (L-G-L) pointsxi(i=0,1,…,M) as the space collocation points and choose tj(j=0,1,…,N), such that z(tj) are L-G-L points, as the time collocation points.Define space collocation matrixM2, such that Mij2=d2/dx2Lj(xi), 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 (M+1)(N+1)×(M+1)(N+1) for the second order derivative and fractional integral of order β separately. Then (101)A=M2⊗IN+1,B=IM+1⊗Mβ.SupposeuN(x,t) is the numerical solution of (98), defining uijN=uN(xi,tj), solution vector U, and initial vector U0, such that (102)U=u0,0N,u0,1N,…,u0,NN,u1,0N,…,u1,NN,…,uM,0N,…,uM,NNT,U0=u0x0,…,u0x0,u0x1,…,u0x1,…,u0xM,…,u0xMT,where, in the definition of U0, each u0(xj) is repeated N+1 times.The matrix form discretized equation of (98) is obtained as (103)U=U0+BAU.In the discretized equation, initial condition is explicitly involved. After boundary condition added properly, numerical solution can be obtained by solving the matrix equation.Example 35.
Assumex∈Ω=(-1,1), t∈(0,0.5], u0(x)=e-10x2-e-10, u(·,t)∂Ω=0. Numerical solutions with M=N=50 are shown in Figures 7–10.Figure 7
Standard Brownian motion:α=1,β=1.Figure 8
Time-fractional diffusion withα=0.6,β=0.6.Figure 9
Fractional Brownian motion withα=1.5,β=1.Figure 10
Fractional Brownian motion withα=0.5,β=1.Erdélyi-Kober diffusion equation characterizes the marginal density function of the processBα,β(t),t≥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 Figures7 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].


## 4.1. Fractional Ordinary Differential Equations

In this subsection, we consider collocation method for the generalized fractional ordinary differential equation of the form(85)D0+,z,w,2αux=λxux+fx,x∈Ω=0,b;(86)u0=u0.Here 0<α<1,w(x)>0,z(x)>0 and z(x) is a monotone function in Ω.We assumeuN(x)∈Pz,wN is a numerical solution of the equation, xi(i=0,1,…,N) are the chosen interpolation points, uNj=uN(xj). The following discretized equation is obtained:(87)wx-1∑j=0N∑k=0NwxjuNjlkjDζ0cζαPk,∗zx=λxuNx+fx.Let (87) hold on collocation points yj(j=1,…,N); the matrix form is obtained:(88)Mα-ΛU=F.Here Λ is a diagonal matrix with Λii=λ(yi), F=(u0,f(y1),…,f(yN))T.Considering the initial condition, we setS=Mα-Λ with the first row replaced by (1,0,…,0), F′=F with its first element replaced by u0. Then the solution U is obtained by solving the matrix equation SU=F′.Example 30.
Consider the following example:(89)D0+,z,w,2αux=1+xux+Γ1.5r+2.5xr-2/3αΓ1.5r+2.5-α-xr1+x,x∈0,1;(90)u0=0.Here z(x)=x2/3, w(x)=x, r is an arbitrary positive number.
The exact solution of the ordinary differential equation isu(x)=xr. Maximum absolute errors of numerical solutions for r=6,7 and α=0.3,0.6,0.9 are shown in Figures 4 and 5.
Whenz(x)=x2/3,w(x)=x, the scaled polynomial space Pz,wN becomes (91)Pz,wN=spanx2n/3-1,n=1,2,…,N.
Forr=7 the error converges exponentially and reaches machine accuracy at N=12. It is faster than any finite difference method, while for r=6 solution converges algebraically as N increases; however, it still reaches machine accuracy at N=23. The major reason for this is that u(x)=x7∈Pz,w12 and u(x)=x6∉Pz,wN for any N∈N.Figure 4
Log-log plot of the maximum error forr=6.Figure 5
Semilog plot of the maximum error forr=7.Remark 31.
Since spacePw,zN is transformed from classical polynomial space with respect to z(x), 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 z(t).


## 4.2. Hadamard-Type Integral Equations

We consider the following Hadamard-type boundary value problem:(92)1Γα∫axsxμlog⁡xsα-1fsdss=gx,x∈Ω=a,b.In [42], Kilbas discussed the existence of the solution of (92). Explicit formulas for the solution f(t) were established in the following theorem.Theorem 32 (see [42]).
Ifxμg(x)∈AC[a,b], then the Hadamard-type integral equation (92) with 0<α<1 is solvable in Xμ(a,b) and its solution may be represented in the form(93)fx=x-μΓ1-αaμgalog⁡xa-α+∫axlog⁡xs-αsμgs′ds.Here AC[a,b] is the set of absolutely continuous functions on [a,b] and Xμ(a,b) is space of those Lebesgue measurable functions f on [a,b] for which xμ-1f(x) is absolutely integrable.Solution of the Hadamard-type integral equation is exactly the same as generalized fractional derivative of Riemann-Liouville type withz(s)=log⁡(s),w(s)=sμ. The relationship between (93) and Caputo type generalized fractional derivative is(94)fx=aμgax-μΓ1-αlog⁡xa-α+Da+,z,w,2αgx.Supposexj,j=0,1,…,N, are interpolation points; the discretized solution of equation (92) is(95)fNx=aμgax-μΓ1-αlog⁡xa-α+x-μ∑j=0N∑k=0Nxjμgxjlkj(DζacζαPk,∗log⁡x.Rewriting (95) in matrix form, we have(96)F¯=Ga¯+MαG¯,where Mα is generalized fractional differential matrix; G¯=(g(x0),g(x1),…,g(xN))T, F=(fN(x0),fN(x1),…,fN(xN))T. G¯a is a vector about the initial condition of the integral equation, defined by (97)Ga¯=aμgaΓ1-αx0-μlog⁡x0a-α,x1-μlog⁡x1a-α,…,xN-μlog⁡xNa-αT.Example 33.
Assumeg(x)=sin⁡(x-1),μ=1/3,Ω=[1,10]. Solutions of (92) for α=0.3,0.6,0.9 are shown in Figure 6.Figure 6
Solution of (92).


## 4.3. Erdélyi-Kober Fractional Diffusion Equation

In this subsection, we consider the following Erdélyi-Kober fractional diffusion equation [12]:(98)ux,t=u0x+1Γβ∫0tαβsα/β-1tα/β-sα/ββ-1∂2∂x2ux,sds.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 processesBα,β(t) 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 Bα,β(t) and stretched time-fractional diffusion equation, the following proposition is presented in [12].Proposition 34.
The marginal probability density functionfα,β(x,t) of the process Bα,β(t), t≥0, is the fundamental solution of the stretched time-fractional diffusion equation:(99)ux,t=u0x+1Γβ∫0tαβsα/β-1tα/β-sα/ββ-1∂2∂x2ux,sds.Recalling the definition of generalized fractional integral and settingz(t)=tα/β,w(t)=1, the equation can be rewritten as(100)ux,t=u0x+I0+,z,wβuxxx,t.We use collocation method for both space and time discretization. We choose Legendre-Gauss-Lobatto (L-G-L) pointsxi(i=0,1,…,M) as the space collocation points and choose tj(j=0,1,…,N), such that z(tj) are L-G-L points, as the time collocation points.Define space collocation matrixM2, such that Mij2=d2/dx2Lj(xi), 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 (M+1)(N+1)×(M+1)(N+1) for the second order derivative and fractional integral of order β separately. Then (101)A=M2⊗IN+1,B=IM+1⊗Mβ.SupposeuN(x,t) is the numerical solution of (98), defining uijN=uN(xi,tj), solution vector U, and initial vector U0, such that (102)U=u0,0N,u0,1N,…,u0,NN,u1,0N,…,u1,NN,…,uM,0N,…,uM,NNT,U0=u0x0,…,u0x0,u0x1,…,u0x1,…,u0xM,…,u0xMT,where, in the definition of U0, each u0(xj) is repeated N+1 times.The matrix form discretized equation of (98) is obtained as (103)U=U0+BAU.In the discretized equation, initial condition is explicitly involved. After boundary condition added properly, numerical solution can be obtained by solving the matrix equation.Example 35.
Assumex∈Ω=(-1,1), t∈(0,0.5], u0(x)=e-10x2-e-10, u(·,t)∂Ω=0. Numerical solutions with M=N=50 are shown in Figures 7–10.Figure 7
Standard Brownian motion:α=1,β=1.Figure 8
Time-fractional diffusion withα=0.6,β=0.6.Figure 9
Fractional Brownian motion withα=1.5,β=1.Figure 10
Fractional Brownian motion withα=0.5,β=1.Erdélyi-Kober diffusion equation characterizes the marginal density function of the processBα,β(t),t≥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 Figures7 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].


## 5. 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.



---

*Source: 3734617-2019-01-01.xml*