Fourth order quasi-compact difference schemes for (tempered) space fractional diffusion equationsThis work was supported by the National Natural Science Foundation of China under Grant Nos. 11271173 and 11471150.

# Fourth order quasi-compact difference schemes for (tempered) space fractional diffusion equations††thanks: This work was supported by the National Natural Science Foundation of China under Grant Nos. 11271173 and 11471150.

Yanyan Yu School of Mathematics and Statistics, Gansu Key Laboratory of Applied Mathematics and Complex Systems, Lanzhou University, Lanzhou 730000, People’s Republic of China (yuyy1029@126.com).    Weihua Deng School of Mathematics and Statistics, Gansu Key Laboratory of Applied Mathematics and Complex Systems, Lanzhou University, Lanzhou 730000, People’s Republic of China (dengwh@lzu.edu.cn).    Yujiang Wu School of Mathematics and Statistics, Gansu Key Laboratory of Applied Mathematics and Complex Systems, Lanzhou University, Lanzhou 730000, People’s Republic of China (myjaw@lzu.edu.cn).
###### Abstract

The continuous time random walk (CTRW) underlies many fundamental processes in non-equilibrium statistical physics. When the jump length of CTRW obeys a power-law distribution, its corresponding Fokker-Planck equation has space fractional derivative, which characterizes Lévy flights. Sometimes the infinite variance of Lévy flight discourages it as a physical approach; exponentially tempering the power-law jump length of CTRW makes it more ‘physical’ and the tempered space fractional diffusion equation appears. This paper provides the basic strategy of deriving the high order quasi-compact discretizations for space fractional derivative and tempered space fractional derivative. The fourth order quasi-compact discretization for space fractional derivative is applied to solve space fractional diffusion equation and the unconditional stability and convergence of the scheme are theoretically proved and numerically verified. Furthermore, the tempered space fractional diffusion equation is effectively solved by its counterpart of the fourth order quasi-compact scheme; and the convergence orders are verified numerically.

Key words. space fractional derivative, tempered space fractional derivative, shifted Grünwald discretization, quasi-compact difference scheme, numerical stability and convergence.

Subject classifications. 65M06, 65M12, 26A33

## 1 Introduction

In recent years, more and more scientific and engineering problems are involved in fractional calculus. They range from relaxation oscillation phenomena [1] to viscoelasticity [2], and from control theory [3] to transport problem [4]. The fractional diffusion equation has been put forward as a more suitable model to describe ion channel gating dynamics [5] and subdiffusive anomalous transport in an external field [6], which are resulted in from the continuous time random walk (CTRW) in the scaling limit. The CTRW is a mathematical formalization of a path that consists of a succession of random steps including the elements of random waiting time and jump length; and it underlies many fundamental stochastic processes in statistical physics. When the first moment of the distribution of waiting time and the second moment of jump length are finite, the probability density function (PDF) of the particle’s location and time satisfies the classical diffusion equation. However, if the jump length obeys the power-law distribution, the PDF of the particle’s location and time is the solution of space fractional diffusion equation; and the corresponding dynamics is called Lévy flight. Sometimes the jumps of the particles are limited by the finite size of the physical system and the infinite variance of Lévy flight discourages it as a physical approach. So the power-law distribution of the jump length is expected to be truncated [7] or exponentially tempered [8]. For the CTRW with the distribution of the tempered jump length [9], the corresponding PDF of the particles satisfies the tempered space fractional diffusion equation [8].

