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

Abstract

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.

1. 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–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 force-free case is usually written in the following way [11]: 𝜕𝑢(𝑥,𝑡)()()𝜕𝑡=−𝑉𝑥,𝑡𝜕𝑢𝑥,𝑡𝜕𝑥+𝐷₊()𝜕𝑥,𝑡^𝛼𝑢(𝑥,𝑡)𝜕₊𝑥^𝛼+𝐷₋()𝜕𝑥,𝑡^𝛼𝑢(𝑥,𝑡)𝜕₋𝑥^𝛼(),𝑢()+𝑠𝑥,𝑡𝐿≤𝑥≤𝑅,0<𝑡≤𝑇,𝑥,𝑡=0=𝑢₀(𝑥)𝑢()(),𝐿≤𝑥≤𝑅,𝐿,𝑡=𝑢𝑅,𝑡=0,0≤𝑡≤𝑇, (1) 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 𝑢₀(𝑥) 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]: 𝜕^𝛼𝑓(𝑥)𝜕₊𝑥^𝛼=1Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛𝑥_𝐿𝑓(𝜉)𝑑𝜉(𝑥−𝜉)^{𝛼−𝑛+1},𝜕^𝛼𝑓(𝑥)𝜕₋𝑥^𝛼=(−1)^𝑛Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛𝑅_𝑥𝑓(𝜉)𝑑𝜉(𝜉−𝑥)^{𝛼−𝑛+1}, (2) 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.

2. Fractional Weighted Average Methods

To present the new finite difference method, we give some notations: Δ𝑡 is the time step, Δ𝑥 is the spatial step, the coordinates of the mesh points are 𝑥_𝑗=𝐿+𝑗Δ𝑥, 𝑗=0,1,2,…,𝑁, 𝑁=(𝑅−𝐿)/Δ𝑥, and 𝑡_𝑚=𝑚Δ𝑡, 𝑚=0,1,2,…,𝑀, 𝑀=𝑇/Δ𝑡, and the values of the solution 𝑢(𝑥,𝑡) at these grid points are 𝑢(𝑥_𝑗,𝑡_𝑚)≡𝑢^𝑚_𝑗≃𝑈^𝑚_𝑗, where we denote by 𝑈^𝑚_𝑗 the numerical estimate of the exact value of 𝑢(𝑥,𝑡) at the point (𝑥_𝑗,𝑡_𝑚). Define 𝑉_𝑗^𝑚+1/2=𝑉(𝑥_𝑗,𝑡_𝑚+1/2), 𝐷_𝑗^𝑚+1/2=𝐷(𝑥_𝑗,𝑡_𝑚+1/2) and 𝑠_𝑗^𝑚+1/2=𝑠(𝑥_𝑗,𝑡_𝑚+1/2).

The centered time difference scheme is [20] 𝜕𝑢𝜕𝑡_{(𝑥𝑗,𝑡𝑚+1/2)}=𝑢_𝑗^𝑚+1−𝑢^𝑚_𝑗()Δ𝑡+𝑂Δ𝑡², (3) and the backward space difference scheme is 𝜕𝑢𝜕𝑥_{(𝑥𝑗,𝑡𝑚)}=𝑢^𝑚_𝑗−𝑢^𝑚_𝑗−1().Δ𝑥+𝑂Δ𝑥 (4)

According to the shifted Grünwald-Letnikov definition [8], the definition (2) can be written as 𝜕^𝛼𝑢(𝑥)𝜕₊𝑥^𝛼_{(𝑥𝑗,𝑡𝑚)}=1ℎ^𝛼𝑗+1_𝑘=0𝑔_𝑘^(𝛼)𝑢^𝑚_{𝑗−𝑘+1}(ℎ),𝜕+𝑂_𝛼𝑢(𝑥)𝜕₋𝑥^𝛼_{(𝑥𝑗,𝑡𝑚)}=1ℎ^{𝛼𝑁−𝑗+1}_𝑘=0𝑔_𝑘^(𝛼)𝑢^𝑚_{𝑗+𝑘−1}+𝑂(ℎ). (5) Here, 𝑔_𝑘^(𝛼)=(−1)^𝑘(^𝛼_𝑘) can be evaluated recursively: 𝑔₀^(𝛼)=1,𝑔_𝑘^(𝛼)=1−𝛼+1𝑘𝑔^(𝛼)_𝑘−1. (6)

