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

Authors: Lijuan Su; Pei Cheng
Journal: Advances in Mathematical Physics
 (2013)

---

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


---

## Body

## 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]:
(1)∂u(x,t)∂t=-V(x,t)∂u(x,t)∂x+D+(x,t)∂αu(x,t)∂+xα+D-(x,t)∂α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)=u(R,t)=0,0≤t≤T,
where V(x,t)>0 is the drift of the process, that is, the mean advective velocity, α is the order of fractional differentiation, D+(x,t)=(1+β)D(x,t)/2, D-(x,t)=(1-β)D(x,t)/2, D(x,t)>0 is the coefficient of dispersion, and -1≤β≤1 indicates the relative weight of forward versus backward transition probability. The function u0(x) is the initial condition, the boundary conditions are zero Dirichlet boundary conditions, and the function s(x,t) is a source/sink term. The ∂αu(x,t)/∂+xα and ∂αu(x,t)/∂-xα 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 f(x) for x∈[L,R] defined by [12]:
(2)∂αf(x)∂+xα=1Γ(n-α)dndxn∫Lx‍f(ξ)dξ(x-ξ)α-n+1,∂αf(x)∂-xα=(-1)nΓ(n-α)dndxn∫xR‍f(ξ)dξ(ξ-x)α-n+1,
where n-1<α≤n (n 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 Section2, 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:Δt is the time step, Δx is the spatial step, the coordinates of the mesh points are xj=L+jΔx, j=0,1,2,…,N, N=(R-L)/Δx, and tm=mΔt, m=0,1,2,…,M, M=T/Δt, and the values of the solution u(x,t) at these grid points are u(xj,tm)≡ujm≃Ujm, where we denote by Ujm the numerical estimate of the exact value of u(x,t) at the point (xj,tm). Define Vjm+1/2=V(xj,tm+1/2), Djm+1/2=D(xj,tm+1/2) and sjm+1/2=s(xj,tm+1/2).The centered time difference scheme is [20]
(3)∂u∂t|(xj,tm+1/2)=ujm+1-ujmΔt+O(Δt)2,
and the backward space difference scheme is
(4)∂u∂x|(xj,tm)=ujm-uj-1mΔx+O(Δx).According to the shifted Grünwald-Letnikov definition [8], the definition (2) can be written as
(5)∂αu(x)∂+xα|(xj,tm)=1hα∑k=0j+1‍gk(α)uj-k+1m+O(h),∂αu(x)∂-xα|(xj,tm)=1hα∑k=0N-j+1‍gk(α)uj+k-1m+O(h).
Here, gk(α)=(-1)k(αk) can be evaluated recursively:
(6)g0(α)=1,gk(α)=(1-α+1k)gk-1(α).In the weighted average method, (1) can be evaluated at the intermediate point of the grid (xj,tm+1/2) by the following formula:
(7)∂u(x,t)∂t|(xj,tm+1/2)=-V(xj,tm+1/2)[λ∂u(x,t)∂x|(xj,tm)-V(xj,tm+1/2)+(1-λ)∂u(x,t)∂x|(xj,tm+1)]+D+(xj,tm+1/2)[λ∂αu(x,t)∂+xα|(xj,tm)+D+(xj,tm+1/2)+(1-λ)∂αu(x,t)∂+xα|(xj,tm+1)]+D-(xj,tm+1/2)[λ∂αu(x,t)∂+xα|(xj,tm)+D-(xj,tm+1/2)+(1-λ)∂αu(x,t)∂+xα|(xj,tm+1)]+s(x,t)|(xj,tm+1/2),
where 0≤λ≤1 is the weighting coefficient.Applying (3)~(5) to (7), letting h=Δx, and neglecting the truncation error, we get the FWA difference scheme
(8)Ujm+1+(1-λ)[∑k=0N-j+1‍gk(α)Uj+k-1m+1-rjm+1/2Uj-1m+1+rjm+1/2Ujm+1Ujm+1+(1-λ)-ξjm+1/2∑k=0j+1‍gk(α)Uj-k+1m+1Ujm+1+(1-λ)-ηjm+1/2∑k=0N-j+1‍gk(α)Uj+k-1m+1]=Ujm-λ[∑k=0N-j+1‍gk(α)Uj+k-1m-rjm+1/2Uj-1m+rjm+1/2Ujm=Ujm-λ-ξjm+1/2∑k=0j+1‍gk(α)Uj-k+1m=Ujm-λ-ηjm+1/2∑k=0N-j+1‍gk(α)Uj+k-1m]+Δtsjm+1/2,j=1,2,…,N-1,m=0,1,2,…,M-1,
where rjm+1/2=Vjm+1/2Δt/h, ξjm+1/2=D+,jm+1/2Δt/hα, ηjm+1/2=D-,jm+1/2Δt/hα, and the initial values are calculated by Uj(0)=u0(xj), j=1,2,…,N-1. Generally, the quantity rjm+1/2 is called the Courant (or CFL) number, the ξjm+1/2 and ηjm+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 coefficientsgk(α) given in (6) with 1<α≤2 satisfy the following properties:
(9)g0(α)=1,g1(α)=-α<0,1≥g2(α)≥g3(α)≥⋯≥0,∑k=0∞‍gk(α)=0,∑k=0m‍gk(α)≤0(m≥1).Theorem 2.
When0≤λ≤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
(10)αΔtDmaxhα+ΔtVmaxh≤12λ-1,
where Vmax=maxL≤x≤R,0≤t≤TV(x,t) and Dmax=maxL≤x≤R,0≤t≤TD(x,t).Proof.
The FWA scheme (8) can be rewritten as [I+(1-λ)A]Um+1=(I-λA)Um, m=0,1,2,…,M-1; here, Um=[U0m,U1m,U2m,…,UNm]T, A=(aij), i,j=0,1,2,…,N.
The matrix entriesai,j for i=1,2,…,N-1 and j=0,1,…,N are defined by
(11)ai,j={rim+1/2-(ξim+1/2+ηim+1/2)g1(α),j=i-rim+1/2-ξim+1/2g2(α)-ηim+1/2g0(α),j=i-1-ξim+1/2g0(α)-ηim+1/2g2(α),j=i+1-ξim+1/2gi-j+1(α),j<i-1-ηim+1/2gj-i+1(α),j>i+1
while a0,j=aN,j=0, for j=0,1,…,N.
According to Lemma1 and the Gerschgorin theorem, the eigenvalues of the matrix A (noted ωi) are in the disks centered at ai,i=rim+1/2-(ξim+1/2+ηim+1/2)g1(α), with radius
(12)Ri=∑j=0,j≠iN‍|ai,j|=rim+1/2+ξim+1/2∑k=0,k≠1i+1‍gk(α)+ηim+1/2∑k=0,k≠1N-i+1‍gk(α)≤rim+1/2-(ξim+1/2+ηim+1/2)g1(α).
Therefore, we have
(13)0≤ωi≤2[rim+1/2-(ξim+1/2+ηim+1/2)g1(α)]=2[Dim+1/2αΔthα+Vim+1/2Δth].
Next, note thatωi is an eigenvalue of A if and only if (1-λωi)/(1+(1-λ)ωi) is an eigenvalue of the matrix (I-λA)[I+(1-λ)A]-1. Because of 0≤ωi and 0≤λ≤1, we get
(14)1-λωi1+(1-λ)ωi=1-ωi1+(1-λ)ωi≤1.
In addition,(1-λωi)/(1+(1-λ)ωi)≥-1 as long as (2λ-1)ωi≤2.
Hence, when0≤λ≤1/2, we can find that -1≤(1-λωi)/(1+(1-λ)ωi)≤1 always holds; that is, |(1-λωi)/(1+(1-λ)ωi)|≤1. Then, the FWA (8) is unconditionally stable. On the other hand, when 1/2<λ≤1, from (13) and ωi≤2/(2λ-1), we get the limited condition (αΔtDmax/hα)+(ΔtVmax/h)≤1/(2λ-1) where Vmax=maxL≤x≤R,0≤t≤TV(x,t) and Dmax=maxL≤x≤R,0≤t≤TD(x,t). Therefore, the FWA (8) is conditionally stable.Let(15)S=αΔtDmaxhα+ΔtVmaxh,
the stability limit S× is S×=1/(2λ-1).In addition, taking into account (3)~(5), for arbitrary Δx and Δt, we derive that this method is consistent with a local truncation error O(Δx+Δt), except for the FCN method, whose accuracy is of (Δt)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.


## 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, V(x,t)=V, D(x,t)=D, and s(x,t)=0; the initial condition is
(16)u0(x)=20π∫0+∞‍cos[(x-0.1V)ξ]e0.1Dcos(πα/2)ξαdξ.In this case, the analytical solution of (1) solved by the Fourier transform methods is [12]
(17)u(x,t)=20π∫0+∞‍cos[(x-V(t+0.1))ξ]eD(t+0.1)cos(πα/2)ξαdξ.In the following numerical experiments, the data are chosen as follow:α=1.9, D=2, V=2, T=2.5, L=-5, and R=15.The numerical solutions are obtained from the FWA scheme (8) discussed above, with different λ, Δt, t, S, and h. From (15), the values of Δt for fixed S and Δx=h are
(18)Δt=S((αDmax/hα)+(Vmax/h))=S((αD/hα)+(V/h)).The computational results are shown in Figures1, 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, Δx=1/40, and small S=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 Δt=2.8×10-4 calculated from (18). Figure 2 has the same assumptions as Figure 1 but for S=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 S×=1/(2λ-1)=1.25.Figure 1
Numerical solutions of (1) by means of the FWA method for λ=0.9, Δx=1/40, and S=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 Figure1 but for S=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 Δx=1/40 and S=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 Δx=1/40 and large S=100 after 10, 50, and 100 time steps. Meanwhile, we can gain the large time step Δt=2.3×10-2 calculated from (18), which is much larger than Δt=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.


## 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 O(Δx+Δt), except for the FCN method, whose accuracy with respect to the time step is of (Δt)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.



---

*Source: 129404-2013-06-20.xml*