It seems that there are less works for the numerical solutions of tempered space fractional diffusion equation [10]. However, for the space fractional diffusion or advection-diffusion equation, much progress has been made for its numerical methods, e.g., [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Transforming the Riemann-Liouville fractional derivative to Caputo fractional derivative, the space fractional Fokker-Planck equation is solved by the method of lines in [11]. Using the superconvergence of Grünwald discretization at a particular point, a second order finite difference scheme is proposed in [17]. Based on the difference discretization and spline approximation to the Riemann-Liouville fractional derivative, a second order scheme is presented for the three dimensional space fractional partial differential equations in [20]. Currently, the most popular discretization scheme for the space Riemann-Liouville fractional derivative seems to be the weighted and shifted Grünwald (WSGD) operator. The first order WSGD operator is firstly presented and detailedly discussed in [12, 13, 14] and the second order convergence is obtained by using extrapolation method [15, 16]. The second order WSGD operator is given in [18]; and the third order compact WSGD (CWSGD) is presented in [19]. Following the idea of weighting and shifting Grünwald operator, this paper provides the basic strategy of deriving the quasi-compact scheme with any desired convergence orders for space fractional diffusion equation; and it can also be extended to solve the tempered space fractional diffusion equation. The fourth order quasi-compact scheme is detailedly discussed in solving space fractional diffusion equation, including stability and convergence analysis and numerical verification of convergence orders. The fourth order quasi-compact scheme for tempered space fractional diffusion equation is also proposed and effectively used to solve the equation; and the convergence orders are numerically verified.

The outline of this paper is as follows. In Sec. 2, the high order quasi-compact discretizations are presented to approximate space Riemann-Liouville fractional derivative. In Sec. 3, following the obtained quasi-compact discretizations, the high order quasi-compact scheme for the one dimensional space fractional diffusion equation is designed and its stability and convergence analysis are performed. Sec. 4 focuses on the quasi-compact scheme and the corresponding stability and convergence analysis in two dimensional case. The high order quasi-compact discretizations is extended to tempered space fractional derivative in Sec. 5 and the corresponding scheme is derived to solve tempered space fractional diffusion equation. In Sec. 6, numerical experiments are performed to testify the efficiency and verify the convergence orders of the schemes. We conclude the paper with some discussions in the last section.

## 2 Quasi-compact discretizations for Riemann-Liouville space fractional derivatives

We first introduce some definitions and lemmas, including Riemann-Liouville fractional derivatives and shifted Grünwald-Letnikov formulations.

###### Definition 2.1

[21] If the function is defined in the interval and regular enough, then the -th order left and right Riemann-Liouville fractional derivatives are, respectively, defined as

and

 xDαbu(x)=(−1)nΓ(n−α)dndxn∫bx(s−x)n−α−1u(s)ds,n−1<α

where can be and can be .

And the standard left and right Grünwald-Letnikov formulations which can be potentially used to approximate the left and right Riemann-Liouville fractional derivatives are, respectively, given as

and

 xDαbu(x)=limh→01hα[b−xh]∑k=0g(α)ku(x+kh), \hb@xt@.01(2.4)

where the Grünwald weights are the coefficients of the power series expansion of . For getting the stable scheme, a shifted Grünwald-Letnikov operator is proposed to approximate the left Riemann-Liouville fractional derivative with first order accuracy [15].

###### Lemma 2.2 ([15])

Let , , and , . For any integer , define the left shifted Grünwald-Letnikov operator by

 Δαpu(x):=1hα∞∑k=0g(α)ku(x−(k−p)h). \hb@xt@.01(2.5)

Then we have

 Δαpu(x)=−∞Dαxu(x)+n−1∑l=1aαp,l−∞Dα+lxu(x)hl+O(hn), \hb@xt@.01(2.6)

uniformly in , where the weights are the coefficients of the power series expansion of the function , and the first four terms of the coefficients are , , , and .

To approximate the right Riemann-Liouville fractional derivative , the right shifted Grünwald-Letnikov operator is given by In the finite interval , the shifted Grünwald-Letnikov fractional derivatives are

 ~Δαpu(x)=1hα[x−ah]+p∑k=0g(α)ku(x−(k−p)h) \hb@xt@.01(2.7)

and

 ~Λαpu(x)=1hα[b−xh]+p∑k=0g(α)ku(x+(k−p)h). \hb@xt@.01(2.8)

In the remaining analysis of the paper, for a function defined in the bounded interval, we suppose that it has been zero extended to whenever the value of outside of the bounded interval is used.

### 2.1 Fourth order quasi-compact approximation to the Riemann-Liouville fractional derivative

According to the definitions of the shifted Grünwald-Letnikov fractional derivatives, we know that can be any integer. In order to ensure that the nodes in (LABEL:eq7) or (LABEL:eq8) are within the bounded interval, we need to choose the integer when approximating non-periodic fractional differential equation in the bounded interval. Inspired by the shifted Grünwald-Letnikov operator and the Taylor expansion, we derive the following fourth order combined quasi-compact approximations.

###### Theorem 2.3

Let and all the derivatives of up to order belong to . Then the following quasi-compact approximation has fourth order accuracy, i.e.,

 Px−∞Dαxu(x)=μ1Δα1u(x)+μ0Δα0u(x)+μ−1Δα−1u(x)+O(h4), \hb@xt@.01(2.9)

where , called CWSGD operator; is the centered difference operator; and the coefficients , , , and are the functions of and

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩μ1=(1+α)(2+α)/12,μ0=−(−2+α)(2+α)/6,μ−1=(−2+α)(−1+α)/12,bα2=(4+α−α2)/24. \hb@xt@.01(2.10)

In fact, under the assumptions of the theorem, we know that for any fixed order and the coefficients , , and the following equalities hold.

 μ1Δα1u(x)+μ0Δα0u(x)+μ−1Δα−1u(x)=−∞Dαxu(x)+bα2−∞Dα+2xu(x)h2+O(h4)=(1+h2bα2∂2∂x2)−∞Dxαu(x)+O(h4)=Px−∞Dxαu(x)+O(h4), \hb@xt@.01(2.11)

where . Then we get (LABEL:2.1). Since , we have for any function ,

 Pxu=(1+h2bα2∂2∂x2)u+O(h4).

In a similar way, we can derive the fourth order quasi-compact approximation for the right Riemann-Liouville fractional derivative:

 PxxDα+∞u(x)=μ1Λα1u(x)+μ0Λα0u(x)+μ−1Λα−1u(x)+O(h4). \hb@xt@.01(2.12)

For defined in a bounded interval, supposing its zero extension to satisfies the assumptions of Theorem LABEL:th2.1, the following approximations hold:

and

 PxxDαbu(x)=μ1~Λα1u(x)+μ0~Λα0u(x)+μ−1~Λα−1u(x)+O(h4). \hb@xt@.01(2.14)

Now using the CWSGD operator, we solve a two-point boundary value problem to numerically verify the above statements.

###### Example 2.4

Consider the steady state fractional diffusion problem

 0Dαxu(x)=720x6−αΓ(7−α),x∈(0,1),

with and the boundary conditions , . Its exact solution is .

Using the quasi-compact scheme (LABEL:fors11) to solve Example LABEL:exforth1 leads to the desired convergence orders; see Table LABEL:tab1.

### 2.2 Fifth order quasi-compact approximation to the Riemann-Liouville fractional derivative

In this subsection, we present a fifth order quasi-compact approximation given as follows.

###### Theorem 2.5

Let and all the derivatives of up to order belong to . Then the following quasi-compact approximation has fifth order accuracy, i.e.,

 P5x−∞Dαxu(x)=μ1Δα1f(x)+μ0Δα0f(x)+μ−1Δα−1f(x)+O(h5), \hb@xt@.01(2.15)

where , called 5-CWSGD operator, and

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩γ1=350+331α−15α2−75α3−15α41724−2α−570α2−30α3+30α4,γ2=566−329α−135α2+105α3−15α41724−2α−570α2−30α3+30α4,μ1=566+329α−135α2−105α3−15α41724−2α−570α2−30α3+30α4,μ0=862+α−285α2+15α3+15α4862−a−285α2−15α3+15α4,μ−1=350−331α−15α2+75α3−15α41724−2α−570α2−30α3+30α4. \hb@xt@.01(2.16)

The way of deriving (LABEL:fif) is similar to the derivation of the fourth order quasi-compact approximation. On one hand, from (LABEL:scheme0), we know for different parameter there exist three equalities

 Δαpu(x)=−∞Dαxu(x)+4∑k=1aαp,k−∞Dα+kxu(x)hk+O(h5),p=1,0,−1. \hb@xt@.01(2.17)

On the other hand, in view of the Taylor expansion we know

 −∞Dαxu(x−h)=−∞Dαxu(x)+(−1)k4∑k=11k!−∞Dα+kxu(x)hk+O(h5),−∞Dαxu(x+h)=−∞Dαxu(x)+4∑k=11k!−∞Dα+kxu(x)hk+O(h5). \hb@xt@.01(2.18)

So in order to get the fifth order quasi-compact approximation, combining (LABEL:2.2.1) and (LABEL:2.2.2), we need to eliminate the low order terms corresponding to (), which can be done by solving the algebraic equation

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩μ1+μ0+μ−1−γ1−γ2=1,μ1aα1,1+μ0aα0,1+μ−1aα−1,1+γ1−γ2=0,μ1aα1,2+μ0aα0,2+μ−1aα−1,2−γ1/2−γ2/2=0,μ1aα1,3+μ0aα0,3+μ−1aα−1,3+γ1/3!−γ2/3!=0,μ1aα1,4+μ0aα0,4+μ−1aα−1,4−γ1/4!−γ2/4!=0. \hb@xt@.01(2.19)

Eq. (LABEL:fifcoef) is the solution of (LABEL:fifcoef5). Then we get Theorem LABEL:th2.2. Next we utilize the 5-CWSGD operator to solve Example LABEL:ex2; and the numerical results are presented in Table LABEL:tab2, from which the accuracy of the 5-CWSGD operator is verified.

###### Example 2.6

We again consider the steady state fractional diffusion problem simulated in Example LABEL:exforth1, i.e.,

 0Dαxu(x)=720x6−αΓ(7−α),x∈(0,1),

with and the boundary conditions , ; and the exact solution .

###### Remark 2.7

As the fifth order quasi-compact scheme is not stable in solving the time-dependent space fractional differential equation, we detailedly discuss the fourth order quasi-compact schemes in Sections 3 and 4.

## 3 Quasi-compact scheme for one dimensional space fractional diffusion equation

Based on the fourth order quasi-compact discretization to the Riemann-Liouville space fractional derivative, we develop the Crank-Nicolson quasi-compact scheme of the two-sided space fractional diffusion equations. Here we consider the initial boundary value problem of the space fractional diffusion equation

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂u(x,t)∂t=K1aDαxu(x,t)+K2xDαbu(x,t)+f(x,t),(x,t)∈(a,b)×(0,T],u(x,0)=u0(x),x∈[a,b],u(a,t)=ϕa(t),u(b,t)=ϕb(t),t∈[0,T], \hb@xt@.01(3.1)

where . The diffusion coefficients and are nonnegative constants and they satisfy . If , then and , then . In the following analysis of the numerical method, we suppose that (LABEL:equ1) has an unique and sufficiently smooth solution.

### 3.1 CN-CWSGD scheme

The time interval is partitioned into a uniform mesh with the step size and the space interval into another uniform mesh with the step sized , where are two positive integers. Then the set of grid points can be denoted by and . Let , , and for . The maximum norm and the discrete norm are defined as

 ∥u∥∞=max1≤j≤M−1|uj|,∥u∥2=hM−1∑j=1u2j. \hb@xt@.01(3.2)

We use the Crank-Nicolson technique for the time discretization of (LABEL:equ1) and get

 \hb@xt@.01(3.3)

In space, the fourth order quasi-compact discretizations are used to approximate the Riemann-Liouville fractional derivatives. This implies that

 Pxun+1j−unjτ=K1τ2LDαhunj+K2τ2RDαhunj+K1τ2LDαhun+1j+K2τ2RDαhun+1j                    +Pxfn+1/2j+Rn+1/2j, \hb@xt@.01(3.4)

where

 LDαhunj=:μ1~Δα1unj+μ0~Δα0unj+μ−1~Δα−1unj=1hαj+1∑k=0w(α)kunj−k+1,
 RDαhunj=:μ1~Λα1unj+μ0~Λα0unj+μ−1~Λα−1unj=1hαM−j+1∑k=0w(α)kunj+k−1,

the coefficients , , and , and . Then the above equation can be rewritten as

 Pxun+1j−K1τ2LDαhun+1j−K2τ2RDαhun+1j=Pxunj+K1τ2LDαhunj+K2τ2RDαhunj+τPxfn+1/2j+τRn+1/2j. \hb@xt@.01(3.5)

Denoting as the numerical approximation of , we obtain the Crank-Nicolson quasi-compact scheme for (LABEL:equ1)

 PxUn+1j−K1τ2LDαhUn+1j−K2τ2RDαhUn+1j=PxUnj+K1τ2LDαhUnj+K2τ2RDαhUnj+τPxfn+1/2j. \hb@xt@.01(3.6)

For convenience, the approximation scheme (LABEL:fourthccd) can be written in matrix form

 (Pα−Bα)Un+1=(Pα+Bα)Un+τFn+Hn, \hb@xt@.01(3.7)

where , , is given by

 Aα=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝w(α)1w(α)0w(α)2w(α)1w(α)0⋮w(α)2w(α)1w(α)M−2⋯⋱⋱w(α)0w(α)M−1w(α)M−2⋯w(α)2w(α)1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, \hb@xt@.01(3.8)
 Pα=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝1−2bα2bα2bα2(1−2bα2)bα2⋯bα21−2bα2bα2bα21−2bα2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

and

 \hb@xt@.01(3.9)

### 3.2 Stability and convergence analysis

In this subsection, we prove that the CN quasi-compact scheme has fourth order accuracy in space and is unconditionally stable. Now we give some important lemmas to be used in the analyses.

###### Lemma 3.1

([23]) Let be a Toeplitz matrix with a generating function . Let and denote the smallest and largest eigenvalues of , respectively. Then we have

 fmin≤λmin(H)≤λmax(H)≤fmax,

where and denote the minimum and maximum values of , respectively. In particular, if and , then is negative definite.

###### Lemma 3.2

([24]) Let be a positive semi-definite matrix. Then there exists a unique -square positive semi-definite matrix such that . Such a matrix is called the square root of , denoted by .

###### Theorem 3.3

The matrix is negative definite, and is also negative definite, where is given by (LABEL:mA) and defined in (LABEL:matrixfourth).

In fact, the generating function [23] of satisfies

 f(α,x)=fAα(x)+fATα(x)=(∞∑k=0w(α)ke−i(k−1)x+∞∑k=0w(α)kei(k−1)x)=μ1(∞∑k=0g(α)ke−i(k−1)σ+∞∑k=0g(α)kei(k−1)σ)+μ0(∞∑k=0g(α)ke−ikσ+∞∑k=0g(α)keikσ)    +μ−1(∞∑k=0g(α)ke−i(k+1)σ+∞∑k=0g(α)kei(k+1)σ)=μ1((1−e−iσ)αeiσ+(1−eiσ)αe−iσ)+μ0((1−e−iσ)α+(1−eiσ)α)    +μ−1((1−e−iσ)αe−iσ+(1−eiσ)αeiσ)=(2sin(σ2))α(μ1(ei(απ2−ασ2+σ)+e−i(απ2−ασ2+σ))+μ0(ei(απ2−ασ2)+e−i(απ2−ασ2))    +μ−1(ei(απ2−ασ2−σ)+e−i(απ2−ασ2−σ)))=2(2sin(x2))α(μ1cos(απ2−αx2+x)+μ0cos(απ2−αx2)+μ−1cos(απ2−αx2−x)), \hb@xt@.01(3.10)

where and denote the generating functions of the matrix and , respectively. Since is a real-valued and even function, it’s reasonable to consider its principal value on .

Together with Fig. LABEL:fig1, we have that for on . Then from Lemma LABEL:sy1, we know the matrix is negative definite. Rewriting as , it can be clearly seen that is negative definite.

###### Theorem 3.4

The difference scheme (LABEL:fourthccd) with is unconditionally stable.

Proof. Define the round-off error as , where is the exact solution of the discretized equation (LABEL:fourthccd) and is the numerical solution of the discretized equation (LABEL:fourthccd) obtained in finite precision arithmetic. Since satisfies the discretized equation exactly, round-off error must also satisfy the discretized equation [22]. Thus we obtain the following error equation

 Pxϵn+1j−K1τ2LDαhϵn+1j−K2τ2RDαhϵn+1j=Pxϵnj+K1τ2LDαhϵnj+K2τ2RDαhϵnj. \hb@xt@.01(3.11)

Since the boundary conditions of error equation (LABEL:errorequation) are , we zero extend the solution of the problem (LABEL:errorequation) to the whole real line . So it’s reasonable to replace the symbols and in error equation (LABEL:errorequation) with . Now we have