In the weighted average method, (1) can be evaluated at the intermediate point of the grid (𝑥_𝑗,𝑡_𝑚+1/2) by the following formula: 𝜕𝑢(𝑥,𝑡)𝜕𝑡_{(𝑥𝑗,𝑡𝑚+1/2)}=−𝑉𝑥_𝑗,𝑡_𝑚+1/2𝜆𝜕𝑢(𝑥,𝑡)𝜕𝑥_{(𝑥𝑗,𝑡𝑚)}+()()1−𝜆𝜕𝑢𝑥,𝑡𝜕𝑥_{(𝑥𝑗,𝑡𝑚+1)}+𝐷₊𝑥_𝑗,𝑡_𝑚+1/2𝜆𝜕^𝛼𝑢(𝑥,𝑡)𝜕₊𝑥^𝛼_{(𝑥𝑗,𝑡𝑚)}+()𝜕1−𝜆^𝛼𝑢(𝑥,𝑡)𝜕₊𝑥^𝛼_{(𝑥𝑗,𝑡𝑚+1)}+𝐷₋𝑥_𝑗,𝑡_𝑚+1/2𝜆𝜕^𝛼𝑢(𝑥,𝑡)𝜕₊𝑥^𝛼_{(𝑥𝑗,𝑡𝑚)}+()𝜕1−𝜆^𝛼𝑢(𝑥,𝑡)𝜕₊𝑥^𝛼_{(𝑥𝑗,𝑡𝑚+1)}+𝑠(𝑥,𝑡)|_{(𝑥𝑗,𝑡𝑚+1/2)}, (7) where 0≤𝜆≤1 is the weighting coefficient.

Applying (3)~(5) to (7), letting ℎ=Δ𝑥, and neglecting the truncation error, we get the FWA difference scheme 𝑈_𝑗^𝑚+1+(1−𝜆)−𝑟_𝑗^𝑚+1/2𝑈^𝑚+1_𝑗−1+𝑟_𝑗^𝑚+1/2𝑈_𝑗^𝑚+1−𝜉_𝑗^{𝑚+1/2𝑗+1}_𝑘=0𝑔_𝑘^(𝛼)𝑈^𝑚+1_{𝑗−𝑘+1}−𝜂_𝑗^{𝑚+1/2𝑁−𝑗+1}_𝑘=0𝑔_𝑘^(𝛼)𝑈^𝑚+1_{𝑗+𝑘−1}=𝑈^𝑚_𝑗−𝜆−𝑟_𝑗^𝑚+1/2𝑈^𝑚_𝑗−1+𝑟_𝑗^𝑚+1/2𝑈^𝑚_𝑗−𝜉_𝑗^{𝑚+1/2𝑗+1}_𝑘=0𝑔_𝑘^(𝛼)𝑈^𝑚_{𝑗−𝑘+1}−𝜂_𝑗^{𝑚+1/2𝑁−𝑗+1}_𝑘=0𝑔_𝑘^(𝛼)𝑈^𝑚_{𝑗+𝑘−1}+Δ𝑡𝑠_𝑗^𝑚+1/2,𝑗=1,2,…,𝑁−1,𝑚=0,1,2,…,𝑀−1, (8) where 𝑟_𝑗^𝑚+1/2=𝑉_𝑗^𝑚+1/2Δ𝑡/ℎ, 𝜉_𝑗^𝑚+1/2=𝐷^𝑚+1/2_+,𝑗Δ𝑡/ℎ^𝛼, 𝜂_𝑗^𝑚+1/2=𝐷^𝑚+1/2_−,𝑗Δ𝑡/ℎ^𝛼, and the initial values are calculated by 𝑈_𝑗⁽⁰⁾=𝑢₀(𝑥_𝑗), 𝑗=1,2,…,𝑁−1. Generally, the quantity 𝑟_𝑗^𝑚+1/2 is called the Courant (or CFL) number, the 𝜉_𝑗^𝑚+1/2 and 𝜂_𝑗^𝑚+1/2 are associated with the diffusion coefficients.

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.

3. 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: 𝑔₀^(𝛼)=1,𝑔₁^(𝛼)=−𝛼<0,1≥𝑔₂^(𝛼)≥𝑔₃^(𝛼)≥⋯≥0,_∞𝑘=0𝑔_𝑘^(𝛼)=0,_𝑚𝑘=0𝑔_𝑘^(𝛼)≤0(𝑚≥1). (9)

