Second order WSGD operators II: A new family of difference schemes for space fractional advection diffusion equation

# Second order WSGD operators II: A new family of difference schemes for space fractional advection diffusion equation

Can Li, Weihua Deng111Corresponding author. E-mail address: dengwh@lzu.edu.cn.
Department of Applied Mathematics, School of Sciences, Xi’an University of Technology,
Xi’an, Shaanxi 710054, P.R.China
School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, P.R.China
###### Abstract

The second order weighted and shifted Grünwald difference (WSGD) operators are developed in [Tian et al., arXiv:1201.5949] to solve space fractional partial differential equations. Along this direction, we further design a new family of second order WSGD operators; by properly choosing the weighted parameters, they can be effectively used to discretize space (Riemann-Liouville) fractional derivatives. Based on the new second order WSGD operators, we derive a family of difference schemes for the space fractional advection diffusion equation. By von Neumann stability analysis, it is proved that the obtained schemes are unconditionally stable. Finally, extensive numerical experiments are performed to demonstrate the performance of the schemes and confirm the convergent orders.

¡¡¡¡

Mathematics Subject Classification(2010): 26A33, 65M06, 65M12

Key words: Riemann-Liouville fractional derivative, WSGD operator, Fractional advection diffusion equation, Finite difference approximation, Stability

## 1 Introduction

As an extension of classical calculus, fractional calculus also has more than three centuries of history. However, it seems that the applications of fractional calculus to physics and engineering attract enough attention only in the last few decades. Nowadays, it has been found that a broad range of non-classical phenomena in the applied sciences and engineering can be well described by some fractional kinetic equations [18]. The fractional calculus is becoming more and more popular, especially in describing the anomalous diffusions [1, 10, 18, 8], which arise in physics, chemistry, biology, and other complex dynamics. With the appearance of various kinds of fractional partial differential equations (PDEs), finding the ways to effectively solve them becomes a natural topic. Based on the integral transformations, sometimes we can find the analytical solutions of linear fractional PDEs with constant coefficients [9, 10, 18]. Even in these cases, most of the time their solutions are expressed by transcendental functions or infinite series. So looking for the numerical solutions of fractional PDEs becomes a realistic expectation in practical applications; and designing efficient and robust numerical schemes for fractional PDEs is the basis to this task.

Over the past decade, the finite difference methods are implemented in simulating the space fractional advection diffusion equations [13, 5, 15, 14, 16, 23]. However, most of the time, the accuracy of finite difference scheme or finite volume method is first order. In recent years, high order discretizations and fast methods for space fractional PDEs with Riemann-Liouville fractional derivatives attract many authors’ interests. Based on the Toeplitz-like structure of the difference matrix, Wang et al [28] solve the algebraic equation corresponding to fractional diffusion equation with a cost; Pang and Sun [21] propose a multigrid method to solve the discretized system of the fractional diffusion equation. By using the linear spline approximation, Sousa and Li provide a second order discretization for the Riemann-Liouville fractional derivatives and establish an unconditionally stable finite difference method for one-dimensional fractional diffusion equation in [24]; the results on two-dimensional two-sided space fractional convection diffusion equation in finite domain can be seen in [3]. Ortigueira [19] gives the “fractional centred derivative” to approximate the Riesz fractional derivative with second order accuracy, and this scheme is used by Çelik and Duman in [11] to approximate fractional diffusion equation with the Riesz fractional derivative in a finite domain.

More recently, Tian et al [27] propose the second order difference approximations, called weighted and shifted Grünwald difference (WSGD) approximations, to the Riemann-Liouville space fractional derivatives. The basic idea of the WSGD approximation is to cancel the low order terms by combining the Grünwald difference operators with different shifts and weights. Along this direction, this paper further introduces a new family of WSGD opertors by providing more to be chosen parameters, called second order WSGD operators II. As specific applications, the second order WSGD operators II are used to discretize the space fractional derivative of the space fractional advection diffusion equation. And the von Neumann stability analyses are performed to the obtained numerical schemes. We theoretically prove and numerically verify that with proper chosen parameters for the WSGD space discretizations, the obtained implicit and Crank-Nicolson schemes are both unconditionally von Neumann stable.

The outline of this paper is organized as follows. In Section 2, we first introduce a new family of second order discretizations of the Riemann-Liouville fractional derivatives. Then two kinds of difference schemes for the one dimensional space fractional advection diffusion equation are presented in Section 3, and the corresponding numerical stabilities are performed. Furthermore, in Section 4, the numerical schemes and the proof of numerical stability are extended to two dimensional cases. Finally, in Section 5, extensive numerical experiments are carried out to demonstrate the performance of the proposed numerical schemes and confirm theoretical analysis.

## 2 Second order WSGD discretizations for the Riemann-Liouville fractional derivatives

