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

Authors: Lijuan Su; Pei Cheng
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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
(1)∂u(x,t)∂t+v(x)∂u(x,t)∂x=d(x)∂αu(x,t)∂xα+s(x,t),L≤x≤R,0<t≤T,u(x,t=0)=u0(x),L≤x≤R,u(L,t)=0,u(R,t)=0,0<t≤T,
where the function v(x)>0 is the drift of the process, α is the fractional order of the spatial derivative, physical considerations restrict 1<α≤2, the functions d(x)>0 may be interpreted as the coefficient of dispersion, and s(x,t) is a source/sink term. Moreover, when (1) is convection dominated diffusion equation, we have v(x)≫d(x).Equation (1) contains a Riemann-Liouville fractional derivative of order α of a function f(x), for x∈[L,R], -∞≤L<R≤∞, defined by [9]
(2)∂αf(x)∂xα=1Γ(n-α)dndxn∫Lx‍f(ξ)(x-ξ)α+1-ndξ,
where n is an integer such that n-1<α≤n, Γ(·) 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Δt. 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 t direction.The outline of the paper is as follows. In Section2, 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ψ(x)=1+v2(x) and denote the characteristic direction associated with the operator ut+v(x)ux by τ=τ(x) (see [22]), where
(3)∂∂τ(x)=1ψ(x)∂∂t+v(x)ψ(x)∂∂x.
Then, the left part of (1) can be rewritten as
(4)∂u(x,t)∂t+v(x)∂u(x,t)∂x=ψ(x)∂u(x,t)∂τ(x).Next, in order to formulate the fractional CFDM of (1), we give some notations. Let Δt be the temporal mesh or time step, tm=mΔt, m=0,1,2,…,M, M=T/Δt, Δx is the spatial mesh, the coordinates of the mesh points are xi=L+iΔx, i=0,1,2,…,N, N=(R-L)/Δx, with M and N are positive integers. The values of the solution u(x,t) on these grid points are u(xi,tm)≡uim≃Uim, and the numerical estimate of the exact value of u(x,t) at the point (xi,tm) is denoted by Uim. Finally, define vi=v(xi),  di=d(xi), and Sim=s(xi,tm).The characteristic derivative is approximated in the following manner:(5)ψ(x)∂u(x,tm+1)∂τ≈ψ(x)u(x,tm+1)-u(x-,tm)[(x-x-)2+Δt2]1/2=u(x,tm+1)-u(x-,tm)Δt,
where x-=x-v(x)Δt.Then, by (5) the differential operators at the point (xi,tm+1) are discretized to be
(6)[∂u(x,t)∂t+v(x)∂u(x,t)∂x](xi,tm+1)=uim+1-u-imΔt+O(Δt+min(Δx,Δx2Δt)),
where x-i=xi-viΔt, u-im=u(x-i,tm) is calculated by the piecewise-linear interpolation.Letri=viΔt/Δx, and let symbol [ri] be the integral part of ri. Let i0=i-[ri]-1 and ri*=ri-[ri], using linear interpolation to u-im, we get [23]
(7)u-im=(1-ri*)ui0+1m+ri*ui0m.
Generally, the quantity ri 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]:
(8)∂αf(x)∂xα=limM→∞1Δxα∑k=0M‍gk(α)f(x-(k-1)Δx),
where M is a positive integer, Δx=(x-L)/M, and gk(α)=(-1)k(αk), which can be calculated by
(9)g0(α)=1,gk(α)=(1-α+1k)gk-1(α).Following (5), we can discretize the Riemann-Liouville operator [12]
(10)∂αf(x)∂xα=1Δxα∑k=0[l]+1‍ωk(α)f(x-(k-1)Δx)+O(Δxp),
where p is a positive integer, l=(x-L)/Δx, and [l] is the integral part of l.Here, according to the definition ofωk(α), we can draw the following conclusion [9, 12]: when p=1, ωk(α)=gk(α); when p=2, ωk(α)=(3α/2α)∑j=0k‍gj(α)gk-j(α)3-j.According to the shifted Grünwald-Letnikov formulae (8) and (10), the definition (2) can be discretized by
(11)∂αu∂xα(xi,tm+1)=1Δxα∑k=0i+1‍gk(α)ui-k+1m+1+O(Δx),
with p=1.Inserting (6), (7), and (11) into (1), neglecting the truncation error, after rearranging the terms, we finally get the fractional characteristic difference scheme
(12)Uim+1-ξi∑k=0i+1‍gk(α)Ui-k+1m+1=(1-ri*)Ui0+1m+ri*Ui0m+ΔtSim+1,
where ξi=diΔt/Δxα is associated with the diffusion coefficient, and the initial values are Ui0=u0(xi).


