# A Fully Discrete Galerkin Method for a Nonlinear Space-Fractional Diffusion Equation

Authors: Yunying Zheng; Zhengang Zhao
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

The spatial transport process in fractal media is generally anomalous. The space-fractional advection-diffusion equation can be used to characterize such a process. In this paper, a fully discrete scheme is given for a type of nonlinear space-fractional anomalous advection-diffusion equation. In the spatial direction, we use the finite element method, and in the temporal direction, we use the modified Crank-Nicolson approximation. Here the fractional derivative indicates the Caputo derivative. The error estimate for the fully discrete scheme is derived. And the numerical examples are also included which are in line with the theoretical analysis.


---

## Body

## 1. Introduction

The normal diffusive motion is modeled to describe the standard Brownian motion. The relation between the flow and the divergence of the particle displacement represents(1.1)J(x,t)=-a∂c∂x+bc,
where J is the diffusive flow. Inserting the above equation into the equation of mass conservation (1.2)∂J∂x=-∂c∂t,
we obtain the standard convection-diffusion equation. From the viewpoint of physics, it means that during the method of time random walkers, the overall particle displacement up to time t can be represented as a sum of independent random steps, in the case that both the mean-squared displacement per step and the mean time needed to perform a step are finite. The measured variance growth in the direction of flow of tracer plumes is typically at a Fickian rate, 〈(c-c¯)2〉~t.The transport process in fractal media cannot be described with the normal diffusion. The process is nonlocal and it does not follow the classical Fickian law. It depicts a particle in spreading tracer cloud which has a standard deviation, and which grows liket2α for some 0<α<1, excluding the Fickian case α=1/2. The description of anomalous diffusion means that the measure variance growth in the direction of flow has a deviation from the Fickian case, it follows the super-Fickian rate 〈(c-c¯)2〉~t2α when α>1/2, or does the subdiffusion rate 〈(c-c¯)2〉~t2α if 0<α<1/2. With the help of the continuous time random walk and the Fourier transform, the governing equation with space fractional derivative can be derived as follows (1.3)∂u∂t=D(a(u)aDxβu)+b(u)Du+f(x,t,u),0<β<1,
where D denotes integer derivative respect to x, and Dβ is fractional derivative. There are some authors studying the spacial anomalous diffusion equation in theoretical analysis and numerical simulations [1–10]. Now the fractional anomalous diffusion becomes a hot topic because of its widely applications in the evolution of various dynamical systems under the influence of stochastic forces. For example, it is a well-suited tool for the description of anomalous transport processes in both absence and presence of external velocities or force fields. Since the groundwater velocities span many orders of magnitude and give rise to diffusion-like dispersion (a term that combines molecular diffusion and hydrodynamic dispersion), the fractional diffusion is an important process in hydrogeology. It can be used to describe the systems with reactions and diffusions across a wide range of applications including nerve cell signaling, animal coat patterns, population dispersal, and chemical waves. In general, fractional anomalous diffusions have numerous applications in statistical physics, biophysics, chemistry, hydrogeology, and biology [4, 11–20].In this paper, we mainly study one kind of typical nonlinear space-fractional partial differential equations by using the finite element method, which reads in the following form:(1.4)∂u∂t=D(a(u)aDxβu)+b(u)Du+f(x,t,u),x∈Ω,t∈(0,T],u∣t=0=φ(x),x∈Ω,u∣∂Ω=g,t∈(0,T],
where Ω is a spacial domain with boundary ∂Ω, Dβ is the βth (0<β<1) order fractional derivative with respect to the space variable x in the Caputo sense (which will be introduced later on),a,b,f are functions ofx,t,u,φ and g are known functions which satisfy the conditions requested by the theorem of error estimations.The rest of this paper is constructed as follows. In Section2 the fractional integral, fractional derivative, and the fractional derivative spaces are introduced. The error estimates of the finite element approximation for (1.4) are studied in Section 3, and in Section 4, numerical examples are taken to verify the theoretical results derived in Section 3.


## 2. Fractional Derivative Space

