Second order WSGD operators II: A new family of difference schemes for space fractional advection diffusion equation
The second order weighted and shifted Grünwald difference (WSGD) operators are developed in [Tian et al., arXiv:1201.5949] to solve space fractional partial differential equations. Along this direction, we further design a new family of second order WSGD operators; by properly choosing the weighted parameters, they can be effectively used to discretize space (Riemann-Liouville) fractional derivatives. Based on the new second order WSGD operators, we derive a family of difference schemes for the space fractional advection diffusion equation. By von Neumann stability analysis, it is proved that the obtained schemes are unconditionally stable. Finally, extensive numerical experiments are performed to demonstrate the performance of the schemes and confirm the convergent orders.
Mathematics Subject Classification(2010): 26A33, 65M06, 65M12
Key words: Riemann-Liouville fractional derivative, WSGD operator, Fractional advection diffusion equation, Finite difference approximation, Stability
As an extension of classical calculus, fractional calculus also has more than three centuries of history. However, it seems that the applications of fractional calculus to physics and engineering attract enough attention only in the last few decades. Nowadays, it has been found that a broad range of non-classical phenomena in the applied sciences and engineering can be well described by some fractional kinetic equations . The fractional calculus is becoming more and more popular, especially in describing the anomalous diffusions [1, 10, 18, 8], which arise in physics, chemistry, biology, and other complex dynamics. With the appearance of various kinds of fractional partial differential equations (PDEs), finding the ways to effectively solve them becomes a natural topic. Based on the integral transformations, sometimes we can find the analytical solutions of linear fractional PDEs with constant coefficients [9, 10, 18]. Even in these cases, most of the time their solutions are expressed by transcendental functions or infinite series. So looking for the numerical solutions of fractional PDEs becomes a realistic expectation in practical applications; and designing efficient and robust numerical schemes for fractional PDEs is the basis to this task.
Over the past decade, the finite difference methods are implemented in simulating the space fractional advection diffusion equations [13, 5, 15, 14, 16, 23]. However, most of the time, the accuracy of finite difference scheme or finite volume method is first order. In recent years, high order discretizations and fast methods for space fractional PDEs with Riemann-Liouville fractional derivatives attract many authors’ interests. Based on the Toeplitz-like structure of the difference matrix, Wang et al  solve the algebraic equation corresponding to fractional diffusion equation with a cost; Pang and Sun  propose a multigrid method to solve the discretized system of the fractional diffusion equation. By using the linear spline approximation, Sousa and Li provide a second order discretization for the Riemann-Liouville fractional derivatives and establish an unconditionally stable finite difference method for one-dimensional fractional diffusion equation in ; the results on two-dimensional two-sided space fractional convection diffusion equation in finite domain can be seen in . Ortigueira  gives the “fractional centred derivative” to approximate the Riesz fractional derivative with second order accuracy, and this scheme is used by Çelik and Duman in  to approximate fractional diffusion equation with the Riesz fractional derivative in a finite domain.
More recently, Tian et al  propose the second order difference approximations, called weighted and shifted Grünwald difference (WSGD) approximations, to the Riemann-Liouville space fractional derivatives. The basic idea of the WSGD approximation is to cancel the low order terms by combining the Grünwald difference operators with different shifts and weights. Along this direction, this paper further introduces a new family of WSGD opertors by providing more to be chosen parameters, called second order WSGD operators II. As specific applications, the second order WSGD operators II are used to discretize the space fractional derivative of the space fractional advection diffusion equation. And the von Neumann stability analyses are performed to the obtained numerical schemes. We theoretically prove and numerically verify that with proper chosen parameters for the WSGD space discretizations, the obtained implicit and Crank-Nicolson schemes are both unconditionally von Neumann stable.
The outline of this paper is organized as follows. In Section 2, we first introduce a new family of second order discretizations of the Riemann-Liouville fractional derivatives. Then two kinds of difference schemes for the one dimensional space fractional advection diffusion equation are presented in Section 3, and the corresponding numerical stabilities are performed. Furthermore, in Section 4, the numerical schemes and the proof of numerical stability are extended to two dimensional cases. Finally, in Section 5, extensive numerical experiments are carried out to demonstrate the performance of the proposed numerical schemes and confirm theoretical analysis.
2 Second order WSGD discretizations for the Riemann-Liouville fractional derivatives
There are some different definitions for the fractional derivatives, such as the Grünwald-Letnikov derivative, the Riemann-Liouville derivative, the Caputo derivative, and other modified definitions for the practical applications . Most of the time, they are not completely equivalent; in particular, the Grünwald-Letnikov and Riemann-Liouville derivatives are equivalent if ignoring the regularity requirements of the functions being performed. Usually, the more popular space fractional derivative is the Riemann-Liouville derivative. Let the function be defined on the finite interval , then the left and right Riemann-Liouville fractional derivatives of order are, respectively, defined by
Based on the definition of the Grünwald-Letnikov derivative and its relationship with the left and right Riemann-Liouville derivatives, we have 
where , is a positive integer, is the gamma function, and is the normalized Grünwald weights. However, the difference scheme based on (3) and/or (4) to discretize space fractional derivatives for the time dependent problems is unconditionally unstable . To overcome this problem, Meerschaert and Tadjeran firstly introduce the shifted Grünwald-Letnikov approximation formulas 
If a function has support on the finite interval , then it can be extended to the whole real line . In this case the symbol (or ) in left (right) Riemann-Liouville fractional derivative (1) ((2)) can be replaced by (or ). And the following properties hold.
 Let , and its Fourier transform belong to . Let and . Define
