A Weighted Average Finite Difference Method for the Fractional Convection-Diffusion Equation

A weighted average finite difference method for solving the two-sided space-fractional convection-diffusion equation is given, which is an extension of the weighted average method for ordinary convection-diffusion equations. Stability, consistency, and convergence of the new method are analyzed. A simple and accurate stability criterion valid for this method, arbitrary weighted factor, and arbitrary fractional derivative is given. Some numerical examples with known exact solutions are provided.


Introduction
The history of the fractional derivatives and integrals can date back to the 17th century. However, only after 124 years later, Lacroix first put forward a result of the simplest fractional calculus. Nowadays, the fractional derivatives and integrals have many important applications in various fields of physics [1][2][3], finance [4,5], hydrology [6], engineering [7], mathematics [8], science, and so forth.
Anomalous diffusion is perhaps the most frequently studied complex problem. Classical (integer-order) partial differential equation of diffusion and wave has been extended to the equation with fractional time and/or space by means of fractional operator [9]. Furthermore, it has been extended to the problems of every kind of nonlinear fractional differential equations, and to present the solutions to the problems of initial and boundary values subject to above equations is another rapidly developing field of fractional operator applications. In general, all of these equations have important background of practice applications, such as dispersion in fractals and porous media [10], semiconductor, turbulence, and condensed matter physics.
As a special case of anomalous diffusion, the two-sided space-fractional convection-diffusion equation for the forcefree case is usually written in the following way [11]: where ( , ) > 0 is the drift of the process, that is, the mean advective velocity, is the order of fractional differentiation, + ( , ) = (1 + ) ( , )/2, − ( , ) = (1 − ) ( , )/2, ( , ) > 0 is the coefficient of dispersion, and −1 ≤ ≤ 1 indicates the relative weight of forward versus backward transition probability. The function 0 ( ) is the initial condition, the boundary conditions are zero Dirichlet boundary conditions, and the function ( , ) is a source/sink term. The ( , )/ + and ( , )/ − in (1) are the Riemann-Liouville fractional derivatives. Equation (1) is a special case of the space-fractional Fokker-Planck equation, which more adequately describes the movement of solute in an aquifer than the traditional second-order Fokker-Planck equation.
The left-sided (+) and the right-sided (−) fractional derivatives in (1) are the Riemann-Liouville fractional derivatives of order of a function ( ) for ∈ [ , ] defined by [12]: Advances in Mathematical Physics where − 1 < ≤ ( is an integer) and Γ(⋅) is the Gamma function.
In some cases, there are some methods to solve fractional partial differential equations and get the analytical solutions [12], such as Fourier transform methods, Laplace transform methods, Mellin transform methods, the method of images, and the method of separation of variables. In this paper, the exact solution of (1) can be obtained by Fourier transform methods. However, as in the cases of integer-order differential equations, there are only very few cases of fractional partial differential equations in which the closed-form analytical solutions are available. Therefore, numerical means have to be used in general.
Many of researches on the numerical methods for solving fractional partial differential equations have been proposed, for example, L2 or L2C methods [13], standard or shifted Grünwald-Letnikov formulae [14], convolution formulae [15], homotopy perturbation method, and so forth. For example, Langlands and Henry [16] use L1 scheme form [8] to discretize the Riemann-Liouville fractional time derivative of order between 1 and 2. Yuste [17] considered a Grünwald-Letnikov approximation for the Riemann-Liouville time fractional derivative and used a weighted average for the second-order space derivative. Lin and Xu [18] proposed the method based on a finite difference scheme in time, Legendre's spectral method in space, and so on.
In this paper, based on shifted Grünwald-Letnikov formula, we consider a fractional weighted average (FWA) finite difference method, which is very close to the classical WA methods for ordinary (nonfractional) partial differential equations. The FWA method has some better properties than the fractional explicit and full implicit methods [19], such as higher-order accuracy in time step when weighting coefficient = 1/2.
The rest of this paper is organized as follows. In Section 2, the FWA finite difference method is developed. The stability and convergence of the method are proved in Section 3. Some numerical examples are given in Section 4. Finally, we draw our conclusions in Section 5.
The centered time difference scheme is [20] ( , +1/2 ) and the backward space difference scheme is According to the shifted Grünwald-Letnikov definition [8], the definition (2) can be written as ( ) Here, ( ) = (−1) ( ) can be evaluated recursively: In the weighted average method, (1) can be evaluated at the intermediate point of the grid ( , +1/2 ) by the following formula: where 0 ≤ ≤ 1 is the weighting coefficient.
Obviously, the scheme is explicit when = 1 and the scheme is fully implicit when = 0. particularly, when = 1/2, the FWA scheme is called the fractional Crank-Nicholson (FCN) scheme.

Stability and Accuracy Analysis
In this section, we study the stability of the FWA method and discuss the truncating error. According to our analysis, we can get a conclusion which is similar to the result of classical WA methods. In fact, the following theorem can be viewed as a generalization of these stability conditions for classical WA methods [20]. Lemma 1. The coefficients ( ) given in (6) with 1 < ≤ 2 satisfy the following properties: where max = max ≤ ≤ ,0≤ ≤ ( , ) and max = max ≤ ≤ ,0≤ ≤ ( , ).

Let
the stability limit × is × = 1/(2 − 1). In addition, taking into account (3)∼(5), for arbitrary Δ and Δ , we derive that this method is consistent with a local truncation error (Δ + Δ ), except for the FCN method, whose accuracy is of (Δ ) 2 with respect to the time step [21]. Therefore, according to Lax's equivalence theorem, the FWA method converges at the same rate, too.
Remark 3. Instead of (4), if forward space difference scheme is used, Theorem 2 still holds, and its proof does not change basically. However, if centered space difference scheme is used, we cannot obtain the same conclusion as Theorem 2.
The computational results are shown in Figures 1, 2, and 3. Figures 1 and 2 show two different cases where the FWA method is stable and unstable according to the theoretical predictions of Theorem 2. Figure 1  which proves that the FWA method is stable. At the moment, we gain the very small time step Δ = 2.8 × 10 −4 calculated from (18). Figure 2 has the same assumptions as Figure 1 but for = 1.3 after 1000 time steps, and the large errors between numerical solutions and exact solutions obviously prove that the FWA method is unstable. In the both figures, because of = 0.9, the stability limit is × = 1/(2 − 1) = 1.25. Next, we consider the special case of = 1/2, under the assumption that the FWA method becomes the FCN method. Figure 3 shows numerical solutions obtained by the FCN method with Δ = 1/40 and large = 100 after 10, 50, and 100 time steps. Meanwhile, we can gain the large time step Δ = 2.3 × 10 −2 calculated from (18), which is much larger than Δ = 2.8 × 10 −4 in Figure 1. The numerical solutions approximate well to the exact solutions, and the FCN method is always stable, so it allows the large time steps to be used.

Conclusions
Based on the shifted Grünwald approximation to the fractional derivative, we propose the FWA method in this paper, which can be viewed as a generalization of the classical WA methods for ordinary diffusion equations [17]. The stability of the FWA method depends on weighting parameter , and its accuracy is of order (Δ + Δ ), except for the FCN method, whose accuracy with respect to the time step is of (Δ ) 2 (see [21]).
Obviously, the FCN method is much better and more convenient than the fractional explicit and fully implicit methods because it is not only unconditionally stable, but also of second-order accuracy in time.