There are some different definitions for the fractional derivatives, such as the Grünwald-Letnikov derivative, the Riemann-Liouville derivative, the Caputo derivative, and other modified definitions for the practical applications [10]. Most of the time, they are not completely equivalent; in particular, the Grünwald-Letnikov and Riemann-Liouville derivatives are equivalent if ignoring the regularity requirements of the functions being performed. Usually, the more popular space fractional derivative is the Riemann-Liouville derivative. Let the function be defined on the finite interval , then the left and right Riemann-Liouville fractional derivatives of order are, respectively, defined by

and

 xDαbu(x)=(−1)nΓ(n−α)∂n∂xn∫bx(s−x)n−α−1u(s)ds. (2)

Based on the definition of the Grünwald-Letnikov derivative and its relationship with the left and right Riemann-Liouville derivatives, we have [10]

 xDαbu(x)=1hα[b−xh]∑k=0w(α)ku(x+kh)+O(h), (4)

where , is a positive integer, is the gamma function, and is the normalized Grünwald weights. However, the difference scheme based on (3) and/or (4) to discretize space fractional derivatives for the time dependent problems is unconditionally unstable [15]. To overcome this problem, Meerschaert and Tadjeran firstly introduce the shifted Grünwald-Letnikov approximation formulas [15]

 xDαbu(x)=1hα[b−xh]+q∑k=0w(α)ku(x+(k−p)h)+O(h). (6)

If a function has support on the finite interval , then it can be extended to the whole real line . In this case the symbol (or ) in left (right) Riemann-Liouville fractional derivative (1) ((2)) can be replaced by (or ). And the following properties hold.

###### Lemma 1.

[15] Let , and its Fourier transform belong to . Let and . Define

 Aαh,pu(x):=1hα∞∑k=0w(α)ku(x−(k−p)h),
 Bαh,pu(x):=1hα∞∑k=0w(α)ku(x+(k−p))h).

Then

 Aαh,pu(x)=−∞Dαxu(x)+O(h),
 Bαh,qu(x)=xDα∞u(x)+O(h),

uniformly in as .

Inspired by the shifted Grünwald difference operators and working along the direction of the ideas given in [27], we derive the following second order approximation for the Riemann-Liouville fractional derivatives.

###### Theorem 2.

Let , and its Fourier transform belong to , and define the left WSGD operator by

 LDα,λ1,λ2,…,λmh,p1,p2,…,pmu(x)=m∑j=1λjAαh,pju(x), (7)

then, for an integer , we have

 LDα,λ1,λ2,…,λmh,p1,p2,…,pmu(x)=−∞Dαxu(x)+O(h2), (8)

uniformly for , where are real numbers and satisfy the following linear system

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩m∑j=1λj=1,m∑j=1λj(pj−α2)=0. (9)
###### Proof.

Using the Fourier transform of the left Riemann-Liouville fractional derivative [10]

 F[−∞Dαxu](ω)=(iω)αF[u](ω), (10)

and taking Fourier transform on the left hand of formula (7), we obtain

 F[LDα,λ1,λ2,…,λmh,p1,p2,…,pmu(x)](ω) =m∑j=1λj(1hα∞∑k=0w(α)ke−(k−pj)h(iω)^u(ω)) (11) =(iω)αm∑j=1λjPh,j(iω)^u(ω),

where

 Ph,j(z)=epjhz(1−e−hzhz)α=1+(pj−α2)hz+O(|z|2h2), z=iω, i=√−1. (12)

Denoting , in view of (12) and (9), we have

 |^ϕ(ω,h)|≤Ch2|iω|α+2|^u(ω)|. (13)

Due to , we have

 |LDα,λ1,λ2,…,λmh,p1,p2,…,pmu−−∞Dαxu|=|ϕ|≤12π∫R|^ϕ(ω,h)|≤C∥F[−∞Dα+2xu](ω)∥L1h2=O(h2). (14)

With the Fourier transform of the right Riemann-Liouville fractional derivative [10]

 F[xDα∞u](ω)=(−iω)αF[u](ω),

by the similar arguments to the above, we get

###### Theorem 3.

Let , and its Fourier transform belong to , and define the right WSGD operator by

 RDα,γ1,γ2,…,γmh,q1,q2,…,qmu(x)=m∑j=1γjBαh,qju(x), (15)

then, for an integer , we have

 RDα,γ1,γ2,…,γmh,q1,q2,…,qmu(x)=xDα∞u(x)+O(h2), (16)

uniformly for , where are real numbers and satisfy

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩m∑j=1γj=1,m∑j=1γj(qj−α2)=0. (17)

Here are determined by solving the above linear system (17).

###### Remark 1.

