High-Order Algorithms for Riesz Derivative and their Applications (III)

1. Introduction

In recent years, fractional calculus has attracted increasing interests due to its applications in physics, mechanics, etc. For more details, see the recent publications [1, 2, 3, 5, 6, 7, 8, 9, 10, 14, 15, 17, 19, 20, 21, 22], and references cited therein. The Riemann-Liouville (R-L) derivative and Caputo derivative are commonly used, respectively defined below:

and

where n − 1 < α < n ∈ Z⁺.

A linear combination of the left R-L derivative (defined above) and the right R-L derivative given as

defines a new derivative, i.e., the Riesz derivative,

The left R-L derivative reflects the dependence on the history, while the right R-L derivative the dependence upon the future. So the Riesz derivative value off(x) at x relies on the whole space (a, b), but the values are different at differentx ∈(a, b). From an angle of application, the case with α ∈ (0, 2) is mostly attracted attention.

From the studies available, the high-order numerical algorithms for RL derivatives were firstly established in [16]. The high-order algorithms for Caputo derivatives were firstly constructed in [12]. Shortly after, some other high-order algorithms for Caputo derivatives were also appeared, for example see [1, 3, 14] and references cited therein. And the high-order algorithms for Riesz derivatives were derived in [5, 6]. In [5], Ding et al. constructed the fourth-order schemes for the Riesz derivative and applied them to the space Riesz fractional diffusion equation. In that paper they established the fourth-order schemes for R-L derivatives [5] from which some similar fourth-order schemes were derived by other people. In [6], Ding et al. continued to establish the sixth-order, eighth-order, tenth-order and twelfth-order schemes for the Riesz derivative by using the Fourier analysis, where the fourth-order compact scheme for R-L derivative was especially highlighted. Then they used the two schemes, among them to the Riesz space fractional reaction-dispersion equation, where the rigorous error analysis was given. The other odd-order (third-order, fifth-order, seventh-order, ninth-order, eleventh-order) schemes can be established from [6] by choosing the different parameters therein. In spite of these, it is absolutely necessary to construct intuitionist and straightforward schemes for the Riesz derivatives. Here we find an interesting and enlightening way to establish high-order schemes (from 2nd-order to 6th-order) by using the corresponding generating functions [16]. Then we use these schemes to solve the Riesz space fractional turbulent diffusion equation.

In the following, we briefly introduce fractional modelling in this respect. From the known first Fick’s law

one obtain the following advection-diffusion equation,

If the advection-diffusion process at any position x ∈(a, b) relies on the whole space (a, b) (i.e., long-range interactions), then the classical Fick’s law does not work well yet. However, the fractional derivative can well characterize such long-range interactions. Now we generalize the typical Fick’s law to a fractional version,

where d₁, d₂, d_α are positive constants. It follows that the Riesz space fractional turbulent diffusion equation (or, Riesz type turbulent diffusion equation) is obtained,

If there is a source term, then one has

where the Riesz partial derivative with order α ∈ (0, 1) is given by

Equation (1.1) is subject to the following initial value condition,

and the boundary value conditions (here choosing the homogeneous condition for brevity),

The remainder of this paper is outlined as follows. In Section 2, five kinds of high-order (2nd-order, 3rd-order, ···, 6th-order) algorithms for the Riesz derivatives are developed. In the next Section 3, we apply the derived schemes to the Riesz-type turbulent diffusion equation (1.1). Here we only use the 2nd-order, 4th-order, 6th-order schemes to (1.1). The convergence orders are 𝒪(τ² + h²), 𝒪(τ² + h⁴), and 𝒪(τ² + h⁶), where τ and h are temporal and spatial stepsizes, respectively. In Section 4, numerical examples are presented which support the theoretical analysis. The last Section5 concludes this article. In Appendices A and B we provide the detailed proofs of the main results from Section 2.

2. High-order numerical schemes

If f^(k)(a+) = 0 (k = 0, 1,...,p − 1), then it follows from [16] that the left R-L derivative has the following approximations

