Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation

# Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation

Zhibo Wang Corresponding author. Email: zhibowangok@gmail.com. Department of Mathematics, University of Macau, Av. Padre Tomás Pereira Taipa, Macau, China    Seakweng Vong Email: swvong@umac.mo. Department of Mathematics, University of Macau, Av. Padre Tomás Pereira Taipa, Macau, China
###### Abstract

In this paper, compact finite difference schemes for the modified anomalous fractional sub-diffusion equation and fractional diffusion-wave equation are studied. Schemes proposed previously can at most achieve temporal accuracy of order which depends on the order of fractional derivatives in the equations and is usually less than two. Based on the idea of weighted and shifted Grünwald difference operator, we establish schemes with temporal and spatial accuracy order equal to two and four respectively.

Keywords: Modified anomalous fractional sub-diffusion equation, Fractional diffusion-wave equation, Compact difference scheme, Weighted and shifted Grünwald difference operator

## 1 Introduction

Since fractional differential equations turn out to model many physical processes more accurately than the classical ones, in the past decades, increasing attentions have been made on these equations. Readers can refer to the books [1, 2] for theoretical results on fractional differential equations. This paper concerns with methods for obtaining accurate numerical approximations to the solutions of fractional sub-diffusion equations and fractional diffusion-wave equations. A fractional sub-diffusion equation is an integro-partial differential equation obtained from the classical diffusion equation by replacing the first-order time derivative by a fractional derivative of order between zero and one. When the time derivative is of order between one and two, we get a fractional diffusion-wave equation. Fractional derivatives of order between zero and one are widely used in describing anomalous diffusion processes [3], while fractional diffusion-wave equations have applications in modeling universal electromagnetic, acoustic, and mechanical responses [4, 5].

Numerical methods for the modified anomalous fractional sub-diffusion equation and the diffusion-wave equation have been considered by many authors, one may refer to [6][17] and the references therein. We point out here that one of the main tasks for developing accurate finite difference scheme of fractional differential equations is to discretize the fractional derivatives. Noticing that fractional derivatives are defined through integrals, one can approximate the derivatives by interpolating polynomials [18]. Using this idea, Sousa and Li developed a second order discretization for the Riemann-Liouville fractional derivative [19].

We note that there are alternative ways to tackle the problem. In [20], a shifted Grünwald formula was proposed by Meerschaert and Tadjeran to approximate fractional derivatives of order for fractional advection-dispersion flow equations. It is also worth to mention that stability of forward-Euler scheme and weighted averaged difference scheme based on Grünwald-Letnikov approximation were analyzed in [21, 22] for time fractional diffusion equations. Very recently, accurate finite difference schemes based on weighted and shifted Grünwald difference operator were developed for solving space fractional diffusion equations in [23, 24].

Inspired by their work on the weighted and shifted Grünwald difference operator, in this paper, we establish high order schemes for the fractional diffusion-wave equation and the modified anomalous fractional sub-diffusion equation, which were proposed and studied in [8, 9] and [25][27], respectively. We remark that the order of temporal accuracy of schemes proposed previously can at most be a fraction depending on the fractional derivatives in the equations and is usually less than two. We show that the schemes proposed in this paper are of order , where and are the temporal and spatial step sizes respectively.

This paper is organized as follows. Some preliminaries will be given in the next section. Compact schemes are proposed and studied in Section 3 and 4 for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation respectively. In the last section, numerical tests are carried out to justify the theoretical results.

## 2 Preliminaries

We first recall that the Caputo fractional derivative of order for a function is defined as

with being the gamma function and, for , the Riemann-Liouville fractional derivative of order for is defined as

Closely related to the fractional derivatives of a function is the Riemann-Liouville fractional integral which is given by

 aIαtf(t)=1Γ(α)∫taf(s)(t−s)1−αds.

In order to develop second order approximation of the Riemann-Liouville fractional derivative, we consider the shifted Grünwald approximation [20] to the Riemann-Liouville fractional derivative given by

 Aατ,rf(t)=τ−α∞∑k=0g(α)kf(t−(k−r)τ),

where for . Inspired by [28], we similarly introduce the shifted operator to the Riemann-Liouville fractional integral defined by

 Bατ,rf(t)=τα∞∑k=0ω(α)kf(t−(k−r)τ),

where for .

The following lemma was given in [29].

###### Lemma 2.1