Note that the desired second order accuracy is obtained by weighting the series expansion of the trigonometric polynomials in frequency space, i.e.,

 Ph,j(z)=epjhz(1−e−hzhz)α=1+(pj−α2)hz+(p2j2−α2pj+3α2+α24)h2z2+O(|z|3h3), z=ik. (18)

For example,

 3∑j=1λjPh,j(z)=3∑j=1λj+3∑j=1(λj(pj−α2))hz+Cαh2z2+O(|z|3h3). (19)

where

 Cα=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩α−α24+3α2+α24, (λ1,λ2,λ3)=(α2,2−α2,0),2−α24+3α2+α24,(λ1,λ2,λ3)=(2+α4,0,2−α4),4−α24+3α2+α24, (λ1,λ2,λ3)=(2+α2,−2−α2,1). (20)

For the similar method for the classical derivatives one can refer to [6]. And following this idea we can get higher order accurate difference approximations for the Riemann-Liouville fractional derivatives. In the following sections, we focus on the second order difference approximations with three parameters .

Supposing that a function has support on the finite domain , then we have

If the function is non-periodic we are more interested in that are integers and satisfy for numerical consideration. For the convenience, we usually take in practical computations. Here are determined by solving the above system of linear equations (9). When , the linear system (9) has the unique solution [27]. However, for , using the knowledge of linear algebra, we know that the system (9) has infinitely many solutions. Now we assume that is given in linear system of equations (9) with , then we have

 λ2=2(p1−p3)λ1+2p3−α2(p3−p2), λ3=2(p2−p1)λ1−2p2+α2(p3−p2), p2≠p3. (22)

If is given in linear algebra system (9), then we have

 λ1=2(p2−p3)λ2+2p3−α2(p3−p1), λ3=2(p1−p2)λ2+2p1+α2(p3−p1), p1≠p3. (23)

And if is known, then the solution of the linear system (9) gives

 λ1=2(p3−p2)λ3−2p2+α2(p1−p2), λ2=2(p3−p2)λ3+2p2−α2(p1−p2), p1≠p2. (24)

Furthermore, if the integers satisfy , then the -WSGD operators in finite domain can be written as the following form

 LDα,λ1,λ2,λ3h,1,0,−1u(x)=λ1hα[x−ah]+1∑k=0w(α)ku(x−(k−1)h)+λ2hα[x−ah]∑k=0w(α)ku(x−kh)+λ3hα[x−ah]−1∑k=0w(α)ku(x−(k+1)h), (25) RDα,λ1,λ2,λ3h,1,0,−1u(x)=λ1hα[b−xh]+1∑k=0w(α)ku(x+(k−1)h)+λ2hα[b−xh]∑k=0w(α)ku(x+kh)+λ1hα[b−xh]−1∑k=0w(α)ku(x+(k+1)h),