in which h is the steplength. Here we only show interests in p = 2, 3, 4, 5, 6.

The convolution (or weight) coefficients

where

By tedious but direct calculations, one has

and

Here

On the other hand, if f_(k)(b−) = 0(k = 0, 1,...,p−1), then one has the approximations below,

where h is the stepsize.

Based on (2.1) and (2.2), if f(x), together with its derivatives, has homogeneous boundary value conditions, one easily gets

Here, we limit our interests in a ∈ (0, 1). The case a ∈ (1, 2) can be similarly studied. When

The properties of convolution coefficients

For the space fractional differential equations, we use the following properties of coefficients

Proof.

We only prove p = 2, the rest cases can be almost similarly shown. Let

f1(α, θ)=∑ℓ=0∞ϖ2,ℓ(α) cos (ℓθ)

which can be expanded as

f1(α, θ)=∑ℓ=0∞ϖ2,ℓ(α) cos (ℓθ)=12∑ℓ=0∞ϖ2,ℓ(α)( exp (iℓθ)+ exp (−iℓθ))                 = 12[(1− exp (iθ))α(32−12 exp (iθ))α+(1− exp (−iθ))α(32 −12 exp (−iθ))α]

Note that f₁(α, θ) is a real-value and even function, so we need only consider θ ∈ [0, π].

Using the following equations

(1− exp (±iθ))α=(2 sin θ2)α exp (±iα(θ−π2))

and

(x−yi)α=(x2+y2)α2 exp (iαϕ), ϕ=− arctan yx,

we can rewrite f₁(α, θ) as

f1(α, θ)=(2 sin θ2)α(λ12(θ)+μ12(θ))α2 cos α(θ−π2+ϕ1)

where

λ1(θ)=3− cos θ, μ1(θ)= sin θ, ϕ1=− arctan μ1(θ)λ1(θ)

Let

z(θ)=θ−π2+ϕ1,0≤θ≤π

Then

z′(θ)=(θ−π2+ϕ1)'=3sin2(θ2)1+3sin2(θ2)≥0

Hence z(θ) is an increasing function in [0 ,π] and

z min (θ)=z(0)=−π2, z max (θ)=z(π)=0.

Evidently α ∈ (0, 1) and θ ∈ [0 ,π] imply

. Furthermore, one has

f1(α, θ)=(2 sin θ2)α(λ12(θ)+μ12(θ))α2 cos α(θ−π2+ϕ1)≥0.

All this ends the proof. □

For p = 4, α can not attain to 1. But we have following theorem.

Theorem 2.2.

If

then the following inequality holds:

∑ℓ=0∞ϖ4,ℓ(α) cos (ℓθ)≥0, θ∈[−π, π].

Proof.

Let

. By almost the same reasoning as those in Theorem 2.1 , we can get

 f2(α, θ)=(2 sin θ2)α(λ22(θ)+μ22(θ))α2 cos α(θ−π2+ϕ2)

where

λ2(θ)=25−23 cos θ+13 cos 2θ−3 cos 3θμ2(θ)=23 sin θ−13 sin 2θ+3 sin 3θ, ϕ2=− arctan μ2(θ)λ2(θ)

Since

λ2(θ)=14( cos (θ)−12)2+24cos2(θ) sin 2(θ2)+172>0

and

μ2(θ)=[12( cos (θ)−1312)2+7124] sin (θ)≥0

it immediately follows that

. We need only consider 0 ≤ θ ≤ π, therefore

.

Obviously, if

, then f₂(α, θ) ≥ 0. A sufficient condition for

is

−π2≤ min α(θ−π2+ϕ2)≤0

i.e.,

0<α≤ min {ππ−θ−2ϕ2}.

Let y(θ) = π − θ − 2 φ₂, then

y′(θ)=1920(5 cos θ−1)sin4(θ2)a22(θ)+b22(θ)

It is clear that

is a unique maximum point of y(θ) when θ ∈ [0, π], i.e.,