## 2.1. The Idea of Characteristic Methods

Letψ(x)=1+v2(x) and denote the characteristic direction associated with the operator ut+v(x)ux by τ=τ(x) (see [22]), where
(3)∂∂τ(x)=1ψ(x)∂∂t+v(x)ψ(x)∂∂x.
Then, the left part of (1) can be rewritten as
(4)∂u(x,t)∂t+v(x)∂u(x,t)∂x=ψ(x)∂u(x,t)∂τ(x).Next, in order to formulate the fractional CFDM of (1), we give some notations. Let Δt be the temporal mesh or time step, tm=mΔt, m=0,1,2,…,M, M=T/Δt, Δx is the spatial mesh, the coordinates of the mesh points are xi=L+iΔx, i=0,1,2,…,N, N=(R-L)/Δx, with M and N are positive integers. The values of the solution u(x,t) on these grid points are u(xi,tm)≡uim≃Uim, and the numerical estimate of the exact value of u(x,t) at the point (xi,tm) is denoted by Uim. Finally, define vi=v(xi),  di=d(xi), and Sim=s(xi,tm).The characteristic derivative is approximated in the following manner:(5)ψ(x)∂u(x,tm+1)∂τ≈ψ(x)u(x,tm+1)-u(x-,tm)[(x-x-)2+Δt2]1/2=u(x,tm+1)-u(x-,tm)Δt,
where x-=x-v(x)Δt.Then, by (5) the differential operators at the point (xi,tm+1) are discretized to be
(6)[∂u(x,t)∂t+v(x)∂u(x,t)∂x](xi,tm+1)=uim+1-u-imΔt+O(Δt+min(Δx,Δx2Δt)),
where x-i=xi-viΔt, u-im=u(x-i,tm) is calculated by the piecewise-linear interpolation.Letri=viΔt/Δx, and let symbol [ri] be the integral part of ri. Let i0=i-[ri]-1 and ri*=ri-[ri], using linear interpolation to u-im, we get [23]
(7)u-im=(1-ri*)ui0+1m+ri*ui0m.
Generally, the quantity ri 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]:
(8)∂αf(x)∂xα=limM→∞1Δxα∑k=0M‍gk(α)f(x-(k-1)Δx),
where M is a positive integer, Δx=(x-L)/M, and gk(α)=(-1)k(αk), which can be calculated by
(9)g0(α)=1,gk(α)=(1-α+1k)gk-1(α).Following (5), we can discretize the Riemann-Liouville operator [12]
(10)∂αf(x)∂xα=1Δxα∑k=0[l]+1‍ωk(α)f(x-(k-1)Δx)+O(Δxp),
where p is a positive integer, l=(x-L)/Δx, and [l] is the integral part of l.Here, according to the definition ofωk(α), we can draw the following conclusion [9, 12]: when p=1, ωk(α)=gk(α); when p=2, ωk(α)=(3α/2α)∑j=0k‍gj(α)gk-j(α)3-j.According to the shifted Grünwald-Letnikov formulae (8) and (10), the definition (2) can be discretized by
(11)∂αu∂xα(xi,tm+1)=1Δxα∑k=0i+1‍gk(α)ui-k+1m+1+O(Δx),
with p=1.Inserting (6), (7), and (11) into (1), neglecting the truncation error, after rearranging the terms, we finally get the fractional characteristic difference scheme
(12)Uim+1-ξi∑k=0i+1‍gk(α)Ui-k+1m+1=(1-ri*)Ui0+1m+ri*Ui0m+ΔtSim+1,
where ξi=diΔt/Δxα is associated with the diffusion coefficient, and the initial values are Ui0=u0(xi).


