High order schemes for the tempered fractional diffusion equations

# High order schemes for the tempered fractional diffusion equations

Department of Applied Mathematics, School of Sciences, Xi’an University of Technology,
Xi’an, Shaanxi 710054, P.R. China.
Beijing Computational Science Research Center, Beijing 100084, P.R. China.
School of Mathematics and Statistics, Gansu Key Laboratory of Applied Mathematics and Complex Systems,
Lanzhou University, Lanzhou 730000, P.R. China.
###### Abstract

Lévy flight models whose jumps have infinite moments are mathematically used to describe the superdiffusion in complex systems. Exponentially tempering Lévy measure of Lévy flights leads to the tempered stable Lévy processes which combine both the -stable and Gaussian trends; and the very large jumps are unlikely and all their moments exist. The probability density functions of the tempered stable Lévy processes solve the tempered fractional diffusion equation. This paper focuses on designing the high order difference schemes for the tempered fractional diffusion equation on bounded domain. The high order difference approximations, called the tempered and weighted and shifted Grünwald difference (tempered-WSGD) operators, in space are obtained by using the properties of the tempered fractional calculus and weighting and shifting their first order Grünwald type difference approximations. And the Crank-Nicolson discretization is used in the time direction. The stability and convergence of the presented numerical schemes are established; and the numerical experiments are performed to confirm the theoretical results and testify the effectiveness of the schemes.

¡¡¡¡

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

Key words: tempered fractional calculus, tempered-WSGD operator, superconvergent, stability and convergence

## 1 Introduction

The probability density function of Lévy flights [20, 23] has a characteristic function of stretched Gaussian form, causing the asymptotic decay as . It produces that the second moment diverges, i.e., . The divergent second moments may not be feasible for some even non-Brownian physical processes of practical interest which take place in bounded domains and involve observables with finite moments. To overcome the divergence of the variance, many techniques are adopted. By simply discarding the very large jumps, Mantegna and Stanley  introduce the truncated Lévy flights and show that the obtained stochastic process ultraslowly converges to a Gaussian. From the point of view of an experimental study, because of the limited time, no expected Gaussian behavior can be observed. Two other modifications to achieve finite second moments are proposed by Sokolov et al , who add a high-order power-law factor, and Chechkin et al. , who add a nonlinear friction term. Exponentially tempering the probability of large jumps of Lévy flights, i.e., making the Lévy density decay as with , to get finite moments seems to be the most popular one; and the corresponding tempered fractional differential equations are derived [1, 4, 5, 6]. In fact, the power-law waiting time can also be exponentially tempered . This paper focused on providing high order schemes for numerically solving tempered fractional diffusion equations, which involve tempered fractional derivatives. In fact, the tempered fractional integral has a long history. Buschman’s earlier work  reports the fractional integration with weak singular and exponential kernels; for more detailed discussions, see Srivastava and Buschman’s book  and the references therein. The definitions of tempered fractional calculus are much similar to the ones of fractional substantial calculus ; but they are introduced from completely different physical backgrounds; e.g., fractional substantial calculus is used to characterize the functional distribution of anomalous diffusion . Mathematically the fractional substantial calculus is time-space coupled operator but the tempered fractional calculus is not; numerically the fractional substantial calculus is discretized in the time direction , but here the tempered fractional calculus is treated as space operator.

Tempered fractional calculus is the generalization of fractional calculus, or fractional calculus is the special/limiting case of tempered fractional calculus. Some important progresses have been made for numerically solving the fractional partial differential equations (PDEs), e.g., the finite difference methods are used to simulate the space fractional advection diffusion equations [17, 22, 18]. Recently, it seems that more efforts of the researchers are put on the high order schemes and fast algorithms. Based on the Toeplitz-like structure of the matrix corresponding the finite difference methods of fractional PDEs, Wang et al  numerically solve the fractional diffusion equations with the computational cost. Later, Pang and Sun  propose a multigrid method to solve the discretized system of the fractional diffusion equation. By introducing the linear spline approximation, Sousa and Li present a second order discretization for the Riemann-Liouville fractional derivatives, and establish an unconditionally stable weighted finite difference method for the one-dimensional fractional diffusion equation in . Ortigueira  gives the “fractional centred derivative” to approximate the Riesz fractional derivative with second order accuracy; and this method is used by Çelik and Duman in  to approximate fractional diffusion equation with the Riesz fractional derivative in a finite domain. More recently, by weighting and shifting the Grünwald discretizations, Tian et al  propose a class of second order difference approximations, called WSGD operators, to the Riemann-Liouville fractional derivatives.

