# A Second-Order Uniformly Stable Explicit Asymmetric Discretization Method for One-Dimensional Fractional Diffusion Equations

Authors: Lin Zhu
Journal: Complexity
 (2019)

---

## Abstract

Using the asymmetric discretization technique, an explicit finite difference scheme is constructed for one-dimensional spatial fractional diffusion equations (FDEs). The spatial fractional derivative is approximated by the weighted and shifted Grünwald difference operator. The scheme can be solved explicitly by calculating unknowns in the different nodal-point sequences at the odd time-step and the even time-step. The uniform stability is proven and the error between the discrete solution and analytical solution is theoretically estimated. Numerical examples are given to verify theoretical analysis.


---

## Body

## 1. Introduction

Fractional operator has been used to solve a certain class of problems arising from the real world [1–7]. Fractional spatial derivatives are suitable to describe anomalous dispersion and diffusion models. In particular, when the second-order derivative is replaced by the Riemann-Liouville fractional derivative of order α∈(1,2) in space, the fractional equation represents a superdiffusive flow model.Up to now, several finite difference methods have been presented for the initial-boundary value problem of the one-dimensional spatial fractional diffusion equations (FDEs) [8–15]. Because of the nonlocal property of fractional differential operators, many numerical methods for FDEs have full coefficient matrices which require storage of O(K2) and computational cost of O(K3) for a problem of size K, where K is the number of discrete points in the spatial coordinate. How to save storage space and reduce the number of calculations is essential. It is well known that implicit scheme is more stable, but more complex in calculation than the explicit scheme.The main purpose of this paper is to construct a semi-implicit finite difference scheme for solving one-dimensional spatial FDE by combining a second-order weighted and shifted Grünwald difference (WSGD) operator [16, 17] with the asymmetric discretization technique introduced by Saulyev [18]. Formally, the obtained finite difference scheme is implicit; however we can obtain the discrete solution of the scheme explicitly by using different nodal-point stencils at odd time-step and even time-step. Instead of solving the linear algebraic system, we compute the unknowns according to the sequences in left-to-right order at the odd time-step and then compute the unknowns according to the sequences in right-to-left order at the even time-step. We prove that the scheme is uniformly stable when α∈[(17-1)/2,2], and the error between the numerical and analytical solutions in discrete l2 norm is of order O(Δt2h-2(α-1)+Δt+h2), where α is the order of spatial fractional derivative, and h and Δt are the space and time mesh sizes. The error estimate here is better than that in [19], where the error in discrete l2 norm is of order O(Δt2h-2(α-1)+Δt+h). Some examples are presented to show that the numerical results match our theoretical analysis.This asymmetric technique has been used to construct parallel algorithms in several research results [20–23]. Within the scope of our knowledge, there are only two papers including this paper and the one Rui and Huang presented [19] to analyze the error for the asymmetric technique and get the error estimates between the discrete and analytical solutions. Former papers presented the stability property, mentioned that the truncation error is of order OΔth-1+Δt+h for parabolic problems [20, 21], or presented the asymmetric technique [22, 23] in numerical calculations.The rest of the paper is organized as follows. In Section2 the semi-implicit finite difference scheme is presented, and in Section 3 it is proven that the scheme proposed is unconditionally stable for the problem with the order of fractional derivative belonging to [(17-1)/2,2]. In Section 4 the error estimation is analyzed, and in Section 5 numerical experiments are carried out to verify our theoretical analysis. Finally some concluding remarks are given in the last section.Throughout the paper, we useC to denote a generic constant, which may have different values in different context.


## 2. The Semi-Implicit Finite Difference Method for FDES

In this section, combining a second-order WSGD operator for the Riemann-Liouville fractional derivatives with the asymmetric discretization technique, the finite difference scheme is constructed for FDEs.Consider the following one-dimensional space-fractional diffusion equation:(1)∂ux,t∂t=cx,t∂αux,t∂xα+sx,t,x∈L,R,t∈0,T,(2)uL,t=0,uR,t=bRt,t∈0,T,(3)ux,0=ϕx,x∈L,R.where α is the fractional order and 1<α≤2 [8, 9].


### 2.1. Second-Order Operator for Riemann-Liouville Fractional Derivatives