In this section, we firstly introduce the fractional integral (or Riemann-Liouville integral), the Caputo fractional derivative, and their corresponding fractional derivative space.Definition 2.1.
Theαth order left and right Riemann-Liouville integrals of function u(x) are defined as follows
(2.1)aIxαu(x)=1Γ(α)∫ax(x-s)α-1u(s)ds,xIbαu(x)=1Γ(α)∫xb(s-x)α-1u(s)ds,
where α>0, and Γ(·) is the Gamma function.Definition 2.2.
Theαth order Caputo derivative of function u(x) is defined as,
(2.2)aDxαu(x)=aIxn-αdnu(x)dxn,n-1<α<n∈Z+,xDbαu(x)=(-1)nxIbn-αdnu(x)dxn,n-1<α<n∈Z+.
The αth order Riemann-Liouville derivative of function u(x) is defined by changing the order of integration and differentiation.Lemma 2.3 (see [8]).
Ifu(0)=u′(0)=⋯=u(n-1)(0)=0, then the Caputo fractional derivative is equal to the Riemann-Liouville derivative.Definition 2.4.
The fractional derivative spaceJα(Ω) is defined as follows:
(2.3)Jα(Ω)={u∈L2(Ω):Daxαu∈L2(Ω),n-1≤α<n},
endowed with the seminorm
(2.4)|u|Jα=‖Daxαu‖L2(Ω),
and the norm
(2.5)‖u‖Jα=(|u|Jα2+∑k≤[α]‖Dku‖2)1/2.LetJ0α(Ω) denote the closure of C0∞(Ω) with respect to the above norm and seminorm.Definition 2.5.
Define the seminorm(2.6)|u|Hα=‖|iw|αF(u)‖L2(Ω),
and the norm
(2.7)‖u‖Hα=(|u|Hα2+∑k≤[α]‖Dku‖2)1/2,
where i is the imaginary unit, and F is the Fourier transform, and which can define another fractional derivative space Hα(Ω).LetH0α(Ω) denote the closure of C0∞(Ω) with respect to the norm and seminorm.Definition 2.6.
The fractional spaceJsα(Ω) is defined below
(2.8)Jsα(Ω)={u∈L2(Ω):Daxαu∈L2(Ω),xDbαu∈L2(Ω),n-1≤α<n},
endowed with the seminorm
(2.9)|u|Jsα=|(Daxαu,Dxbαu)1/2|L2(Ω),
and the norm
(2.10)‖u‖Jsα=(∑k≤[α]‖Dku‖2+|u|Jsα2)1/2.Theorem 2.7 (see [3, 6]).
Jsα, Jα, and Hα are equal with equivalent seminorm and norm.
The following are some useful results.Lemma 2.8 (see [3]).
Foru∈J0α(Ω), 0<β<α, then
(2.11)Daxαu(x)=Daxα-βDaxβu.Lemma 2.9 (see [2]).
Foru∈H0α(Ω), one has
(2.12)‖u‖L2(Ω)≤c|u|H0α.
For 0<β<α,
(2.13)|u|H0β(Ω)≤c|u|H0α.SinceJsα, Jα, and Hα are equal with equivalent seminorm and norm, the norms with each space which will be used following are without distinction, and the notations are used seminorm |·|α and norm ∥·∥α.


## 3. Finite Element Approximation

