# High-Order Algorithms for Riesz Derivative and Their Applications(I)

Authors: Hengfei Ding; Changpin Li; YangQuan Chen
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We firstly develop the high-order numerical algorithms for
the left and right Riemann-Liouville derivatives. Using these derived schemes,
we can get high-order algorithms for the Riesz fractional derivative. Based on
the approximate algorithm, we construct the numerical scheme for the space
Riesz fractional diffusion equation, where a fourth-order scheme is proposed
for the spacial Riesz derivative, and where a compact difference scheme is
applied to approximating the first-order time derivative. It is shown that the
difference scheme is unconditionally stable and convergent. Finally, numerical
examples are provided which are in line with the theoretical analysis.


---

## Body

## 1. Introduction

In recent decades, fractional differential equations are undergoing rapid development due to their wide application in physics, engineering, economics, and some other research realms [1–3]. Several methods have been introduced for some special fractional differential equations to seek their analytical solutions, such as the integral transformation method (Laplace's transform, Fourier's transform, and Mellin's transform), Adomian decomposition method, the method of separating variables, and so on [4–7]. However, the exact solutions of most fractional differential equations cannot be obtained. So it becomes important to develop numerical methods for these equations. It is worth mentioning that the fractional linear multistep (also high order) methods for the Riemann-Liouville integrals and derivatives were firstly proposed in [8]. High-order numerical methods for Caputo derivatives were firstly constructed in [9]. Other fractional differential equations include Caputo or (and) Riemann-Liouville derivatives; there also existed some numerical approaches [10–13]. So, in this paper, our main aim is to construct high-order numerical methods for Riesz fractional derivative and Riesz fractional differential equations.The Riesz fractional derivative was derived from the kinetics of chaotic dynamics [14, 15]. For the Riesz fractional differential equations, there have existed several analytical and numerical methods. Zhang and Liu [16] studied the analytical solutions of space Riesz and time Caputo fractional partial differential equations. Later on, they [17] used Galerkin finite element approximation and a backward difference technique in Riesz fractional advection-dispersion equation. By using the Laplace and Fourier transform, Chen et al. [18] obtained the analytical solution for the space Riesz fractional reaction dispersion equation, and they also constructed an explicit finite difference scheme for it. Shen et al. [19] presented the explicit and implicit difference schemes for the Riesz fractional advection-dispersion equation and the space Riesz time Caputo fractional advection-dispersion equation. Recently, they [20] proposed a novel numerical method for the space Riesz fractional advection-dispersion equation based on fractional central difference operator. Özdemir et al. [21] were concerned with the numerical solutions of a two-dimensional problem. Yang et al. [22] obtained some numerical solutions for the space Riesz fractional diffusion and advection-dispersion equations on a finite domain by using three numerical methods. Çelik and Duman [23] numerically approximated the fractional diffusion equation with the space Riesz fractional derivative on a finite domain according to the fractional central difference operator. Very recently, Wang et al. [24] consider the coupled nonlinear Schrödinger equations by using similar methods as in [23]. In the above mentioned numerical methods, the accuracy of the numerical methods for the Riesz fractional derivative is not more than 2; this motivates us to investigative high-order numerical schemes for it.In the present paper, we study the following space Riesz fractional diffusion equation:(1)∂u(x,t)∂t=K∂αu(x,t)∂|x|α+f(x,t),mmmmmma<x<b,0<t≤T,
with the corresponding initial and boundary value conditions:
(2)u(x,0)=ψ(x),a≤x≤b,RLDa,xα-2u(x,t)|x=a=ϕ(t),0≤t≤T,RLDx,bα-2u(x,t)|x=b=φ(t),0≤t≤T,
where K>0 is a diffusion coefficient and f(x,t), ϕ(t), φ(t), and ψ(x) are sufficiently smooth functions. ∂α/∂|x|α is the Riesz fractional derivative with order α(1<α<2) defined as follows [22]:
(3)∂αu(x,t)∂|x|α=-Ψα(RLDa,xα+RLDx,bα)u(x,t),
where Ψα=(1/2)sec(πα/2) and RLDa,xα and RLDx,bα are the left- and right-side Riemann-Liouville fractional derivatives represented in the following form:
(4)RLDa,xαu(x,t)=1Γ(2-α)∂2∂x2∫ax‍u(ξ,t)(x-ξ)α-1dξ,mmmmmmmmmmmmmmmmmmx∈[a,b],RLDx,bαu(x,t)=1Γ(2-α)∂2∂x2∫xb‍u(ξ,t)(ξ-x)α-1dξ,mmmmmmmmmmmmmmmmmmx∈[a,b];RLDa,xα-2 and RLDx,bα-2 are left- and right-side Riemann-Liouville fractional integral operators and defined by (1<α<2):
(5)RLDa,xα-2u(x,t)=1Γ(2-α)∫ax‍(x-ξ)1-αu(ξ,t)dξ,RLDx,bα-2u(x,t)=1Γ(2-α)∫xb‍(ξ-x)1-αu(ξ,t)dξ,
in which Γ(·) indicates Euler's Gamma function.The structure of the paper is outlined as follows. In Section2, we list existing numerical approximation of the Riesz fractional derivative just for reference. Then, we focus on constructing one second-order scheme and two fourth-order numerical schemes. In Section 3, we only present one fourth-order numerical algorithm for the space Riesz fractional diffusion equation on a finite domain. In Sections 4 and 5, we discuss the stability and convergence of the difference scheme, respectively. In Section 6, some numerical examples are shown which support the theoretical results derived in the above sections. Finally, the last section concludes this paper.


## 2. The Numerical Approximation of the Riesz Derivative

Firstly, we take the mesh pointsxm=a+mh, m=0,1,…,M, and tn=nτ, n=0,1,…,N, where h=(b-a)/M and τ=T/N; that is, h and τ are the uniform spatial stepsize and temporal stepsize.As we all know, for everyα(1<α<2) the left and right Riemann-Liouville derivatives exist and coincide with the left and right Grünwald-Letnikov derivatives under suitable conditions, respectively, where the left and right Grünwald-Letnikov derivative definitions with order α are given below [6]:
(6)GLDa,xαu(xm,t)=1hα∑k=0m‍ϖk(α)u(xm-k,t)+O(h),GLDx,bαu(xm,t)=1hα∑k=0M-m‍ϖk(α)u(xm+k,t)+O(h),
in which ϖk(α)=(-1)k(αk)=((-1)kΓ(1+α))/(Γ(1+k)Γ(1+α-k)).So, the Riesz derivative with orderα∈(1,2) can be discretized in the following ways.(I) By the Standard Grünwald-Letnikov Formula. Based on the above assumption and (3), we can obtain the first-order approximation formula
(7)∂αu(xm,t)∂|x|α=-Ψαhα(∑k=0m‍ϖk(α)u(xm-k,t)mmmmmmmmmmm+∑k=0M-m‍ϖk(α)u(xm+k,t))+O(h).(II) By the Shifted Grünwald-Letnikov Formula. In [25], Meerschaert and Tadjeran show that the above standard Grünwald-Letnikov formula is often unstable for time dependent problems. Hence, they proposed the following shifted Grünwald-Letnikov formulas for the left and right Riemann-Liouville derivatives in order to overcome the instability:
(8)RLDa,xαu(xm,t)=∑k=0m+1‍ϖk(α)u(xm-k+1,t)+O(h),RLDx,bαu(xm,t)=∑k=0M-m+1‍ϖk(α)u(xm+k-1,t)+O(h).Therefore, the modified first-order approximation scheme is constructed as follows:(9)∂αu(xm,t)∂|x|α=-Ψαhα(∑k=0m+1‍ϖk(α)u(xm-k+1,t)mmmmmmmmmmm+∑k=0M-m+1‍ϖk(α)u(xm+k-1,t))+O(h).(III) By the L2 Approximation Method. Note that the left and right Riemann-Liouville derivatives can be rewritten as (1<α<2):
(10)RLDa,xαu(x,t)=∑k=01‍xk-αΓ(k+1-α)∂ku(a,t)∂xk+1Γ(2-α)∫ax‍∂2u(ξ,t)∂ξ2(x-ξ)1-αdξ,RLDx,bαu(x,t)=∑k=01‍(b-x)k-αΓ(k+1-α)∂ku(b,t)∂xk+1Γ(2-α)∫xb‍∂2u(ξ,t)∂ξ2(ξ-x)1-αdξ.Hence, we can obtain a first-order scheme for the left and right Riemann-Liouville fractional derivatives [22]:
(11)RLDa,xαu(xm,t)=1Γ(3-α)hα{(1-α)(2-α)u(x0,t)mαmmmmmmmmm+(2-α)[u(x1,t)-u(x0,t)]mα-1mmmmmmmmm+∑k=0m-1‍dk(α)[u(xm-k+1,t)mmmmmmmmmmmmmm-2u(xm-k,t)mmmmmmmmmmmmmm+u(xm-k-1,t)](1-α)(2-α)u(x0,t)mα}+O(h),(12)RLDx,bαu(xm,t)=1Γ(3-α)hα{(1-α)(2-α)u(xM,t)(M-m)αmmmmmmmmm+(2-α)[u(xM,t)-u(xM-1,t)](M-m)α-1mmmmmmmmm+∑k=0M-m-1‍dk(α)[u(xm+k-1,t)mmmmmmmmmmmmmmmm-2u(xm+k,t)mmmmmmmmmmmmmmmm+u(xm+k+1,t)](1-α)(2-α)u(xM,t)(M-m)α}+O(h),
where dk(α)=(k+1)2-α-k2-α, k=0,1,…,m-1, or k=0,1,…,M-m-1.Therefore, applying the above two formulas and (3) gives
(13)∂αu(xm,t)∂|x|α=-ΨαΓ(3-α)hα×{(1-α)(2-α)u(x0,t)mαmmmmmmmm+(2-α)[u(x1,t)-u(x0,t)]mα-1mmmmmmmm+∑k=0m-1‍dk(α)[u(xm-k+1,t)-2u(xm-k,t)mmmmmmmmmmmmmm+u(xm-k-1,t)]mmmmmmmm+(1-α)(2-α)u(xM,t)(M-m)αmmmmmmmm+(2-α)[u(xM,t)-u(xM-1,t)](M-m)α-1mmmmmmmm+∑k=0M-m-1‍dk(α)[u(xm+k-1,t)-2u(xm+k,t)mmmmmmmmmmmmmmmm+u(xm+k+1,t)](1-α)(2-α)u(x0,t)mα}+O(h),
in which dk(α) is defined as above.(IV) By the Spline Interpolation Method. In [26], Sousa proposed a second-order scheme by linear spline interpolation method for the left and right Riemann-Liouville derivatives,
(14)RLDa,xαu(xm,t)=1Γ(4-α)hα∑k=0m+1‍z-m,k(α)u(xk,t)+O(h2),
where
(15)z-m,k(α)={c-m-1,k-2c-m,k+c-m+1,k,k≤m-1,-2c-m,k+c-m+1,k,k=m,c-m+1,k,k=m+1,0,k>m+1,
in which
(16)c-j,k={(j-1)3-α-j2-α(j-3+α),k=0,(j-k+1)3-α-2(j-k)3-α+(j-k-1)3-α,1≤k≤j-1,1,k=j,RLDx,bαu(xm,t)=1Γ(4-α)hα∑k=m-1M‍z~m,k(α)u(xk,t)+O(h2),
where
(17)z~m,k(α)={0,k<m-1,c~m-1,m-1,k=m-1,-2c~m,m+c~m-1,m,k=m,c~m-1,k-2c~m,k+c~m+1,k,m+1≤k≤M,
in which
(18)c~j,k={1,k=j,(k-j+1)3-α-2(k-j)3-α+(k-j-1)3-α,j+1≤k≤M-1,(3-α-M+j)(M-j)2-α+(M-j-1)3-α,k=M,
with j=m-1,m,m+1.Combining (14), (16), and (3) gives
(19)∂αu(xm,t)∂|x|α=-ΨαΓ(4-α)hα∑k=0M‍zm,k(α)u(xk,t)+O(h2),
where
(20)zm,k(α)={z-m,k(α),k<m-1,z-m,m-1(α)+z~m,m-1(α),k=m-1,z-m,m(α)+z~m,m(α),k=m,z-m,m+1(α)+z~m,m+1(α),k=m+1,z~m,k(α),k>m+1.(V) By the Fractional Centered Difference Method. In [27], Ortigueira introduced a symmetrical fractional centered difference operator as follows:
(21)Δhαu(x,t)=∑k=-∞∞‍(-1)kΓ(α+1)Γ((α/2)-k+1)Γ((α/2)+k+1)×u(x-kh,t).Later on, Çelik and Duman [23] proved that the above symmetrical fractional centered difference operator for the Riesz fractional derivative has the following relation:
(22)∂αu(xm,t)∂|x|α=-1hαΔhαu(xm,t)+O(h2).(VI) By the Weighted and Shifted Grünwald-Lentikov Formulas. In [28], Tian et al. proposed the second-order and third order numerical schemes for the left and right Riemann-Liouville derivatives:
(23)RLDa,xαu(xm,t)=ν1hα∑k=0m+l1‍ϖk(α)u(xm-k+l1,t)+ν2hα∑k=0m+l2‍ϖk(α)u(xm-k+l2,t)+O(h2),(24)RLDx,bαu(xm,t)=ν1hα∑k=0M-m+l1‍ϖk(α)u(xm+k-l1,t)+ν2hα∑k=0M-m+l2‍ϖk(α)u(xm+k-l2,t)+O(h2),
where l1 and l2 are two arbitrary integers and l1-l2≠0, ν1=(α-2l2)/2(l1-l2), and ν2=(2l1-α)/2(l1-l2).And(25)RLDa,xαu(xm,t)=κ1hα∑k=0m+l1‍ϖk(α)u(xm-k+l1,t)+κ2hα∑k=0m+l2‍ϖk(α)u(xm-k+l2,t)+κ3hα∑k=0m+l3‍ϖk(α)u(xm-k+l3,t)+O(h3),RLDx,bαu(xm,t)=κ1hα∑k=0M-m+l1‍ϖk(α)u(xm+k-l1,t)+κ2hα∑k=0M-m+l2‍ϖk(α)u(xm+k-l2,t)+κ3hα∑k=0M-m+l3‍ϖk(α)u(xm+k-l3,t)+O(h3),
in which l1, l2, and l3 are three arbitrary integers and (l1-l2)(l2-l3)(l1-l3)≠0, κ1=(12l2l3-(6l2+6l3+1)α+3α2)/12(l2l3-l1l2-l1l3+l12), κ2=(12l1l3-(6l1+6l3+1)α+3α2)/12(l1l3-l1l2-l2l3+l22), and κ3=(12l1l2-(6l1+6l2+1)α+3α2)/12(l1l2-l1l3-l2l3+l32).Naturally, we can obtain the following second-order and third order numerical formulas for the Riesz fractional derivative:(26)∂αu(xm,t)∂|x|α=-Ψαhα(ν1∑k=0m+l1‍ϖk(α)u(xm-k+l1,t)mmmmm+ν2∑k=0m+l2‍ϖk(α)u(xm-k+l2,t)mmmmm+ν1∑k=0M-m+l1‍ϖk(α)u(xm+k-l1,t)mmmmm+ν2∑k=0M-m+l2‍ϖk(α)u(xm+k-l2,t))+O(h2),(27)∂αu(xm,t)∂|x|α=-Ψαhα(κ1∑k=0m+l1‍ϖk(α)u(xm-k+l1,t)mmmmm+κ2∑k=0m+l2‍ϖk(α)u(xm-k+l2,t)mmmmm+κ3∑k=0m+l3‍ϖk(α)u(xm-k+l3,t)mmmmm+κ1∑k=0M-m+l1‍ϖk(α)u(xm+k-l1,t)mmmmm+κ2∑k=0M-m+l2‍ϖk(α)u(xm+k-l2,t)mmmmm+κ3∑k=0M-m+l3‍ϖk(α)u(xm+k-l3,t))+O(h3),
respectively.Here, we construct another second-order scheme and two fourth-order numerical schemes for the Riesz fractional derivative. In order to construct the new computational schemes, we introduce the following theorem.Lemma 1 (see [29]).
Letα>0, u(x,t)∈C0∞(R) with respect to x. The Fourier transforms of the left and right Riemann-Liouville fractional derivatives with respect to x are
(28)Fx(RLD-∞,xαu(x,t))=(iω)αu^(ω,t),Fx(RLDx,∞αu(x,t))=(-iω)αu^(ω,t),
where u^(ω,t) denotes the Fourier transform of the function u(x,t) with respect to x; that is,
(29)u^(ω,t)=∫R‍exp⁡(-iωx)u(x,t)dx.In [30], Tuan and Gorenflo introduced the following left fractional central difference operator:
(30)CΔ-hαu(x,t)=∑k=0∞‍ϖk(α)u(x-(k-α2)h,t).Similarly, we define the following right fractional central difference operator:(31)CΔ+hαu(x,t)=∑k=0∞‍ϖk(α)u(x+(k-α2)h,t).Analogous to the integer-order finite difference formula, we define the following fractional average operator:(32)μ±hαu(x-sh,t)=u(x±(s-(α/2))h,t)+u(x±(s+(α/2))h,t)2.Then we can get the following fractional left and right average central difference operators based on (30), (31), and (32), respectively:
(33)AC⁡Δ-hαu(x,t)=μ-hα(CΔ-hαu(x,t))=∑j=0∞‍(-1)jmmm×(αj)μ-hα(u(x-(j-α2)h,t))=12∑j=0∞(-1)jmmmm×(αj)(u(x-jh,t)mmmmmmmmmm+u(x-(j-α)h,t)),(34)AC⁡Δ+hαu(x,t)=μ+hα(CΔ+hαu(x,t))=∑j=0∞‍(-1)jmmm×(αj)μ+hα(u(x+(j-α2)h,t))=12∑j=0∞(-1)jmmmm×(αj)(u(x+jh,t)mmmmmmmmmm+u(x-(j-α)h,t)).Here, we always assume thatμ±hα can commute with the infinite summation.For the fractional left and right average central difference operators defined in (33) and (34), one has the following result.Theorem 2.
Letu(x,t) and the Fourier transform of RLD-∞,xα+2u(x,t) and RLDx,+∞α+2u(x,t) with respect to x both be in L1(R); then
(35)RLD-∞,xαu(x,t)=ACΔ-hαu(x,t)hα+O(h2),RLDx,+∞αu(x,t)=ACΔ+hαu(x,t)hα+O(h2)
uniformly holds for x∈R.Proof.
Here, we only prove (35). As u(x,t) with respect to x belongs to L1(R), then the Fourier transform of the fractional average central difference operator (33) exists and has the following form:
(36)Fx{AC⁡Δ-hαu(x,t)hα;ω}=12hα∑j=0∞‍(-1)jmmmmmm×(αj)(exp⁡(-iωjh)mmmmmmmmmmmm+exp⁡(-iω(j-α)h))u^(ω,t)=1hα(∑j=0∞‍(-1)j(αj)exp⁡(-ijωh))×(1+exp⁡(iωαh)2)u^(ω,t)=(iω)α(1-exp⁡(-iωh)iωh)α(1+exp⁡(iωαh)2)u^(ω,t).
Note that the function((1-exp⁡(-iωh))/iωh)α((1+exp⁡(iωαh))/2) has the following Taylor expansion:
(37)(1-exp⁡(-iωh)iωh)α(1+exp⁡(iωαh)2)=1+α(3α+1)24(iωh)2+O(|iωh|)4.
If we denote
(38)ϕ^(ω,h)=Fx{AC⁡Δ-hαu(x,t)hα;ω}-Fx(RLD-∞,xαu(x,t)),
then from (36), (37), and Lemma 1, we have
(39)|ϕ^(ω,h)|≤C1h2|(iω)α+2u^(ω,t)|.
In light of the conditionF(RLD-∞,xα+2u(x,t))∈L1(R), that is,
(40)∫R‍|F(RLD-∞,xα+2u(x,t))|dω<C2,
we obtain
(41)|AC⁡Δ-hαu(x,t)hα-RLD-∞,xαu(x,t)|=|ϕ(ω,h)|=12π|∫R‍exp⁡(iωh)ϕ^(ω,h)dω|≤12π∫R‍|ϕ^(ω,h)|dω≤C12π(∫R‍|(iω)α+2u^(ω,t)|dω)h2=C12π(∫R‍|F(RLD-∞,xα+2u(x,t))|dω)h2≤Ch2=O(h2),
where C=C1C2/2π. This finishes the proof.Next, we construct fourth-order difference scheme for the left and right Riemann-Liouville derivatives based on (33) and (34) by the following theorem.Theorem 3.
Letu(x,t) and the Fourier transform of RLD-∞,xα+4u(x,t) and RLDx,+∞α+4u(x,t) with respect to x both be in L1(R); then
(42)RLD-∞,xαu(x,t)=1hαACΔ-hαu(x,t)1+(α(3α+1)/24)δx2+O(h4),RLDx,+∞αu(x,t)=1hαACΔ+hαu(x,t)1+(α(3α+1)/24)δx2+O(h4)
uniformly holds for x∈R, where δx2 denotes second-order central difference operator with respect to x and is defined by δx2u(xj,t)=u(xj+1,t)-2u(xj,t)+u(xj-1,t).Proof.
It is similar to Theorem2.
Combining (3) and Theorems 2 and 3, we can get the following high-order difference schemes for the Riesz derivative:
(43)∂αu(x,t)∂|x|α=-Ψα2hα[∑j=0∞(-1)j(αj)mmmmmm×(u(x-jh,t)mmmmmmmm+u(x-(j-α)h,t))mmmmm+∑j=0∞(-1)j(αj)mmmmmmmm×(u(x+jh,t)mmmmmmmmmm+u(x+(j-α)h,t))∑j=0∞(-1)j]+O(h2),(44)∂αu(x,t)∂|x|α=-Ψα2hα[∑j=0∞(-1)j(αj)mmmmmm×u(x-jh,t)+u(x-(j-α)h,t)1+(α(3α+1)/24)δx2+∑j=0∞(-1)j(αj)mmmm×u(x+jh,t)+u(x+(j-α)h,t)1+(α(3α+1)/24)δx2]+O(h4).
Moreover, let(45)u~(x,t)={u(x,t),x∈[a,b],0,x∉[a,b];
then formulas (43) and (44) change into
(46)∂αu(x,t)∂|x|α=-Ψα2hα[∑j=0[(x-α)/h](-1)j(αj)mmmmmmmmm×(u(x-jh,t)mmmmmmmmmmm+u(x-(j-α)h,t))+∑j=0[(b-x)/h](-1)j(αj)×(u(x+jh,t)+u(x+(j-α)h,t))∑j=0[(x-α)/h](-1)j]+O(h2),(47)∂αu(x,t)∂|x|α=-Ψα2hα[∑j=0[(x-α)/h](-1)j(αj)×u(x-jh,t)+u(x-(j-α)h,t)1+(α(3α+1)/24)δx2+∑j=0[(b-x)/h](-1)j(αj)×u(x+jh,t)+u(x+(j-α)h,t)1+(α(3α+1)/24)δx2∑j=0[(x-α)/h](-1)j]+O(h4).
Finally, we derive another fourth-order numerical method for the Riesz fractional derivative which is presented in the following theorem.Theorem 4.
Letu(x,t) lie in C7(R) whose partial derivatives up to order seven with respect to x belong to L1(R). Set
(48)Lθu(x,t)=∑k=-∞∞‍gk(α)u(x-(k+θ)h,t),θ=-1,0,1,
in which
(49)gk(α)=(-1)kΓ(α+1)Γ((α/2)-k+1)Γ((α/2)+k+1);
then one has
(50)∂αu(x,t)∂|x|α=1hα[α24L-1u(x,t)-(1+α12)L0u(x,t)+α24L1u(x,t)]+O(h4).Proof.
Here, we use the Fourier transform method to prove it [23]. From [27], we know that the generating function with coefficients gk(α) satisfies
(51)|2sin(x2)|α=∑k=-∞∞‍gk(α)exp⁡(ikx).
From (3) and Lemma 1, we get the Fourier transform of the Riesz fractional derivative as follows:
(52)Fx{∂αu(x,t)∂|x|α;ω}=-Ψα[(iω)α+(-iω)α]u^(ω,t)=-|ω|αu^(ω,t).
Applying the Fourier transform to the difference operator(53)1hα[α24L-1u(x,t)-(1+α12)L0u(x,t)+α24L1u(x,t)]
with respect to x and using (51) give
(54)Fx{1hα[α24L-1u(x,t)-(1+α12)L0u(x,t)+α24L1u(x,t)];ω}=1hα[α24∑k=-∞∞‍gk(α)exp⁡(-i(k-1)ωh)u^(ω,t)-(1+α12)∑k=-∞∞‍gk(α)exp⁡(-ikωh)u^(ω,t)+α24∑k=-∞∞‍gk(α)exp⁡(-i(k+1)ωh)u^(ω,t)]=-1hα[1+α12(1-cos⁡(ωh))]|2sin(ωh2)|αu^(ω,t).
Set
(55)|ω|αu^(ω,t)=C^(h,ω)-Fx{1hα[α24L-1u(x,t)mmmmmmmmmmmm-(1+α12)L0u(x,t)mmmmmmmmmmmm+α24L1u(x,t)]1hα;ω};
then
(56)C^(h,ω)=|ω|α{|2sin(ωh/2)ωh|α}1-[1+α12(1-cos⁡(ωh))]×|2sin(ωh/2)ωh|α}u^(ω,t)=|ω|α{(11920+α-11152)1-[1+α24(ωh)2-α288(ωh)4+O(ωh)6]·[1-α24(ωh)2+α(11920+α-11152)(ωh)4+O(ωh)6(11920+α-11152)]}u^(ω,t)=-|ω|α{α(α1152+112880)(ωh)4-O(ωh)6(α1152+112880)}u^(ω,t).
Sinceu(x,t)∈C7(R) and its partial derivatives up to order seven with respect to x belong to L1(R), there exists a positive constant C~1 such that
(57)|u^(ω,t)|≤C~1(1+|ω|)-7.
So, using (56) and (57) leads to
(58)|C^(h,ω)|≤C~2h4|ω|4+α|u^(ω,t)|≤C~2h4(1+|ω|)4+α|u^(ω,t)|≤C~3h4(1+|ω|)α-3,
where C~3=C~1C~2.
At this moment, taking the inverse Fourier transformation in both sides of (55) and noting (52) give
(59)∂αu(x,t)∂|x|α=1hα[α24L-1u(x,t)-(1+α12)L0u(x,t)+α24L1u(x,t)]-C(h,ω).
In view of (58), we have
(60)|C(h,x)|=12π|∫R‍C^(h,ω)exp⁡⁡(iωx)dω|≤12π∫R‍|C^(h,ω)|dω≤C~32π(∫R‍(1+|ω|)α-3dω)h4=C~h4,
where C~=C~3/(2-α)π; that is to say,
(61)∂αu(x,t)∂|x|α=1hα[α24L-1u(x,t)-(1+α12)L0u(x,t)+α24L1u(x,t)]+O(h4).
This finishes the proof.Furthermore, (61) can be rewritten as
(62)∂αu(x,t)∂|x|α=α24hα∑k=-∞∞‍gk(α)u(x-(k+1)h,t)-(1+α12)1hα∑k=-∞∞‍gk(α)u(x-kh,t)+α24hα∑k=-∞∞‍gk(α)u(x-(k-1)h,t)+O(h4).Combining (45), one can get
(63)∂αu(xm,t)∂|x|α=α24hα∑k=-M+m+1m-1‍gk(α)u(xm-(k+1),t)+α24hα∑k=-M+m+1m-1‍gk(α)u(xm-(k-1),t)-(1+α12)1hα∑k=-M+m+1m-1‍gk(α)u(xm-k,t)+O(h4).


## 3. The Numerical Method for the Space Riesz Fractional Diffusion Equation

In this section, we only develop the fourth-order numerical method for the Riesz fractional diffusion equation (1). The second-order methods derived in the preceding section for (1) are easy so they are omitted here.By Taylor expansion, one has(64)u(xm,tn-1)=u(xm,tn)-τ∂u(xm,tn)∂t+τ22∂2u(xm,tn)∂t2+O(τ3),(65)∂u(xm,tn-1)∂t=∂u(xm,tn)∂t-τ∂2u(xm,tn)∂t2+O(τ2).Furthermore, from (65), one gets
(66)∂2u(xm,tn)∂t2=1τ(∂u(xm,tn)∂t-∂u(xm,tn-1)∂t)+O(τ).Substituting (66) into (64) yields
(67)τ2(∂u(xm,tn)∂t+∂u(xm,tn-1)∂t)=u(xm,tn)-u(xm,tn-1)+O(τ3).From (1) and (67), we have
(68)2u(xm,tn)-τK∂αu(xm,tn)∂|x|α=2u(xm,tn-1)+τK∂αu(xm,tn-1)∂|x|α+τf(xm,tn)+τf(xm,tn-1)+O(τ3).Letumn be the approximation solution of u(xm,tn). Substituting (63) into (68) and removing the high-order terms, one has the following finite difference scheme for (1):
(69)umn-μ1∑k=-M+m+1m-1‍gk(α)um-k-1n-μ1∑k=-M+m+1m-1‍gk(α)um-k+1n+μ2∑k=-M+m+1m-1‍gk(α)um-kn=umn-1+μ1∑k=-M+m+1m-1‍gk(α)um-k-1n-1+μ1∑k=-M+m+1m-1‍gk(α)um-k+1n-1-μ2∑k=-M+m+1m-1‍gk(α)um-kn-1+τ2fmn+τ2fmn-1,
where μ1=ατK/48hα and μ2=(τK/2hα)(1+(α/12)).Denoting(70)Un=(u1n,u2n,…,uM-1n)T,Fn=(f1n,f2n,…,fM-1n)T,
then system (69) can be written in the following matrix form:
(71)(I+H)Un=(I-H)Un-1+τ2Fn+τ2Fn-1,
where I is a (M-1)×(M-1) identity matrix, H=μ2G-μ1G+-μ1G-,
(72)G=(g0(α)g-1(α)g-2(α)⋯g4-M(α)g3-M(α)g2-M(α)g1(α)g0(α)g-1(α)g-2(α)⋯g4-M(α)g3-M(α)g2(α)g1(α)g0(α)g-1(α)g-2(α)⋯g4-M(α)⋮⋱⋱⋱⋱⋱⋮gM-4(α)⋯g2(α)g1(α)g0(α)g-1(α)g-2(α)gM-3(α)gM-4(α)⋯g2(α)g1(α)g0(α)g-1(α)gM-2(α)gM-3(α)gM-4(α)⋯g2(α)g1(α)g0(α)),G+=(0g0(α)g-1(α)⋯g5-M(α)g4-M(α)g3-M(α)0g1(α)g0(α)g-1(α)⋯g5-M(α)g4-M(α)0g2(α)g1(α)g0(α)g-1(α)⋯g5-M(α)⋮⋱⋱⋱⋱⋱⋮0⋯g3(α)g2(α)g1(α)g0(α)g-1(α)0gM-3(α)⋯g3(α)g2(α)g1(α)g0(α)0gM-2(α)gM-3(α)⋯g3(α)g2(α)g1(α)),G-=(g-1(α)g-2(α)g-3(α)⋯g3-M(α)g2-M(α)0g0(α)g-1(α)g-2(α)g-3(α)⋯g3-M(α)0g1(α)g0(α)g-1(α)g-2(α)g-3(α)⋯0⋮⋱⋱⋱⋱⋱⋮gM-5(α)⋯g1(α)g0(α)g-1(α)g-2(α)0gM-4(α)gM-5(α)⋯g1(α)g0(α)g-1(α)0gM-3(α)gM-4(α)gM-5(α)⋯g1(α)g0(α)0),(73)gk(α)=(-1)kΓ(α+1)Γ((α/2)-k+1)Γ((α/2)+k+1),k=2-M,3-M,…,M-3,M-2.


## 4. Stability Analysis

Now we perform the detailed stability analysis for the difference scheme (71). Firstly, we introduce some lemmas and a definition for the following discussion.Lemma 5 (see [31]).
LetA be an M-1 order positive definite matrix. Then for any parameter σ≥0, the following two inequalities:
(74)∥(I+σA)-1∥∞≤1,∥(I+σA)-1(I-σA)∥∞≤1,
hold.Definition 6 (see [32]).
Let Toeplitz matrixQM have the following form:
(75)QM=(ϱ0ϱ-1⋯ϱ2-Mϱ1-Mϱ1ϱ0ϱ-1⋯ϱ2-M⋮ϱ1ϱ0⋱⋮ϱM-2⋯⋱⋱ϱ-1ϱM-1ϱM-2⋯ϱ1ϱ0).
If the diagonals {ϱk}k=-M+1M-1 are the Fourier coefficients of function U(x,t), that is,
(76)ϱk=12π∫-ππ‍U(x,t)e-ikxdx,
then the function U(x,t) is called the generating function of QM.Lemma 7 (see [33]).
For the above Toeplitz matrixQM, let U(x,t) be a 2π-periodic continuous real-valued function defined on [-π,π]. Denote λmin⁡(QM) and λmax⁡(QM) as the smallest and largest eigenvalues of QM, respectively. Then one has
(77)Fmin⁡≤λmin⁡(QM)≤λmax⁡(QM)≤Fmax⁡,
where Fmin⁡ and Fmax⁡ are the minimum and maximum values of U(x,t). Moreover, if Fmin⁡<Fmax⁡, then all eigenvalues of QM satisfy
(78)Fmin⁡<λ(QM)<Fmax⁡,
for all M>0. And furthermore if Fmin⁡≥0, then QM is positive definite.Theorem 8.
The matrixH is positive definite.Proof.
From Definition6, we know that the generating functions of the matrices G, G+, and G- are
(79)FG(x)=∑k=-∞∞‍gk(α)eikx,FG+(x)=∑k=-∞∞‍gk(α)ei(k+1)x,FG-(x)=∑k=-∞∞‍gk(α)ei(k-1)x,
respectively.
Therefore the generating function of the matrixH is
(80)FH(x)=μ2FG(x)-μ1FG+(x)-μ1FG-(x);
that is,
(81)FH(x)=μ2∑k=-∞∞‍gk(α)eikx-μ1∑k=-∞∞‍gk(α)ei(k+1)x-μ1∑k=-∞∞‍gk(α)ei(k-1)x=(μ2-2μ1cos⁡x)∑k=-∞∞‍gk(α)eikx=Kτ12hα[6+αsin2(x2)]|2sin(x2)|α≥0.
It immediately follows thatH is a positive definite matrix in view of Lemma 7.Theorem 9.
The difference scheme (71) is unconditionally stable.Proof.
LetUn and U~n be the exact and numerical solutions of difference equation (71). Since the matrix (I+H)-1 is invertible, then we can obtain the following error equation:
(82)En=QEn-1,
where En=U~n-Un and Q=(I+H)-1(I-H).
By (82), we have
(83)En=Qn-1E0.
Using (83) and Lemma 5, one has
(84)∥En∥∞≤∥Qn-1∥∞∥E0∥∞≤∥Q∥∞n-1∥E0∥∞≤∥E0∥∞,
which means difference scheme (71) is unconditionally stable, and so is the difference scheme (69).


## 5. Convergence Analysis

In this section, we give the local truncation error analysis of difference scheme (69).Theorem 10.
The local truncation error of numerical scheme (69) for approximating the space Riesz fractional diffusion equation (1) at the point (xm,tn) is O(τ2+h4).Proof.
We define the local truncation error at the point(xm,tn) as
(85)Rmn=2[u(xm,tn)-u(xm,tn-1)]τ+2τ[μ2∑k=-M+m+1m-1‍gk(α)u(xm-k,tn)-μ1∑k=-M+m+1m-1‍gk(α)u(xm-k-1,tn)-μ1∑k=-M+m+1m-1‍gk(α)u(xm-k+1,tn)]+2τ[μ2∑k=-M+m+1m-1‍gk(α)u(xm-k,tn-1)-μ1∑k=-M+m+1m-1‍gk(α)u(xm-k-1,tn-1)-μ1∑k=-M+m+1m-1‍gk(α)u(xm-k+1,tn-1)]-f(xm,tn)-f(xm,tn-1).
Utilizing the Taylor expansion, we have(86)2[u(xm,tn)-u(xm,tn-1)]τ=∂u(xm,tn)∂t+∂u(xm,tn-1)∂t+O(τ2).
Based on (63) and (86), we obtain
(87)Rmn=[∂u(xm,tn)∂t-K∂αu(xm,tn)∂|x|α-f(xm,tn)]+[∂u(xm,tn-1)∂t-K∂αu(xm,tn)∂|x|α-f(xm,tn)∂u(xm,tn-1)∂t]+O(τ2+h4).
Substituting (1) into (87), we get the following local truncation error:
(88)Rmn=O(τ2+h4).
Next, we study the convergence of difference scheme (69).Theorem 11.
The order of the convergence of difference scheme (69) is O(τ2+h4); that is to say, there exists a positive constant C such that
(89)|u(xm,tk)-umk|≤C(τ2+h4).Proof.
Letemn=u(xm,tk)-umk. Then it follows from (69) and (85) that
(90)emn-μ1∑k=-M+m+1m-1‍gkem-k-1n-μ1∑k=-M+m+1m-1‍gkem-k+1n+μ2∑k=-M+m+1m-1‍gkem-kn=emn-1+μ1∑k=-M+m+1m-1‍gkem-k-1n-1+μ1∑k=-M+m+1m-1‍gkem-k+1n-1-μ2∑k=-M+m+1m-1‍gkem-kn-1+τ2Rmn.
Denote(91)En=(e1n,e2n,…,eM-1n)T,Rn=(R1n,R2n,…,RM-1n)T;
then (90) can be written in a matrix form
(92)(I+H)En=(I-H)En-1+τ2Rn.
Furthermore, we can obtain
(93)En=QEn-1+S=Q(QEn-2+S)+S=⋯=(Qn-1+Qn-2+⋯+I)S,
where S=(τ/2)(I+H)-1Rn.
According to Lemma5 and Theorem 10, we get
(94)∥S∥∞=∥τ2(I+H)-1Rn∥∞≤τ2∥(I+H)-1∥∞∥Rn∥∞≤τC1(τ2+h4).
Combining (93) and (94) gives
(95)∥En∥∞=∥(Qn-1+Qn-2+⋯+I)S∥∞≤∥(Qn-1+Qn-2+⋯+I)∥∞∥S∥∞≤(∥Qn-1∥∞+∥Qn-2∥∞+⋯+∥I∥∞)∥S∥∞≤(∥Q∥∞n-1+∥Q∥∞n-2+⋯+∥I∥∞)∥S∥∞≤nτC1(τ2+h4).
Noting thatnτ≤Nτ=T, we easily get
(96)∥En∥∞≤TC1(τ2+h4)=C(τ2+h4),
where C=TC1. This completes the proof.


## 6. Numerical Examples

In the present section, some numerical examples are presented to demonstrate the theoretical analysis.Example 1.
Consider the functionu(x)=x2(1-x)2, x∈[0,1].
The Riesz fractional derivative of the functionu(x) is
(97)∂αu(x)∂|x|α=-sec(π2α){1Γ(3-α)[x2-α+(1-x)2-α]-6Γ(4-α)[x3-α+(1-x)3-α]+12Γ(5-α)[x4-α+(1-x)4-α]}.
Table1 lists the absolute errors and convergence orders at x=0.5 by numerical scheme (46) for function u(x) with different α.Table 1
The absolute errors and convergence orders of Example1 by numerical scheme (46).
α
h
The absolute error
The convergence order
0.2
1
10
3.630966
e
-
003
—
1
20
9.120270
e
-
004
1.9932
1
40
2.285315
e
-
004
1.9967
1
80
5.719787
e
-
005
1.9984
0.4
1
10
5.124542
e
-
003
—
1
20
1.289681
e
-
003
1.9904
1
40
3.234889
e
-
004
1.9952
1
80
8.100606
e
-
005
1.9976
0.6
1
10
4.629914
e
-
003
—
1
20
1.164707
e
-
003
1.9910
1
40
2.920982
e
-
004
1.9954
1
80
7.314118
e
-
005
1.9977
0.8
1
10
2.652282
e
-
003
—
1
20
6.653815
e
-
004
1.9950
1
40
1.666617
e
-
004
1.9973
1
80
4.170691
e
-
005
1.9986
1.2
1
10
2.371107
e
-
003
—
1
20
5.878877
e
-
004
2.0119
1
40
1.464631
e
-
004
2.0050
1
80
3.655831
e
-
005
2.0023
1.4
1
10
3.614913
e
-
003
—
1
20
8.891952
e
-
004
2.0234
1
40
2.207491
e
-
004
2.0101
1
80
5.500869
e
-
005
2.0047
1.6
1
10
3.311900
e
-
003
—
1
20
8.088550
e
-
004
2.0337
1
40
2.001866
e
-
004
2.0145
1
80
4.981273
e
-
005
2.0068
1.8
1
10
1.748268
e
-
003
—
1
20
4.265789
e
-
004
2.0350
1
40
1.055454
e
-
004
2.0149
1
80
2.625954
e
-
005
2.0069Example 2.
Consider the functionu(x)=x6(1-x)6, x∈[0,1].
The Riesz fractional derivative of the functionu(x) is
(98)∂αu(x)∂|x|α=12etsec(π2α){Γ(7)Γ(7-α)[x6-α+(1-x)6-α]-6Γ(8)Γ(8-α)[x7-α+(1-x)7-α]+15Γ(9)Γ(9-α)[x8-α+(1-x)8-α]-20Γ(10)Γ(10-α)[x9-α+(1-x)9-α]+15Γ(8)Γ(11-α)[x10-α+(1-x)10-α]-6Γ(12)Γ(12-α)[x11-α+(1-x)11-α]+Γ(13)Γ(13-α)[x12-α+(1-x)12-α]}.
Table2 lists the absolute errors and convergence orders at x=0.5 by numerical scheme (47) for function u(x) with different α. These computational results confirm the second-order and fourth-order of the numerical formulas (46) and (47), respectively.Table 2
The absolute errors and convergence orders of Example2 by numerical scheme (47).
α
h
The absolute error
The convergence order
0.2
1
20
1.571776
e
-
008
—
1
40
9.923695
e
-
010
3.9854
1
80
6.218021
e
-
011
3.9963
1
160
3.888658
e
-
012
3.9991
0.4
1
20
6.255711
e
-
008
—
1
40
3.955450
e
-
009
3.9833
1
80
2.479333
e
-
010
3.9958
1
160
1.550707
e
-
011
3.9990
0.6
1
20
1.660570
e
-
007
—
1
40
1.051535
e
-
008
3.9811
1
80
6.593628
e
-
010
3.9953
1
160
4.124394
e
-
011
3.9988
0.8
1
20
3.501076
e
-
007
—
1
40
2.219986
e
-
008
3.9792
1
80
1.392502
e
-
009
3.9948
1
160
8.711096
e
-
011
3.9987
1.2
1
20
7.029869
e
-
007
—
1
40
4.446255
e
-
008
3.9828
1
80
2.787163
e
-
009
3.9957
1
160
1.743236
e
-
010
3.9990
1.4
1
20
4.372243
e
-
008
—
1
40
3.985260
e
-
009
3.4556
1
80
2.690303
e
-
010
3.8888
1
160
1.712778
e
-
011
3.9734
1.6
1
20
3.751536
e
-
006
—
1
40
2.440706
e
-
007
3.9421
1
80
1.540687
e
-
008
3.9857
1
160
9.653215
e
-
010
3.9964
1.8
1
20
1.592788
e
-
005
—
1
40
1.034107
e
-
006
3.9451
1
80
6.524548
e
-
008
3.9864
1
160
4.087469
e
-
009
3.9966Example 3.
We consider the equation(99)∂u(x,t)∂t=K∂αu(x,t)∂|x|α+f(x,t),mmmmmm0<x<1,0<t≤1,
where
(100)f(x,t)=(2α+1)t2αx4(1-x)4+t2α+1sec(π2α)×{12Γ(5-α)[x4-α+(1-x)4-α]-240Γ(6-α)[x5-α+(1-x)5-α]+2160Γ(7-α)[x6-α+(1-x)6-α]-10080Γ(8-α)[x7-α+(1-x)7-α]+20160Γ(9-α)[x8-α+(1-x)8-α]},
together with the initial condition u(x,0)=0 and homogeneous boundary value conditions.
The analytical solution of (99) is u(x,t)=t2α+1x4(1-x)4.
Table3 lists the maximum error and time and space convergence orders for different α. Numerical results show that the convergence order of the difference scheme (69) is O(τ2+h4), which are in line with our theoretical analysis.
Figures1 and 2 show the comparison of the analytical and numerical solutions with α=1.1 at t=0.2 and α=1.9 at t=0.8, respectively. It can be seen that the numerical solutions are in excellent agreement with the analytical solutions.
Figures3, 4, 5, and 6 display the numerical and analytical solutions surface with α=1.2 and α=1.8. By Figures 3–6, we find that the corresponding two groups of figures are almost the same, respectively.Table 3
The maximum errors and temporal and spatial convergence orders of Example3 by difference scheme (66).
α
The maximum errors
Temporal convergence orders
Spatial convergence orders
1.2
h
=
1
4
,
τ
=
1
2
3.070974
e
-
004
—
—
h
=
1
8
,
τ
=
1
8
2.038855
e
-
005
1.9564
3.9129
h
=
1
16
,
τ
=
1
32
1.318397
e
-
006
1.9755
3.9509
h
=
1
32
,
τ
=
1
128
8.371766
e
-
008
1.9886
3.9771
h
=
1
64
,
τ
=
1
512
5.279883
e
-
009
1.9935
3.9870
h
=
1
128
,
τ
=
1
2048
3.486707
e
-
010
1.9603
3.9206
1.4
h
=
1
4
,
τ
=
1
2
3.715656
e
-
004
—
—
h
=
1
8
,
τ
=
1
8
2.465865
e
-
005
1.9567
3.9135
h
=
1
16
,
τ
=
1
32
1.603083
e
-
006
1.9716
3.9432
h
=
1
32
,
τ
=
1
128
1.022389
e
-
007
1.9854
3.9708
h
=
1
64
,
τ
=
1
512
6.472320
e
-
009
1.9908
3.9815
h
=
1
128
,
τ
=
1
2048
4.083361
e
-
010
1.9932
3.9865
1.6
h
=
1
4
,
τ
=
1
2
4.193160
e
-
004
—
—
h
=
1
8
,
τ
=
1
8
2.775395
e
-
005
1.9586
3.9173
h
=
1
16
,
τ
=
1
32
1.810618
e
-
006
1.9691
3.9381
h
=
1
32
,
τ
=
1
128
1.158839
e
-
007
1.9829
3.9657
h
=
1
64
,
τ
=
1
512
7.363046
e
-
009
1.9881
3.9762
h
=
1
128
,
τ
=
1
2048
4.662698
e
-
010
1.9905
3.9811
1.8
h
=
1
4
,
τ
=
1
2
4.417615
e
-
004
—
—
h
=
1
8
,
τ
=
1
8
2.895409
e
-
005
1.9657
3.9314
h
=
1
16
,
τ
=
1
32
1.884230
e
-
006
1.9709
3.9417
h
=
1
32
,
τ
=
1
128
1.204909
e
-
007
1.9835
3.9670
h
=
1
64
,
τ
=
1
512
7.658088
e
-
009
1.9879
3.9758
h
=
1
128
,
τ
=
1
2048
4.855499
e
-
010
1.9896
3.9793Figure 1
Comparison between the analytical solution and the numerical solution att=0.2 with α=1.1 in Example 3  (τ=1/50,h=1/50).Figure 2
Comparison between the analytical solution and the numerical solution att=0.8 with α=1.9 in Example 3  (τ=1/100,h=1/80).Figure 3
The numerical solution surface whenα=1.2 in Example 3  (τ=1/500,h=1/100).Figure 4
The analytical solution surface whenα=1.2 in Example 3.Figure 5
The numerical solution surface whenα=1.8 in Example 3  (τ=1/50,h=1/80).Figure 6
The analytical solution surface whenα=1.8 in Example 3.


## 7. Conclusions

In this paper, we construct one second-order method and two fourth-order methods for the Riesz fractional derivative. Later on, a fourth-order numerical formula is used for the Riesz space fractional diffusion equation. The stability and convergence of this numerical method are analyzed by the matrix method in detail. Finally, some numerical results are given to demonstrate the effectiveness of numerical method. From the above results, it is possible to claim that the methods and techniques discussed in this paper are useful for solving some other fractional differential equations with Riesz derivatives.



---

*Source: 653797-2014-05-22.xml*