Suppose . The Fourier transform of the Riemann-Liouville fractional integral satisfy the following:

 F[−∞Iαtf(t)]=(iω)−α^f(ω),

where denotes the Fourier transform of .

Lemma 2.1 can be used to obtain

###### Lemma 2.2

(i) Let and its Fourier transform belong to , and define the weighted and shifted Grünwald difference operator by

 Dατ,p,qf(t)=α−2q2(p−q)Aατ,pf(t)+2p−α2(p−q)Aατ,qf(t), (1)

then we have

 Dατ,p,qf(t)=−∞Dαtf(t)+O(τ2)

for , where and are integers and .

(ii) Let and belong to . Define the weighted and shifted difference operator by

 Iατ,p,qf(t)=2q+α2(q−p)Bατ,pf(t)+2p+α2(p−q)Bατ,qf(t),

then we have

 Iατ,p,qf(t)=−∞Iαtf(t)+O(τ2)

for , where and are integers and .

Proof. The proof of (i) can be found in [23]. The proof of (ii) is similar to that of (i) but we give it here for the completeness of our presentation.

Referring to the definition of , we let

 Iατ,p,qf(t)=τα[μ1∞∑k=0ω(α)kf(t−(k−p)τ)+μ2∞∑k=0ω(α)kf(t−(k−q)τ)]. (2)

Taking Fourior transform on (2), we get

 F[Iατ,p,qf](ω)=τα∞∑k=0ω(α)k[μ1e−iω(k−p)τ+μ2e−iω(k−q)τ]F[f](ω)=τα[μ1(1−e−iωτ)−αeiωτp+μ2(1−e−iωτ)−αeiωτq]F[f](ω)=(iω)−α[μ1Wp(iωτ)+μ2Wq(iωτ)]F[f](ω), (3)

where

 (4)