## 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 ofgk(α) in (9), when 1<α≤2, the coefficients gk(α) have the following properties:(1)
g0(α)=1, g1(α)=-α<0 and gk(α)≥0, k≥2,(2)
∑k=0∞‍gk(α)=0 and ∑k=0,k≠1∞‍gk(α)=-g1(α)=α,(3)
∑k=0m‍gk(α)≤0, m≥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 AUm+1=BUm+Sm+1, where Um=[U1m,U2m,…,UN-1m]T, Sm+1=[ΔtS1m+1,ΔtS2m+1,…,ΔtSN-1m+1]T, A=[ai,j] and B=[bi,j] are two (N-1)×(N-1) matrices.
Wheni=1,2,…,N-1 and j=1,…,N-1, we have
(13)ai,j={-ξigi-j+1(α),j≤i-1,1-ξig1(α),j=i,-ξig0(α),j=i+1,0,j>i+1,bi,j={1-ri*,j=i-[ri],ri*,j=i-[ri]-1,0,otherwise,
where 0≤ri*≤1, so 0≤1-ri*≤1, too. Thus, B≥0 holds.
For matrixB, we can easily get
(14)∑j=1N-1bi,j={0,i<[ri],1-ri*,i=[ri],1,i>[ri].
Therefore,
(15)∥B∥∞=max1≤i≤N-1∑j=1N-1‍|bi,j|=max1≤i≤N-1∑j=1N-1‍bi,j≤1.
For matrixA, we can easily see that ai,j≤0 for all i≠j, and
(16)ai,i=1-ξig1(α)=1+ξiα>0.
At the same time, according to Lemma 1, we have
(17)∑j=1,j≠iN-1‍|ai,j|=ξi∑k=0,k≠1i‍gk(α)≤ξiα,
for i=1,2,…,N-1, we can obtain
(18)|ai,i|-∑j=1,j≠iN-1‍|ai,j|≥1.
Therefore, A is diagonally dominant by rows. According to the matrix theory [24, 25],
(19)∥A-1∥∞≤1min1≤i≤N-1{|ai,i|-∑j≠i|ai,j|}≤1.
Then, from (15) and (19), we can get
(20)∥A-1B∥∞≤∥A-1∥∞∥B∥∞≤1.
Thus, the characteristic finite difference scheme (12) is unconditionally stable.In addition, taking into account (6) and (11), we have
(21)T(x,t)=O(Δt+min(Δx,Δx2Δt))=O(Δx+Δt),
so this method is consistent with a local truncation error which is O(Δx+Δt).


### 3.2. Convergence Analysis

Letu(xi,tm)(i=1,2,…,N-1; m=0,1,2,…,M) be the exact solution of (1) at mesh point (xi,tm). Define eim=u(xi,tm)-Uim, i=1,2,…,N-1; m=0,1,2,…,M and Em=[e1m,e2m,…,eN-1m]T. Using E0=0, substitution into (12) leads to
(22)AEm+1=BEm+Rm+1,
where Rm=[R1m,R2m,…,RN-1m]T, and Rim=ΔtT(xi,tm), i=1,2,…,N-1; m=0,1,2,…,M.From (6), (11), and (21), we have the following.Lemma 3.
There is a positive constantC, such that
(23)|Rim|≤CΔt(Δx+Δt),
where i=1,2,…,N-1; m=0,1,2,…,M.Let∥Em∥∞=max1≤i≤N-1|eim|, because mΔt≤T is finite, we can obtain the following result.Theorem 4.
LetUim be the numerical solution computed by the CFDM (12). Then there is a positive constant C1, such that
(24)|Uim-u(xi,tm)|≤C1(Δx+Δt),i=1,2,…,N-1;m=0,1,2,…,M.Proof.
From (22) and E0=0, we can get
(25)Em+1=A-1(BEm+Rm+1)=A-1[BA-1(BEm-1+Rm)+Rm+1]=(A-1B)mA-1R1+(A-1B)m-1A-1R2+⋯+A-1Rm+1.
According to (19) and Lemma 3, using (25), we finally get
(26)∥Em+1∥∞≤∥(A-1B)m∥∞∥A-1∥∞∥R1∥∞+∥(A-1B)m-1∥∞∥A-1∥∞∥R2∥∞+⋯+∥A-1∥∞∥Rm+1∥∞≤CΔt(m+1)(Δx+Δt).
BecausemΔt≤T, from (26), we obtain the error estimate ∥Em∥∞≤CT(Δx+Δt). Therefore, we get the result |Uim-u(xi,tm)|≤C1(Δx+Δt), where C1=CT is a positive constant.Finally, according to Theorem4, we draw the following conclusion: when Δt→0 and Δx→0, |eim|→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.


## 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 AUm+1=BUm+Sm+1, where Um=[U1m,U2m,…,UN-1m]T, Sm+1=[ΔtS1m+1,ΔtS2m+1,…,ΔtSN-1m+1]T, A=[ai,j] and B=[bi,j] are two (N-1)×(N-1) matrices.
Wheni=1,2,…,N-1 and j=1,…,N-1, we have
(13)ai,j={-ξigi-j+1(α),j≤i-1,1-ξig1(α),j=i,-ξig0(α),j=i+1,0,j>i+1,bi,j={1-ri*,j=i-[ri],ri*,j=i-[ri]-1,0,otherwise,
where 0≤ri*≤1, so 0≤1-ri*≤1, too. Thus, B≥0 holds.
For matrixB, we can easily get
(14)∑j=1N-1bi,j={0,i<[ri],1-ri*,i=[ri],1,i>[ri].
Therefore,
(15)∥B∥∞=max1≤i≤N-1∑j=1N-1‍|bi,j|=max1≤i≤N-1∑j=1N-1‍bi,j≤1.
For matrixA, we can easily see that ai,j≤0 for all i≠j, and
(16)ai,i=1-ξig1(α)=1+ξiα>0.
At the same time, according to Lemma 1, we have
(17)∑j=1,j≠iN-1‍|ai,j|=ξi∑k=0,k≠1i‍gk(α)≤ξiα,
for i=1,2,…,N-1, we can obtain
(18)|ai,i|-∑j=1,j≠iN-1‍|ai,j|≥1.
Therefore, A is diagonally dominant by rows. According to the matrix theory [24, 25],
(19)∥A-1∥∞≤1min1≤i≤N-1{|ai,i|-∑j≠i|ai,j|}≤1.
Then, from (15) and (19), we can get
(20)∥A-1B∥∞≤∥A-1∥∞∥B∥∞≤1.
Thus, the characteristic finite difference scheme (12) is unconditionally stable.In addition, taking into account (6) and (11), we have
(21)T(x,t)=O(Δt+min(Δx,Δx2Δt))=O(Δx+Δt),
so this method is consistent with a local truncation error which is O(Δx+Δt).


## 3.2. Convergence Analysis

Letu(xi,tm)(i=1,2,…,N-1; m=0,1,2,…,M) be the exact solution of (1) at mesh point (xi,tm). Define eim=u(xi,tm)-Uim, i=1,2,…,N-1; m=0,1,2,…,M and Em=[e1m,e2m,…,eN-1m]T. Using E0=0, substitution into (12) leads to
(22)AEm+1=BEm+Rm+1,
where Rm=[R1m,R2m,…,RN-1m]T, and Rim=ΔtT(xi,tm), i=1,2,…,N-1; m=0,1,2,…,M.From (6), (11), and (21), we have the following.Lemma 3.
There is a positive constantC, such that
(23)|Rim|≤CΔt(Δx+Δt),
where i=1,2,…,N-1; m=0,1,2,…,M.Let∥Em∥∞=max1≤i≤N-1|eim|, because mΔt≤T is finite, we can obtain the following result.Theorem 4.
LetUim be the numerical solution computed by the CFDM (12). Then there is a positive constant C1, such that
(24)|Uim-u(xi,tm)|≤C1(Δx+Δt),i=1,2,…,N-1;m=0,1,2,…,M.Proof.
From (22) and E0=0, we can get
(25)Em+1=A-1(BEm+Rm+1)=A-1[BA-1(BEm-1+Rm)+Rm+1]=(A-1B)mA-1R1+(A-1B)m-1A-1R2+⋯+A-1Rm+1.
According to (19) and Lemma 3, using (25), we finally get
(26)∥Em+1∥∞≤∥(A-1B)m∥∞∥A-1∥∞∥R1∥∞+∥(A-1B)m-1∥∞∥A-1∥∞∥R2∥∞+⋯+∥A-1∥∞∥Rm+1∥∞≤CΔt(m+1)(Δx+Δt).
BecausemΔt≤T, from (26), we obtain the error estimate ∥Em∥∞≤CT(Δx+Δt). Therefore, we get the result |Uim-u(xi,tm)|≤C1(Δx+Δt), where C1=CT is a positive constant.Finally, according to Theorem4, we draw the following conclusion: when Δt→0 and Δx→0, |eim|→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 v(x)=v and d(x)=d are two constant functions, the function s(x,t)=0, and the initial condition u0(x)=δ(x) 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]:
(27)u(x,t)=1π∫0+∞‍cos[(x-vt)ξ+dtsin(πα2)ξα]×edtcos(πα/2)ξαdξ.In the numerical experiments, the data are chosen as follows:v=2, T=1.5, L=-2, and R=4. For simplicity, we assume that the initial time t0=0.1, and the initial values U0=u(x,t0) are calculated from (27). The numerical solutions are obtained from the CFDM scheme (12) discussed above. The errors of numerical solutions are measured in L∞ and L2 norms which are defined as follows:
(28)E∞,Δt=maxi|U(x(i),t)-u(x(i),t)|,E2,Δt=h∑i=0N‍(U(x(i),t)-u(x(i),t))2,
where U(x(i),t) and u(x(i),t) are the numerical solution and analytical solution at (x(i),t), respectively. h is the space step, and Δt is the time step. The ratios of convergence in time are calculated by
(29)log(El,Δt1El,Δt2)[log(Δt1Δt2)]-1,l=2,∞,
where Δt1 and Δt2 are two time steps, and very small spatial step size Δx 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 d=0.1, the space step Δx=1/90, and the time step Δt=1/30 at t=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 withd=0.1, Δx=1/90, and Δt=1/30 at t=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. Table1 shows L∞ and L2 errors and ratios in space with the diffusion coefficient d=0.1 and the small time step Δt=1/400 at time t=1.5, for two different fractional orders α=1.9 and 1.6, respectively. The second and fourth columns show the L∞ and L2 errors, respectively. Column three shows the ratio of the CFDM in L∞ norm, and column five shows the ratio of the CFDM in L2 norm. It is clearly shown that the CFDM method is of first-order accuracy in space in both L∞ and L2 norms.Table 1
The errors and ratios in space calculated by CFDM (12) with d=0.1 and Δt=1/400 at time t=1.5, for two different fractional orders α=1.9 and 1.6, respectively.
Δ
x
∥
E
∥
∞-CFDM
Ratio
∥
E
∥
2-CFDM
Ratio
α
=
1.9
1/25
8.9731
E
-
2
—
7.4998
E
-
2
—
1/30
7.5702
E
-
2
0.9325
6.2978
E
-
2
0.9581
1/35
6.5100
E
-
2
0.9788
5.3975
E
-
2
1.0001
1/40
5.6805
E
-
2
1.0207
4.6976
E
-
2
1.0401
1/45
5.0149
E
-
2
1.0581
4.1379
E
-
2
1.0771
α
=
1.6
1/60
1.0621
E
-
1
—
7.1299
E
-
2
—
1/65
9.7561
E
-
2
1.0612
6.5310
E
-
2
1.0961
1/70
8.9922
E
-
2
1.1002
6.0045
E
-
2
1.1342
1/75
8.3116
E
-
2
1.1408
5.5381
E
-
2
1.1720
1/80
7.7015
E
-
2
1.1813
5.1218
E
-
2
1.2108Meanwhile, Table2 shows L∞ and L2 errors and ratios in time with d=0.02 and Δx=1/300 at time t=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 L∞ and L2 norms, too. Furthermore, from the table, we can see that even with the very large time steps size of Δt=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 O(Δx+Δt), which confirms the result of Theorem 4.Table 2
The errors and ratios in time calculated by CFDM (12) with d=0.02 and Δx=1/300 at time t=1.5, for two different fractional orders α=1.9 and 1.6, respectively.
Δ
t
∥
E
∥
∞-CFDM
Ratio
∥
E
∥
2-CFDM
Ratio
α
=
1.9
1/5
9.4691
E
-
2
—
4.7080
E
-
2
—
1/10
4.5550
E
-
2
1.0558
2.3455
E
-
2
1.0052
1/15
2.9714
E
-
2
1.0536
1.5526
E
-
2
1.0175
1/20
2.1900
E
-
2
1.0607
1.1555
E
-
2
1.0268
1/30
1.4163
E
-
2
1.0749
7.5847
E
-
3
1.0383
α
=
1.6
1/10
1.3447
E
-
1
—
5.1437
E
-
2
—
1/15
9.6286
E
-
2
0.8238
3.6205
E
-
2
0.8661
1/18
8.4253
E
-
2
0.7322
3.1293
E
-
2
0.7997
1/20
7.6641
E
-
2
0.8987
2.8581
E
-
2
0.8604
1/25
6.4439
E
-
2
0.7771
2.4270
E
-
2
0.7327


### 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. Table3 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 L2 errors calculated by the CFDM (12) and the implicit SFDM for two different diffusion coefficients d=0.1 and d=0.02, with some different space steps Δx and time steps Δt for α=1.9 at time t=1.5. From the table, we can easily see that the L2 errors generated by the CFDM is obviously much smaller than the L2 errors generated by the implicit SFDM, which has one- even two-order higher accuracy than the latter.Table 3
Comparison ofL2 errors calculated by the CFDM (12) and the implicit SFDM with α=1.9, different space steps Δx and time steps Δt at time t=1.5, for two different diffusion coefficients d=0.1 and d=0.02.
d
Δ
x
Δ
t
∥
E
∥
2-CFDM
∥
E
∥
2-im-SFDM
0.1
1/60
1/30
4.0456
E
-
3
1.1639
E
-
1
1/50
1/40
3.3588
E
-
3
1.0382
E
-
1
1/80
1/50
2.1298
E
-
3
8.0739
E
-
2
0.02
1/80
1/20
7.9441
E
-
3
5.6804
E
-
1
1/90
1/30
6.4445
E
-
3
4.8836
E
-
1
1/100
1/45
5.3891
E
-
3
4.1255
E
-
1Figure 2
TheL∞ errors between exact solutions and numerical solutions calculated by the CFDM (12) and the implicit SFDM at different time levels with α=1.9, Δx=1/50, and Δt=1/40 for d=0.1.Moreover, Figure2 shows that the L∞ errors between exact solutions and numerical solutions calculated by the CFDM and the SFDM at different time levels with α=1.9, the diffusion coefficient d=0.1, the space step Δx=1/50, and the time step Δt=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 d=0.02 at time t=1.5 in Figure 3. To keep similar L∞ error, for the CFDM, the grid sizes are chosen as Δt=1/30 and Δx=1/100, and it generates the L∞ error of 1.2168×10-2. However, for the implicit SFDM, we have to choose very small grid sizes, which are at least Δt=1/5000 and Δx=1/1000. In this case, the implicit SFDM generates the L∞ error of 4.6630×10-2, 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 att=1.5 with d=0.02. The solid line corresponds to the exact solution and the dotted line corresponds to numerical solution of the CFDM with α=1.9, Δx=1/100, and Δt=1/30. The stared line corresponds to numerical solution of the implicit SFDM with Δx=1/1000 and Δt=1/5000.In Section4 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. Table4 shows that the L2 errors of three different explicit schemes with different Δx and Δt for d=0.02 and α=1.9 at time t=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
TheL2 errors calculated by the upwind, central and Lax-Wendroff explicit difference schemes with different space steps Δx and time steps Δt for d=0.02 and α=1.9 at time t=1.5.
Δ
x
Δ
t
∥
E
∥
2-Upwind
∥
E
∥
2-Central
∥
E
∥
2-Lax-Wendroff
1/50
1/1000
2.5462
E
-
1
3.3689
E
-
2
1.1582
E
-
2
1/80
1/1500
1.7972
E
-
1
1.9049
E
-
2
2.5064
E
-
3
1/100
1/2000
1.5122
E
-
1
1.4560
E
-
2
1.0258
E
-
3
1/120
1/3000
1.3178
E
-
1
9.2617
E
-
3
5.9643
E
-
4


## 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 d=0.1, the space step Δx=1/90, and the time step Δt=1/30 at t=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 withd=0.1, Δx=1/90, and Δt=1/30 at t=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. Table1 shows L∞ and L2 errors and ratios in space with the diffusion coefficient d=0.1 and the small time step Δt=1/400 at time t=1.5, for two different fractional orders α=1.9 and 1.6, respectively. The second and fourth columns show the L∞ and L2 errors, respectively. Column three shows the ratio of the CFDM in L∞ norm, and column five shows the ratio of the CFDM in L2 norm. It is clearly shown that the CFDM method is of first-order accuracy in space in both L∞ and L2 norms.Table 1
The errors and ratios in space calculated by CFDM (12) with d=0.1 and Δt=1/400 at time t=1.5, for two different fractional orders α=1.9 and 1.6, respectively.
Δ
x
∥
E
∥
∞-CFDM
Ratio
∥
E
∥
2-CFDM
Ratio
α
=
1.9
1/25
8.9731
E
-
2
—
7.4998
E
-
2
—
1/30
7.5702
E
-
2
0.9325
6.2978
E
-
2
0.9581
1/35
6.5100
E
-
2
0.9788
5.3975
E
-
2
1.0001
1/40
5.6805
E
-
2
1.0207
4.6976
E
-
2
1.0401
1/45
5.0149
E
-
2
1.0581
4.1379
E
-
2
1.0771
α
=
1.6
1/60
1.0621
E
-
1
—
7.1299
E
-
2
—
1/65
9.7561
E
-
2
1.0612
6.5310
E
-
2
1.0961
1/70
8.9922
E
-
2
1.1002
6.0045
E
-
2
1.1342
1/75
8.3116
E
-
2
1.1408
5.5381
E
-
2
1.1720
1/80
7.7015
E
-
2
1.1813
5.1218
E
-
2
1.2108Meanwhile, Table2 shows L∞ and L2 errors and ratios in time with d=0.02 and Δx=1/300 at time t=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 L∞ and L2 norms, too. Furthermore, from the table, we can see that even with the very large time steps size of Δt=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 O(Δx+Δt), which confirms the result of Theorem 4.Table 2
The errors and ratios in time calculated by CFDM (12) with d=0.02 and Δx=1/300 at time t=1.5, for two different fractional orders α=1.9 and 1.6, respectively.
Δ
t
∥
E
∥
∞-CFDM
Ratio
∥
E
∥
2-CFDM
Ratio
α
=
1.9
1/5
9.4691
E
-
2
—
4.7080
E
-
2
—
1/10
4.5550
E
-
2
1.0558
2.3455
E
-
2
1.0052
1/15
2.9714
E
-
2
1.0536
1.5526
E
-
2
1.0175
1/20
2.1900
E
-
2
1.0607
1.1555
E
-
2
1.0268
1/30
1.4163
E
-
2
1.0749
7.5847
E
-
3
1.0383
α
=
1.6
1/10
1.3447
E
-
1
—
5.1437
E
-
2
—
1/15
9.6286
E
-
2
0.8238
3.6205
E
-
2
0.8661
1/18
8.4253
E
-
2
0.7322
3.1293
E
-
2
0.7997
1/20
7.6641
E
-
2
0.8987
2.8581
E
-
2
0.8604
1/25
6.4439
E
-
2
0.7771
2.4270
E
-
2
0.7327


## 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. Table3 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 L2 errors calculated by the CFDM (12) and the implicit SFDM for two different diffusion coefficients d=0.1 and d=0.02, with some different space steps Δx and time steps Δt for α=1.9 at time t=1.5. From the table, we can easily see that the L2 errors generated by the CFDM is obviously much smaller than the L2 errors generated by the implicit SFDM, which has one- even two-order higher accuracy than the latter.Table 3
Comparison ofL2 errors calculated by the CFDM (12) and the implicit SFDM with α=1.9, different space steps Δx and time steps Δt at time t=1.5, for two different diffusion coefficients d=0.1 and d=0.02.
d
Δ
x
Δ
t
∥
E
∥
2-CFDM
∥
E
∥
2-im-SFDM
0.1
1/60
1/30
4.0456
E
-
3
1.1639
E
-
1
1/50
1/40
3.3588
E
-
3
1.0382
E
-
1
1/80
1/50
2.1298
E
-
3
8.0739
E
-
2
0.02
1/80
1/20
7.9441
E
-
3
5.6804
E
-
1
1/90
1/30
6.4445
E
-
3
4.8836
E
-
1
1/100
1/45
5.3891
E
-
3
4.1255
E
-
1Figure 2
TheL∞ errors between exact solutions and numerical solutions calculated by the CFDM (12) and the implicit SFDM at different time levels with α=1.9, Δx=1/50, and Δt=1/40 for d=0.1.Moreover, Figure2 shows that the L∞ errors between exact solutions and numerical solutions calculated by the CFDM and the SFDM at different time levels with α=1.9, the diffusion coefficient d=0.1, the space step Δx=1/50, and the time step Δt=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 d=0.02 at time t=1.5 in Figure 3. To keep similar L∞ error, for the CFDM, the grid sizes are chosen as Δt=1/30 and Δx=1/100, and it generates the L∞ error of 1.2168×10-2. However, for the implicit SFDM, we have to choose very small grid sizes, which are at least Δt=1/5000 and Δx=1/1000. In this case, the implicit SFDM generates the L∞ error of 4.6630×10-2, 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 att=1.5 with d=0.02. The solid line corresponds to the exact solution and the dotted line corresponds to numerical solution of the CFDM with α=1.9, Δx=1/100, and Δt=1/30. The stared line corresponds to numerical solution of the implicit SFDM with Δx=1/1000 and Δt=1/5000.In Section4 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. Table4 shows that the L2 errors of three different explicit schemes with different Δx and Δt for d=0.02 and α=1.9 at time t=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
TheL2 errors calculated by the upwind, central and Lax-Wendroff explicit difference schemes with different space steps Δx and time steps Δt for d=0.02 and α=1.9 at time t=1.5.
Δ
x
Δ
t
∥
E
∥
2-Upwind
∥
E
∥
2-Central
∥
E
∥
2-Lax-Wendroff
1/50
1/1000
2.5462
E
-
1
3.3689
E
-
2
1.1582
E
-
2
1/80
1/1500
1.7972
E
-
1
1.9049
E
-
2
2.5064
E
-
3
1/100
1/2000
1.5122
E
-
1
1.4560
E
-
2
1.0258
E
-
3
1/120
1/3000
1.3178
E
-
1
9.2617
E
-
3
5.9643
E
-
4


## 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.



---

*Source: 648595-2013-03-24.xml*