where the parameters satisfy the following linear system

 {λ1+λ2+λ3=1,λ1−λ3=α2. (26)

With the help of the knowledge of linear algebra, the solutions of the system of the linear algebraic equations (26) can be collected by the following three sets

 S1={λ2=2+α2−2λ1, λ3=λ1−α2, λ1 is~{}given}, (27)
 S2={λ1=2+α4−λ22, λ3=2−α4−λ22, λ2 is~{}given}, (28)

and

 S3={λ1=α2+λ3, λ2=2−α2−2λ3, λ3 is~{}given}. (29)

It produces a second order approximation if we take any choice of in . Particularly, if we take for in , it recovers the second order approximations presented in [27]. After rearranging the weights , the second order approximations for the Riemann-Liouville fractional derivatives at point give

where the weights are

 g(α)0=λ1w(α)0, g(α)1=λ1w(α)1+λ2w(α)0,g(α)k=λ1w(α)k+λ2w(α)k−1+λ3w(α)k−2, k≥2. (31)

## 3 Second order WSGD schemes: one dimensional case

We first consider the initial boundary value problem of the following one dimensional advection diffusion equation

where can stand for the concentration of a solute at a point at time [1], is the source term, and are, respectively, the advection and diffusion coefficients, and and are the left and right Riemann-Liouville fractional derivatives of order , respectively, defined by (1) and (2) with .

For the numerical approximations, we define , is the equidistant grid size in space, for , so that . Denote In the space discretizations, we choose the WSGD operators and to approximate the Riemann-Liouville fractional derivatives and , respectively. And the first order space derivative is approximated by the standard central difference. If the time derivative is approximated by the implicit Euler discretization, we have

 un+1j−unjτ+vxun+1j+1−un+1j−12h=dx(LDα,λ1,λ2,λ3h,1,0,−1+RDα,λ1,λ2,λ3h,1,0,−1)un+1j+fn+1j+O(τ+h2),

Dropping the truncation error term, we get the implicit WSGD scheme

 Un+1j=Unj−τ⋅vxUn+1j+1−Un+1j−12h+dxτhαj+1∑k=0g(α)kUn+1j−k+1+dxτhαN−j+1∑k=0g(α)kUn+1j+k−1+τfn+1j. (33)

To improve the accuracy of the numerical scheme in time direction, using the Crank-Nicolson time discretization, we get the following Crank-Nicolson-WSGD scheme

 Un+1j−Unjτ+vxUnj+1−Unj−14h−(dx2hαj+1∑k=0g(α)kUnj−k+1+dx2hαN−j+1∑k=0g(α)kUnj+k−1) (34)

### 3.1 Stability analysis

To study the stability of the new scheme, we use the von Neumann linear stability analysis [2, 22]; and assume that the solution of the problem (32) can be zero extended to the whole real line . Letting be the perturbation error and reasonably replacing the symbols and in schemes (33) and (34) with (since should be zero beyond the considered domain), we obtain the following error equations

 ϵn+1j=ϵnj−vxτϵn+1j+1−ϵn+1j−12h+dxτ2hα∞∑k=0g(α)kϵn+1j−k+1+dxτ2hα∞∑k=0g(α)kϵn+1j+k−1; (35)
 ϵn+1j=ϵnj−vxτϵnj+1−ϵnj−14h+(dxτ2hα∞∑k=0g(α)kϵnj−k+1+dxτ2hα∞∑k=0g(α)kϵnj+k−1) (36) −vxτϵn+1j+1−ϵn+1j−14h+dxτ2hα∞∑k=0g(α)kϵn+1j−k+1+dxτ2hα∞∑k=0g(α)kϵn+1j+k−1.
###### Theorem 4.

Let be chosen in set or or , and assume that the trigonometric polynomial , for all , , then the implicit WSGD scheme (33) is unconditionally von Neumann stable, where is defined by (40).

###### Proof.

Let be the solution of (35), where , is the amplitude at time level and is the phase angle with wavelength . We need to show that the amplification factor satisfies the relation for all in . In fact, by substituting the expressions of and into (35), we get the amplification factor of the implicit WSGD difference scheme

 ρ(θ)=11+ν2(eiθ−e−iθ)−λ2Q(θ,α;λ1,λ2,λ3), (37)

where and

 (38)

In light of the relation (31), we get

 Q(θ,α;λ1,λ2,λ3) =∞∑k=0g(α)kei(1−k)θ+∞∑k=0g(α)ke−i(1−k)θ +λ1e−iθ∞∑k=0w(α)keikθ+λ2∞∑k=0w(α)keikθ+λ3eiθ∞∑k=0w(α)keikθ.

By a straightforward calculation, the trigonometric polynomial gives

 Q(θ,α;λ1,λ2,λ3) =λ1(eiθ(1−e−iθ)α+e−iθ(1−eiθ)α)+λ2((1−e−iθ)α+(1−eiθ)α) +λ3(e−iθ(1−e−iθ)α+eiθ(1−eiθ)α).

In view of is a real-valued and even function for given , we just consider its principal value on . Furthermore, using the relation

 eiϕ−eiφ=2isin(ϕ−φ2)eiϕ+φ2, (39)

after a straightforward calculation yields

 Q(θ,α;λ1,λ2,λ3)=(2sin(θ2))α (λ1cos(α2(θ−π)−θ)+λ2cos(α2(θ−π))+λ3cos(α2(θ−π)+θ)). (40)

Then the corresponding amplification factor is obtained as

 ρ(θ)=11+iνsin(θ)−λ2Q(θ,α;λ1,λ2,λ3).

Due to the functions is non-positive, we get

 |ρ(θ)|=1√(1−λ2Q(θ,α;λ1,λ2,λ3))2+(νsin(θ))2≤1,

which means that the implicit WSGD difference scheme is unconditionally stable. ∎

###### Theorem 5.

Let be chosen in set or or , and assume that the trigonometric polynomial for all , then the Crank-Nicolson-WSGD scheme (34) is unconditionally von Neumann stable, where is defined by (40).

###### Proof.

Inserting into (36), we get the amplification factor of Crank-Nicolson-WSGD difference scheme

 ρ(θ)=1−ν2(eiθ−e−iθ)+λ2(∑∞k=0g(α)kei(1−k)θ+∑∞k=0g(α)ke−i(1−k)θ)1+ν2(eiθ−e−iθ)−λ2(∑∞k=0g(α)kei(1−k)θ+∑∞k=0g(α)ke−i(1−k)θ), (41)

with .

By direct calculation, we have

 ρ(θ)=1−iνsin(θ)+λ2Q(θ,α;λ1,λ2,λ3)1+iνsin(θ)−λ2Q(θ,α;λ1,λ2,λ3), (42)

where is given by (40). Furthermore, using the inequality

 (1−a)2+b2(1+a)2+b2≤1−2a1