High-order numerical methods for the Riesz space fractional advection-dispersion equations

by   Libo Feng, et al.
Xiamen University

In this paper, we propose high-order numerical methods for the Riesz space fractional advection-dispersion equations (RSFADE) on a finite domain. The RSFADE is obtained from the standard advection-dispersion equation by replacing the first-order and second-order space derivative with the Riesz fractional derivatives of order α∈(0,1) and β∈(1,2], respectively. Firstly, we utilize the weighted and shifted Grünwald difference operators to approximate the Riesz fractional derivative and present the finite difference method for the RSFADE. Specifically, we discuss the Crank-Nicolson scheme and solve it in matrix form. Secondly, we prove that the scheme is unconditionally stable and convergent with the accuracy of O(τ^2+h^2). Thirdly, we use the Richardson extrapolation method (REM) to improve the convergence order which can be O(τ^4+h^4). Finally, some numerical examples are given to show the effectiveness of the numerical method, and the results are excellent with the theoretical analysis.



There are no comments yet.


page 1

page 2

page 3

page 4


A high-order and fast scheme with variable time steps for the time-fractional Black-Scholes equation

In this paper, a high-order and fast numerical method is investigated fo...

On high-order schemes for tempered fractional partial differential equations

In this paper, we propose third-order semi-discretized schemes in space ...

On Theoretical and Numerical Aspect of Fractional Differential Equations with Purely Integral Conditions

In this paper, we are interested in the study of a problem with fraction...

Numerical solution of the general high-dimensional multi-term time-space-fractional diffusion equations

In this article, an advanced differential quadrature (DQ) approach is pr...

Finite Difference Weerakoon-Fernando Method to solve nonlinear equations without using derivatives

This research was mainly conducted to explore the possibility of formula...

Finite difference and numerical differentiation: General formulae from deferred corrections

This paper provides a new approach to derive various arbitrary high orde...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

There have been increasing interests in the description of the physical and chemical processes by means of equations involving fractional derivatives over the last decades. And, fractional derivatives have been successfully applied into many sciences, such as physics b1 , biology b2 , chemistry b3 , hydrology b4 ; b5 ; b6 ; b7 , and even finance b8 . In groundwater hydrology the fractional advection-dispersion equation (FADE) is utilized to model the transport of passive tracers carried by fluid flow in a porous medium b6 ; b9 .

Considerable numerical methods for solving the FADE have been proposed. Kilbas et al. Kilbas06 introduced the theory and applications of fractional differential equations. Meerschaert and Tadjeran b10 developed practical numerical methods to solve the one-dimensional space FADE with variable coefficients on a finite domain. Liu et al. b6

transformed the space fractional Fokker-Planck equation into a system of ordinary differential equations (method of lines), which was then solved using backward differentiation formulas. Momani and Odibat

b9 developed two reliable algorithms, the Adomian decomposition method and variational iteration method, to construct numerical solutions of the space-time FADE in the form of a rapidly convergent series with easily computable components. Zhuang et al. b11 discussed a variable-order fractional advection-diffusion equation with a nonlinear source term on a finite domain. Liu et al. b12 proposed an approximation of the Lévy-Feller advection-dispersion process by employing a random walk and finite difference methods. In addition, other finite difference methods b13 , finite element method b23 , finite volume method b24 , homotopy perturbation method b25 and spectral method b26 ; b27 are also employed to approximate the FADE.

In this paper, we consider the following RFADE:


subject to the initial condition:


and the zero Dirichlet boundary conditions:


where , , and represent the average fluid velocity and the dispersion coefficient. The Riesz space fractional operators and on a finite domain are defined respectively as

where , , and

where represents the Euler gamma function.

The fractional kinetic equation (1) possesses a physical meaning (see b21 ; b22 for further details). Physical considerations of a fractional advection-dispersion transport model restrict , , and we assume and so that the flow is from left to right. In the case of and , Eq.(1) reduces to the classical advection-dispersion equation (ADE). In this paper, we only consider the fractional cases: when , Eq.(1) reduces to the Riesz fractional diffusion equation (RFDE) b21 and when , the Riesz fractional advection-dispersion equation (RFADE) is obtained b22 .

For the RFADE (1), Anh and Leonenko b18 presented a spectral representation of the mean-square solution without the non-homogeneous part and for some range of values of the parameters. Later, Shen et al. b19 derived the fundamental solution of Eq.(1) and discussed the numerical approximation of Eq.(1) using finite difference method with first convergence order. Another method based on the spectral approach and the weak solution formulation was given in Leonenko and Phillips b20 . In addition, Zhang et al. b14 use the Galerkin finite element method to approximate the RFADE. Besides, Yang et al. b28 applied the -approximation method, the standard/shifted Grünwald method, and the matrix transform method (MTM) to solve the RFADE. And, Ding et al. b29 also consider the numerical solution of the RFADE by using improved matrix transform method and the (2, 2) Pade approximation. Most of the numerical methods proposed by these authors are low order or lack stability analysis.