Now from the shifted Grünwald finite difference formula, a second-order WSGD operator is introduced for Riemann-Liouville derivatives. First, some preliminary knowledge involving the Riemann-Liouville derivatives is given [24].Definition 1.
Theα(n-1<α<n) order left and right Riemann-Liouville fractional derivatives of the function u(x) on (a,b) are, respectively, defined as follows:
(a) The left Riemann-Liouville fractional derivative:(4)Daxαux=1Γn-αdndxn∫axuξx-ξα-n+1dξ
(b) The right Riemann-Liouville fractional derivative:(5)Dxbαux=-1nΓn-αdndxn∫xbuξξ-xα-n+1dξ,if α=n, Daxαux=dn/dxnux and Dxbαux=-1ndn/dxnux.The fractional derivative can be approximated by the standard Grüwald formula or the shifted Grüwald formula [8, 9] as follows:(6)∂αux,t∂xα=limh→0⁡1hα∑k=0Mgkux-kh,t,(7)∂αux,t∂xα=limh→0⁡1hα∑k=0Mgkux-k-1h,t,where h=(x-L)/M and the coefficients are defined as(8)gk=Γk-αΓ-αΓk+1=-1kαk.Lemma 2 (see [8, 9, 24]).
The coefficients{gk} satisfy the following properties for 1<α⩽2,(9)g0=1,g1=-α<0,1≥g2≥g3≥⋯≥0;∑k=0mgk≤0,m∈Z+;∑k=0∞gk=0.LetΔt=T/(2N) and h=R-L/K be the time increment and spacial mesh size, respectively. Here K and N are positive integers. For n=0,1,…,T/2∆t, and i=0,…,K, we denote(10)tn=nΔt,xi=L+ih,uin=uxi,tn,cin=cxi,tn,sin=sxi,tn.Based on the weighted and shifted Grünwald difference (WSGD) operators, the following second-order operator is introduced to approximate the Riemann-Liouville fractional derivative [16] at the grid point (xl,tn). For l=1,…,K-1,(11)∂αuxl,tn∂xα=1hα1-α2∑k=0lgkuxl-kh,tn+α2∑k=0l+1gkuxl-k-1h,tn+Oh2.The simplified form of the approximate formula (11) for the Riemann-Liouville fractional derivative is(12)∂αuxl,tn∂xα=1hα∑k=0l+1wkuxl-k-1h,tn+Oh2,where(13)w0=α2,wk=α2gk+1-α2gk-1,k⩾1.From Lemma2 and some straightforward calculations, we obtain the following properties of the coefficient {wk} in (12) [16].Lemma 3 (see [16]).
The coefficients{wk} satisfy the following properties for α∈[(17-1)/2,2], (14)1⩾w0⩾w2⩾w3⩾w4⩾⋯⩾0,w1<0;∑k=0Mwk⩽0,M∈Z+;∑k=0∞wk=0.Combining (12) with the implicit Euler discretization, the finite difference scheme for solving equations (1)-(3) is introduced.(15)vln+1-vlnΔt=cln+1hα∑k=0l+1wkvl-k+1n+1+sln+1+OΔt+h2,where l=1,…,K-1.Now, the following semi-implicit finite difference scheme is constructed using the asymmetric technique [18, 19].For1≤l≤K-1 and 0≤n≤N=T/2Δt, one goal is to find {vl2n+2} such that(16)avl2n+1-vl2nΔt=cl2n+1hα∑k=2l+1wkvl-k+12n+1+w1+α2vl2n+1-α2vl2n+α2vl+12n+sl2n+1,bvl2n+2-vl2n+1Δt=cl2n+2hα∑k=2l+1wkvl-k+12n+1+w1+α2vl2n+1-α2vl2n+2+α2vl+12n+2+sl2n+2.The initial and boundary conditions are as follows:(17)vl0=ϕxl,l=0,1,…,K;v0n=0,vKn=bRR,tn.The scheme is an implicit scheme formally. In order to solve it explicitly, we calculate the unknowns according to the following nodal-point sequences. Here different nodal-point stencils are given at levels 2n+1 and 2n+2 (see Figures 1-2).Figure 1
The finite difference stencil for point (l,2n+1).Figure 2
The finite difference stencil for point (l,2n+2).The specific calculation processes are as follows. Let(18)pln=clnΔthα.From (16),(19)1-pl2n+1w1+α2vl2n+1=∑k=2l+1pl2n+1wkvl-k+12n+1+1-pl2n+1α2vl2n+pl2n+1α2vl+12n+Δtsl2n+1,and(20)1+pl2n+2α2vl2n+2=∑k=2l+1pl2n+2wkvl-k+12n+1+1+pl2n+2w1+α2vl2n+1+pl2n+2α2vl+12n+2+Δtsl2n+2.(1) Solve (19) according to the following sequence from time-step t2n to time-step t2n+1: (2n+1,1)→(2n+1,2)→⋯→(2n+1,l-1)→(2n+1,l)→⋯ (cf. Figure 3).Figure 3Once we findvl-k+12n+1 for k≤l-1, combining with vl2n and vl+12n, vl2n+1 can be obtained explicitly.(2) Solve (20) according to the following sequence from time t2n+1 to time t2n+2: (2n+2,K-1)→(2n+2,K-2)→⋯→(2n+2,l)→(2n+2,l-1)→⋯ (cf. Figure 4).Figure 4Similarly,vl2n+2 can be derived explicitly.In conclusion, the discrete solutions of (19) and (20) can be obtained explicitly one by one instead of solving the linear algebraic system.


## 2.1. Second-Order Operator for Riemann-Liouville Fractional Derivatives

