A High-Accuracy MOC/FD Method for Solving Fractional Advection-Diffusion Equations

Abstract

Fractional-order diffusion equations are viewed as generalizations of classical diffusion equations, treating super-diffusive flow processes. In this paper, in order to solve the fractional advection-diffusion equation, the fractional characteristic finite difference method is presented, which is based on the method of characteristics (MOC) and fractional finite difference (FD) procedures. The stability, consistency, convergence, and error estimate of the method are obtained. An example is also given to illustrate the applicability of theoretical results.

1. Introduction

The history of fractional calculus is almost as long as integer-order calculus. However, because of lack of application background, fractional calculus was developed very slowly. Up to now, fractional calculus has been applied to almost every field of hydrology [1], physics [2, 3], engineering [4], mathematics [5], and science [6, 7], such as diffusion [8], oscillation, and viscoellastic dynamics. The fractional derivatives in space are often used to model anomalous diffusion because particles spread faster than the classical models predict. The main physical purpose for investigating diffusion equation of fractional order is to describe phenomena of anomalous diffusion in transport progresses through complex or disordered systems including fractal media, and fractional kinetic equations have been proved especially useful in the context of anomalous slow diffusion.

There are several different methods to solve differential equations of fractional order [9], for instance, the method of images, the Fourier transform methods, the method of separation of variables, the Mellin transform methods, and the Laplace transform methods. Wyss [10] considered the time fractional diffusion equation and the solution is given in closed form in terms of the Fox functions. In the last decade, different numerical methods have been developed to solve the fractional diffusion equations [11, 12]. Several explicit finite difference schemes were proposed and analyzed in [13]. Liu et al. [14] used L2 scheme form [5] to discretize the fractional Fokker-Planck equation. Deng [15] proposed finite element method for the space and time fractional Fokker-Planck equation. Meerschaert and Tadjeran [16] considered a Grünwald-Letnikov approximation for the two-sided space-fractional partial differential equations.

Many difficult problems arise in the numerical simulation of fluid flow within porous media in petroleum reservoir simulation and in subsurface contaminant transport and remediation [17]. For example, their mathematical models are basically advection dominated. Fractional derivatives play an important role in modelling particle transport in anomalous diffusion [18]. The space-fractional advection-dominated diffusion equation [19] is 𝜕𝑢(𝑥,𝑡)(𝑥)()𝜕𝑡+𝑣𝜕𝑢𝑥,𝑡(𝑥)𝜕𝜕𝑥=𝑑^𝛼𝑢(𝑥,𝑡)𝜕𝑥^𝛼()𝑢()+𝑠𝑥,𝑡,𝐿≤𝑥≤𝑅,0<𝑡≤𝑇,𝑥,𝑡=0=𝑢₀(𝑥)𝑢()(),𝐿≤𝑥≤𝑅,𝐿,𝑡=0,𝑢𝑅,𝑡=0,0<𝑡≤𝑇, (1) where the function 𝑣(𝑥)>0 is the drift of the process, 𝛼 is the fractional order of the spatial derivative, physical considerations restrict 1<𝛼≤2, the functions 𝑑(𝑥)>0 may be interpreted as the coefficient of dispersion, and 𝑠(𝑥,𝑡) is a source/sink term. Moreover, when (1) is convection dominated diffusion equation, we have 𝑣(𝑥)≫𝑑(𝑥).

Equation (1) contains a Riemann-Liouville fractional derivative of order 𝛼 of a function 𝑓(𝑥), for 𝑥∈[𝐿,𝑅], −∞≤𝐿<𝑅≤∞, defined by [9] 𝜕^𝛼𝑓(𝑥)𝜕𝑥^𝛼=1Γ()𝑑𝑛−𝛼𝑑𝑥^𝑛𝑛𝑥_𝐿𝑓(𝜉)(𝑥−𝜉)^{𝛼+1−𝑛}𝑑𝜉, (2) where 𝑛 is an integer such that 𝑛−1<𝛼≤𝑛, Γ(⋅) is the gamma function.