In order to achieve second order accuracy, we let the coefficients and satisfy the following system:

 {μ1+μ2=1,(p+α2)μ1+(q+α2)μ2=0,

which implies that and .

Denote , then by (3), (4) and Lemma 2.1, we have

 |^g(ω,τ)|≤Cτ2|iω|2−α|F[f](ω)|.

Thus

 |Iατ,p,qf−−∞Iαtf|=|g|≤12π∫R|^g(ω,τ)|dω≤C∥(iω)2−αF[f](ω)∥L1τ2=O(τ2).□

We are now ready to establish our high order compact schemes in the next two sections.

## 3 A compact scheme for the fractional sub-diffusion equation

Consider the following modified anomalous fractional sub-diffusion equation

 ∂u(x,t)∂t=(κ1∂α∂tα+κ2∂β∂tβ)[∂2u(x,t)∂x2]+f(x,t),0≤x≤L,0

subject to

 u(x,0)=0,0≤x≤L,u(0,t)=φ1(t),u(L,t)=φ2(t),0

where , . We have used and to denote the Riemann-Liouville fractional operators and with respect to the time variable .

###### Remark 3.1

(i) We note that, in [25, 26, 27], the fractional derivatives in the equation are of order and with . We have changed the notations here in order to match the presentation for the two types of equations discussed in this paper.

(ii) Without loss of generality, we have assumed the initial condition . If , one may consider the equation for instead.

To develop a finite difference scheme for the problem (5)–(6), we let and be the spatial and temporal step sizes respectively, where and are some given positive integers. For and , denote . For any grid function , we introduce the following notations:

 δxuki−12=1h(uki−uki−1),  δ2xuki=1h(δxuki+12−δxuki−12),
 Hui=⎧⎪⎨⎪⎩(1+h212δ2x)ui=112(ui−1+10ui+ui+1),1≤i≤M−1,ui,i=0 or M.

With this notations, we study the problem under the following inner product and norms:

 ⟨u,v⟩=hM−1∑i=1uivi,  ⟨δxu,δxv⟩=hM−1∑i=0(δxui+12)(δxvi+12),
 ∥u∥2=⟨u,u⟩,  ∥u∥∞=max0≤i≤M|ui|.

It is easy to check that, if , the following identity holds

 ⟨δ2xu,v⟩=−⟨δxu,δxv⟩,

which plays important role in our analysis.

Note that one can continuously extend the solution to be zero for . By choosing in (i) of Lemma 2.2, we get in (1), which gives

 ∂α∂tα[uxx(xi,tn+1)]=τ−α(2+α2n+1∑k=0g(α)kδ2xun+1−ki−α2n∑k=0g(α)kδ2xun−ki)+O(τ2+h2)=τ−αn+1∑k=0λ(α)kδ2xun+1−ki+O(τ2+h2),

where

 λ(α)0=2+α2g(α)0,  λ(α)k=2+α2g(α)k−α2g(α)k−1,  k≥1. (7)

We can therefore consider an weighted Crank-Nicolson type discretization for equation (5) given by

 un+1i−uniτ=κ1τ−α2(n+1∑k=0λ(α)kδ2xun+1−ki+n∑k=0λ(α)kδ2xun−ki)+κ2τ−β2(n+1∑k=0λ(β)kδ2xun+1−ki+n∑k=0λ(β)kδ2xun−ki)+12(fni+fn+1i).

In order to raise the accuracy in the spatial direction, we need the following lemma:

###### Lemma 3.1

([30]) Denote . If , then it holds that

 112[f′′(xi−1)+10f′′(xi)+f′′(xi+1)]=1h2[f(xi−1)−2f(xi)+f(xi+1)]+h4360∫10[f(6)(xi−sh)+f(6)(xi+sh)]ζ(s)ds.

Based on Lemma 3.1, we therefore propose the following compact scheme:

 H(un+1i−uni)=κ1τ1−α2(n+1∑k=0λ(α)kδ2xun+1−ki+n∑k=0λ(α)kδ2xun−ki) 0≤n≤N−1,  1≤i≤M−1, (8) un0=φn1,unM=φn2,1≤n≤N, u0i=0,0≤i≤M. (9)

It is easy to see that at each time level, the difference scheme (3)–(9) is a linear tridiagonal system with strictly diagonal dominant coefficient matrix, thus the difference scheme has a unique solution.

The following lemmas are critical for establishing the convergence of the proposed scheme.

###### Lemma 3.2

Let be defined as in (7), then for any positive integer and real vector , it holds that

 k−1∑n=0(n∑p=0λ(α)pvn+1−p)vn+1≥0.

Proof. For simplicity of presentation, in this proof, we denote without ambiguity. One can easily check that, to prove the above quadratic form is nonnegative is equivalent to proving the symmetric Toeplitz matrix T is positive semi-definite, where

 T=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝λ0λ12λ22⋯λk−12λ12λ0λ12⋱⋮λ22λ12⋱⋱λ22⋮⋱⋱λ0λ12λk−12⋯λ22λ12λ0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

Notice that the generating function (see [31]) of T is given by

 f(α,x)=λ0+12∞∑k=1λkeikx+12∞∑k=1λke−ikx=2+α2g0+12∞∑k=1(2+α2gk−α2gk−1)eikx+12∞∑k=1(2+α2gk−α2gk−1)e−ikx=2+α4(1−eix)α+2+α4(1−e−ix)α−α4eix(1−eix)α−α4e−ix(1−e−ix)α.

As mentioned in [23], we only need to consider the principal value of on which gives

 f(α,x)=2+α4[2isin(−x2)eix2]α+2+α4[2isin(x2)e−ix2]α−α4eix[2isin(−x2)eix2]α−α4e−ix[2isin(x2)e−ix2]α=[2sin(x2)]α{2+α2cos[α2(π−x)]−α2cos[α2(π−x)−x]}.

Let , then one can easily check that

 hx(α,x)=α(2+α)2sin(x2)cos[α2(π−x)−x2]≥0.

Therefore is nondecreasing with respect to and , which implies that . The lemma now follows as a result of the Grenander-Szegö Theorem [31].

###### Remark 3.2

We note here that the function can not be proved to be nonnegative by considering differentiation with respect to as in [23].

###### Lemma 3.3

(Grownall’s inequality [32]) Assume that and are nonnegative sequences, and the sequence satisfies

 ϕ0≤g0,   ϕn≤g0+n−1∑l=0pl+n−1∑l=0klϕl,   n≥1,

where . Then the sequence satisfies

 ϕn≤(g0+n−1∑l=0pl)exp(n−1∑l=0kl),   n≥1.

With all the preparation, we can now show the convergence and stability of the compact finite difference scheme (3)–(9).

###### Theorem 3.1

Assume that is the solution of (5)–(6) and is the solution of the finite difference scheme (3)–(9), respectively. Denote

 eki=u(xi,tk)−uki,0≤i≤M,0≤k≤N.

Then there exists a positive constant such that

 ∥ek∥≤~c1(τ2+h4),0≤k≤N.

Proof. We can easily get the error equation:

 H(ek+1i−eki)=κ1τ1−α2k∑l=0λ(α)lδ2x(ek+1−li+ek−li)+κ2τ1−β2k∑l=0λ(β)lδ2x(ek+1−li+ek−li)+τRk+1i0≤k≤N−1,  1≤i≤M−1,ek0=ekM=0,1≤k≤Ne0i=0,0≤i≤M, (10)

where .

Denote , multiplying (10) by and summing in , we obtain

 h(ek+1+ek)TC(ek+1−ek)=κ1τ1−α2k∑l=0λ(α)l⟨δ2x(ek+1−l+ek−l),ek+1+ek⟩+κ2τ1−β2k∑l=0λ(β)l⟨δ2x(ek+1−l+ek−l),ek+1+ek⟩+τh(ek+1+ek)TRk+1=−κ1τ1−α2k∑l=0λ(α)l⟨δx(ek+1−l+ek−l),δx(ek+1+ek)⟩−κ2τ1−β2k∑l=0λ(β)l⟨δx(ek+1−l+ek−l),δx(ek+1+ek)⟩+τh(ek+1+ek)TRk+1. (11)

Summing up for and noticing that

 h(ek+1+ek)TC(ek+1−ek)=h(ek+1TCek+1−ekTCek),
 henTCen≥23∥en∥2,

and

 τh(ek+1+ek)TRk+1≤τ3(∥ek+1∥2+∥ek∥2)+3τ2∥Rk+1∥2,

we get, by Lemma 3.2, that

 23∥en∥2≤−κ1τ1−α2n−1∑k=0k∑l=0λ(α)l⟨δx(ek+1−l+ek−l),δx(ek+1+ek)⟩−κ2τ1−β2n−1∑k=0k∑l=0λ(β)l⟨δx(ek+1−l+ek−l),δx(ek+1+ek)⟩+τ3n−2∑k=0(∥ek+1∥2+∥ek∥2)+3τ2n−2∑k=0∥Rk+1∥2+τh(en+en−1)TRn≤13∥en∥2+τ3∥en−1∥2+3τ24∥Rn∥2+3τ4∥Rn∥2+τ3n−1∑k=1∥ek∥2+τ3n−2∑k=1∥ek∥2+3τ2n−2∑k=0∥Rk+1∥2≤13∥en∥2+3τ24∥Rn∥2+2τ3n−1∑k=1∥ek∥2+3τ2n−1∑k=0∥Rk+1∥2,

which gives

 ∥en∥2≤9τ24∥Rn∥2+2τn−1∑k=1∥ek∥2+9τ2n−1∑k=0∥Rk+1∥2≤2τn−1∑k=1∥ek∥2+c(τ2+h4)2.

then the desired result follows by Lemma 3.3.

###### Remark 3.3

By using similar techniques, we can show that the compact scheme (3)–(9) is stable for . In fact, suppose that is the solution of

 H(vk+1i−vki)=κ1τ1−α2(k+1∑l=0λ(α)lδ2xvk+1−li+k∑l=0λ(α)lδ2xvk−li) 0≤k≤N−1,  1≤i≤M−1, (12) vk0=φk1,vkM=φk2,1≤k≤N, v0i=ρi,0≤i≤M, (13)

where . Then, by subtracting (3)–(9) from (3.3)–(13), we obtain the following equations for ,

 +κ2τ1−β2(k+1∑l=0λ(β)lδ2xεk+1−li+k∑l=0λ(β)lδ2xεk−li)+κ2τ1−β2(k+1∑l=0λ(β)l+k∑l=0λ(β)l)δ2xρi, (14) 0≤k≤N−1,  1≤i≤M−1, εk0=εkM=0,1≤k≤N, ε0i=0,0≤i≤M.

Notice that [13, 14], and are less than 0 for and

 τ−αk∑l=0g(α)l=1Γ(1−α)+O(τ), τ−βk∑l=0g(β)l=1Γ(1−β)+O(τ).

We can therefore obtain

 τ−α2(k+1∑l=0λ(α)l+k∑l=0λ(α)l) =τ−α2[(1+α2)k+1∑l=0g(α)l+k∑l=0g(α)l−α2k−1∑l=0g(α)l] =1Γ(1−α)+O(τ),

and

 τ−β2(k+1∑l=0λ(β)l+k∑l=0λ(β)l)=1Γ(1−β)+O(τ).