y max (θ)=y max ( arccos 15)=π− arccos 15+2 arctan 1916317

it follows that

 min {ππ−θ−2φ2}=ππ− arccos 15+2 arctan 1916317≈0.8439,

i.e.,

0<α≤0.8439

This finishes the proof. □

Theorems2.1 and2.2 are very suitable for numerically analyzing R-L space fractional partial differential equations. But for numerically analyzing R-L time fractional partial differential equations, the monotonicity of the coefficients

Theorem 2.3.

The second-order coefficients

satisfy:

1. ϖ2,0(α)=(32)α>0 ,ϖ2,1(α)=−4α3(32)α<0ϖ2,2(α)=α(8α−3)9(32)α ϖ2,3(α)=−4α(α−1)(8α−7)81(32)αϖ2,4(α)=α(α−1)(64α2−176α+123)486(32)α

2. When 0 < α < 1,

   

   and

   

   hold for ≥ 4,

3. When 1 < α < 2,

   

   and

   

   hold for ≥ 5,

P r o o f. See Appendix A. □

From the results of the above theorem, one can see that: If α = 1, then

Besides their monotonicity, studying bounds of these coefficients is also of importance, which can be used to analyze the stability and convergence for time fractional differential equations. In [4], Dimitrov gave the bounds for first-order

Theorem 2.4.

The first-order coefficients

(0 < α < 1) satisfy:

1. 

   , ℓ ≥ 3

   where

   

   

2. 

   where

   

In this paper, we can give tighter estimates for the lower bounds. See the following theorem.

Theorem 2.5.

The first-order coefficients

(0 < α < 1) satisfy:

1. 

   , ℓ ≥ 3

   where

   

2. 

   ℓ ≥ 3.

   where

   

Next, we list the comparison theorem for

In the following, we give the bounds for the corresponding coefficients.

Theorem 2.7.

The first-order coefficients

, (0 < α < 1) satisfy:

1. 

   , ℓ ≥ 4, where

   B¯(1+α,ℓ)=(1−α)α(1+α)6(3ℓ)2(2+α)B¯(1+α,ℓ)=α(1+α)2(3ℓ+1)2+α

2. 

   where

   S¯(1+α,ℓ)=(1−α)α(1+α)2(3+2α)(3ℓ)3+2α,S¯(1+α,ℓ)=3α2(3ℓ)1+α

Now we start to show the bounds of the second-order coefficients.

Theorem 2.8.

The second-order coefficients

, and

for 0 < α < 1 satisfy:

1. 

   where

   B2L(α,ℓ)=(32)α[(1+(13)ℓ)α(1−α)2(2ℓ)2(1+α)−(1−(13)ℓ−1)α222α+11+(α+1)ℓ],

   B2L(α,ℓ)=(32)α[(1+(13)ℓ)α2α+1(ℓ+1)α+1−α2(1−α)242α+12(1−(13)ℓ−1)(2ℓ)4(α+1)]

2. 

   where

   B¯(1+α,ℓ)=(32)1+α[(1+(13)ℓ)(1−α)α(1+α)6(3ℓ)2(2+α)+(1−α)2α2(1+α)2216(1−(13)ℓ−3)(6ℓ)4(2+α)−α(1+α)22(13+(13)ℓ−1)(3ℓ)2+α],B¯(1+α,ℓ)=(32)1+α[(1+(13)ℓ)α(α+1)3α+22(ℓ+1)α+2+(1−(13)ℓ−3)α2(1+α)232(2+α)24(1+(2+α)ℓ)−(1−α)α(1+α)26(13+(13)ℓ−1)(3ℓ−1)2(2+α)]

All proofs for Theorems 2.4-2.8 are given in Appendix B.

3. Numerical methods for the Riesz-type turbulent diffusion equation

Define t_k = kτ, k = 0, 1 ,...,N, and

3.1. The 2nd-order scheme in space

Firstly, using the Crank-Nicolson method for the Riesz space fractional turbulent diffusion equation (1.1) in time direction, we obtain