Theorem 2. When 0≤𝜆≤1/2, the FWA (8) is unconditionally stable, based on the shifted Grünwald approximation (5) to the fractional equation (1) with 1<𝛼≤2. When 1/2<𝜆≤1, the FWA (8) is conditionally stable if 𝛼Δ𝑡𝐷_maxℎ^𝛼+Δ𝑡𝑉_maxℎ≤1,2𝜆−1 (10) where 𝑉_max=max_{𝐿≤𝑥≤𝑅,0≤𝑡≤𝑇}𝑉(𝑥,𝑡) and 𝐷_max=max_{𝐿≤𝑥≤𝑅,0≤𝑡≤𝑇}𝐷(𝑥,𝑡).

Proof. The FWA scheme (8) can be rewritten as [𝐼+(1−𝜆)𝐴]𝑈^𝑚+1=(𝐼−𝜆𝐴)𝑈^𝑚, 𝑚=0,1,2,…,𝑀−1; here, 𝑈^𝑚=[𝑈^𝑚₀,𝑈^𝑚₁,𝑈^𝑚₂,…,𝑈^𝑚_𝑁]^𝑇, 𝐴=(𝑎_𝑖𝑗), 𝑖,𝑗=0,1,2,…,𝑁.
The matrix entries 𝑎_𝑖,𝑗 for 𝑖=1,2,…,𝑁−1 and 𝑗=0,1,…,𝑁 are defined by 𝑎_𝑖,𝑗=𝑟_𝑖^𝑚+1/2−𝜉_𝑖^𝑚+1/2+𝜂_𝑖^𝑚+1/2𝑔₁^(𝛼),𝑗=𝑖−𝑟_𝑖^𝑚+1/2−𝜉_𝑖^𝑚+1/2𝑔₂^(𝛼)−𝜂_𝑖^𝑚+1/2𝑔₀^(𝛼),𝑗=𝑖−1−𝜉_𝑖^𝑚+1/2𝑔₀^(𝛼)−𝜂_𝑖^𝑚+1/2𝑔₂^(𝛼),𝑗=𝑖+1−𝜉_𝑖^𝑚+1/2𝑔^(𝛼)_{𝑖−𝑗+1},𝑗<𝑖−1−𝜂_𝑖^𝑚+1/2𝑔^(𝛼)_{𝑗−𝑖+1},𝑗>𝑖+1 (11) while 𝑎_0,𝑗=𝑎_𝑁,𝑗=0, for 𝑗=0,1,…,𝑁.
According to Lemma 1 and the Gerschgorin theorem, the eigenvalues of the matrix 𝐴 (noted 𝜔_𝑖) are in the disks centered at 𝑎_𝑖,𝑖=𝑟_𝑖^𝑚+1/2−(𝜉_𝑖^𝑚+1/2+𝜂_𝑖^𝑚+1/2)𝑔₁^(𝛼), with radius 𝑅_𝑖=^𝑁_{𝑗=0,𝑗≠𝑖}𝑎_𝑖,𝑗=𝑟_𝑖^𝑚+1/2+𝜉_𝑖^{𝑚+1/2𝑖+1}_{𝑘=0,𝑘≠1}𝑔_𝑘^(𝛼)+𝜂_𝑖^{𝑚+1/2𝑁−𝑖+1}_{𝑘=0,𝑘≠1}𝑔_𝑘^(𝛼)≤𝑟_𝑖^𝑚+1/2−𝜉_𝑖^𝑚+1/2+𝜂_𝑖^𝑚+1/2𝑔₁^(𝛼). (12) Therefore, we have 0≤𝜔_𝑖≤2𝑟_𝑖^𝑚+1/2−𝜉_𝑖^𝑚+1/2+𝜂_𝑖^𝑚+1/2𝑔₁^(𝛼)𝐷=2^𝑖𝑚+1/2𝛼Δ𝑡ℎ^𝛼+𝑉_𝑖^𝑚+1/2Δ𝑡ℎ. (13)
Next, note that 𝜔_𝑖 is an eigenvalue of 𝐴 if and only if (1−𝜆𝜔_𝑖)/(1+(1−𝜆)𝜔_𝑖) is an eigenvalue of the matrix (𝐼−𝜆𝐴)[𝐼+(1−𝜆)𝐴]⁻¹. Because of 0≤𝜔_𝑖 and 0≤𝜆≤1, we get 1−𝜆𝜔_𝑖1+(1−𝜆)𝜔_𝑖𝜔=1−^𝑖1+(1−𝜆)𝜔_𝑖≤1. (14)
In addition, (1−𝜆𝜔_i)/(1+(1−𝜆)𝜔_i)≥−1 as long as (2𝜆−1)𝜔_𝑖≤2.
Hence, when 0≤𝜆≤1/2, we can find that −1≤(1−𝜆𝜔_𝑖)/(1+(1−𝜆)𝜔_𝑖)≤1 always holds; that is, |(1−𝜆𝜔_𝑖)/(1+(1−𝜆)𝜔_𝑖)|≤1. Then, the FWA (8) is unconditionally stable. On the other hand, when 1/2<𝜆≤1, from (13) and 𝜔_𝑖≤2/(2𝜆−1), we get the limited condition (𝛼Δ𝑡𝐷_max/ℎ^𝛼)+(Δ𝑡𝑉_max/ℎ)≤1/(2𝜆−1) where 𝑉_max=max_{𝐿≤𝑥≤𝑅,0≤𝑡≤𝑇}𝑉(𝑥,𝑡) and 𝐷_max=max_{𝐿≤𝑥≤𝑅,0≤𝑡≤𝑇}𝐷(𝑥,𝑡). Therefore, the FWA (8) is conditionally stable.

