High-order algorithms for Riesz derivative and their applications (II)

Abstract

In this paper, we firstly develop two high-order approximate formulas for the Riesz fractional derivative. Secondly, we propose a temporal second order numerical method for a fractional reaction-dispersion equation, where we discretize the Riesz fractional derivative by using two numerical schemes. We prove that the numerical methods for a spatial Riesz fractional reaction dispersion equation are both unconditionally stable and convergent, and the orders of convergence are $\mathcal{O}({\tau }^2+h^6)$ and $\mathcal{O}({\tau }^2+h^8)$, in which τ and h are spatial and temporal step sizes, respectively. Finally, we test our numerical schemes and observe that the numerical results are in good agreement with the theoretical analysis.

Introduction

Recently, fractional differential equations have attracted increasing interests mainly because they have been proved to be successful in modeling many natural phenomena [1], [19], [20], [27], [30], [31]. However, most of fractional differential equations cannot be easily analytically solved. So, more and more researchers are interested in studying the numerical methods for the fractional differential equations [8], [11], [15], [24], [25], [38]. Meerschaert and Tadjeran [23] considered the finite difference method for the space fractional advection–dispersion equation with Riemann–Liouville derivative, where a shifted Grünwald–Letnikov scheme was proposed. In [21], [22], Murillo and Yuste constructed the difference schemes for the fractional diffusion equation by using the so-called L1 formula and for the fractional diffusion-wave equation by using the L2 formula, respectively. In [17], Liu et al. used L1 scheme to discrete the time Caputo derivative and the shifted Grünwald–Letnikov formulae to discrete the left- and right-sided Riemann–Liouville derivatives. Yang et al. [35] obtained numerical solutions for the space Riesz fractional diffusion and advection–dispersion equations on a finite domain by using the L1, L2 methods. Shen et al. [32] used a second-order numerical scheme for the space Riesz fractional advection–dispersion equation based on fractional central difference operator. In addition, Podlubny firstly proposed the matrix approach to discrete fractional calculus and developed the idea to approximate fractional partial differential equations in [28]. Li et al. [18] firstly developed a finite element method for solving the nonlinear time–space fractional subdiffusion and superdiffusion equations, respectively. They analyzed the stability and error estimates for the semi-discrete and fully discrete schemes. Jiang and Ma [13] constructed a high-order finite element method for solving the time Caputo fractional partial differential equations. A space fractional reaction–diffusion equation has been considered by Burrage et al. [2]. In [6], Deng and Hesthaven proposed discontinuous Galerkin methods for fractional diffusion equations. Wang et al. [16] proposed a difference method for the fractional diffusion equation. The Petrov–Galerkin spectral method for fractional initial-value problems and fractional final-value problems was developed by Zayernouri and Karniadakis [36]. Later on, they developed an exponentially accurate fractional spectral collocation method for solving steady-state and time-dependent fractional partial differential equations [37]. It is worth noting that these difference schemes have high-accuracy in space and the order of accuracy in time is second-order. Very recently, Ding et al. [7] constructed two different fourth-order numerical schemes for the Riesz derivative and they applied one of them to the spatial Riesz fractional diffusion equation.

As far as we know, there are no numerical formulas more than sixth-order of the accuracy for the Riesz fractional derivative. The aim of this paper is to construct several high-order approximate schemes (the orders more than 6) for the Riesz fractional derivative and use sixth- and eighth-order schemes to solve the following space Riesz fractional reaction-dispersion equation

$$\frac{\partial u(x,t)}{\partial t}=-u(x,t)+K_{\alpha }\frac{{\partial }^{\alpha }u(x,t)}{\partial {|x|}^{\alpha }}+f(x,t),1<\alpha <2,a<x<b,0<t\le T.$$

Suppose that the above equation (1) is subject to the initial and boundary value conditions as below

$$u(x,0)=\psi (x),a\le x\le b,$$

$${}_{\mathit{RL}}D_{a,x}^{\alpha -2}u(x,t){|}_{x=a}=0,0\le t\le T,$$

$${}_{\mathit{RL}}D_{x,b}^{\alpha -2}u(x,t){|}_{x=b}=0,0\le t\le T,$$

where $K_{\alpha }>0$ is the dispersion coefficient. $\frac{{\partial }^{\alpha }u(x,t)}{\partial {|x|}^{\alpha }}$ denotes the Riesz fractional derivative with order α $(1<\alpha <2)$ defined by [14]

$$\frac{{\partial }^{\alpha }u(x,t)}{\partial {|x|}^{\alpha }}=-{\Psi }_{\alpha }({}_{\mathit{RL}}D_{a,x}^{\alpha }+{}_{\mathit{RL}}D_{x,b}^{\alpha })u(x,t),1<\alpha <2,$$

where ${\Psi }_{\alpha }=\frac{1}{2}\sec (\frac{\pi \alpha }{2})$, ${}_{\mathit{RL}}D_{a,x}^{\alpha }$ and ${}_{\mathit{RL}}D_{x,b}^{\alpha }$ are the left- and right-sided Riemann–Liouville derivatives represented in the following forms

$${}_{\mathit{RL}}D_{a,x}^{\alpha }u(x,t)=\frac{1}{\Gamma (2-\alpha )}\frac{{\partial }^2}{\partial x^2}\underset{a}{\overset{x}{\int }}\frac{u(\eta ,t)}{{(x-\eta )}^{\alpha -1}}d\eta ,x\in [a,b],$$

and

$${}_{\mathit{RL}}D_{x,b}^{\alpha }u(x,t)=\frac{1}{\Gamma (2-\alpha )}\frac{{\partial }^2}{\partial x^2}\underset{x}{\overset{b}{\int }}\frac{u(\eta ,t)}{{(\eta -x)}^{\alpha -1}}d\eta ,x\in [a,b].$$

