# Numerical Computation and Stability Analysis for the Fractional Subdiffusions with Spatial Variable Coefficients

Authors: Lei Ren
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this paper, we propose an efficient compact finite difference method for a class of time-fractional subdiffusion equations with spatially variable coefficients. Based on theL2-1σ approximation formula of the time-fractional derivative and a fourth-order compact finite difference approximation to the spatial derivative, an efficient compact finite difference method is developed. The local truncation error and the solvability of the developed method are discussed in detail. The unconditional stability of the resulting scheme and also its convergence of second-order in time and fourth-order in space are rigorously proved using a discrete energy analysis method. Numerical examples are provided to demonstrate the accuracy and the theoretical results.


---

## Body

## 1. Introduction

The fractional differential equations are widely used to describe plenty of nature phenomena in physics, biology, and chemistry [1–10]. Fractional subdiffusion is the evolution equation for the probability density function that describes particles diffusing with mean square displacement x2t∼tα0<α<1 [8].In this paper, we consider a class of time-fractional subdiffusion equations with variable coefficients as follows:(1)∂u∂tx,t+D0Ctαux,t=∂∂xφx∂u∂x+fx,t,x,t∈0,L×0,T,with its initial boundary conditions(2)u0,t=ϕ0t,uL,t=ϕLt,t∈0,T,(3)ux,0=γx,x∈0,L,where D0Ctαux,t is the Caputo fractional derivative(4)D0Ctαux,t=1Γ1−α∫0t∂u∂sx,st−s−αds,0<α<1.In practice, the motion time and the status of particles are described by the time drift term∂u/∂t. Equation (1) has many important applications in various fields [11–13]. Generally, various numerical methods are used to obtain the numerical solutions of fractional differential, such as the finite difference method [14–18], the finite element method [19], the spectral method [20, 21], and so on.In the construction of high-order numerical methods for equation (1), it is crucial to obtain a stable high-order approximation to the Caputo time-fractional derivative D0Ctαvx,t0<α<1. A common technique to design a second-order or more higher-order approximation to the fractional derivative is to use the weighted and shifted Grünwald difference operator. The operator is often used to handle the Riemann–Liouville-type fractional derivative (see [22–25]). Some of its applications to the Caputo-type fractional derivative D0Ctαvx,t0<α<1 can be found in [24, 26, 27]. Another alternative method was given in [28], where a modified L1 approximation formula, called L1-2 formula, was used to discretize the Caputo time-fractional derivative D0Ctαvx,t0<α<1, but the rigorous convergence analysis for the corresponding difference scheme has not been available. In [29], a new approximation to the Caputo time-fractional derivative D0Ctαvx,t0<α<1 with the convergence order 3−α, called L2-1σ formula, was constructed.Until now, there is a little discussion for time-fractional subdiffusion equations with spatial variable coefficients [30–32]. Zhao and Xu [33] proposed an efficient compact difference scheme for solving the fractional subdiffusion equation with the spatial variable coefficient. The stability and the convergence of the scheme are proved. In [34], Vong proposed a compact finite difference scheme for spatially variable coefficient fractional subdiffusion with Neumann boundary conditions. The analysis of stability and convergence is studied by the energy method. Unfortunately, the temporal accuracy of the aforementioned methods was still less than two order. In general, the weighted and shifted Grünwald difference operator [22–27] is a usual approximation method to the Caputo time-fractional derivative D0Ctαvx,t0<α<1. Recently, an efficient L1-2 method [28, 29] was used to discretize the Caputo time-fractional derivative. In this paper, we propose some new methods to construct a high-order compact finite difference method for equation (1).


## 2. Compact Finite Difference Method