LetΩ=[a,b], and 0≤β<1. Define α=(1+β)/2. In this section, we will formulate a fully discrete Galerkin finite element method for a type of nonlinear anomalous diffusion equation as follows.Problem 1 (Nonlinear spacial anomalous diffusion equation).
We  consider equations of the form(3.1)∂u∂t=D(a(u)aDxβu)+b(u)Du+f(x,t,u),(x,t)∈Ω×(0,T],u(x,t)=ϕ(x,t),x∈∂Ω×(0,T],u(x,0)=g(x),x∈Ω¯.
We always assume that
(3.2)0<m<a(u)<M,0<m<b(u)<M,0<m<f(u)<M.The algorithm and analysis in this paper are applicable for a large class of linear and nonlinear functions (including polynomials and exponentials) in the unknown variables. Throughout the paper, we assume the following mild Lipschitz continuity conditions ona,b,  and  f: there exist positive constants L and c such that for x∈Ω,t∈(0,T], and s,t∈R, (3.3)|a(x,t,s)-a(x,t,r)|≤L|s-r|,(3.4)|b(x,t,s)-b(x,t,r)|≤L|s-r|,(3.5)|f(x,t,s)-f(x,t,r)|≤L|s-r|.In order to derive a variational form of Problem1, we suppose that u is a sufficiently smooth solution of Problem 1. Multiplying an arbitrary v∈H0α(Ω) in both sides yields (3.6)∫Ω∂u∂tvdx=∫ΩD(a(u)aDxβu)vdx+∫Ωb(u)Duvdx+∫Ωf(x,t,u)vdx.Rewriting the above expression yields(3.7)∫Ω∂u∂tvdx+∫Ωa(u)aDxβuDvdx-∫Ωb(u)Duvdx=∫Ωf(x,t,u)vdx.We define the associated bilinear formA:J0α(Ω)×J0α(Ω)→R as (3.8)A(u,v)=(a(u)aDxβu,Dv)-(b(u)Du,v),
where (·,·) denotes the inner product on L2(Ω) and J0α(Ω).For givenf∈J-α(Ω), we define the associated function F:J0α(Ω)→R as (3.9)F(v)=〈f,v〉.Definition 3.1.
A functionu∈J0α(Ω) is a variational solution of Problem 1 provided that
(3.10)(∂u∂t,v)+A(u,v)=F(v),∀v∈J0α(Ω).Now we are ready to describe a fully discrete Galerkin finite element method to solve nonlinear Problem1. In our new scheme, the finite element trial and test spaces for Problem 1 are chosen to be same.For a positive integerN, let ∏t={tn}n=0N be a uniform partition of the time interval (0,T] such that tn=nτ, where τ=T/N, and let tn-1/2=tn-τ/2. Throughout the paper, we use the following notation for a function ϕ: (3.11)ϕn=ϕ(tn),∂¯tϕn=ϕn-ϕn-1τ,ϕ¯n=ϕn+ϕn-12,ϕ̃n=3ϕn-1-ϕn-22.Let𝒦h={K} be a partition of spatial domain Ω. Define hk as the diameter of the element K and h=maxK∈𝒦hhK. And let Sh be a finite element space (3.12)Sh={v∈H0α(Ω):v∣K∈Pr-1(K),K∈Kh},
where Pr-1(K) is the set of polynomials of degree r-1 on a given domain K. And the functions in Sh are continuous on Ω. Our fully discrete quadrature scheme to solve Problem 1 is to find uh: for v∈Sh such that (3.13)(∂¯tuhn,v)+(a(ũhn)aDxβu¯hn,Dv)-(b(ũhn)Du¯hn,v)=〈f(ũhn),v〉.The linear systems in the above equation requires selecting the value ofuh0 and uh1. Given uh0 depending on the initial data g(x), we select uh1 by solving the following predictor-corrector linear systems: (3.14)(uh1,0-uh0τ,v)+(a(uh0)aDxβuh1,0+uh02,Dv)-(b(uh0)Duh1,0+uh02,v)=〈f(uh0),v〉,(uh1-uh0τ,v)+(a(uh1,0+uh02)Daxβuh1+uh02,Dv)-(b(uh1,0+uh02)Duh1+uh02,v)=〈f(uh1,0+uh02),v〉.Lemma 3.2.
Foru,v,w∈Js,0α(Ω),0<m≤a(u)≤M,α=(1+β)/2, there exist constants γ1,γ2 such that
(3.15)(a(u)aDxβu,Dv)≤γ1‖u‖α⋅‖v‖α,(a(w)aDxβv,Dv)≥γ2‖v‖α2.Proof.
With the assumption ofa(u) in (3.3) and the property of dual space
(3.16)(a(w)aDxβu,Dv)≤‖a(w)aDxβu‖1-α⋅‖Dv‖-(1-α)≤Mc‖u‖1-α+β⋅‖v‖-(1-α)+1≤γ1‖u‖α⋅‖v‖α,(a(w)aDxβv,Dv)=-(Da(w)aDxβv,v)=-(Dax(1-β)/2a(w)aDxβv,Dxb(1+β)/2v)≥m|v|Jsα2≥γ2‖v‖α2.Lemma 3.3 (see [2]).
ForΩ⊂Rn,α>n/4,v,w∈H0α(Ω),ε>0, one has
(3.17)(vb(w),∇v)≤c0(qε)-p/qp‖∇b(w)‖p⋅‖v‖2+ε‖v‖α2,
where p=4α/(4α-n),q=4α/n.Theorem 3.4.
Letuhn be bounded, then for a sufficiently small step τ, there exists a unique solution uhn∈Sh satisfying scheme (3.13).Proof.
As scheme represents a finite system of problem, the continuity and coercivity of(u¯hn,ω¯hn)/τ+A(u¯hn,ω¯hn) is the sufficient and essential condition for the existence and uniqueness of uhn. Let v=u¯hn,w=ω¯hn, then
(3.18)(v,v)τ+A(v,v)=(v,v)τ+(a(w)aDxβv,Dv)-(b(w)Dv,v)≥‖v‖2τ+γ2‖v‖α-c0‖Db(w)‖2‖v‖2-ε‖v‖α2=(γ2-ε)‖v‖α2+(τ-1-c0‖Db(w)‖2)‖v‖2≥c‖v‖α2.
For the chosen sufficiently small τ, the above inequality holds.
(3.19)(v,w)τ+A(v,w)=(v,w)τ+(a(u)aDxβv,Dw)+(Db(u)v,Dw)≤‖u‖⋅‖w‖τ+γ1‖v‖α‖w‖α+‖v‖⋅‖D(b(u)w)‖≤‖u‖⋅‖w‖τ+γ1‖v‖α‖w‖α+M‖v‖⋅‖w‖h≤c‖v‖α‖w‖α.
Hence, the scheme (3.13) is uniquely solvable for uhn.
Letρn=Phun-un, and θn=uhn-Phun, then
(3.20)uhn-un=uhn-Phun+Phun-un=θn+ρn,
where Phun is a Rits-Galerkin projection operator defined as follows:
(3.21)(a(w)aDxβ(un-Phun),Dv)=0,(a(u0)aDxβ(un-Phun),Dv)=0.Lemma 3.5.
Leta(u),b(u) be smooth functions on Ω, 0<m≤a(u),b(u)≤M, and Phun is defined as above, then
(3.22)‖Daxα(un-Phun)‖≤chk+1-α‖u‖k+1,∥(Phun-un)∥≤chk+1‖u‖k+1.Proof.
Using the definition ofPhun, one gets
(3.23)‖Daxα(Phun-un)‖2=|(Daxα(Phun-un),Daxα(Phun-un))|≤c‖Daxα(Phun-un)‖⋅‖Daxα(χ-un)‖,
where χ∈Sh. Utilizing the interpolation of Ihun leads to
(3.24)‖Daxα(Phun-un)‖≤inf⁡χ∈Shc‖χ-u‖α≤c‖Ihun-un‖α≤chk+1-α‖u‖k+1.
Next we estimate ∥Phun-un∥. For allϕ∈L2(Ω), w is the solution of the following equation:
(3.25)-aDx2αw=ϕ,w∈Ω,w=0,w∈∂Ω.
So we have(3.26)‖w‖2α≤γ3‖ϕ‖.
For allχ∈Sh, with the help of approximation properties of Sh and the weak form, we can obtain
(3.27)(Phun-un,ϕ)=-(Phun-un,Dax2αw)=-(Dxbα(Phun-un),Daxαw)=-(Dxbα(Phun-un),Daxα(w-χ))≤‖Phun-un‖α‖w-χ‖α≤‖Phun-un‖αinf⁡χ∈Sh‖w-χ‖α≤chr-α‖u‖rhα‖w‖2α=chr‖u‖r‖ϕ‖,‖Phun-un‖=sup⁡0≠ϕ∈L2(Ω)(Phun-un,ϕ)‖ϕ‖≤chr‖u‖r.Lemma 3.6 (see [21]).
LetTh,0<h≤1, denote a quasiuniform family of subdivisions of a polyhedral domain Ω⊂Rd. Let (K′,P,N) be a reference finite element such that P⊂Wl,p(K′)∩Wm,q(K′) is a finite-dimensional space of functions on K′,N is a basis for P′, where 1≤p≤∞,1≤p≤∞, and 0≤m≤l. For K∈Th, let (K,PK,NK) be the affine equivalent element, and Vh=v:v is measurable and v|K∈PK,forallK∈Th. Then there exists a constant C=C(l,p,q) such that
(3.28)[∑k∈Th‖v‖Wl,p(K)2]1/p≤Chm-l+min⁡(0,d/p-d/q)⋅[∑k∈Th‖v‖Wm,q(K)q]1/q.
The following Gronwall’s lemma is useful for the error analysis later on.Lemma 3.7 (see [2]).
LetΔt,H and an,bn,cn,γn (for integer n≥0) be nonnegative numbers such that
(3.29)aN+Δt∑n=0Nbn≤Δt∑n=0Nγnan+Δt∑n=0Ncn+H,
for N≥0. Suppose that Δtγn<1, for all n, and set σn=(1-Δtγn)-1. Then
(3.30)aN+Δt∑n=0Nbn≤exp⁡(Δt∑n=0Nσnγn){Δt∑n=0Ncn+H},
for N≥0.
The following norms are also used in the analysis:(3.31)‖|v|‖∞,k=max⁡0≤n≤N‖vn‖k,‖|v|‖0,k=[∑n=0Nτ‖vn‖k2]1/2.Theorem 3.8.
Assume that Problem1 has a solution u satisfying utt,uttt∈L2(0,T,L2(Ω)) with u,ut∈L2(0,T,Hk+1). If Δt≤ch, then the finite element approximation is convergent to the solution of Problem 1 on the interval (0,T], as Δt,h→0. The approximation uh also satisfies the following error estimates
(3.32)‖u-uh‖0,α≤C(hk+1‖ut‖0,k+1+hk+1-α‖u‖0,k+1+τ2‖utt‖0,0+τhk+1-α‖utt‖0,k+1+τ2‖uttt‖0,0),(3.33)‖u-uh‖∞,0≤C(hk+1‖ut‖0,k+1+hk+1-α‖u‖0,k+1+τ2‖uttt‖0,0+τhk+1-α‖utt‖0,k+1+τ2‖utt‖0,0+hk+1‖u‖∞,k+12).Proof.
Fort=tn-τ/2=tn-1/2,n=0,1,…,N, find un-1/2 such that
(3.34)(∂tun-1/2,v)+(a(un-1/2)Daxun-1/2,Dv)-(b(un-1/2)Dun-1/2,v)=〈f(un-1/2),v〉.
Subtracting the above equation from the fully discrete scheme (3.13), and substituting uhn-un=(uhn-Phun)+(Phun-un)=θn+ρn into it, we obtain the following error formulation relating to θn and ρn:
(3.35)(∂¯tθn,v)+(a(ũhn)Daxβθ¯n,Dv)-(b(ũhn)Dθn¯,v)=(a(ũhn)DaxβIhu¯n,Dv)+(b(ũhn)DaxβIhu¯n,v)+(∂tun-1/2,v)-(∂¯tIhun,v)+(a(un-1/2)Daxβun-1/2,Dv)-(b(un-1/2)Dun-1/2,v)+(f(ũhn),v)-(f(un-1/2),v)=-(a(ũhn)Daxβρ¯n,Dv)+{(a(un-1/2)Daxβun-1/2-a(ũhn)Daxβun¯,Dv)}+(b(ũhn)Dρ¯n,v)+{(b(ũhn)Du¯n-(b(un-1/2)Dun-1/2,Dv)}+{(f(ũhn)-f(un-1/2),v)}+{(∂tun-1/2-∂¯tIhun,v)}=R1(v)+R2(v)+R3(v)+R4(v)+R5(v)+R6(v).
Setting v=θn¯, we obtain
(3.36)(∂¯tθn,θn¯)+(a(ũn)aDxβθ¯n,Dθ¯n)-(b(ũn)Dθ¯n,θ¯n)=R1(θ¯n)+R2(θ¯n)+R3(θ¯n)+R4(θ¯n)+R5(θ¯n)+R6(θ¯n).
Note that
(3.37)(∂¯tθn,θn¯)=(θn-θn-1τ,θn+θn-12)=12τ(‖θn‖2-‖θn-1‖2).
According to (3.2) and Lemma 3.2, we have
(3.38)(a(ũn)aDxβθ¯n,Dθ¯n)≥m|θ¯n|α2≥c(|θn|α2+|θn-1|α2).
From Lemma 3.3, the following inequality can be derived:
(3.39)(b(ũn)θ¯n,Dθ¯n)≤c0ε2-c1‖Db(ũn)‖c2‖θ¯n‖2+ε3‖θ¯n‖α2=c0ε2-c1‖Db(ũn)‖c2‖θn+θn-12‖2+ε3‖θn+θn-12‖α2≤c3ε2-c1‖Db(ũn)‖c2(‖θn‖2+‖θn-1‖2)+c4ε3(‖θn‖α2+‖θn-1‖α2).
Substituting (3.37)–(3.39) into (3.36) then multiplying (3.36) by 2τ, summing from n=1 to N, we have
(3.40)‖θn‖2-‖θ2‖2+τ∑n=1N(2mc-2c4ε3)(‖θn‖α2+‖θn-1‖α2)≤2τ∑n=1Nc3ε2-c1‖Db(ũhn)‖c2(‖θn‖2+‖θn-1‖2)+2τ∑n=3N[R1(θ¯n)+R2(θ¯n)+R3(θ¯n)+R4(θ¯n)+R5(θ¯n)+R6(θ¯n)].
We now estimate R1 to R6 in the right hand of (3.40),
(3.41)R1(θ¯n)=(Dax1-α(a(ũhn)Daxβρ¯n),Daxαθ¯n)≤M(Daxαρ¯n,Daxαθ¯n)≤M‖Daxαρ¯n‖‖Daxαθ¯n‖≤ε4‖θ¯n‖α2+c524ε4‖ρ¯n‖α2=ε42‖ρn+ρn-1‖α2+c5216ε4‖θn+θn-1‖α2≤ε4c6(‖θn‖α2+‖θn-1‖α2)+c7ε4(‖ρn‖α2+‖ρn-1‖α2).
Secondly, we deduce the estimation ofR2,
(3.42)R2(θ¯n)=(-a(ũhn)Daxβu¯n,Dθ¯n)+(a(un-1/2)Daxβun-1/2,Dθ¯n)=((a(un-1/2)-a(ũhn))Daxβu¯n,Dθ¯n)+(a(un-1/2)(Daxβun-1/2-aDxβu¯n)),Dθ¯n)=R21+R22,
where
(3.43)R21=([a(un-1/2)-a(ũn)]Daxβũn),Dθ¯n)≤c84ε5‖[a(un-1/2)-a(ũn)]Daxβũn‖1-α2+ε5‖Dθ¯n‖α-12≤c9‖a(un-1/2)-a(ũn)‖‖Daxβũn‖1-α2+ε5‖θ¯n‖α2≤c9L‖un-1/2-ũn‖+ε5‖θ¯n‖α2,R22=(a(un-1/2)(Daxβun-1/2-aDxβu¯n),Dθ¯n)≤c104ε6‖a(un-1/2)[Daxβun-1/2-aDxβu¯n]‖1-α2+ε6‖Dθ¯n‖-1+α2≤c10‖a(un-1/2)‖2‖Daxβun-1/2-aDxβu¯n‖1-α2+ε6‖θ¯n‖α2≤c10M2‖un-1/2-u¯n‖α2+ε6c(‖θn‖α2+‖θn-1‖α2).
The estimations of ∥ũn-un-1/2∥ and ∥u¯n-un-1/2∥α can be derived as follows:
(3.44)‖ũn-un-1/2‖=‖32[un-1/2-τ2utn-1/2+uttn-1/22!(τ2)2+O(τ3)]-12[un-1/2-3τ2utn-1/2+uttn-1/22!(3τ2)2+O(τ3)]-un-1/2‖≤c11τ2‖utt(tn-1/2)‖≤c11τ2∫tn-1tn‖utt(⋅,s)‖ds,‖u¯n-un-1/2‖α=‖τ-1{∫tn-1/2tn(s-tn)2utt(s)ds+∫tn-1tn-1/2(s-tn-1)2utt(s)ds}‖α≤c12τ‖∫tn-1tnutt(s)ds‖α≤c12τ∫tn-1tn‖utt(s)‖αds≤c12τhk+1-α∫tn-1tn‖utt(s)‖k+1ds.
Thirdly, it is turn to consider R3,
(3.45)R3(θ¯n)=(b(ũhn)Dρ¯n,θ¯n)≤‖b(ũhn)Dρ¯n‖-α‖θ¯n‖α≤c134ε7‖b(ũhn)‖2‖ρ¯n‖1-α2+ε7‖θ¯n‖α2≤c144ε7(‖ρn‖1-α2+‖ρn-1‖1-α2)+ε7c15(‖θn‖α2+‖θn-1‖α2).
Next,
(3.46)R4(θ¯n)=(b(ũhn)Du¯n,θ¯n)-(b(un-1/2)Dun-1/2,θ¯n)=((b(ũhn)Du¯n-b(un-1/2)Du¯n),θ¯n)+(b(un-1/2)(Du¯n-Dun-1/2),θ¯n)=R41+R42,
where
(3.47)R41≤c164ε8‖{b(ũhn)-b(un-1/2)}Du¯n‖1-α2+ε8‖θ¯n‖α2≤c16L4ε8‖ũhn-un-1/2‖2‖Du¯n‖1-α2+ε8‖θ¯n‖α2=c16L4ε8‖ũhn-ũn+ũn-un-1/2‖2‖Du¯n‖1-α2+ε8‖θ¯n‖α2≤c17‖ũhn-ũn‖2+c17‖ũn-un-1/2‖2+ε8‖θ¯n‖α2≤c17‖θ̃n+ρ̃n‖2+c17‖ũn-un-1/2‖2+ε8‖θ¯n‖α2≤c18(‖θ̃n‖2+‖ρ̃n‖2)+c17‖ũn-un-1/2‖2+ε8‖θ¯n‖α2.
Rewriting R42 by the aid of (3.20), we have
(3.48)R42≤c194ε9‖u¯n-un-1/2‖2+ε9‖θ¯n‖α2.
The estimation of R5 is deduced as follows:
(3.49)R5(θn¯)≤‖f(ũhn)-f(un-1/2)‖‖θ¯n‖≤L‖ũhn-un-1/2‖‖θ¯n‖≤Lc204ε10‖ũhn-un-1/2‖2+ε10‖θ¯n‖2≤Lc21(‖θ̃n+ρ̃n‖2+‖ũn-un-1/2‖2)+ε10‖θ¯n‖2≤c22(‖θ̃n‖2+‖ρ̃n‖2)+Lc21‖ũn-un-1/2‖2+ε10‖θ¯n‖2.
Last, we estimate R6,
(3.50)R6(θ¯n)=(∂tun-1/2,θ¯n)-(∂¯tPhun,θ¯n)=(∂tun-1/2-∂¯un,θ¯n)+(∂¯tun-∂¯tPhun,θ¯n)=(∂tun-1/2-∂¯tun,θ¯n)+(∂¯tρn,θ¯n)≤‖∂tun-1/2-∂¯tun‖‖θn‖+‖∂¯tρn‖‖θn‖,
where
(3.51)‖∂tun-1/2-∂¯un‖=(2τ)-1c23‖∫tn-1/2tn(s-tn)2uttt(s)ds+∫tn-1tn-1/2(s-tn-1)2uttt(s)ds‖≤c23τ‖∫tn-1tnuttt(s)ds‖≤c23τ∫tn-1tn‖uttt(s)‖ds,‖∂¯tρn‖=‖ρn-ρn-1τ‖≤τ-1‖∫tn-1tnρtn(s)ds‖≤τ-1∫tn-1tn‖ut(s)‖ds≤τ-1∫tn-1tn1ds∫tn-1tn‖ut(s)‖ds=∫tn-1tn‖ut(s)‖≤hk+1∫tn-1tn‖ut(s)‖k+1ds.
The∥θ2∥ should be estimated with (3.14). Let n=1 then subtracting (3.34) from the two equations of (3.14), respectively, one gets
(3.52)(∂¯tθ1,0,v)+(a(uh0)Daxβθ¯1,0,Dv)-b((uh0)Dθ¯1,0,v)=-(a(uh0)Daxβρ¯1,0,Dv)-{(a(u1/2)Daxβu1/2-a(uh0)Daxβu¯1,0,Dv)}+(b(uh0)Dρ¯1,0,v)+{(b(uh0)Du¯1,0-b(u1/2)Du1/2,Dv)}+{(f(uh0)-f(u1/2),v)}+{(∂tu1/2-∂¯tIhu1,0,v)}=R1(v)+R2(v)+R3(v)+R4(v)+R5(v)+R6(v),(∂¯tθ1,v)+(a(uh0)aDxβθ¯1,Dv)-b((uh0)Dθ¯1,v)=-(a(uh0+uh1,02)Daxβρ¯1,Dv)-{(a(u1/2)Daxβu1/2-a(uh0+uh1,02)Daxβu¯1,Dv)}+(b(uh0+uh1,02)Dρ¯1,v)+{(b(uh0+uh1,02)Du¯1-(b(u1/2)Du1/2,Dv)}+{(f(uh0+uh1,02)-f(u1/2),v)}+{(∂tu1/2-∂¯tIhu1,v)}=R1(v)+R2(v)+R3(v)+R4(v)+R5(v)+R6(v).
Settingv=θ1,0¯, and using the similar estimation (see (3.40)), one has
(3.53)‖θ1,0‖2≤c{τ2∫t0t1‖utt(s)‖α2ds+h2(k+1-α)‖u‖k+12+τ2∫t0t1‖uttt(s)‖2ds+ch2(k+1)∫t0t1‖ut(s)‖k+12ds}.
Lettingv=θ1¯, applying the above result of θ1,0, and using the similar estimation (see (3.53)), we get
(3.54)‖θ1‖2≤c{τ2∫t0t1‖utt(s)‖α2ds+h2(k+1-α)‖u‖k+12+τ2∫t0t1‖uttt(s)‖2ds+ch2(k+1)∫t0t1‖ut(s)‖k+12ds}.
UsingT=Nτ and Gronwall’s lemma, we get
(3.55)‖|θ|‖0,α2=∑n=0Nτ‖θ‖α2.
Hence, using the interpolation property and(3.56)‖|u-uh|‖0,α≤‖|θ|‖0,α+‖|ρ|‖0,α,
the estimate (3.32) holds.
Also using the interpolation property, Gronwall’s lemma, and the approximation properties, we get(3.57)‖|u-uh|‖∞,0≤‖|θ|‖∞,0+‖|ρ|‖∞,0≤max⁡0≤n≤N‖|θn|‖2+h2k+1‖|u|‖∞,k+12,
which is just the estimate (3.33).


## 4. Numerical Examples

In this section, we present the numerical results which confirm the theoretical analysis in Section3.LetK denote a uniform partition on [0,a],  and Sh the space of continuous piecewise linear functions on K, that is, k=1.. In order to implement the Galerkin finite element approximation, we adapt finite element discrete along the space axis, and finite difference scheme along the time axis. We associate shape function of space Xh with the standard basis of hat functions on the uniform grid of size h=1/n. We have the predicted rates of convergence if the condition Δt=ch of (4.1)‖u-uh‖0,α~O(h2-α),‖u-uh‖∞,0~O(h2-α),
provided that the initial value φ(x) is smooth enough.Example 4.1.
The following equation(4.2)∂u∂t=D(u2D0x0.5u(x,t))-2x(x-1)(2x1.5Γ(2.5)-x0.5Γ(1.5))e-2tDu-u(x,t)-u2e-t(2x0.5Γ(1.5)-x-0.5Γ(0.5)),0≤x≤1,0≤t≤1,u(x,0)=x(x-1),0≤x≤1,u(0,t)=u(1,t)=0,0≤t≤1,
has a unique solution u(x,t)=e-tx(x-1).
If we selectΔt=ch and note that the initial value u0 is smooth enough, then we have
(4.3)‖u-uh‖0,0.75~O(h1.25),‖u-uh‖∞,0~O(h1.25).Table1 includes numerical calculations over a regular partition of [0,1]. We can observe the experimental rates of convergence agree with the theoretical rates for the numerical solution.Table 1
Numerical error result for Example4.1.
h∥u-uh∥∞,0cvge. rate∥u-uh∥0,0.75cvge. rate1/52.2216E−003—1.0213E−003—1/101.3551E−0030.71326.0779E−0040.748751/205.5865E−0041.27842.3188E−0041.39011/403.0515E−0040.87241.0545E−0041.13671/801.2423E−0041.29643.9883E−0051.40271/1605.1033E−0051.28352.1310E−0050.9042Example 4.2.
The functionu(x,t)=cos(t)x2(2-x)2 solves the equation in the following form:
(4.4)∂u∂t=Dx1.70u(x,t)+b(u)Du-u[4(1-x)+tant]+f(x,t),x∈(0,2),t∈[0,1),u(x,0)=x2(2-x)2,0≤x≤2,u(0,t)=0,u(2,t)=0,0≤t≤1,
where
(4.5)b(u)=ucos⁡t,f(x,t)=cos⁡tcos⁡(0.85π)[24(x2.3+(2-x)2.3)Γ(3.3)-24(x1.3+(2-x)1.3)Γ(2.3)-8(x0.3+(2-x)0.3)Γ(1.3)].
If we selectΔt=ch, then
(4.6)‖u-uh‖0,0.85~O(h1.15),‖u-uh‖∞,0~O(h1.15).Table2 shows the error results at different size of space grid. We can observe that the experimental rates of convergence still support the theoretical rates.Table 2
Numerical error result for Example4.2.
h∥u-uh∥∞,0cvge. rate∥u-uh∥0,0.85cvge. rate1/51.3010E−001—3.2223E−002—1/104.6402E−0021.48781.4133E−0021.18901/201.6843E−0021.46206.2946E−0031.16691/406.6019E−0031.68432.8571E−0031.13951/802.7979E−0031.23861.3137E−0031.12091/1601.2665E−0031.14346.0848E−0041.1103Example 4.3.
Consider the following space-fractional differential equation with the nonhomogeneous boundary conditions,(4.7)∂u∂t=D0x1.7u(x,t)-3x∫0xudx-2x0.3e-tΓ(1.3),0≤x≤1,0≤t≤1,u(x,0)=x2,0≤x≤1,u(0,t)=0,u(1,t)=e-t,0≤t≤1,
whose exact solution is u(x,t)=e-tx2.
We still chooseΔt=ch, then get the convergence rates
(4.8)‖u-uh‖0,0.85~O(h1.15),‖u-uh‖∞,0~O(h1.15).The numerical results are presented in Table3 which are in line with the theoretical analysis.Table 3
Numerical error result for Example4.3.
h∥u-uh∥∞,0cvge. rate∥u-uh∥0,0.85cvge. rate1/58.3052E−002—2.8009E−002—1/103.6038E−0021.20451.0086E−0021.47351/201.3839E−0021.38073.2327E−0031.64141/405.0631E−0031.45071.1789E−0031.45541/801.8920E−0031.42015.9555E−0040.98511/1609.9899E−0040.92143.0034E−0040.9876


## 5. Conclusion

In this paper, we propose a fully discrete Galerkin finite element method to solve a type of fractional advection-diffusion equation numerically. In the temporal direction we use the modified Crank-Nicolson method, and in the spatial direction we use the finite element method. The error analysis is derived on the basis of fractional derivative space. The numerical results agree with the theoretical error estimates, demonstrating that our algorithm is feasible.



---

*Source: 171620-2011-10-05.xml*