uniformly in as .
Inspired by the shifted Grünwald difference operators and working along the direction of the ideas given in , we derive the following second order approximation for the Riemann-Liouville fractional derivatives.
Let , and its Fourier transform belong to , and define the left WSGD operator by
then, for an integer , we have
uniformly for , where are real numbers and satisfy the following linear system
With the Fourier transform of the right Riemann-Liouville fractional derivative 
by the similar arguments to the above, we get
Let , and its Fourier transform belong to , and define the right WSGD operator by
then, for an integer , we have
uniformly for , where are real numbers and satisfy
Here are determined by solving the above linear system (17).
Note that the desired second order accuracy is obtained by weighting the series expansion of the trigonometric polynomials in frequency space, i.e.,
For the similar method for the classical derivatives one can refer to . And following this idea we can get higher order accurate difference approximations for the Riemann-Liouville fractional derivatives. In the following sections, we focus on the second order difference approximations with three parameters .
Supposing that a function has support on the finite domain , then we have
If the function is non-periodic we are more interested in that are integers and satisfy for numerical consideration. For the convenience, we usually take in practical computations. Here are determined by solving the above system of linear equations (9). When , the linear system (9) has the unique solution . However, for , using the knowledge of linear algebra, we know that the system (9) has infinitely many solutions. Now we assume that is given in linear system of equations (9) with , then we have
If is given in linear algebra system (9), then we have
And if is known, then the solution of the linear system (9) gives
Furthermore, if the integers satisfy , then the -WSGD operators in finite domain can be written as the following form
where the parameters satisfy the following linear system
With the help of the knowledge of linear algebra, the solutions of the system of the linear algebraic equations (26) can be collected by the following three sets
It produces a second order approximation if we take any choice of in . Particularly, if we take for in , it recovers the second order approximations presented in . After rearranging the weights , the second order approximations for the Riemann-Liouville fractional derivatives at point give
where the weights are
3 Second order WSGD schemes: one dimensional case
We first consider the initial boundary value problem of the following one dimensional advection diffusion equation
where can stand for the concentration of a solute at a point at time , is the source term, and are, respectively, the advection and diffusion coefficients, and and are the left and right Riemann-Liouville fractional derivatives of order , respectively, defined by (1) and (2) with .
For the numerical approximations, we define , is the equidistant grid size in space, for , so that . Denote In the space discretizations, we choose the WSGD operators and to approximate the Riemann-Liouville fractional derivatives and , respectively. And the first order space derivative is approximated by the standard central difference. If the time derivative is approximated by the implicit Euler discretization, we have
Dropping the truncation error term, we get the implicit WSGD scheme
To improve the accuracy of the numerical scheme in time direction, using the Crank-Nicolson time discretization, we get the following Crank-Nicolson-WSGD scheme
3.1 Stability analysis
To study the stability of the new scheme, we use the von Neumann linear stability analysis [2, 22]; and assume that the solution of the problem (32) can be zero extended to the whole real line . Letting be the perturbation error and reasonably replacing the symbols and in schemes (33) and (34) with (since should be zero beyond the considered domain), we obtain the following error equations
Let be the solution of (35), where , is the amplitude at time level and is the phase angle with wavelength . We need to show that the amplification factor satisfies the relation for all in . In fact, by substituting the expressions of and into (35), we get the amplification factor of the implicit WSGD difference scheme
In light of the relation (31), we get
By a straightforward calculation, the trigonometric polynomial gives
In view of is a real-valued and even function for given , we just consider its principal value on . Furthermore, using the relation
after a straightforward calculation yields
Then the corresponding amplification factor is obtained as
Due to the functions is non-positive, we get
which means that the implicit WSGD difference scheme is unconditionally stable. ∎