The standard finite difference methods (SFDM) for solving (1) were discussed in [16, 19, 20], and their stability was given in those papers, too. However, these methods often generate numerical solutions with severe nonphysical oscillations. Because of the nature of hyperbolic types [17, 21], the convection-diffusion problems involve moving sharp fronts or boundary layers. While upwind finite difference methods eliminate nonphysical oscillations, they tend to generate numerical solutions with excessive numerical diffusion and smear the moving steep fronts. The efficient grid refinements are required if the standard methods are used in computations. But for high dimensional large-scale problems, it may lead to very large linear systems and need long computation time.

In this paper, we will consider a new numerical difference scheme proposed by the characteristic method, which is named the fractional characteristic finite difference method (CFDM). It is proved that the new method is unconditionally stable, consistent, and convergent; therefore, there is no stability limitation on the size of Δ𝑡. Moreover, the CFDM allows large time steps to be used without the loss of accuracy, which is much better than the SFDM and other numerical methods. The reason is that, for problems with significant convection, the solution changes much less rapidly in the characteristic 𝜏 direction than in the 𝑡 direction.

The outline of the paper is as follows. In Section 2, we describe the fractional characteristic finite difference schemes for solving (1). In Section 3, the stability, consistency, and convergence are proved, and the maximum error estimate is derived, too. The numerical experiments are given in Section 4. The conclusions of this paper are given in Section 5.

2. Fractional Characteristic Difference Method

In this section, we consider the combination of characteristic methods with finite difference techniques [22] and seek a new numerical method of (1), which is the convection-dominated problem and reflects the almost hyperbolic nature.

2.1. The Idea of Characteristic Methods

Let 𝜓(𝑥)=1+𝑣²(𝑥) and denote the characteristic direction associated with the operator 𝑢_𝑡+𝑣(𝑥)𝑢_𝑥 by 𝜏=𝜏(𝑥) (see [22]), where 𝜕(𝑥)=1𝜕𝜏𝜓(𝑥)𝜕+𝑣(𝑥)𝜕𝑡𝜓(𝑥)𝜕.𝜕𝑥 (3) Then, the left part of (1) can be rewritten as 𝜕𝑢(𝑥,𝑡)(𝑥)()𝜕𝑡+𝑣𝜕𝑢𝑥,𝑡(𝑥)()𝜕𝑥=𝜓𝜕𝑢𝑥,𝑡(𝑥).𝜕𝜏 (4)

Next, in order to formulate the fractional CFDM of (1), we give some notations. Let Δ𝑡 be the temporal mesh or time step, 𝑡^𝑚=𝑚Δ𝑡, 𝑚=0,1,2,…,𝑀, 𝑀=𝑇/Δ𝑡, Δ𝑥 is the spatial mesh, the coordinates of the mesh points are 𝑥_𝑖=𝐿+𝑖Δ𝑥, 𝑖=0,1,2,…,𝑁, 𝑁=(𝑅−𝐿)/Δ𝑥, with 𝑀 and 𝑁 are positive integers. The values of the solution 𝑢(𝑥,𝑡) on these grid points are 𝑢(𝑥_𝑖,𝑡^𝑚)≡𝑢^𝑚_𝑖≃𝑈^𝑚_𝑖, and the numerical estimate of the exact value of 𝑢(𝑥,𝑡) at the point (𝑥_𝑖,𝑡^𝑚) is denoted by 𝑈^𝑚_𝑖. Finally, define 𝑣_𝑖=𝑣(𝑥_𝑖),  𝑑_𝑖=𝑑(𝑥_𝑖), and 𝑆^𝑚_𝑖=𝑠(𝑥_𝑖,𝑡^𝑚).

The characteristic derivative is approximated in the following manner: 𝜓(𝑥)𝜕𝑢𝑥,𝑡^𝑚+1(𝑥)𝑢𝜕𝜏≈𝜓𝑥,𝑡^𝑚+1−𝑢𝑥,𝑡^𝑚(𝑥−𝑥)²+Δ𝑡^21/2=𝑢𝑥,𝑡^𝑚+1−𝑢𝑥,𝑡^𝑚,Δ𝑡 (5) where 𝑥=𝑥−𝑣(𝑥)Δ𝑡.