For positive integersMandN, we let τ=T/N be the time step and h=L/M be the spatial step. Let tn=nτ0≤n≤N and tn−α/2=n−α/2τ1≤n≤N. We partition 0,L into a mesh by the mesh points xi=ih0≤i≤M. Denote xi−1/2=i−1/2h1≤i≤M. Define spatial difference operators as follows:(5)δxui−1/2=1hui−ui−1,δx2ui=1h2ui+1−2ui+ui−1,ℋxui=ui+h212δx2ui−h212δxφ′φui.Let(6)ck,nα=a0α,k=0,n=1,akα+bk+1α−bkα,0≤k≤n−2,n≥2,an−1α−bn−1α,k=n−1,where(7)σ=1−α2ηk=k+σk≥0,a0α=η01−α,b0α=0,akα=ηk1−α−ηk−11−α,bkα=12−αηk2−α−ηk−12−α−12ηk1−α+ηk−11−αk≥1.Lemma 1 (see [35]).
It holds(8)0<bkα<α22−αk+σ1−α−k−1+σ1−αk≥1,(9)1−α2k+σ−α<akα−bkα<k+σ1−α−k−1+σ1−αk≥1,(10)ck−1,nα>ck,nα>0,n≥2,1≤k≤n−1,(11)1−αc0,nα−σc1,nα>0,n≥2.Lemma 2.
(assumeut∈C30,T). When n≥1, we have(12)D0Ctαutn−α/2=1μ∑k=1ncn−k,nαutk−utk−1+Oτ3−α,where μ=ταΓ2−α.Proof.
The proof follows from Lemma 2 of [29].Lemma 3 (see [35]).
Assumeut∈C30,T. We obtain(13)dudtt1−α/2=2−ατut1−ut0−1−αdudtt0+Oτ2,and for n≥2,(14)dudttn−α/2=12τ3−αutn−4−2αutn−1+1−αutn−2+Oτ2.Lemma 4 (see [35]).
Supposeut∈C20,T. When n≥1, we have(15)utn−α/2=α2utn−1+1−α2utn+Oτ2.Lemma 5.
Supposeux∈C6,30,L×0,T. It holds that(16)ℋx∂∂xφx∂u∂xxi,tn=δxψδxuxi,tn+Rxi,tn,1≤i≤M−1,1≤n≤N,where(17)ψ=φ−h212φ′2φ−12φ″,Rxi,tn=Oh4.Proof.
The proof follows from (31)–(33) in [32].
Next, we discretize (1) into a compact finite difference scheme. We define(18)δtun−1/2=1τun−un−11≤n≤N,Dt^un=12τ3−αun−4−2αun−1+1−αun−22≤n≤N,un,α/2=α2un−1+1−α2un1≤n≤N.
Supposeux,t be the solution of (1), we define the grid functions as follows: (19)Uin=uxi,tn,Win=∂u∂txi,tn,Zin=∂∂xφx∂u∂xxi,tn,φi=φxi,fin=fxi,tn,fin−α/2=fxi,tn−α/2,ϕ0n=ϕ0tn,ϕLn=ϕLtn,γi=γxi.
We consider equation (1) at the point xi,tn−α/2, and(20)Win−α/2+D0Ctαuxi,tn−α/2=dZin−α/2+fin−α/2.
From Lemmas 2–4, we have(21)Dt^Uin+1μ∑k=1ncn−k,nαUik−Uik−1=Zin,α/2+fin−α/2+Rtαin,1≤i≤M−1,2≤n≤N,(22)2−αδtUi1/2+1μc0,1αUi1−Ui0=Zi1,α/2+1−αWi0+fi1−α/2+Rtαi1,1≤i≤M−1,where(23)Rtαin=Oτ2+Oτ3−α,1≤i≤M−1,1≤n≤N.
Using Lemma 5, we obtain(24)ℋxDt^Uin+1μ∑k=1ncn−k,nαℋxUik−Uik−1=δxψδxUin,α/2+ℋxfin−α/2+Rtxαin,1≤i≤M−1,2≤n≤N,(25)2−αℋxδtUi1/2+1μc0,1αℋxUi1−Ui0=δxψδxUi1,α/2+1−αℋxWi0+ℋxfi1−α/2+Rtxαi1,1≤i≤M−1,where(26)Rtxαin=ℋxRtαin+dRxin,1≤i≤M−1,1≤n≤N.
SinceD0Ctαuxi,0=0 [4], we have Wi0=fi∗,0, where(27)fi∗,0=∂2γ∂x2xi+fxi,0,1≤i≤M−1.
We obtain the following compact finite difference scheme:(28)ℋxDt^uin+1μ∑k=1ncn−k,nαℋxuik−uik−1=δxψδxuin,α/2+ℋxfin−α/2,1≤i≤M−1,2≤n≤N,2−αℋxδtui1/2+1μc0,1αℋxui1−ui0=δxψδxui1,α/2+1−αℋxfi∗,0+ℋxfi1−α/2,1≤i≤M−1,u0n=ϕ0n,uMn=ϕLn,1≤n≤N,ui0=γi,0≤i≤M.


## 3. Truncation Error and Solvability