So far, there are limited works addressing the finite difference schemes for the tempered fractional diffusion equations. Baeumera and Meerschaert  provide finite difference and particle tracking methods for solving the tempered fractional diffusion equation with the second order accuracy. The stability and convergence of the provided schemes are discussed. Cartea and del-Castillo-Negrete  derive a general finite difference scheme to numerically solve a Black-Merton-Scholes model with tempered fractional derivatives. Recently, Marom and Momoniat  compare the numerical solutions of three kinds of fractional Black-Merton-Scholes equations with tempered fractional derivatives. And the stability and convergence of the presented schemes are not given. To the best of our knowledge, there is no published work to the high order difference schemes for the tempered fractional diffusion equation. In this paper, with the similar method presented in [22, 1], we first propose the first order shifted Grünwald type approximation for the tempered fractional calculus; then motivated by the idea in , we design a series of high order schemes, called the tempered-WSGD operators, by weighting and shifting the first order Grünwald type approximations to the tempered fractional calculus. The obtained high order schemes are applied to solve the tempered fractional diffusion equation and the Crank-Nicolson discretization is used in the time direction. The unconditionally numerical stability and convergence are detailedly discussed; and the corresponding numerical experiments are carried out to illustrate the effectiveness of the schemes.

The remainder of the paper is organized as follows. In Sec. 2, we introduce the definitions of the tempered fractional calculus and derive their first order shifted Grünwald type approximations and the high order discretizations, the tempered-WSGD operators. In Sec. 3, the tempered fractional diffusion equation is numerically solved by using the tempered-WSGD operators to approximate the space derivative and the Crank-Nicolson discretization to the time derivative; and the numerical stability and convergence are discussed. The effectiveness and convergence orders of the presented schemes are numerically verified in Sec. 4. And the concluding remarks are given in the last section.

## 2 Definitions of the tempered fractional calculus and the derivation of the tempered-WSGD operators

We first introduce the definitions of the tempered fractional integral and derivative then focus on deriving their high order discretizations, the tempered-WSGD operators.

### 2.1 Definitions and Fourier transforms of the tempered fractional calculus

We introduce the definitions of the tempered fractional calculus and perform their Fourier transforms.

###### Definition 1 ([2, 5]).

Let be piecewise continuous on or corresponding to the right integral and integrable on any finite subinterval of or corresponding to the right integral, , . Then

• the left Riemann-Liouville tempered fractional integral of order is defined to be

• the right Riemann-Liouville tempered fractional integral of order is defined to be

 xD−σ,λbu(x)=1Γ(σ)∫bxe−λ(ξ−x)(ξ−x)σ−1u(ξ)dξ.
###### Definition 2 ([13, 15, 28]).

For , let be -times continuously differentiable on or corresponding to the right derivative and its -times derivative be integrable on any subinterval of or corresponding to the right derivative. Then

• the left Riemann-Liouville fractional derivative:

• the right Riemann-Liouville fractional derivative:

 xDαbu(x)=(−1)nΓ(n−α)dndxn∫bxu(ξ)(ξ−x)α−n+1dξ.
###### Definition 3 ([1, 5]).

For , let be -times continuously differentiable on or corresponding to the right derivative and its -times derivative be integrable on any subinterval of or corresponding to the right derivative, . Then

• the left Riemann-Liouville tempered fractional derivative:

• the right Riemann-Liouville tempered fractional derivative:

 xDα,λbu(x)=eλxxDαb(e−λxu(x))=(−1)neλxΓ(n−α)dndxn∫bxe−λξu(ξ)(ξ−x)α−n+1dξ.