Now from the shifted Grünwald finite difference formula, a second-order WSGD operator is introduced for Riemann-Liouville derivatives. First, some preliminary knowledge involving the Riemann-Liouville derivatives is given [24].Definition 1.
Theα(n-1<α<n) order left and right Riemann-Liouville fractional derivatives of the function u(x) on (a,b) are, respectively, defined as follows:
(a) The left Riemann-Liouville fractional derivative:(4)Daxαux=1Γn-αdndxn∫axuξx-ξα-n+1dξ
(b) The right Riemann-Liouville fractional derivative:(5)Dxbαux=-1nΓn-αdndxn∫xbuξξ-xα-n+1dξ,if α=n, Daxαux=dn/dxnux and Dxbαux=-1ndn/dxnux.The fractional derivative can be approximated by the standard Grüwald formula or the shifted Grüwald formula [8, 9] as follows:(6)∂αux,t∂xα=limh→0⁡1hα∑k=0Mgkux-kh,t,(7)∂αux,t∂xα=limh→0⁡1hα∑k=0Mgkux-k-1h,t,where h=(x-L)/M and the coefficients are defined as(8)gk=Γk-αΓ-αΓk+1=-1kαk.Lemma 2 (see [8, 9, 24]).
The coefficients{gk} satisfy the following properties for 1<α⩽2,(9)g0=1,g1=-α<0,1≥g2≥g3≥⋯≥0;∑k=0mgk≤0,m∈Z+;∑k=0∞gk=0.LetΔt=T/(2N) and h=R-L/K be the time increment and spacial mesh size, respectively. Here K and N are positive integers. For n=0,1,…,T/2∆t, and i=0,…,K, we denote(10)tn=nΔt,xi=L+ih,uin=uxi,tn,cin=cxi,tn,sin=sxi,tn.Based on the weighted and shifted Grünwald difference (WSGD) operators, the following second-order operator is introduced to approximate the Riemann-Liouville fractional derivative [16] at the grid point (xl,tn). For l=1,…,K-1,(11)∂αuxl,tn∂xα=1hα1-α2∑k=0lgkuxl-kh,tn+α2∑k=0l+1gkuxl-k-1h,tn+Oh2.The simplified form of the approximate formula (11) for the Riemann-Liouville fractional derivative is(12)∂αuxl,tn∂xα=1hα∑k=0l+1wkuxl-k-1h,tn+Oh2,where(13)w0=α2,wk=α2gk+1-α2gk-1,k⩾1.From Lemma2 and some straightforward calculations, we obtain the following properties of the coefficient {wk} in (12) [16].Lemma 3 (see [16]).
The coefficients{wk} satisfy the following properties for α∈[(17-1)/2,2], (14)1⩾w0⩾w2⩾w3⩾w4⩾⋯⩾0,w1<0;∑k=0Mwk⩽0,M∈Z+;∑k=0∞wk=0.Combining (12) with the implicit Euler discretization, the finite difference scheme for solving equations (1)-(3) is introduced.(15)vln+1-vlnΔt=cln+1hα∑k=0l+1wkvl-k+1n+1+sln+1+OΔt+h2,where l=1,…,K-1.Now, the following semi-implicit finite difference scheme is constructed using the asymmetric technique [18, 19].For1≤l≤K-1 and 0≤n≤N=T/2Δt, one goal is to find {vl2n+2} such that(16)avl2n+1-vl2nΔt=cl2n+1hα∑k=2l+1wkvl-k+12n+1+w1+α2vl2n+1-α2vl2n+α2vl+12n+sl2n+1,bvl2n+2-vl2n+1Δt=cl2n+2hα∑k=2l+1wkvl-k+12n+1+w1+α2vl2n+1-α2vl2n+2+α2vl+12n+2+sl2n+2.The initial and boundary conditions are as follows:(17)vl0=ϕxl,l=0,1,…,K;v0n=0,vKn=bRR,tn.The scheme is an implicit scheme formally. In order to solve it explicitly, we calculate the unknowns according to the following nodal-point sequences. Here different nodal-point stencils are given at levels 2n+1 and 2n+2 (see Figures 1-2).Figure 1
The finite difference stencil for point (l,2n+1).Figure 2
The finite difference stencil for point (l,2n+2).The specific calculation processes are as follows. Let(18)pln=clnΔthα.From (16),(19)1-pl2n+1w1+α2vl2n+1=∑k=2l+1pl2n+1wkvl-k+12n+1+1-pl2n+1α2vl2n+pl2n+1α2vl+12n+Δtsl2n+1,and(20)1+pl2n+2α2vl2n+2=∑k=2l+1pl2n+2wkvl-k+12n+1+1+pl2n+2w1+α2vl2n+1+pl2n+2α2vl+12n+2+Δtsl2n+2.(1) Solve (19) according to the following sequence from time-step t2n to time-step t2n+1: (2n+1,1)→(2n+1,2)→⋯→(2n+1,l-1)→(2n+1,l)→⋯ (cf. Figure 3).Figure 3Once we findvl-k+12n+1 for k≤l-1, combining with vl2n and vl+12n, vl2n+1 can be obtained explicitly.(2) Solve (20) according to the following sequence from time t2n+1 to time t2n+2: (2n+2,K-1)→(2n+2,K-2)→⋯→(2n+2,l)→(2n+2,l-1)→⋯ (cf. Figure 4).Figure 4Similarly,vl2n+2 can be derived explicitly.In conclusion, the discrete solutions of (19) and (20) can be obtained explicitly one by one instead of solving the linear algebraic system.