Assumeux,t∈C6,30,L×0,T, and there exists a positive constant C∗ such that(29)Rtαin≤C∗τ2,Rxin≤C∗h4.Applying (29) to (26), we have the following theorem:Theorem 1.
Supposeux,t∈C6,30,L×0,T. The truncation error Rtxαin satisfies(30)Rtxαin≤C1τ2+h4,1≤i≤M−1,1≤n≤N,where C1 is a positive constant.
For convenience, we construct the matrix form of the compact difference scheme (28). Let(31)wk,nα=c0,1α,k=0,n=1,−c0,1α,k=1,n=1,c0,nα,k=0,n≥2,ck,nα−ck−1,nα,1≤k≤n−1,n≥2,−cn−1,nα,k=n,n≥2.
Define the column vectors as follows:(32)un=u1n,u2n,…,uM−1nT,un−2,∗=u1n−2,∗,u2n−2,∗,…,uM−1n−2,∗T,fn−α/2=f1n−α/2,f2n−α/2,…,fM−1n−α/2T,f∗,0=f1∗,0,g2∗,0,…,fM−1∗,0T,φ′φn=ζn=ζ1n,ζ2n,…,ζM−1nT,ψn=ψ1n,ψ2n,…,ψM−1nT,where(33)uin−2,∗=∑k=2nwk,nαuin−k,1≤i≤M−1,2≤n≤N.
Define the followingM−1-order tridiagonal or diagonal matrices:(34)A=1h2ψx1/2+ψx3/2−ψx3/2−ψx3/2ψx3/2+ψx5/2−ψx5/2⋱⋱⋱−ψxM−5/2ψxM−5/2+ψxM−3/2−ψxM−3/2−ψxM−3/2ψxM−3/2+ψxM−1/2,B1=112tridiag1,10,1,B2=h24tridiag1,0,−1,Q=diagζ1,ζ2,…,ζM−1,Q1=2−ατ+w0,1αμB1+B2Q+2−α2h2A,Qn=3−α2τ+w0,nαμB1+B2Q+2−α2h2A,2≤n≤N,Q˜n=2−ατ−w1,nαμB1+B2Q−α2h2A,1≤n≤N.
According to some simple calculations, the scheme (28) can be expressed in the matrix form as (35)Qnun=Q˜nun−1−B1+B2Q1−α2τun−2+1μun−2,∗−fn−α/2+rn,2≤n≤N,Q1u1=Q˜1u0+B1+B2Q1−αf∗,0+f1−α/2+r1,where rn absorbs the boundary values of the solution vector and the source term.Theorem 2.
The compact difference scheme (28) is uniquely solvable if and only if the matrices Qn 1≤n≤N are all nonsingular.
Assume(36)c1≤ψ≤c2,φ′φ=ζ≤c.
A sufficient condition for the matricesQn 1≤n≤N to be all nonsingular is given by(37)hc<8.Corollary 1.
The compact difference scheme (28) is uniquely solvable if condition (37) holds true.Proof.
At first, we prove the matrixQn2≤n≤N is nonsingular. Since w0,nα=a0α+b1αn≥2, we have Qn=tridiagpi−1∗,qi∗,pi+1∗, where(38)pi−1∗=112+h24ζi3−α2τ+a0α+b1αμ−2−α2h2ψxi−1/2,i≠1,pi+1∗=112+h24ζi3−α2τ+a0α+b1αμ−2−α2h2ψxi+1/2,(39)qi∗=563−α2τ+a0α+b1αμ+2−α2h2ψxi−1/2+ψxi+1/2.
Sinceqi∗>01≤i≤M−1, we have(40)pi−1∗+pi+1∗≤16+h24ζi−1+ζi+13−α2τ+a0α+b1αμ+2−α2h2ψxi−1/2+ψxi+1/2≤16+h12c3−α2τ+a0α+b1αμ+2−α2h2ψxi−1/2+ψxi+1/2<563−α2τ+a0α+b1αμ+2−α2h2ψxi−1/2+ψxi+1/2=qi∗.
Similarly,(41)p2∗≤112+h24ζ23−α2τ+a0α+b1αμ+2−α2h2ψx3/2<563−α2τ+a0α+b1αμ+2−α2h2ψx3/2+ψx1/2=q1∗,(42)pM−2∗≤112+h24ζM−23−α2τ+a0α+b1αμ+2−α2h2ψxM−3/2<563−α2τ+a0α+b1αμ+2−α2h2ψxM−3/2+ψxM−1/2=qM−1∗.
Therefore,Qn 1≤n≤N are strictly diagonally dominant and thus nonsingular [36].


## 4. Stability and Convergence of the Proposed Difference Scheme