Let 𝑆=𝛼Δ𝑡𝐷_maxℎ^𝛼+Δ𝑡𝑉_maxℎ, (15) 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 (Δ𝑡)² 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.

4. Numerical Simulations

In this section, we apply the FWA scheme (8) to solve the two-sided space-fractional convection-diffusion equation (1) with 𝛽=0, 𝑉(𝑥,𝑡)=𝑉, 𝐷(𝑥,𝑡)=𝐷, and 𝑠(𝑥,𝑡)=0; the initial condition is 𝑢₀(𝑥)=20𝜋₀^+∞cos(𝑥−0.1𝑉)𝜉𝑒^{0.1𝐷cos(𝜋𝛼/2)𝜉𝛼}𝑑𝜉. (16)

In this case, the analytical solution of (1) solved by the Fourier transform methods is [12] 𝑢()=𝑥,𝑡20𝜋₀^+∞cos(𝑥−𝑉(𝑡+0.1))𝜉𝑒^{𝐷(𝑡+0.1)cos(𝜋𝛼/2)𝜉𝛼}𝑑𝜉. (17)

In the following numerical experiments, the data are chosen as follow: 𝛼=1.9, 𝐷=2, 𝑉=2, 𝑇=2.5, 𝐿=−5, and 𝑅=15.

The numerical solutions are obtained from the FWA scheme (8) discussed above, with different 𝜆, Δ𝑡, 𝑡, 𝑆, and ℎ. From (15), the values of Δ𝑡 for fixed 𝑆 and Δ𝑥=ℎ are 𝑆Δ𝑡=𝛼𝐷_max/ℎ^𝛼+𝑉_max=𝑆/ℎ𝛼𝐷/ℎ^𝛼+().𝑉/ℎ (18)

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 shows numerical solutions obtained by the FWA method (8) with 𝜆=0.9, Δ𝑥=1/40, and small 𝑆=1.2 after 100, 1000, and 4000 time steps. The numerical solutions compare well to the exact solutions, which proves that the FWA method is stable. At the moment, we gain the very small time step Δ𝑡=2.8×10⁻⁴ 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.

Figure 1 

Numerical solutions of (1) by means of the FWA method for 𝜆=0.9, Δ𝑥=1/40, and 𝑆=1.2. The numerical solutions are shown after 100 (dots), 1000 (stars), and 4000 (circles) time steps. The lines correspond to the exact solutions.

Figure 2 

The same as Figure 1 but for 𝑆=1.3. The errors between numerical solution and exact solution after 1000 time steps are shown by line.

Figure 3 

Numerical solutions of (1) by means of the FCN method for Δ𝑥=1/40 and 𝑆=100. The numerical solutions are shown after 10 (dots), 50 (stars), and 100 (circles) time steps. The lines correspond to the exact solutions.

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⁻² calculated from (18), which is much larger than Δ𝑡=2.8×10⁻⁴ 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.

5. 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 (Δ𝑡)² (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.

Acknowledgments

This research was supported by the National Natural Science Foundations of China (Grants nos. 11126179 and 11226247), the 211 Project of Anhui University (nos. 02303319 and 12333010266), the Scientific Research Award for Excellent Middle-Aged and Young Scientists of Shandong Province (no. BS2010HZ012), and the Nature Science Foundation of Anhui Provincial (no. 1308085QA15). The authors acknowledge the anonymous reviewers for their helpful comments.