Here $\Gamma (\cdot )$ indicates the Euler's Gamma function. In particular, when $\alpha =2$, the Riesz derivative $\frac{{\partial }^{\alpha }u(x,t)}{\partial {|x|}^{\alpha }}$ reduces to the general second-order partial derivative $\frac{{\partial }^{2}u(x,t)}{\partial x^2}$.

In addition, operators ${}_{\mathit{RL}}D_{a,x}^{\alpha -2}$ and ${}_{\mathit{RL}}D_{x,b}^{\alpha -2}$ denote the left- and right-sided Riemann–Liouville integral operators and are defined as follows $(1<\alpha <2)$:

$${}_{\mathit{RL}}D_{a,x}^{\alpha -2}u(x,t)=\frac{1}{\Gamma (2-\alpha )}\underset{a}{\overset{x}{\int }}{(x-\eta )}^{1-\alpha }u(\eta ,t)d\eta ,$$

and

$${}_{\mathit{RL}}D_{x,b}^{\alpha -2}u(x,t)=\frac{1}{\Gamma (2-\alpha )}\underset{x}{\overset{b}{\int }}{(\eta -x)}^{1-\alpha }u(\eta ,t)d\eta .$$

The Dirichlet boundary value conditions (or the first kind of boundary value conditions) (1.2) and (1.3) for Eq. (1) mean that,

$${}_{\mathit{RL}}D_{a,x}^{\alpha -2}u(x,t){|}_{x=a}=\underset{x\to a+}{\lim }(\frac{1}{\Gamma (2-\alpha )}\underset{a}{\overset{x}{\int }}{(x-\xi )}^{1-\alpha }u(\xi ,t)d\xi ),$$

and

$${}_{\mathit{RL}}D_{x,b}^{\alpha -2}u(x,t){|}_{x=b}=\underset{x\to b-}{\lim }(\frac{1}{\Gamma (2-\alpha )}\underset{x}{\overset{b}{\int }}{(\xi -x)}^{1-\alpha }u(\xi ,t)d\xi ).$$

If the integral operators and the limits, ${}_{\mathit{RL}}D_{a,x}^{2-\alpha }$ and ${\lim }_{x\to a+}$, ${}_{\mathit{RL}}D_{x,b}^{2-\alpha }$ and ${\lim }_{x\to b-}$, can be exchanged, then it follows from Eqs. (1.2) and (1.3) that

$$0={}_{\mathit{RL}}D_{a,x}^{2-\alpha }\{\underset{x\to a+}{\lim }(\frac{1}{\Gamma (2-\alpha )}\underset{a}{\overset{x}{\int }}{(x-\xi )}^{1-\alpha }u(\xi ,t)d\xi )\}=\underset{x\to a+}{\lim }{}_{\mathit{RL}}D_{a,x}^{2-\alpha }({}_{\mathit{RL}}D_{a,x}^{\alpha -2}u(x,t))=\underset{x\to a+}{\lim }u(x,t)=u(a,t),$$

and

$$0={}_{\mathit{RL}}D_{x,b}^{2-\alpha }\{\underset{x\to b-}{\lim }(\frac{1}{\Gamma (2-\alpha )}\underset{x}{\overset{b}{\int }}{(\xi -x)}^{1-\alpha }u(\xi ,t)d\xi )\}=\underset{x\to b-}{\lim }{}_{\mathit{RL}}D_{x,b}^{2-\alpha }({}_{\mathit{RL}}D_{x,b}^{\alpha -2}u(x,t))=\underset{x\to b-}{\lim }u(x,t)=u(b,t).$$

So the boundary value conditions (1.2) and (1.3) can be converted into

$$u(a,t)=0,0\le t\le T,$$

$$u(b,t)=0,0\le t\le T,$$

respectively. This is why the boundary value conditions like (1.4) and (1.5) were chosen in all the past studies, where the assumption of the exchangeability of the integral operators and the limits was used. Throughout this paper, such exchangeability of the integral operators and the limits is always assumed.

Generally speaking, definite conditions of fractional differential equations with Riemann–Liouville or Riesz derivatives are quite different from those of typical differential equations. So to propose definite conditions for Riemann–Liouville type differential equations must be very careful. For more details, see Appendix A.

We remark that there exists another equivalent definition for the Riesz derivative by utilizing the Fourier transform on an infinite domain under homogeneous Dirichlet boundary conditions [35]

$$\frac{{\partial }^{\alpha }u(x,t)}{\partial {|x|}^{\alpha }}=-{(-\mathrm{\Delta })}^{\frac{\alpha }{2}}u(x,t),$$

where the fractional power of the Laplace operator is defined as follows:

$$-{(-\mathrm{\Delta })}^{\frac{\alpha }{2}}u(x,t)=-{\mathcal{F}}_{\omega }^{-1}{|\omega |}^{\alpha }{\mathcal{F}}_{x}u(x,t),$$

in which ${\mathcal{F}}_x$ and ${\mathcal{F}}_{\omega }^{-1}$ denote the Fourier and inverse Fourier transforms, respectively.

The plan of this paper is listed as follows. In Section 2, we construct several high-order numerical formulas for the Riesz fractional derivative. In Section 3, two new finite difference schemes are proposed for solving a Riesz fractional reaction-dispersion equation. In Section 4, we discuss the stability and convergence of the proposed numerical methods in details. Some numerical examples are taken to show the effectiveness and accuracy of the numerical schemes in Section 5. At last, some remarks are concluded in Section 6.