If , then the left and right Riemann-Liouville tempered fractional derivatives and reduce to the left and right Riemann-Liouville fractional derivatives and defined in Definition 2.

###### Definition 4.

The variants of the left and right Riemann-Liouville tempered fractional derivatives are defined as [1, 5, 24]

and

 xDα,λbu(x)={xDα,λbu(x)−λαu(x),0<α<1,xDα,λbu(x)+αλα−1∂xu(x)−λαu(x),1<α<2, (2.2)

where denotes the classic first order derivative .

###### Remark 1.

In the above definitions, the ‘a’ can be extended to ‘’ and ‘b’ to ‘’. In the following analysis, we assume that is defined on and whenever necessary can be smoothly zero extended to or or even . Then ; ; ; and .

###### Lemma 1 ([2, 16, 1]).

Let and its -times derivative belong to , . Then the Fourier transforms of the left and right Riemann-Liouville tempered fractional integrals are

 F(−∞D−σ,λxu(x))=(λ+iω)−σ^u(ω); (2.3)

and

 F(xD−σ,λ+∞u(x))=(λ−iω)−σ^u(ω) (2.4)

and the Fourier transforms of the left and right Riemann-Liouville tempered fractional derivatives are

 F(−∞Dα,λxu(x))=(λ+iω)α^u(ω); (2.5)

and

 F(xDα,λ+∞u(x))=(λ−iω)α^u(ω) (2.6)

and the Fourier transforms of the variants of the left and right Riemann-Liouville tempered fractional derivatives give

 F(−∞Dα,λxu(x))={(λ+iω)α^u(ω)−λα^u(ω),0<α<1,(λ+iω)α^u(ω)−αiωλα−1^u(ω)−λα^u(ω),1<α<2; (2.7)

and

 F(xDα,λ+∞u(x))={(λ−iω)α^u(ω)−λα^u(ω),0<α<1,(λ−iω)α^u(ω)+αiωλα−1^u(ω)−λα^u(ω),1<α<2, (2.8)

where the Fourier transform of is defined by

 F(u(x))(ω)=∫Re−iωxu(x)dx,i2=−1.
###### Remark 2 ([3, 10]).

The left and right Riemann-Liouville tempered fractional derivatives can be, respectively, rewritten as

 −∞Dα,λxu(x)=1Γ(n−α)(ddx+λ)n∫x−∞e−λ(x−ξ)u(ξ)(x−ξ)α−n+1dξ;

and

 xDα,λ+∞u(x)=(−1)nΓ(n−α)(ddx−λ)n∫+∞xe−λ(ξ−x)u(ξ)(ξ−x)α−n+1dξ.

### 2.2 Discretizations of the tempered fractional calculus

In this subsection, we derive the Grünwald type discretizations for the tempered fractional calculus. The standard Grünwald discretization generally yields an unstable finite difference scheme when it is used to solve the time dependent fractional PDEs . To remedy this defect, Meerschaert et al introduce a shifted Grünwald formula. The similar numerical unstability also happens for the time dependent tempered fractional PDEs; so the shift for the Grünwald type discretizations of the tempered fractional derivative is also necessary.

###### Lemma 2.

Let , and its Fourier transform belong to ; and . Defining the shifted Grünwald type difference operator

 Aα,λh,pu(x):=1hα+∞∑k=0w(α)ke−(k−p)hλu(x−(k−p)h)−1hα(ephλ(1−e−hλ)α)u(x), (2.9)

then

 Aα,λh,pu(x)=−∞Dα,λxu(x)−λαu(x)+O(h), (2.10)

where denotes the normalized Grünwald weights.

###### Remark 3.

The point is the superconvergent point of the approximation to , i.e., with (the deriving process is similar to the one given in ).

###### Remark 4.

Under the assumption given in Lemma 2, for tempered fractional derivatives defined in (2.7), we have 

 Aα,λh,pu(x)={−∞Dα,λxu(x)+O(h),0<α<1,−∞Dα,λxu(x)+αλα−1∂xu(x)+O(h),1<α<2. (2.11)
###### Proof.

The proof is similar to the one given in . Taking Fourier transform on both sides of (2.9), we obtain

 F[Aα,λh,pu](ω) =1hα+∞∑k=0w(α)ke−(k−p)h(λ+iω)^u(ω)−1hα(ephλ(1−e−hλ)α)^u(ω) (2.12) =eph(λ+iω)(1−e−h(λ+iω)h)α^u(ω)−ephλ(1−e−hλh)α^u(ω) =[(λ+iω)αPh(λ+iω)−λαPh(λ)]^u(ω),

where

 Ph(z)=ephz(1−e−hzhz)α=1+(p−α2)hz+O(|z|2),withz=λ+iωorλ. (2.13)

Denoting

 ^ϕ(ω,h)=F[Aα,λh,pu](ω)−F[−∞Dα,λxu−λαu](ω)=[(λ+iω)α(Ph(λ+iω)−1)−λα(Ph(λ)−1)]^u(ω),

from (2.12) and (2.5) there exists

 |^ϕ(ω,h)|≤C[h|(λ+iω)|α+1+h|λ|α+1]|^u(ω)|.

With the condition , and using the Riemann-Lebesgue Lemma, it yields

 |Aα,λh,pu(x)− −∞Dα,λxu(x)+λαu(x)| =|ϕ|≤12π∫R|^ϕ(ω,h)|dω ≤C∥F[−∞Dα+1,λxu+λα+1u(x)](ω)∥L1h=O(h),

where the property of the Fourier transforms for the left Riemann-Liouville tempered fractional derivatives (2.5) is used. ∎

###### Lemma 3.

Let , and its Fourier transform belong to ; and . Define the tempered shifted Grünwald type difference operator

 Bα,λh,pu(x):=1hα+∞∑k=0w(α)ke−(k−p)hλu(x+(k−p)h)−1hα(ephλ(1−e−hλ)α)u(x). (2.14)

Then

 Bα,λh,qu(x)=xDα,λ+∞u(x)−λαu(x)+O(h). (2.15)
###### Remark 5.

The point is the superconvergent point of the approximation to , i.e., with (the deriving process is similar to the one given in ).

###### Remark 6.

Under the assumption given in Lemma 3, for tempered fractional derivatives defined in (2.8), we have

 Bα,λh,qu(x)={xDα,λ+∞u(x)+O(h),0<α<1,xDα,λ+∞u(x)−αλα−1∂xu(x)+O(h),1<α<2. (2.16)
###### Proof.

Taking Fourier transform on both sides of (2.14), we obtain

 F[Bα,λh,pu](ω) =1hα+∞∑k=0w(α)ke−(k−p)h(λ−iω)^u(ω)−1hα(ephλ(1−e−hλ)α)^u(ω) =eph(λ−iω)(1−e−h(λ−iω)h)α^u(ω)−ephλ(1−e−hλh)α^u(ω) =[(λ−iω)αPh(λ−iω)−λαPh(λ)]^u(ω),

where is defined by (2.13) with or . Denoting , then with the similar method used in the proof of Lemma 2, and using the Fourier transform of the right Riemann-Liouville tempered fractional derivative (2.6), we obtain

 |Bα,λh,pu(x)−xDα,λ+∞u(x)+λαu(x)| =|ϕ|≤12π∫R|^ϕ(ω,h)|dω ≤C∥F[xDα+1,λ+∞u+λα+1u(x)](ω)∥L1h=O(h).

The approximation accuracy of the classic difference operator can be improved by adding the band of discretization stencils . And then the computational cost increases accordingly. However, because of the nonlocal property of the fractional operator, even for the first order discretizations, the stencil covers the whole interval. Without introducing new computational cost, we can improve the approximation accuracy of the discretized fractional operators by modifying the Grünwald type weights. The improved discretized tempered fractional operators are called tempered weighted and shifted Grünwald difference (tempered-WSGD) operators.

###### Theorem 4.

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

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

where and are determined by (2.21)-(2.24). Then, for any integer , there exists

 LDα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)=−∞Dα,λxu(x)−λαu(x)+O(hℓ), (2.18)

uniformly for .

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

 RDα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)=m∑j=1γjBα,λh,pju(x), (2.19)

where and are determined by (2.21)-(2.24). Then, for any integer , there is

 RDα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)=xDα,λ+∞u(x)−λαu(x)+O(hℓ), (2.20)

uniformly for .

• For , are real numbers and satisfy the linear system

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩m∑j=1γj=1,m∑j=1γj[pj−α2]=0. (2.21)
• For , are real numbers and satisfy

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩m∑j=1γj=1,m∑j=1γj[pj−α2]=0,m∑j=1γj[p2j2−αpj2+α6+α(α−1)8]=0. (2.22)
• For , are real numbers and the following hold

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩m∑j=1γj=1,m∑j=1γj[pj−α2]=0,m∑j=1γj[p2j2−αpj2+α6+α(α−1)8]=0,m∑j=1γj[p3j6−αp2j4+(α6+α(α−1)8)pj−α24−α(α−1)12−α(α−1)(α−2)48]=0. (2.23)
• For , are real numbers and the following hold

 (2.24)
###### Proof.

The standard Fourier transforms are again used here. Performing the Fourier transform on the left hand of (2.17), we obtain

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

where or , By a simple Taylor’s expansion, we get

 epjhz(1−e−hzhz)α= (2.26) +[p3j6−αp2j4+(α6+α(α−1)8)pj−α24−α(α−1)12−α(α−1)(α−2)48](hz)3 +[p4j24−αp3j4+12(α6+α(α−1)8))p2j+(−α24−α(α−1)12−α(α−1)(α−2)48)pj +O(|zh|5).

Denoting , in view of (2.26), (2.5), and (2.21)-(2.24), we have

 |^ϕ(ω,h)|≤Chl[|λ+iω|α+ℓ+|λ|α+ℓ]|^u(ω)|. (2.27)

Due to , there exists

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

By the similar arguments we can prove (2.20). ∎

###### Remark 7.

Under the assumptions given by Theorem 4, for the tempered fractional derivatives defined in (2.7) and (2.8), we deduce that

 LDα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)={−∞Dα,λxu(x)+O(hℓ),,0<α<1,−∞Dα,λxu(x)+αλα−1∂xu(x)+O(hℓ),1<α<2; (2.28)

and

 RDα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)={xDα,λ+∞u(x)+O(hℓ),0<α<1,xDα,λ+∞u(x)−αλα−1∂xu(x)+O(hℓ),1<α<2. (2.29)
###### Remark 8.

If , and its Fourier transform belong to ; and . Defining the shifted Grünwald type difference operator

 ~Aα,λh,pu(x):=1hα+∞∑k=0w(α)ke−(k−p)hλu(x−(k−p)h), (2.30)

then

 ~Aα,λh,pu(x)=−∞Dα,λxu(x)+O(h), (2.31)

where denotes the normalized Grünwald weights.

If , and its Fourier transform belong to ; and . Define the shifted Grünwald type difference operator

 ~Bα,λh,pu(x):=1hα+∞∑k=0w(α)ke−(k−p)hλu(x+(k−p)h). (2.32)

Then

 ~Bα,λh,qu(x)=xDα,λ+∞u(x)+O(h). (2.33)

Moreover, if , and its Fourier transform belong to ; and define the left tempered-WSGD operator by

 L~Dα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)=m∑j=1γj~Aα,λh,pju(x), (2.34)

where and are determined by (2.21)-(2.24). Then, for any integer , there exists

 L~Dα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)=−∞Dα,λxu(x)+O(hℓ), (2.35)

uniformly for .

If , and its Fourier transform belong to ; and define the right tempered-WSGD operator by

 R~Dα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)=m∑j=1γj~Bα,λh,pju(x), (2.36)

where and are determined by (2.21)-(2.24). Then, for any integer , there is

 R~Dα,γ1,γ2,…,γmh,p1,p2,…,pmu(x)=xDα,λ+∞u(x)+O(hℓ), (2.37)

uniformly for .

###### Remark 9.

To get the discretizations, including the first and high orders, of the left and right Riemann-Liouville tempered fractional integrals of order :