## 3. Stability of the Semi-Implicit Finite Difference Scheme

In this section, it is proven that the presented scheme (16) is uniformly stable as α∈[(17-1)/2,2].Suppose the coefficient of the diffusion termc(x,t)=c is a constant. Then pin=p is a constant.From (19) and (20), (21)∑k=2l+1-pwkvl-k+12n+1+1-pw1+α2vl2n+1=1-α2pvl2n+α2pvl+12n+Δtsl2n+1,1+α2pvl2n+2-α2pvl+12n+2=∑k=2l+1pwkvl-k+12n+1+1+pw1+α2vl2n+1+Δtsl2n+2.Equation (21) can be written as the following linear system:(22)I+AV2n+1=I-BV2n+ΔtS2n+1,I+BV2n+2=I-AV2n+1+ΔtS2n+2,where for nonnegative integer m, (23)Vm=v1m,v2m,…,vK-1mT,Sm=s1m,s2m,…,sK-1mT.The matrix entries Aij and Bij for i=1,…,K-1 and j=1,…,K-1 are defined by(24)Aij=0j≥i+1-pw1+α2j=i-pwi-j+1j≤i-1,(25)Bij=0j≥i+2-α2pj=i+1α2pj=i0j≤i-1.Using the Fourier method, the stability of the presented scheme (16) with the initial value is proven. In general, we suppose s=0 and bR=0, then (21) becomes(26)∑k=2l+1-pwkvl-k+12n+1+1-pw1+α2vl2n+1=1-α2pvl2n+α2pvl+12n,1+α2pvl2n+2-α2pvl+12n+2=∑k=2l+1pwkvl-k+12n+1+1+pw1+α2vl2n+1.Let (27)vl2n+1=Z2n+1eiβlhwhere β is a nonnegative integer.So(28)∑k=2l+1-pwke-iβk-1h+1-pw1+α2Z2n+1eiβlh=1-α2p+α2peiβhZ2neiβlh,and(29)1+α2p-α2peiβhZ2n+2eiβlh=∑k=2l+1pwke-iβk-1h+1+pw1+α2Z2n+1eiβlh.Then(30)Z2n+2=G2n+2Z2n,where(31)G2n+2=1-α/2p+α/2peiβh1+α/2p-α/2peiβh∑k=2l+1pwke-iβk-1h+1+pw1+α/2∑k=2l+1-pwke-iβk-1h+1-pw1+α/2=1-α/2p+α/2peiβh1+α/2p-α/2peiβh∑k=1lpwk+1e-iβkh+1+pw1+α/2∑k=1l-pwk+1e-iβkh+1-pw1+α/2≡G1G2.Here(32)G1=a1+iαp/2sin⁡βha2-iαp/2sin⁡βh,G2=b1-icb2+ic.According to Lemma3, (33)a1=1-α2p+α2pcos⁡βh,a2=1+α2p-α2pcos⁡βh,b1=1+p2-α22+∑k=1lpwk+1cos⁡βkh,b2=1-p2-α22-∑k=1lpwk+1cos⁡βkh,c=∑k=1lpwk+1sin⁡βkh.Lemma 4.
Ifs=0 and bR=0 in (21) with α∈[(17-1)/2,2], (34)V2n+2≤V2n.Proof.
First, inequalityG2≤1 is proved. It is obvious that b1-ic/b2+ic=b1-ic/b2-ic, so (35)G22-1=b12-b22b22+c2=4p2-α2/2+∑k=1lwk+1cos⁡βkhb22+c2.From Lemmas 2 and 3, we have(36)∑k=1lwk+1cos⁡βkh≤∑k=1lwk+1=∑k=1lα2gk+1+1-α2gk=α2∑k=1lgk+1+1-α2∑k=1lgk≤α2α-1-1-α2=α22-1.So G22≤1. Furthermore(37)G2≤1.In the same way, G1≤1 is obtained. On combining with (37) and (31), we have(38)G2n+2≤1.This completes the proof.Define(39)V2n+1=V^2n+1+V~2n+1,V2n+2=V^2n+2+V~2n+2.Then the linear system (22) can be written as follows:(40)I+AV^2n+1=I-BV2n,I+AV~2n+1=ΔtS2n+1,I+BV^2n+2=I-AV^2n+1,I+BV~2n+2=I-AV~2n+1+ΔtS2n+2.Through simple computations,(41)I+AV^2n+1=I-BV2n,I+BV^2n+2=I-AV^2n+1.Then we have(42)I+BV~2n+2=I-AI+A-1ΔtS2n+1+ΔtS2n+2.Lemma 5.
The eigenvalues of matricesA and B are nonnegative.Proof.
Letλ be an eigenvalue of matrix A. According to Gerschgorim theorem [25], for every eigenvalue λ, it holds(43)All-R≤λ≤All+R,where R is the radius and (44)R=∑j=1,j≠lK-1Alj=∑j=1l-1pwl-j+1=p∑j=2lwj.By inequality (36), it is obvious that(45)All-Rl=-pw1+α2-p∑j=2lwj=pα22-1-p∑j=2lwj⩾0.From (43), we derive that λ⩾0. The eigenvalues of matrix B which are nonnegative can be proven in the same way.The above two lemmas are sewn together to yield the following theorem that guarantees the stability of the presented scheme (16).Theorem 6.
Supposec(x,t)=c is a constant function. If bR=0, the scheme is uniformly stable about the initial value and right-hand side term with α∈[(17-1)/2,2] and (46)V2n+2≤V0+∑m=12n+2ΔtSm.Proof.
Using Lemma5, we obtain (47)I+B≥max⁡λI+B≥1,I+B-1=I+B-1≤1,and(48)I-AI+A-1≤1.By (42),(49)V~2n+2≤I+B-1I-AI+A-1ΔtS2n+1+ΔtS2n+2≤ΔtS2n+1+ΔtS2n+2.From (41) and Lemma 4,(50)V^2n+2≤V2n.Combining (49), (50) with (39),(51)V2n+2=V^2n+2+V~2n+2≤V2n+ΔtS2n+1+ΔtS2n+2≤V2n-2+ΔtS2n-1+ΔtS2n+ΔtS2n+1+ΔtS2n+2⋯≤V0+∑m=12n+2ΔtSm.This completes the proof.