Then, by (5) the differential operators at the point (𝑥_𝑖,𝑡^𝑚+1) are discretized to be 𝜕𝑢(𝑥,𝑡)(𝑥)()𝜕𝑡+𝑣𝜕𝑢𝑥,𝑡𝜕𝑥_{(𝑥𝑖,𝑡}^𝑚+1₎=𝑢_𝑖^𝑚+1−𝑢^𝑚_𝑖Δ𝑡+𝑂Δ𝑡+minΔ𝑥,Δ𝑥²,Δ𝑡 (6) where 𝑥_𝑖=𝑥_𝑖−𝑣_𝑖Δ𝑡, 𝑢^𝑚_𝑖=𝑢(𝑥_𝑖,𝑡^𝑚) is calculated by the piecewise-linear interpolation.

Let 𝑟_𝑖=𝑣_𝑖Δ𝑡/Δ𝑥, and let symbol [𝑟_𝑖] be the integral part of 𝑟_𝑖. Let 𝑖₀=𝑖−[𝑟_𝑖]−1 and 𝑟^∗_𝑖=𝑟_𝑖−[𝑟_𝑖], using linear interpolation to 𝑢^𝑚_𝑖, we get [23] 𝑢^𝑚_𝑖=1−𝑟^∗_𝑖𝑢^𝑚_𝑖0+1+𝑟^∗_𝑖𝑢^𝑚_𝑖0. (7) Generally, the quantity 𝑟_𝑖 is called the Courant (or CFL) number.

2.2. The Discrete Versions of the Fractional Derivative

First, we develop the right-shifted Grünwald discrete approximation to the fractional derivatives by means of improving the standard form of the Grünwald-Letnikov definition [16]: 𝜕^𝛼𝑓(𝑥)𝜕𝑥^𝛼=lim_𝑀→∞1Δ𝑥^𝛼𝑀_𝑘=0𝑔_𝑘^(𝛼)𝑓(𝑥−(𝑘−1)Δ𝑥), (8) where 𝑀 is a positive integer, Δ𝑥=(𝑥−𝐿)/𝑀, and 𝑔_𝑘^(𝛼)=(−1)^𝑘(^𝛼_𝑘), which can be calculated by 𝑔₀^(𝛼)=1,𝑔_𝑘^(𝛼)=1−𝛼+1𝑘𝑔^(𝛼)_𝑘−1. (9)

Following (5), we can discretize the Riemann-Liouville operator [12] 𝜕^𝛼𝑓(𝑥)𝜕𝑥^𝛼=1Δ𝑥^𝛼[𝑙]+1_𝑘=0𝜔_𝑘^(𝛼)𝑓(𝑥−(𝑘−1)Δ𝑥)+𝑂Δ𝑥^𝑝, (10) where 𝑝 is a positive integer, 𝑙=(𝑥−𝐿)/Δ𝑥, and [𝑙] is the integral part of 𝑙.

Here, according to the definition of 𝜔_𝑘^(𝛼), we can draw the following conclusion [9, 12]: when 𝑝=1, 𝜔_𝑘^(𝛼)=𝑔_𝑘^(𝛼); when 𝑝=2, 𝜔_𝑘^(𝛼)=(3^𝛼/2^𝛼)∑^𝑘_𝑗=0𝑔_𝑗^(𝛼)𝑔^(𝛼)_𝑘−𝑗3^−𝑗.

According to the shifted Grünwald-Letnikov formulae (8) and (10), the definition (2) can be discretized by 𝜕^𝛼𝑢𝜕𝑥^𝛼𝑥_𝑖,𝑡^𝑚+1=1Δ𝑥^𝛼𝑖+1_𝑘=0𝑔_𝑘^(𝛼)𝑢^𝑚+1_{𝑖−𝑘+1}+𝑂(Δ𝑥), (11) with 𝑝=1.

Inserting (6), (7), and (11) into (1), neglecting the truncation error, after rearranging the terms, we finally get the fractional characteristic difference scheme 𝑈_𝑖^𝑚+1−𝜉_𝑖^𝑖+1_𝑘=0𝑔_𝑘^(𝛼)𝑈^𝑚+1_{𝑖−𝑘+1}=1−𝑟^∗_𝑖𝑈^𝑚_𝑖0+1+𝑟^∗_𝑖𝑈^𝑚_𝑖0+Δ𝑡𝑆_𝑖^𝑚+1, (12) where 𝜉_𝑖=𝑑_𝑖Δ𝑡/Δ𝑥^𝛼 is associated with the diffusion coefficient, and the initial values are 𝑈⁰_𝑖=𝑢₀(𝑥_𝑖).

3. Stability and Convergence Analysis

In this section, before studying the stability and convergence of the CFDM, we give the following lemma, which was proved previously [16].

Lemma 1. According to the definition of 𝑔_𝑘^(𝛼) in (9), when 1<𝛼≤2, the coefficients 𝑔_𝑘^(𝛼) have the following properties:(1)𝑔₀^(𝛼)=1, 𝑔₁^(𝛼)=−𝛼<0 and 𝑔_𝑘^(𝛼)≥0, 𝑘≥2,(2)∑^∞_𝑘=0𝑔_𝑘^(𝛼)=0 and ∑^∞_{𝑘=0,𝑘≠1}𝑔_𝑘^(𝛼)=−𝑔₁^(𝛼)=𝛼,(3)∑^𝑚_𝑘=0𝑔_𝑘^(𝛼)≤0, 𝑚≥1.

In the following section, we study the stability and convergence of the fractional CFDM for (1) with the given initial and boundary conditions.

3.1. Stability Analysis

Theorem 2. The fractional characteristic finite difference scheme (12) of (1) with 1<𝛼≤2, based on the shifted Grünwald approximation (8) to the space fractional derivative, is unconditionally stable.

Proof. Incorporating the boundary conditions, (12) is rewritten as 𝐴𝑈^𝑚+1=𝐵𝑈^𝑚+𝑆^𝑚+1, where 𝑈^𝑚=[𝑈^𝑚₁,𝑈^𝑚₂,…,𝑈^𝑚_𝑁−1]^𝑇, 𝑆^𝑚+1=[Δ𝑡𝑆₁^𝑚+1,Δ𝑡𝑆₂^𝑚+1,…,Δ𝑡𝑆^𝑚+1_𝑁−1]^𝑇, 𝐴=[𝑎_𝑖,𝑗] and 𝐵=[𝑏_𝑖,𝑗] are two (𝑁−1)×(𝑁−1) matrices.
When 𝑖=1,2,…,𝑁−1 and 𝑗=1,…,𝑁−1, we have 𝑎_𝑖,𝑗=−𝜉_𝑖𝑔^(𝛼)_{𝑖−𝑗+1},𝑗≤𝑖−1,1−𝜉_𝑖𝑔₁^(𝛼),𝑗=𝑖,−𝜉_𝑖𝑔₀^(𝛼)𝑏,𝑗=𝑖+1,0,𝑗>𝑖+1,_𝑖,𝑗=1−𝑟^∗_𝑖,𝑗=𝑖−𝑟_𝑖,𝑟^∗_𝑖,𝑗=𝑖−𝑟_𝑖−1,0,otherwise, (13) where 0≤𝑟^∗_𝑖≤1, so 0≤1−𝑟^∗_𝑖≤1, too. Thus, 𝐵≥0 holds.
For matrix 𝐵, we can easily get ^𝑁−1_𝑗=1𝑏_𝑖,𝑗=0,𝑖<𝑟_𝑖,1−𝑟^∗_𝑖,𝑖=𝑟_𝑖,1,𝑖>𝑟_𝑖. (14) Therefore, ‖𝐵‖_∞=max_{1≤𝑖≤𝑁−1}^𝑁−1_𝑗=1𝑏_𝑖,𝑗=max_{1≤𝑖≤𝑁−1}^𝑁−1_𝑗=1𝑏_𝑖,𝑗≤1. (15)
For matrix 𝐴, we can easily see that 𝑎_𝑖,𝑗≤0 for all 𝑖≠𝑗, and 𝑎_𝑖,𝑖=1−𝜉_𝑖𝑔₁^(𝛼)=1+𝜉_𝑖𝛼>0. (16) At the same time, according to Lemma 1, we have ^𝑁−1_{𝑗=1,𝑗≠𝑖}𝑎_𝑖,𝑗=𝜉_𝑖^𝑖_{𝑘=0,𝑘≠1}𝑔_𝑘^(𝛼)≤𝜉_𝑖𝛼, (17) for 𝑖=1,2,…,𝑁−1, we can obtain 𝑎_𝑖,𝑖−^𝑁−1_{𝑗=1,𝑗≠𝑖}𝑎_𝑖,𝑗≥1. (18) Therefore, 𝐴 is diagonally dominant by rows. According to the matrix theory [24, 25], 𝐴⁻¹_∞≤1min_{1≤𝑖≤𝑁−1}𝑎_𝑖,𝑖−∑_𝑗≠𝑖𝑎_𝑖,𝑗≤1. (19)
Then, from (15) and (19), we can get 𝐴⁻¹𝐵_∞≤𝐴⁻¹_∞‖𝐵‖_∞≤1. (20) Thus, the characteristic finite difference scheme (12) is unconditionally stable.

In addition, taking into account (6) and (11), we have 𝑇(𝑥,𝑡)=𝑂Δ𝑡+minΔ𝑥,Δ𝑥²(),Δ𝑡=𝑂Δ𝑥+Δ𝑡 (21) so this method is consistent with a local truncation error which is 𝑂(Δ𝑥+Δ𝑡).

3.2. Convergence Analysis

Let 𝑢(𝑥_𝑖,𝑡^𝑚)(𝑖=1,2,…,𝑁−1; 𝑚=0,1,2,…,𝑀) be the exact solution of (1) at mesh point (𝑥_𝑖,𝑡^𝑚). Define 𝑒^𝑚_𝑖=𝑢(𝑥_𝑖,𝑡^𝑚)−𝑈^𝑚_𝑖, 𝑖=1,2,…,𝑁−1; 𝑚=0,1,2,…,𝑀 and 𝐸^𝑚=[𝑒^𝑚₁,𝑒^𝑚₂,…,𝑒^𝑚_𝑁−1]^𝑇. Using 𝐸⁰=0, substitution into (12) leads to 𝐴𝐸^𝑚+1=𝐵𝐸^𝑚+𝑅^𝑚+1, (22) where 𝑅^𝑚=[𝑅^𝑚₁,𝑅^𝑚₂,…,𝑅^𝑚_𝑁−1]^𝑇, and 𝑅^𝑚_𝑖=Δ𝑡𝑇(𝑥_𝑖,𝑡^𝑚), 𝑖=1,2,…,𝑁−1; 𝑚=0,1,2,…,𝑀.

From (6), (11), and (21), we have the following.

Lemma 3. There is a positive constant 𝐶, such that 𝑅^𝑚_𝑖≤𝐶Δ𝑡(Δ𝑥+Δ𝑡), (23) where 𝑖=1,2,…,𝑁−1; 𝑚=0,1,2,…,𝑀.

Let ‖𝐸^𝑚‖_∞=max_{1≤𝑖≤𝑁−1}|𝑒^𝑚_𝑖|, because 𝑚Δ𝑡≤𝑇 is finite, we can obtain the following result.

Theorem 4. Let 𝑈^𝑚_𝑖 be the numerical solution computed by the CFDM (12). Then there is a positive constant 𝐶₁, such that 𝑈^𝑚_𝑖−𝑢𝑥_𝑖,𝑡^𝑚≤𝐶₁(Δ𝑥+Δ𝑡),𝑖=1,2,…,𝑁−1;𝑚=0,1,2,…,𝑀. (24)

Proof. From (22) and 𝐸⁰=0, we can get 𝐸^𝑚+1=𝐴⁻¹𝐵𝐸^𝑚+𝑅^𝑚+1=𝐴⁻¹𝐵𝐴⁻¹𝐵𝐸^𝑚−1+𝑅^𝑚+𝑅^𝑚+1=𝐴⁻¹𝐵^𝑚𝐴⁻¹𝑅¹+𝐴⁻¹𝐵^𝑚−1𝐴⁻¹𝑅²+⋯+𝐴⁻¹𝑅^𝑚+1. (25)
According to (19) and Lemma 3, using (25), we finally get 𝐸^𝑚+1_∞≤𝐴⁻¹𝐵^𝑚_∞𝐴⁻¹_∞𝑅¹_∞+𝐴⁻¹𝐵^𝑚−1_∞𝐴⁻¹_∞𝑅²_∞+⋯+𝐴⁻¹_∞𝑅^𝑚+1_∞≤𝐶Δ𝑡(𝑚+1)(Δ𝑥+Δ𝑡). (26)
Because 𝑚Δ𝑡≤𝑇, from (26), we obtain the error estimate ‖𝐸^𝑚‖_∞≤𝐶𝑇(Δ𝑥+Δ𝑡). Therefore, we get the result |𝑈^𝑚_𝑖−𝑢(𝑥_𝑖,𝑡^𝑚)|≤𝐶₁(Δ𝑥+Δ𝑡), where 𝐶₁=𝐶𝑇 is a positive constant.

Finally, according to Theorem 4, we draw the following conclusion: when Δ𝑡→0 and Δ𝑥→0, |𝑒^𝑚_𝑖|→0 holds; that is, the characteristic finite difference scheme (12) is convergent. In summary, the fractional CFDM for (1) is unconditionally stable, consistent, and convergent.

4. Numerical Simulations

In this section, we give some numerical results for a special case of (1) to confirm our theoretical analysis discussed above. In addition, we also compare our new CFDM with the SFDM and point out the advantages of the new method.

For (1), if we assume that 𝑣(𝑥)=𝑣 and 𝑑(𝑥)=𝑑 are two constant functions, the function 𝑠(𝑥,𝑡)=0, and the initial condition 𝑢₀(𝑥)=𝛿(𝑥) is the dirac delta function. Then, under the zero Dirichlet boundary conditions, the exact solution of (1) could be found by the Fourier transform methods [9]: 𝑢()=1𝑥,𝑡𝜋₀^+∞cos(𝑥−𝑣𝑡)𝜉+𝑑𝑡sin𝜋𝛼2𝜉^𝛼×𝑒^{𝑑𝑡cos(𝜋𝛼/2)𝜉𝛼}𝑑𝜉. (27)

In the numerical experiments, the data are chosen as follows: 𝑣=√2, 𝑇=1.5, 𝐿=−2, and 𝑅=4. For simplicity, we assume that the initial time 𝑡₀=0.1, and the initial values 𝑈⁰=𝑢(𝑥,𝑡₀) are calculated from (27). The numerical solutions are obtained from the CFDM scheme (12) discussed above. The errors of numerical solutions are measured in 𝐿_∞ and 𝐿₂ norms which are defined as follows: 𝐸_∞,Δ𝑡=max_𝑖|𝑈(𝑥(𝑖))(𝑥(𝑖),𝐸,𝑡−𝑢,𝑡)|_2,Δ𝑡=ℎ^𝑁_𝑖=0(𝑈(𝑥(𝑖),𝑡)−𝑢(𝑥(𝑖),𝑡))², (28) where 𝑈(𝑥(𝑖),𝑡) and 𝑢(𝑥(𝑖),𝑡) are the numerical solution and analytical solution at (𝑥(𝑖),𝑡), respectively. ℎ is the space step, and Δ𝑡 is the time step. The ratios of convergence in time are calculated by 𝐸log^𝑙,Δ𝑡1𝐸_𝑙,Δ𝑡2logΔ𝑡₁Δ𝑡₂⁻¹,𝑙=2,∞, (29) where Δ𝑡₁ and Δ𝑡₂ are two time steps, and very small spatial step size Δ𝑥 is taken in computation. Similarly, the ratios of convergence in space can be calculated, too.

4.1. The Properties of the CFDM

In this subsection, we use the CFDM (12) to compute the numerical solutions to (1), and list the computational results in the figures and tables. From Figure 1, we can see that the numerical solutions of the CFDM coincide with the analytic solutions, with the diffusion coefficient 𝑑=0.1, the space step Δ𝑥=1/90, and the time step Δ𝑡=1/30 at 𝑡=1.5, for 𝛼=1.9 and 1.6. It is clearly shown that the CFDM is stable for (1), which allows the large time step to be used and gains a very good numerical approximation to the exact solution.

Figure 1 

Exact solutions and numerical solutions obtained by the CFDM with 𝑑=0.1, Δ𝑥=1/90, and Δ𝑡=1/30 at 𝑡=1.5. The solid lines correspond to the exact solutions. The starred and dotted lines correspond to numerical solutions for 𝛼=1.9 and 1.6, respectively.

Next, we examine the ratio of convergence in space for the CFDM. Table 1 shows 𝐿_∞ and 𝐿₂ errors and ratios in space with the diffusion coefficient 𝑑=0.1 and the small time step Δ𝑡=1/400 at time 𝑡=1.5, for two different fractional orders 𝛼=1.9 and 1.6, respectively. The second and fourth columns show the 𝐿_∞ and 𝐿₂ errors, respectively. Column three shows the ratio of the CFDM in 𝐿_∞ norm, and column five shows the ratio of the CFDM in 𝐿₂ norm. It is clearly shown that the CFDM method is of first-order accuracy in space in both 𝐿_∞ and 𝐿₂ norms.

Table 1 

The errors and ratios in space calculated by CFDM (12) with 𝑑=0.1 and Δ𝑡=1/400 at time 𝑡=1.5, for two different fractional orders 𝛼=1.9 and 1.6, respectively.

Meanwhile, Table 2 shows 𝐿_∞ and 𝐿₂ errors and ratios in time with 𝑑=0.02 and Δ𝑥=1/300 at time 𝑡=1.5, for two different fractional orders 𝛼=1.9 and 1.6, respectively. The table shows that the CFDM method is of first-order accuracy in time in both 𝐿_∞ and 𝐿₂ norms, too. Furthermore, from the table, we can see that even with the very large time steps size of Δ𝑡=1/5 to 1/30, which corresponds to the Courant numbers of 14.14 to 84.85, the CFDM can still generate the accurate numerical solutions without noticeable numerical errors. In conclusion, according to Tables 1 and 2, the accuracy of the CFDM method is of order 𝑂(Δ𝑥+Δ𝑡), which confirms the result of Theorem 4.

Table 2 

The errors and ratios in time calculated by CFDM (12) with 𝑑=0.02 and Δ𝑥=1/300 at time 𝑡=1.5, for two different fractional orders 𝛼=1.9 and 1.6, respectively.

4.2. Comparison with the SFDM

To gain a better understanding of the CFDM, we use the well-developed and well-received SFDM to simulate the same model problem for comparison, which usually contains the implicit and explicit finite difference schemes. We only consider the case of 𝛼=1.9 in this subsection.

For the implicit SFDM, it is also unconditionally stable. Table 3 and Figure 2 both show that the CFDM is much better than the implicit SFDM, if they have the same assumptions. Firstly, Table 3 shows the comparison of 𝐿₂ errors calculated by the CFDM (12) and the implicit SFDM for two different diffusion coefficients 𝑑=0.1 and 𝑑=0.02, with some different space steps Δ𝑥 and time steps Δ𝑡 for 𝛼=1.9 at time 𝑡=1.5. From the table, we can easily see that the 𝐿₂ errors generated by the CFDM is obviously much smaller than the 𝐿₂ errors generated by the implicit SFDM, which has one- even two-order higher accuracy than the latter.

Table 3 

Comparison of 𝐿₂ errors calculated by the CFDM (12) and the implicit SFDM with 𝛼=1.9, different space steps Δ𝑥 and time steps Δ𝑡 at time 𝑡=1.5, for two different diffusion coefficients 𝑑=0.1 and 𝑑=0.02.

Figure 2 

The 𝐿∞ errors between exact solutions and numerical solutions calculated by the CFDM (12) and the implicit SFDM at different time levels with 𝛼=1.9, Δ𝑥=1/50, and Δ𝑡=1/40 for 𝑑=0.1.

Moreover, Figure 2 shows that the 𝐿_∞ errors between exact solutions and numerical solutions calculated by the CFDM and the SFDM at different time levels with 𝛼=1.9, the diffusion coefficient 𝑑=0.1, the space step Δ𝑥=1/50, and the time step Δ𝑡=1/40. It is easily to see that the CFDM can generate accurate numerical solution even if relatively large time step is used, and this significantly improves efficiency [17, 26]. However, the implicit SFDM does not have this kind of property.

Next, we compare the CFDM and the implicit SFDM numerical solutions with 𝛼=1.9 and 𝑑=0.02 at time 𝑡=1.5 in Figure 3. To keep similar 𝐿_∞ error, for the CFDM, the grid sizes are chosen as Δ𝑡=1/30 and Δ𝑥=1/100, and it generates the 𝐿_∞ error of 1.2168×10⁻². However, for the implicit SFDM, we have to choose very small grid sizes, which are at least Δ𝑡=1/5000 and Δ𝑥=1/1000. In this case, the implicit SFDM generates the 𝐿_∞ error of 4.6630×10⁻², which is still larger than that of the CFDM. As a result, it is evidently shown that the CFDM can allow large time steps and space steps without the loss of accuracy.

Figure 3 

Numerical solutions and exact solution at 𝑡=1.5 with 𝑑=0.02. The solid line corresponds to the exact solution and the dotted line corresponds to numerical solution of the CFDM with 𝛼=1.9, Δ𝑥=1/100, and Δ𝑡=1/30. The stared line corresponds to numerical solution of the implicit SFDM with Δ𝑥=1/1000 and Δ𝑡=1/5000.

In Section 4 discussed above, we only compare the numerical results of the CFDM with the implicit SFDM. In fact, the SFDM for (1) also contains several other explicit schemes which are conditionally stable. There are the upwind explicit scheme, central explicit scheme, and the Lax-Wendroff explicit scheme, and so forth. Of course, our new fractional CFDM is obviously better than all of them for the space-fractional convection-dominated diffusion problems because the explicit SFDM schemes are conditionally stable and they does not allow large time steps.

Finally, we compare the CFDM with the explicit SFDM schemes simply. Table 4 shows that the 𝐿₂ errors of three different explicit schemes with different Δ𝑥 and Δ𝑡 for 𝑑=0.02 and 𝛼=1.9 at time 𝑡=1.5. As a result, the much refined time steps and space steps have to be used to improve the accuracy of the numerical solution in Table 4, in order to obtain a numerical solution with comparable accuracy to the numerical solution by the CFDM with the relatively coarse spatial grids and time steps in Tables 2 or 3. Moreover, under the same assumptions of Table 3, the three explicit SFDM schemes are all unstable.

Table 4 

The 𝐿₂ errors calculated by the upwind, central and Lax-Wendroff explicit difference schemes with different space steps Δ𝑥 and time steps Δ𝑡 for 𝑑=0.02 and 𝛼=1.9 at time 𝑡=1.5.

5. Conclusions

In actual applications, we often use the space-fractional derivatives to model convection-diffusion problems, and in many diffusion processes arising in physical problems, convection essentially dominates diffusion. In this paper, we propose a new scheme called the fractional CFDM and prove that it is unconditionally stable, consistent, and convergent. For this new method, even though the time steps are very coarse, the resulting CFDM solutions are more accurate than the SFDM solutions with much finer temporal grid sizes. Thus, the CFDM is more efficient and superior to the SFDM, especially for the convection-dominated problems.

Acknowledgments

This research was supported by the National Natural Science Foundations of China Grant (no. 11126179 and no. 11226247), the 211 project of Anhui University (no. 02303319 and no. KJTD002B), 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.