LetSh=uu=u0,u1,…,uM,u0=uM=0, and for any grid functions u,v∈Sh, we define(43)δxu,δxv=h∑i=0M−1δxui+1/2δxvi+1/2,u1=δxu,δxu1/2,u1=u2+u121/2,δxu,δxvψ=h∑i=0M−1ψxi+1/2δxui+1/2δxvi+1/2,δxuψ=δxu,δxuψ1/2.For any grid functionu,v∈Sh, we have (44)δx2u,v=−δxu,δxv,hδx2u≤2u1,hu1≤2u.Using Lemma 2.6 and Theorem 2.1 in [33], we obtain ℋu,u≥1/4u2. Next, we define the discrete inner product as follows:(45)u,v=ℋxu,v,uϵ=u,u1/2.Lemma 6 (see [35]).
For any grid functionu=uin0≤i≤M,0≤n≤N∈Sh, we have(46)ℋxDt^un,un,α/2≥14τEn−En−1,2≤n≤N,where for 1≤n≤N,(47)En=3−αunε2−1−αun−1ε2+4−α1−α2un−un−1ε2≥22−αunε2.Lemma 7.
For any grid functionuin∈Sh, we obtain(48)∑k=1ncn−k,nαℋxuk−uk−1,un,α/2≥12∑k=1ncn−k,nαukε2−uk−1ε2,1≤n≤N.Proof.
Using the definition of the inner product⋅,⋅,(49)∑k=1ncn−k,nαℋxuk−uk−1,un,α/2=∑k=1ncn−k,nαuk−uk−1,un,α/2.
Applying Corollary 1 [29], we have(50)∑k=1ncn−k,nαuk−uk−1,un,α/2≥12∑k=1ncn−k,nαukϵ2−uk−1ϵ2.
The desired result (48) follows immediately.Lemma 8 (discrete Gronwall lemma [37]).
SupposeKn and Sn are nonnegative sequences, and λn satisfies(51)λ0≤g0,λn≤g0+∑l=0n−1Sl+∑l=0n−1Klλl,n≥1,where g0≥0. Then, the sequence λn satisfies(52)λn≤g0+∑l=0n−1Slexp∑l=0n−1Kl,n≥1.Theorem 3.
Supposeun=u0n,u1n,…,uMn be the solution of the compact difference scheme (28) with u0n=uMn=0. Then, when τ≤1/42−α3, it holds(53)unε2≤exp42−α2−α2+α2TGn∗,1≤n≤N,where(54)Gn∗=2−α1+22−α2T1−αΓ2−αG0+8τGn,with(55)G0=3−α+2T1−αΓ2−α+4τα2φε2,(56)Gn=1+22−α2T1−αΓ2−α1−αℋxf∗,0+ℋxf1−α/2ε2+∑k=2nℋxfk−α/2ε2.Proof.
Applying the inner product to equation (28) with un,α/2, we obtain(57)ℋxDt˜un,un,α/2+1μ∑k=1ncn−k,nαℋxuk−uk−1,un,α/2=−δxun,α/2ψ2+ℋxfn−α/2,un,α/2,2≤n≤N.
By Lemmas 6 and 7,(58)14τEn−En−1+12μ∑k=1ncn−k,nαukε2−uk−1ε2≤ℋxfn−α/2,un,α/2,2≤n≤N,where(59)En=3−αunε2−1−αun−1ε2+4−α1−α2un−un−1ε2≥22−αunε2.
Using the definition ofcn−k,nα, we have(60)∑k=1ncn−k,nαukε2−uk−1ε2=bn−1α−an−1αu0ε2−bnαu1ε2+∑k=0n−1akα+bk+1α−bkαun−kε2−∑k=0n−2akα+bk+1α−bkαun−1−kε2.
By the Cauchy–Schwarz inequality and Lemma 6, we obtain(61)ℋxfn−α/2,un,α/2≤142−α2un2+α2un−12+12ℋxfn−α/22.
Let(62)Fn=14τEn+12μ∑k=0n−1akα+bk+1α−bkαun−kε2,n≥1.
Substituting (60) and (61) into (58), we have(63)Fn≤Fn−1+142−α2un2+α2un−12+12μbnαu1ε2+an−1α−bn−1αu0ε2+12ℋxfn−α/22,2≤n≤N,or equivalently,(64)Fn≤F1+142−α2∑k=2nuk2+α2∑k=1n−1uk2+12μ∑k=2nbkαu1ε2+ak−1α−bk−1αu0ε2+12∑k=2nℋxfk−α/22,2≤n≤N.
Taking the inner product of equation (28) with u1,α/2, we obtain(65)2−αℋxδtu1/2,u1,α/2+1μc0,1αℋxu1−u0,u1,α/2=−δxu1,α/2ψ2+1−αℋxf∗,0+ℋxf1−α/2,u1,α/2.
It is easy to see that(66)ℋxδtu1/2,u1,α/2=12τu1ε2−u0ε2+1−αu1−u0ε2.
Applying Lemma 7, we have(67)c0,1αℋxu1−u0,u1,α/2≥12a0αu1ε2−u0ε2.
Substituting (66) and (67) into (65), we obtain(68)12τ2−αu1ε2+1−αu1−u0ε2+12μa0αu1ε2≤121τ2−α+1μa0αu0ε2+142−α2u12+α2u02+121−αℋxf∗,0+ℋxf1−α/22.
Since3−α<22−α and 4−α<42−α, we have(69)F1<12τ2−αu1ε2+1−αu1−u0ε2+12μa0αu1ε2−14τ1−αu0ε2+12μb1αu1ε2.Thus,(70)F1<1212τ3−α+1μa0αu0ε2+12μb1αu1ε2+142−α2u12+α2u02+121−αℋxf∗,0+ℋxf1−α/22.
Using (63), we have(71)Fn≤14τ3−αu0ε2+142−α2∑k=1nuk2+α2∑k=0n−1uk2+12μ∑k=1nbkαu1ε2+ak−1α−bk−1αu0ε2+121−αℋxf∗,0+ℋxf1−α/22+∑k=2nℋxfk−α/22,2≤n≤N.
By the definition ofFn, equation (60), and akα+bk+1α−bkα>0 (see (10)), we have(72)122−ατunε2−142−α2un2≤14τ3−α+12μ∑k=1nak−1α−bk−1αu0ε2+α24u02+12μ∑k=1nbkαu1ε2+121−αℋxf∗,0+ℋxf1−α/22+∑k=2nℋxfk−α/22+142−α2+α2∑k=1n−1uk2,2≤n≤N.
Applying Lemma 6, we obtain(73)122−ατ−2−α2unε2≤14τ3−α+12μ∑k=1nak−1α−bk−1αu0ε2+α2u0ε2+12μ∑k=1nbkαu1ε2+21−αℋxf∗,0+ℋxf1−α/2ε2+∑k=2nℋxfk−α/2ε2+2−α2+α2∑k=1n−1ukε2,2≤n≤N.
Whenτ≤1/42−α3, we have (74)unε2≤2−α3−α+2τμ∑k=1nak−1α−bk−1αu0ε2+4α2τu0ε2+2τμ∑k=1nbkαu1ε2+8τ1−αℋxf∗,0+ℋxf1−α/2ε2+∑k=2nℋxfk−α/2ε2+4τ2−α2−α2+α2∑k=1n−1ukε2,2≤n≤N.
By Lemma 1,(75)τμ∑k=1nak−1α−bk−1α<τ1−αΓ2−αn1−α≤T1−αΓ2−α,τμ∑k=1nbkα<τ1−αΓ2−αn+11−α≤2T1−αΓ2−α.
Therefore, we have(76)unε2≤2−αG0+22T1−αΓ2−αu1ε2+82−ατ1−αℋxf∗,0+ℋxf1−α/2ε2+∑k=2nℋxgk−α/2ε2+4τ2−α2−α2+α2∑k=1n−1ukε2,2≤n≤N,where G0 is defined by (55). Similarly, using (70), we have(77)u1ε2≤2−αG0+8τ1−αℋxf∗,0+ℋxf1−α/22.
Substituting (77) into (76), we get(78)unε2≤2−α1+22−α2T1−αΓ2−αG0+8τGn+4τ2−α2−α2+α2∑k=1n−1ukε2,2≤n≤N,where Gn is defined by (56). From Lemma 3, this completes the proof.
Letein=Uin−uin. From (24), (25), and (28), we get the following error equation:(79)ℋxDt^ein+1μ∑k=1ncn−k,nαℋxeik−eik−1=δx2ein,α/2+Rtxαin,1≤i≤M−1,2≤n≤N,2−αℋxδtei1/2+1μc0,1αℋxei1−ei0=δx2ei1,α/2+Rtxαi1,1≤i≤M−1,e0n=eMn=0,1≤n≤N,ei0=0,0≤i≤M.Theorem 4.
SupposeUin denotes the value of the solution ux,t of (1) at the mesh point xi,tn and let un=u0n,u1n,…,uMn be the solution of the compact difference scheme (28), assume that the condition in Theorem 1 is satisfied. Then, when τ≤1/42−α3, it holds(80)Un−un∗≤C2τ2+h4,1≤n≤N,where(81)C2=C2∗exp42−α2−α2+α2T1/2,C2∗=82−αLTCR21+22−α2T1−αΓ2−α.Proof.
Applying (79) and Theorem 3, we have(82)en∗2≤2−αGRexp42−α2−α2+α2T,1≤n≤N,where(83)GR=81+22−α2T1−αΓ2−ατRtxα1∗2+8τ∑k=2nRtxαk∗2.
Using Theorem 1, we obtain(84)en∗2≤C12τ2+h42.
The proof is completed.