In this paper, based on the weighted and shifted Grünwald difference (WSGD) operators to approximate the Riesz space fractional derivative, we obtain the second order approximation of the RFADE. Furthermore, We propose the finite difference method for the RFADE and obtain the Crank-Nicolson scheme. Moreover, we prove that the Crank-Nicolson scheme is unconditionally stable and convergent with the accuracy of and improve the convergence order to by applying the Richardson extrapolation method.

The outline of the paper is as follows. In Section 2, the WSGD operators and some lemmas are given. In Section 3, we first present the finite difference method for the RFADE, and then derive the Crank-Nicolson scheme. We proceed with the proof of the stability and convergence of the Crank-Nicolson scheme in Section 4. Besides, we further improve the convergence order by applying the Richardson extrapolation method. In order to verify the effectiveness of our theoretical analysis, some numerical examples are carried out and the results are compared with the exact solution in Section 5. Finally, the conclusions are drawn.

2 The approximation for the Riemann-Liouville fractional derivative

First, in the interval , we take the mesh points , , and , , where , , i.e., and are the uniform spatial step size and temporal step size. Now, we give the definition of the Riemann-Liouville fractional derivative.

Definition 1 (b15 ).

The order left and right Riemann-Liouville fractional derivatives of the function on , are given by

  • left Riemann-Liouville fractional derivative:

  • right Riemann-Liouville fractional derivative:

Generally, the standard Grünwald-Letnikov difference formula is applied to approximate the Riemann-Liouville fractional derivative. Meerschaert and Tadjeran b10 showed that the standard Grünwald-Letnikov difference formula was often unstable for time dependent problems and they proposed the shifted Grünwald difference operators

whose accuracies are first order, i.e.,

where are integers and . In fact, the coefficients are the coefficients of the power series of the function ,

for all , and they can be evaluated recursively

Lemma 1 (b13 ).

Suppose that , then the coefficients satisfy

Lemma 2 (b13 ; b16 ).

Suppose that , then the coefficients satisfy

Inspired by the shifted Grünwald difference operators and multi-step method, Tian et al. b16 derive the WSGD operators:

Lemma 3 (b16 ).

Supposing that , let , and

and their Fourier transforms belong to

, then the WSGD operators satisfy

uniformly for , where , are integers and .

Tian et al. prove Lemma 3 under the additional conditions that . In fact, their proof also holds for the case. The proof proceeds the same as b16 , hence we will not repeat it here.

Remark 1.

Considering a well defined function on the bounded interval , if or , the function can be zero extended for or . And then the order left and right Riemann-Liouville fractional derivatives of at each point can be approximated by the WSGD operators with second order accuracy

where and .

When , and , the discrete approximations for the Riemann-Liouville fractional derivatives on the domain are




Now, we discuss the properties of the coefficients and .

Lemma 4.

Suppose that , then the coefficients satisfy


Combining the definition of and the property of , it is easy to derive the value of , and . When , by the definition of

we have as for and . Moreover,

Recalling Lemma 1, when , we obtain when , i.e., . For the sum , we have

Since when , for . ∎

Lemma 5 (b16 ).

Suppose that , then the coefficients satisfy

3 The finite difference method for the RFADE

In this section, we utilize the Eqs.(4)-(7) to approximate the Riesz space fractional derivative and derive the Crank-Nicolson scheme of the equation. We define , , let be a finite domain, setting be a uniform partition of , which is given by for , where and are the time and space steps, respectively. First, we present the semi-discrete form of Eq.(1),


Let be the approximation solution of . Substituting (4)-(7) into (10), we obtain


where and . Denote




Thus, Eq. (3) can be simplified as


The boundary and initial conditions are discretized as

where .

4 Theoretical analysis of the finite difference method

4.1 Stability

Before giving the proof, we start with some useful lemmas.

Lemma 6 (b17 ).

Let be an order positive define matrix, then for any parameter , the following two inequalities


Now, we discuss the property of matrix .

Theorem 1.

Suppose that and , , and are defined as (12)-(14), then the coefficients satisfy

i.e., is strictly diagonally dominant.


It is easy to obtain

where and . First, we consider the signs of . According to Lemma 4, when , , thus . According to Lemma 5, when , , thus . Therefore, when or . For the items and , we have

Since and , then

According to Lemmas 4 and Lemma 5, we have and , thus

as and . Now, for a given , we consider the sum


Thus, the proof is completed. ∎

According to the theorem, it is easy to conclude the following corollaries.

Corollary 1.

The matrix is strictly diagonally dominant as well. Therefore, is invertible and Eq. (15) is solvable.

Corollary 2.

The matrix is symmetric positive definite.


In view of (14), the symmetry of is evident. Let