## 4. Error Estimate for the Finite Difference Scheme

In this section, the error estimate is considered for the presented scheme (16). Suppose c(x,t)=c is a constant.Letuln be the exact solution of (1)~(3) at (xl,tn), then we have(52)ul2n+1-ul2nΔt-sl2n+1-chα∑k=2l+1wkul-k+12n+1+w1+α2ul2n+1-α2ul2n+α2ul+12n=ul2n+1-ul2nΔt-∂u∂tl2n+1-chα∑k=0l+1wkul-k+12n+1+c∂αu∂xαl2n+1+chαα2ul+12n+1-ul2n+1-ul+12n-ul2n.Using Taylor expansion,(53)1hul+12n+1-ul2n+1-ul+12n-ul2n=∂u∂xl2n+1-∂u∂xl2n+h2∂2u∂x2l2n+1-∂2u∂x2l2n+h26∂3u∂x3l2n+1-∂3u∂x3l2n+Oh3=Δt∂2u∂x∂tl2n+1+OΔt2+Δth+Δth2+h3at time-step t2n+1.So, (52) can be written as follows:(54)ul2n+1-ul2nΔt-sl2n+1-chα∑k=2l+1wkul-k+12n+1+w1+α2ul2n+1-α2ul2n+α2ul+12n=OΔt+h2+cΔthα-1α2∂2u∂x∂tl2n+1+OΔt2h1-α+Δth2-α+Δth3-α+h4-α.Since1<α<2,(55)ul2n+1-ul2nΔt-sl2n+1-chα∑k=2l+1wkul-k+12n+1+w1+α2ul2n+1-α2ul2n+α2ul+12n=cΔthα-1α2∂2u∂x∂tl2n+1+OΔt2h1-α+Δt+h2.In the same way,(56)ul2n+2-ul2n+1Δt-sl2n+2-chα∑k=2l+1wkul-k+12n+1+w1+α2ul2n+1-α2ul2n+2+α2ul+12n+2=-cΔthα-1α2∂2u∂x∂tl2n+1+OΔt2h1-α+Δt+h2at time-step t2n+2.Let(57)Um=u1m,u2m,…,uK-1m.Theorem 7.
Letu(xl,tn) be the exact solution of (1)~(3) and vln be the solution of (16). Suppose c(x,t) is a constant function. When h and Δt are sufficiently small, there exists a positive constant C independent of Δt and h such that for m≤T/Δt, (58)Vm-Um≤CΔt2h-2α-1+Δt+h2.Proof.
We defineel2n+1 as follows at the odd time-step t2n+1(59)e02n+1=v02n+1-u02n+1=0,el2n+1=vl2n+1-ul2n+1+cΔt2hα-1α2∂2u∂x∂tl2n+1,1≤l≤K-1,eK2n+1=vK2n+1-uK2n+1=0.Similarly, we define el2n+2 the even time-step t2n+2(60)e02n+2=v02n+2-u02n+2=0,el2n+2=vl2n+2-ul2n+2,1≤l≤K-1,eK2n+2=vK2n+2-uK2n+2=0.According to (56), (55), and (16),(61)el2n+1-el2nΔt=chα∑k=2l+1wkel-k+12n+1+w1+α2el2n+1-α2el2n+α2el+12n-cΔt2h2α-1α2∑k=2l+1wkc∂2u∂x∂tl-k+12n+1-c∂2u∂x∂tl2n+1+OΔt2h1-α+Δt+h2,el2n+2-el2n+1Δt=chα∑k=2l+1wkel-k+12n+1+w1+α2el2n+1-α2el2n+2+α2el+12n+2-cΔt2h2α-1α2∑k=2l+1wkc∂2u∂x∂tl-k+12n+1-c∂2u∂x∂tl2n+1+OΔt2h1-α+Δt+h2.
Using Taylor expansion,(62)c∂2u∂x∂tl-k+12n+1-c∂2u∂x∂tl2n+1=Ok-1h.Through the estimates for the gk in [9], there exist two constants K1 and K2 such that(63)K1kα+1≤gk≤K2kα+1.Combining (12) with (63), we know(64)C1kα+1≤wk≤C2kα+1.Here C1 and C2 are constant.
Using (62) and (64), we have(65)wkc∂2u∂x∂tl-k+12n+1-c∂2u∂x∂tl2n+1≤1kα+1Ck-1h≤C3kαh.Here C and C3 are constant.
Since1<α<2, there exists a positive constant K independent of h such that(66)∑k=2l+1wkc∂2u∂x∂tl-k+12n+1-c∂2u∂x∂tl2n+1≤∑k=2l+1C3kαh≤Kh.So (61) can be written as follows:(67)el2n+1-el2nΔt=chα∑k=2l+1wkel-k+12n+1+w1+α2el2n+1-α2el2n+α2el+12n+OΔt2h-2α-1+Δt+h2,el2n+2-el2n+1Δt=chα∑k=2l+1wkel-k+12n+1+w1+α2el2n+1-α2el2n+2+α2el+12n+2+OΔt2h-2α-1+Δt+h2.
It is obvious that(68)el0=0,0≤l≤K.Set (69)Em=e1m,e2m,…,eK-1m.Using Theorem 6,(70)Em≤CΔt2h-2α-1+Δt+h2.Remark 8.
From (55) the truncation error of the scheme (16) is O∆th-∝-1+∆t+h2. Here it is proven that the error between the discrete and the analytical solutions is bounded by C∆t2h-2∝-1+∆t+h2 by using the asymmetric discretisation technique.