## 5. Applications and Numerical Results

In this section, we apply the proposed compact finite difference method to two model problems in the form (1)–(3). Let Vin=vxi,tn be the value of the solution vx,t of problem (1)–(3) at the mesh point xi,tn. We have from (30) that(85)Vn−vnν≤C3τ2+h4,1≤n≤N,ν=1,2,∞,where C3 is a positive constant. In this paper, we compute norm errors as follows:(86)E2τ,h=max0≤n≤NVn−vn,Eντ,h=max0≤n≤NVn−vnν,ν=1,∞.The temporal and spatial convergence orders are computed, respectively, by(87)Oνtτ,h=log2Eν2τ,hEντ,h,Oνsτ,h=log2Eντ,2hEντ,h,ν=1,2,∞.Example 1.
We consider a problem which is governed by equation (1) in 0,π×0,1, and(88)fx,t=3t2cosx1+2t1−αΓ4−α+t3exsinx+cosx.
The boundary and initial conditions are given by (2) and (3) with(89)ϕ0t=3t21+2t1−αΓ4−α+t3ex,ϕLt=−3t21+2t1−αΓ4−α−t3ex,φx=ex.
It is easy to check thatvx,t=t3cosx is the solution of this problem.
At first, we let the spatial steph=π/400. Table 1 gives the errors Eντ,hν=1,2,∞ and the temporal convergence orders Oνtτ,hν=1,2,∞. Obviously, the computed solution vin has the second-order temporal accuracy, and this coincides well with the analysis. For comparison, in Table 2, we list the corresponding errors Eντ,hν=1,2,∞ and temporal convergence order Oνtτ,hν=1,2,∞ of the compact scheme (6) in [33]. This demonstrates that the scheme (6) in [33] has only the temporal accuracy of order 2−α.
Next, we compute the spatial convergence order of the compact difference scheme (28). Table 3 presents the errors Eντ,hν=1,2,∞ and the spatial convergence orders Oνsτ,hν=1,2,∞ with the time step τ=1/10000. As expected, the compact difference scheme (28) has the fourth-order spatial accuracy. Figure 1 shows that the difference scheme (28) is more accurate with small α.Table 1
The errors and the temporal convergence orders of the compact difference scheme (28) for Example 1h=π/400.
α
τ
E
1
τ
,
h
O
1
t
τ
,
h
E
2
τ
,
h
O
2
t
τ
,
h
E
∞
τ
,
h
O
∞
t
τ
,
h
1/4
1/10
2.022691e − 03
1.145804e −03
1.099575e −03
1/20
5.143506e − 04
1.975452
2.917805e −04
1.973404
2.797617e −04
1.974676
1/40
1.312214e − 04
1.970749
7.354273e −05
1.988229
7.049402e −05
1.988626
1/80
3.315323e − 05
1.984782
1.845280e −05
1.994743
1.768694e −05
1.994817
1/160
8.338128e − 06
1.991354
4.620517e −06
1.997713
4.428856e −06
1.997678
1/320
2.091440e − 06
1.995227
1.155860e −06
1.999088
1.107964e −06
1.999024
1/640
5.238263e − 07
1.997337
2.890005e −07
1.999823
2.770399e −07
1.999745
1/2
1/10
4.513663e − 03
2.896981e −03
2.571810e −03
1/20
1.136538e −03
1.989653
7.299431e −04
1.988694
6.480926e −04
1.988512
1/40
2.846515e −04
1.997376
1.829297e −04
1.996494
1.624428e −04
1.996268
1/80
7.114390e −05
2.000385
4.574314e −05
1.999663
4.062614e −05
1.999452
1/160
1.776905e −05
2.001374
1.142932e −05
2.000816
1.015216e −05
2.000621
1/320
4.437564e −06
2.001526
2.855128e −06
2.001112
2.536350e −06
2.000961
1/640
1.108323e −06
2.001389
7.132375e −07
2.001101
6.336576e −07
2.000979
3/4
1/10
6.789704e −03
4.178798e −03
3.709467e −03
1/20
1.698782e −03
1.998848
1.046757e −03
1.997162
9.289960e −04
1.997468
1/40
4.239228e −04
2.002628
2.614880e −04
2.001110
2.320296e −04
2.001363
1/80
1.056923e −04
2.003932
6.525416e −05
2.002603
5.789376e −05
2.002829
1/160
2.634707e −05
2.004155
1.627953e −05
2.003011
1.444159e −05
2.003179
1/320
6.568883e −06
2.003922
4.061570e −06
2.002949
3.602684e −06
2.003085
1/640
1.638209e −06
2.003529
1.013481e −06
2.002718
8.989079e −07
2.002827Table 2
The errors and the temporal convergence orders of the compact difference scheme (6) in [33] for Example 1h=π/400.
α
τ
E
1
τ
,
h
O
1
t
τ
,
h
E
2
τ
,
h
O
2
t
τ
,
h
E
∞
τ
,
h
O
∞
t
τ
,
h
1/4
1/10
2.2330679e −03
9.799808e −04
1.065741e −03
1/20
7.260649e −04
1.620857
3.186951e −04
1.620577
3.465575e −04
1.620690
1/40
2.310135e −04
1.652121
1.014135e −04
1.651926
1.102737e −04
1.652006
1/80
7.241011e −05
1.673714
3.179070e −05
1.673574
3.456677e −05
1.673631
1/160
2.245119e −05
1.689400
9.857621e −06
1.689293
1.071810e −05
1.689337
1/320
6.904356e −06
1.701213
3.031675e −06
1.701125
3.296226e −06
1.701162
1/640
2.109890e −06
1.710339
9.265034e −07
1.710247
1.007324e −06
1.710288
1/2
1/10
8.500612e −03
3.720403e −03
4.050738e −03
1/20
3.162060e −03
1.426702
1.384433e −03
1.426164
1.507136e −03
1.426376
1/40
1.155523e −03
1.452318
5.060407e −04
1.451970
5.508382e −04
1.452109
1/80
4.176415e −04
1.468209
1.829273e −04
1.467983
1.991085e −04
1.468074
1/160
1.498823e −04
1.478435
6.565537e −05
1.478285
7.146007e −05
1.478345
1/320
5.353812e −05
1.485192
2.345381e −05
1.485091
2.552668e −05
1.485132
1/640
1.906363e −05
1.489743
8.351745e −06
1.489673
9.089699e −06
1.489701
3/4
1/10
2.469894e −02
1.078660e −02
1.175650e −02
1/20
1.069987e −02
1.206856
4.674932e −03
1.206223
5.094444e −03
1.206462
1/40
4.574317e −03
1.225965
1.999121e −03
1.225580
2.178288e −03
1.225730
1/80
1.941433e −03
1.236434
8.486006e −04
1.236208
9.245974e −04
1.236298
1/160
8.206577e −04
1.242270
3.587423e −04
1.242138
3.908554e −04
1.242190
1/320
3.461070e −04
1.245563
1.513052e −04
1.245486
1.648459e −04
1.245517
1/640
1.457787e −04
1.247438
6.373102e −05
1.247394
6.943361e −05
1.247412Table 3
The errors and the spatial convergence orders of the compact difference scheme (28) for Example 1 τ=1/10000.
α
h
E
1
τ
,
h
O
1
s
τ
,
h
E
2
τ
,
h
O
2
s
τ
,
h
E
∞
τ
,
h
O
∞
s
τ
,
h
1/4
π
/
4
4.604088e – 03
4.604088e − 03
3.680647e − 03
π
/
8
2.841961e −04
3.990161
2.841961e − 04
4.017957
2.401639e − 04
3.937868
π
/
16
1.769697e −05
3.997913
1.769697e − 05
4.005312
1.496024e − 05
4.004814
π
/
32
1.103983e −06
4.000813
1.103983e − 06
4.002713
9.380629e − 07
3.995305
π
/
64
6.794511e −07
4.021394
6.794511e − 08
4.022204
5.769242e − 08
4.023231
1/2
π
/
4
6.370119e −03
4.408914e − 03
3.511266e − 03
π
/
8
4.015891e − 04
3.987528
2.722962e − 04
4.017175
2.319704e − 04
3.919979
π
/
16
2.514532e − 05
3.997359
1.695588e − 05
4.005319
1.444969e − 05
4.004830
π
/
32
1.569810e − 06
4.001628
1.056533e − 06
4.004376
9.029468e − 07
4.000254
π
/
64
9.567678e − 08
4.036277
6.381205e − 08
4.049365
5.493287e − 08
4.038900
3/4
π
/
4
6.074741e − 03
4.175553e − 03
3.306275e − 03
π
/
8
3.839021e − 04
3.984013
2.580707e − 04
4.016129
2.220191e − 04
3.896451
π
/
16
2.405107e − 05
3.996565
1.607074e − 05
4.005258
1.382999e − 05
4.004812
π
/
32
1.501251e − 06
4.001864
1.000621e − 06
4.005469
8.619179e − 07
4.004106
π
/
64
9.116310e − 08
4.041571
5.969859e − 08
4.067055
5.211927e − 08
4.047662Figure 1
Error planes forτ=1/20andh=π/100 for different α. (a) α=0.3. (b) α=0.5. (c) α=0.7. (d) α=0.9.
(a)
(b)
(c)
(d)Example 2.
Consider equation (1) in the domain 0,π×0,1, and(90)fx,t=3t21+x21+2t1−αΓ4−α−t36x2+4x+2.
The boundary and initial conditions are given by (2) and (3) with(91)ϕ0t=3t21+2t1−αΓ4−α−2t3,ϕLt=12t21+2t1−αΓ4−α−12t3,φx=1+x2.
It is clear thatvx,t=t31+x2 is the exact analytical solution of this problem.
In Table4, we give the errors Eντ,hν=1,2,∞ and the temporal convergence orders Oνtτ,hν=1,2,∞ with the spatial step h=1/200. Obviously, the computed solution vin has the second-order temporal accuracy.
Table5 presents the errors Eντ,hν=1,2,∞ and the spatial convergence orders Oνsτ,hν=1,2,∞ with the time step τ=1/15000. These results show that the compact difference scheme (28) has the fourth-order spatial accuracy. Figure 2 presents the probability density function ux,t and describes the anomalous subdiffusive particle versus the anomalous diffusion exponent α.Table 4
The errors and the temporal convergence orders of the compact difference scheme (28) for Example 2h=π/400.
α
τ
E
1
τ
,
h
O
1
t
τ
,
h
E
2
τ
,
h
O
2
t
τ
,
h
E
∞
τ
,
h
O
∞
t
τ
,
h
1/4
1/10
6.204449e − 02
4.366942e − 02
3.382734e − 02
1/20
1.574638e − 02
1.978283
1.108138e − 02
1.978486
8.583688e − 03
1.978520
1/40
3.968730e − 03
1.988271
2.792763e − 03
1.988372
2.163275e − 03
1.988380
1/80
9.966436e − 04
1.993528
7.013047e − 04
1.993580
5.432299e − 04
1.993582
1/160
2.497870e − 04
1.996379
1.757634e − 04
1.996407
1.361461e − 04
1.996408
1/320
6.253527e − 05
1.997956
4.400270e − 05
1.997971
3.408443e − 05
1.997970
1/640
1.564637e − 05
1.998842
1.100945e − 05
1.998850
8.527917e − 06
1.998848
1/2
1/10
6.018456e − 02
4.245709e − 02
3.295475e − 02
1/20
1.534351e − 02
1.971765
1.082348e − 02
1.971842
8.398932e − 03
1.972209
1/40
3.881225e − 03
1.983045
2.737771e − 03
1.983090
2.124237e − 03
1.983261
1/80
9.775209e − 04
1.989313
6.895183e − 04
1.989341
5.349612e − 04
1.989439
1/160
2.455614e − 04
1.993044
1.732106e − 04
1.993062
1.343795e − 04
1.993121
1/320
6.158773e − 05
1.995369
4.344152e − 05
1.995381
3.370175e − 05
1.995418
1/640
1.543038e − 05
1.996869
1.088391e − 05
1.996877
8.443565e − 06
1.996899
3/4
1/10
5.263247e − 02
3.710102e − 02
2.900228e − 02
1/20
1.340236e − 02
1.973466
9.448768e − 03
1.973261
7.382180e − 03
1.974047
1/40
3.394668e − 03
1.981145
2.393477e − 03
1.981018
1.869338e − 03
1.981520
1/80
8.571578e − 04
1.985638
6.043931e − 04
1.985549
4.719259e − 04
1.985895
1/160
2.159868e − 04
1.988618
1.523021e − 04
1.988551
1.188994e − 04
1.988818
1/320
5.434321e − 05
1.990771
3.832125e − 05
1.990719
2.991224e − 05
1.990934
1/640
1.365739e − 05
1.992418
9.631070e − 06
1.992377
7.516772e − 06
1.992551Table 5
The errors and the spatial convergence orders of the compact difference scheme (28) for Example 2τ=1/15000.
α
h
E
1
τ
,
h
O
1
s
τ
,
h
E
2
τ
,
h
O
2
s
τ
,
h
E
∞
τ
,
h
O
∞
s
τ
,
h
1/4
π/4
1.661627e − 02
8.933632e − 03
9.781141e − 03
π/8
1.203317e − 03
3.787508
5.667278e − 04
3.978519
5.872186e − 04
4.058033
π/16
7.447394e − 05
4.014137
3.441625e − 05
4.041494
4.064209e − 05
3.852851
π/32
4.673315e − 06
3.994218
2.125617e − 06
4.017136
2.512751e − 06
4.015635
π/64
2.797367e − 07
4.062305
1.199961e − 07
4.146822
1.472003e − 07
4.093415
1/2
π/4
1.571731e − 02
8.392759e − 03
9.222668e − 03
π/8
1.142756e − 03
3.781765
5.330938e − 04
3.976684
5.533294e − 04
4.058974
π/16
7.080945e − 05
4.012432
3.235842e − 05
4.042176
3.857001e − 05
3.842587
π/32
4.448701e − 06
3.992486
1.998328e − 06
4.017276
2.384765e − 06
4.015561
π/64
2.675409e − 07
4.055553
1.129747e − 07
4.144721
1.399255e − 07
4.091116
3/4
π/4
1.464867e − 02
7.752890e − 03
8.557413e − 03
π/8
1.070539e − 03
3.774361
4.932762e − 04
3.974266
5.164638e − 04
4.050436
π/16
6.645084e − 05
4.009905
2.992427e − 05
4.043008
3.609436e − 05
3.838822
π/32
4.182715e − 06
3.989776
1.848879e − 06
4.016593
2.232553e − 06
4.015007
π/64
2.537283e − 07
4.043083
1.056609e − 07
4.129138
1.319669e − 07
4.080447Figure 2
Numerical solutions for (28) when τ=1/20andh=π/40 for different α. (a) α=0.3. (b) α=0.5. (c) α=0.7. (d) α=0.9.
(a)
(b)
(c)
(d)


## 6. Conclusions

In this paper, theL2-1σ approximation formula is applied for the Caputo time-fractional derivative and a fourth-order compact finite difference approximation is used for the spatial derivative. The proposed method extends the known method in [29] to a class of more general time-fractional subdiffusion equations with spatially variable coefficients. We have proved that the proposed scheme is uniquely solvable, unconditionally stable, and convergent. We have also provided the optimal error estimates in the discrete H1, L2, and L∞ norms. The error estimates show that the proposed method has the second-order temporal accuracy and the fourth-order spatial accuracy, and thus it improves the temporal accuracy of the method given in [33]. Numerical results confirm our analysis and demonstrate the efficiency of our method. Applying the proposed method to the time multiterm fractional subdiffusion equations will be our further work.



---

*Source: 4582020-2019-10-20.xml*