(3.1)u(xj,tk+1)−u(xj,tk)τ=12(d2∂2u(xj,tk+1)∂x2−d1∂u(xj,tk+1)∂x

(3.2)+dα∂αu(xj,tk+1)∂|x|α+d2∂2u(xj,tk)∂x2−d1∂u(xj,tk)∂x+dα∂αu(xj,tk)∂|x|α)+s(xj,tk+12)+O(τ2)

Secondly, for the first- and second-order derivatives, we use the following approximations, respectively

(3.3)∂u(xj,tk)∂x=μxδxu(xj,tk)+O(h2)

and

(3.4)∂2u(xj,tk)∂x2=δx2u(xj,tk)+O(h2)

where μ_x and

are defined by

μxδxu(xj,tk)=u(xj+1,tk)−u(xj−1,tk)2h

and

δx2u(xj,tk)=u(xj+1,tk)−2u(xj,tk)+u(xj−1,tk)h2.

Next, we choose the second-order formula to approximate Riesz derivative,

(3.5)∂αu(x,t)∂|x|α=−12cos(πα2)hα(∑ℓ=0∞ϖ2,ℓ(α)u(x−ℓh,t)+∑ℓ=0∞ϖ2,ℓ(α)u(x+ℓh,t))+O(h2)

Substituting (3.3), (3.4) and (3.5) into (3.2) and removing the high-order term yield

(3.6)(d2h2+d12h)uj−1k+1−(2τ+2d2h2)ujk+1+(d2h2−d12h)uj+1k+1−ν[∑ℓ=0jϖ2,ℓ(α)uj−lk+1+∑ℓ=0M−jϖ2,ℓ(α)uj+lk+1]=−(d2h2+d12h)uj−1k2τ−2d2h2)ujk−(d2h2−d12h)uj+1k+l↗[∑ℓ=0jϖ2,ℓ(α)uj−lk+∑ℓ=0M−jϖ2,ℓ(α)uj+lk]−2sjk+12,j=1,…,M−1,k=0,1,…,N−1,

where

Next we discuss the stability and convergence of scheme (3.6).

Theorem 3.1.

The numerical scheme (3.6) is unconditionally stable and convergent with order 𝒪(τ² + h²).

Proof.

Assume that the solution of equation (1.1) can be zero extended to the whole real line R. Suppose

is the approximation solution of equation (3.6) and let

, then the error equation reads as

(3.7)(d2h2+d12h)ℰj−1k+1−(2τ+2d2h2)ℰjk+1+(d2h2−d12h)ℰj+1k+1−v[∑ℓ=0∞ϖ2,ℓ(γ)ℰj−lk+1+∑ℓ=0∞ϖ2,ℓ(γ)ℰj+lk+1]=−(d2h2+d12h)ℰj−1k−(2τ−2d2h2)ℰjk−(d2h2−d12h)ℰj+1k+U[∑ℓ=0∞ϖ2,ℓ(γ)ℰj−lk+∑ℓ=0∞ϖ2,ℓ(γ)ℰj+lk],j=1,…,M−1,k=0,1,…,N−1.

Let

be the solution of equation (3.7),

, where θ ∈ [−π, π] is called the phase angle. The stability condition of scheme (3.6) is |ξ (θ)| ≤ 1 for all θ ∈ [−π, π].

Substituting

into (3.7) gives

ξ(θ)= (2hτ−4d2hsin2(θ2)−2νh∑ℓ=0∞ϖ2,ℓ(γ)cos(ℓθ))−id1sin(θ)(2hτ+4d2hsin2(θ2)+2νh∑ℓ=0∞ϖ2,ℓ(γ)cos(ℓθ))+id1sin(θ).

According to Theorem 2.1, we easily obtain |ξ (θ)| ≤ 1. So scheme (3.6) is unconditionally stable. It is easy to show the convergence order of scheme (3.6) is 𝒪 (τ² +h²). In effect, the proof is almost the same as that of [5]. □