## 5. Numerical Experiments

In this section, numerical examples of one-dimensional cases are presented to show the effectiveness and convergence orders of the scheme (16).We consider the following space-fractional diffusion problem:(71)∂ux,t∂t=cx∂αux,t∂xα+sx,twith 0<x<1 and 0<t≤1.In order to describe the numerical results, the discretel2-norms, the discrete l∞-norms, and the corresponding convergence rates are defined as follows: (72)enl2=ehnl2=∑j=1K-1ejn2h1/2,enl∞=max1≤j≤K-1⁡ejn,Rate=log⁡eh1nv/eh2nvlog⁡h1/h2,v=l2,l∞,where ejn=u(xj,tn)-vjn, n=0,1,…,N-1, and h1, h2 are space mesh sizes.The discretizing mesh sizes are given according to the following rules.Case 1.
SetΔt=Chα, where C is a nonnegative integer. By Theorem 7, the convergence order is O(hα). We get different values of C to verify the theoretical analysis.Case 2.
SetΔt=5h2 and take the convergence order is O(h2).Two numerical experiments are given with different values ofc(x), α, and the source function s(x,t). The errors, corresponding convergence rates, and CPU times are listed in the Tables 1–8. The results in Tables 1–8 are in agreement with the theoretical analysis. In addition, these tables show that, even with very small time-steps, the CPU execution times stay very small as well.Table 1
Error and convergence rates for Example1.
α = 1.6 ,   Δ t = 4 h 1.6 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,10) 8.1450 × 1 0 - 2 - 2.5887 × 1 0 - 2 - 0.001  (20,32) 2.3726 × 1 0 - 2 1.78 5.3215 × 1 0 - 3 2.28 0.001  (40,92) 8.0849 × 1 0 - 3 1.55 1.2802 × 1 0 - 3 2.06 0.003  (80,278) 2.6561 × 1 0 - 3 1.61 3.0010 × 1 0 - 4 2.09 0.026  (160,840) 8.7693 × 1 0 - 4 1.60 7.7019 × 1 0 - 5 1.96 0.12  (320,2548) 2.8887 × 1 0 - 4 1.60 2.3741 × 1 0 - 5 1.70 0.9Table 2
Error and convergence rates for Example1.
α = 1.8 ,   Δ t = 5 h 1.8 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,14) 5.6493 × 1 0 - 2 - 2.3870 × 1 0 - 2 - 0.001  (20,44) 1.7108 × 1 0 - 2 1.72 6.6279 × 1 0 - 3 1.85 0.002  (40,154) 4.8088 × 1 0 - 3 1.83 1.5575 × 1 0 - 3 2.09 0.007  (80,534) 1.3804 × 1 0 - 3 1.82 3.7624 × 1 0 - 4 2.05 0.034  (160,1856) 3.9664 × 1 0 - 4 1.80 9.0186 × 1 0 - 5 2.06 0.389  (320,6462) 1.1384 × 1 0 - 4 1.80 2.1685 × 1 0 - 5 2.06 2.861Table 3
Error and convergence rates for Example1.
α = 1.6 ,   Δ t = 5 h 2 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,20) 3.8690 × 1 0 - 2 - 1.2379 × 1 0 - 2 - 0.001 (20,80) 9.3129 × 1 0 - 3 2.05 2.1854 × 1 0 - 3 2.50 0.003 (40,320) 2.3064 × 1 0 - 3 2.01 4.2095 × 1 0 - 4 2.38 0.017 (80,1280) 5.7525 × 1 0 - 4 2.00 8.9045 × 1 0 - 5 2.24 0.048 (160,5120) 1.4371 × 1 0 - 4 2.00 2.0923 × 1 0 - 5 2.09 0.887 (320,20480) 3.5912 × 1 0 - 5 2.00 5.8087 × 1 0 - 6 1.85 8.326Table 4
Error and convergence rates for Example1.
α = 1.8 ,   Δ t = 5 h 2 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,20) 3.8690 × 1 0 - 2 - 1.3752 × 1 0 - 2 - 0.001  (20,80) 9.3129 × 1 0 - 3 2.05 2.3787 × 1 0 - 3 2.53 0.002  (40,320) 2.3064 × 1 0 - 3 2.01 4.1125 × 1 0 - 4 2.53 0.007  (80,1280) 5.7525 × 1 0 - 4 2.00 6.9437 × 1 0 - 5 2.57 0.053  (160,5120) 1.4371 × 1 0 - 4 2.00 1.1636 × 1 0 - 5 2.58 0.723  (320,20480) 3.5912 × 1 0 - 5 2.00 2.0268 × 1 0 - 6 2.52 7.104Table 5
Error and convergence rates for Example2.
α = 1.6 ,   Δ t = 4 h 1.6 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,10) 8.1450 × 1 0 - 2 - 2.5828 × 1 0 - 2 - 0.001  (20,32) 2.3726 × 1 0 - 2 1.78 5.3327 × 1 0 - 3 2.28 0.001  (40,92) 8.0849 × 1 0 - 3 1.55 1.2920 × 1 0 - 3 2.01 0.003  (80,278) 2.6561 × 1 0 - 3 1.62 3.1102 × 1 0 - 4 2.05 0.039  (160,840) 8.7693 × 1 0 - 4 1.60 8.2125 × 1 0 - 5 1.92 0.131  (320,2548) 2.8887 × 1 0 - 4 1.60 2.5244 × 1 0 - 5 1.70 2.193Table 6
Error and convergence rates for Example2.
α = 1.8 ,   Δ t = 5 h 1.8 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,14) 5.6493 × 1 0 - 2 - 1.9797 × 1 0 - 2 - 0.001  (20,44) 1.7108 × 1 0 - 2 1.72 5.1683 × 1 0 - 3 1.94 0.001  (40,154) 4.8088 × 1 0 - 3 1.83 1.2170 × 1 0 - 3 2.09 0.006  (80,534) 1.3804 × 1 0 - 3 1.82 2.9246 × 1 0 - 4 2.06 0.058  (160,1856) 3.9664 × 1 0 - 4 1.80 6.9376 × 1 0 - 5 2.08 0.456  (320,6462) 1.1384 × 1 0 - 4 1.80 1.6405 × 1 0 - 5 2.08 3.441Table 7
Error and convergence rates for Example2.
α = 1.6 ,   Δ t = 5 h 2 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,20) 3.8690 × 1 0 - 2 - 1.2540 × 1 0 - 2 - 0.001  (20,80) 9.3129 × 1 0 - 3 2.05 2.2198 × 1 0 - 3 2.50 0.002  (40,320) 2.3064 × 1 0 - 3 2.01 4.2601 × 1 0 - 4 2.38 0.011  (80,1280) 5.7525 × 1 0 - 4 2.00 8.9084 × 1 0 - 5 2.26 0.137  (160,5120) 1.4372 × 1 0 - 4 2.00 2.0611 × 1 0 - 5 2.11 0.837  (320,20480) 3.5921 × 1 0 - 5 2.00 5.6286 × 1 0 - 6 1.87 7.769Table 8
Error and convergence rates for Example2.
α = 1.8 ,   Δ t = 5 h 2 1/h,1/Δt e l ∞ Rate e l 2 Rate CPU(s) (10,20) 3.8690 × 1 0 - 2 - 1.2598 × 1 0 - 2 - 0.001  (20,80) 9.3129 × 1 0 - 3 2.05 2.1751 × 1 0 - 3 2.53 0.002  (40,320) 2.3064 × 1 0 - 3 2.01 3.7881 × 1 0 - 4 2.52 0.006  (80,1280) 5.7525 × 1 0 - 4 2.00 6.5455 × 1 0 - 5 2.53 0.19  (160,5120) 1.4372 × 1 0 - 4 2.00 1.1493 × 1 0 - 5 2.51 0.834  (320,20480) 3.5921 × 1 0 - 5 2.00 2.0426 × 1 0 - 6 2.49 10.407The results in Tables9–12 are listed by using the finite difference scheme developed by Rui and Huang [19] for Examples 1 and 2. We set Δt=hα-0.5 and have the same spatial grid sizes h from 1/10 to 1/320, since their error estimates between the numerical and analytical solutions are optimal under the condition Δt=hα-0.5 for α∈(1,2). Comparing our results in Tables 1–8 with those in Tables 9–12, one can see that the error estimate produced by our method is better than that in [19] in the sense that their convergence rate is O(h). Furthermore, the time-steps in Tables 1–8 are smaller than that in Tables 9–12 with the spatial grid step staying the same. However the CPU execution times increase by a few seconds, which can prove the validity of our presented scheme.Table 9
Error and convergence rates for Example1 by [19].
α = 1.6 ,   Δ t = h 1.1 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,14) 5.6493 × 1 0 - 2 - 1.7916 × 1 0 - 2 - 0.001  (20,28) 2.7238 × 1 0 - 2 1.05 6.2433 × 1 0 - 3 1.53 0.001  (40,58) 1.2907 × 1 0 - 2 1.08 2.5313 × 1 0 - 3 1.30 0.001  (80,124) 5.9817 × 1 0 - 3 1.11 1.1917 × 1 0 - 3 1.09 0.015  (160,266) 2.7764 × 1 0 - 3 1.11 6.1784 × 1 0 - 4 0.95 0.064  (320,570) 1.2931 × 1 0 - 3 1.10 3.2979 × 1 0 - 4 0.91 0.186Table 10
Error and convergence rates for Example1 by [19].
α = 1.8 ,   Δ t = h 1.3 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,20) 3.8690 × 1 0 - 2 - 1.3867 × 1 0 - 2 - 0.001  (20,50) 1.5013 × 1 0 - 2 1.35 5.2908 × 1 0 - 3 1.39 0.001  (40,122) 6.0805 × 1 0 - 3 1.30 2.6276 × 1 0 - 3 1.01 0.002  (80,298) 2.4773 × 1 0 - 3 1.30 1.3979 × 1 0 - 3 0.91 0.019  (160,734) 1.0038 × 1 0 - 3 1.30 7.3133 × 1 0 - 4 0.93 0.074  (320,1806) 5.1822 × 1 0 - 4 0.94 3.8051 × 1 0 - 4 0.94 1.295Table 11
Error and convergence rates for Example2 by [19].
α = 1.6 ,   Δ t = h 1.1 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,14) 5.6493 × 1 0 - 2 - 1.8277 × 1 0 - 2 - 0.001  (20,28) 2.7238 × 1 0 - 2 1.05 6.1801 × 1 0 - 3 1.56 0.001  (40,58) 1.2907 × 1 0 - 2 1.08 2.2795 × 1 0 - 3 1.44 0.001  (80,124) 5.9817 × 1 0 - 3 1.11 1.0000 × 1 0 - 3 1.19 0.006  (160,266) 2.7764 × 1 0 - 3 1.11 5.1210 × 1 0 - 4 0.97 0.071  (320,570) 1.2931 × 1 0 - 3 1.10 2.7682 × 1 0 - 4 0.89 0.199Table 12
Error and convergence rates for Example2 by [19].
α = 1.8 ,   Δ t = h 1.3 1 / h , 1 / Δ t e l ∞ Rate e l 2 Rate CPU(s) (10,20) 3.8690 × 1 0 - 2 - 1.2610 × 1 0 - 2 - 0.001  (20,50) 1.5013 × 1 0 - 2 1.37 4.2105 × 1 0 - 3 1.29 0.001  (40,122) 6.0805 × 1 0 - 3 1.30 1.9787 × 1 0 - 3 1.09 0.003  (80,298) 2.4773 × 1 0 - 3 1.30 1.0719 × 1 0 - 3 0.88 0.019  (160,734) 1.0038 × 1 0 - 3 1.30 5.7552 × 1 0 - 4 1.57 0.086  (320,1806) 4.4583 × 1 0 - 4 0.95 3.0530 × 1 0 - 4 0.91 0.633Example 1.
The diffusion coefficientc(x,t), the analytical solution u, and source coefficient s(x,t) are as follows:(73)cx,t=1.0,u=e-tx1+α,s=-e-tx1+α+Γ2+αx,u0,t=0,u1,t=e-t.The numerical results are listed in Tables1–4.Example 2.
The diffusion coefficientc(x,t), the analytical solution u, and source coefficient s(x,t) are as follows:(74)cx,t=xα,u=e-tx1+α,s=-e-tx1+α1+Γ2+α,u0,t=0,u1,t=e-t.The numerical results are listed in Tables5–8.


## 6. Conclusion

In this paper, a uniformly stable explicit finite difference scheme is presented to solve one-dimensional space-fractional diffusion equation. The scheme is constructed combining a WSGD operator with an asymmetric technique. It can be solved explicitly by using different nodal-point stencils at odd and even time-steps. Results of two numerical experiments show that the error estimates are optimal whenΔt=O(Ch2) and Δt=O(Chα), where C is a nonnegative integer. There may exist super-convergence of the proposed semi-implicit finite difference method, which will be studied in the future.



---

*Source: 4238420-2019-05-28.xml*