# A Fully Discrete Discontinuous Galerkin Method for Nonlinear Fractional Fokker-Planck Equation

Authors: Yunying Zheng; Changpin Li; Zhengang Zhao
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

The fractional Fokker-Planck equation is often used to characterize anomalous diffusion. In this paper, a fully discrete approximation
for the nonlinear spatial fractional Fokker-Planck equation is given, where the discontinuous Galerkin finite element approach is utilized in time domain and the Galerkin finite element approach is utilized in spatial domain. The priori error estimate is derived in detail. Numerical examples are presented which are inline with the theoretical convergence rate.


---

## Body

## 1. Introduction

Many models in physics, chemistry are successfully described by the Langevin equation, which has been introduced almost 100 years before. And for some particular cases, say diffusion, the original Langevin equation can be transformed into the Fokker-Planck equation. Hänggi and Thomas [1] associated a Gaussian distribution of the increments of the noise generating process with the classical Fokker-Planck equation. Sun et al. [2] discussed the fractional model for anomalous diffusion. Metzler et al. [3] and Dubkov and Spagnolo [4] derived the fractional Fokker-Planck equation from different anomalous diffusion procedures. Metzler and Klafter [5] discussed fractional kinetic equation and its relation to the fractional Fokker-Planck equation. Dubkov et al. [6] introduced Fokker-Planck equation for Lévy flights. Now the Fokker-Planck equation is one of the best tools for characterizing anomalous diffusion, especially sub-/super-diffusion. Meanwhile the fractional Fokker-Planck equation has been found to be used in relatively wide field of applied sciences, such as plasma physics, population dynamics, biophysics, engineering, neuroscience, nonlinear hydrodynamics, and marking; see [7–13].The Fokker-Planck equation describes the changes of a random function in space and in time. So different assumptions on probability density function lead to a variety of space-time equations. In this paper, we mainly study the model described by the following fractional Fokker-Planck equation, which is a special case in [12]:(1.1)∂u∂t=(DC0,xβb(x,t,u)+DC0,x2βa(x,t,u))u,0<β<1,(x,t)∈Ω×(0,T],
where DC0,xβ denotes left fractional derivative with order β in the sense of Caputo.There are some numerical methods to find the approximate solutions of the fractional differential equations [14–21]. But the discontinuous Galerkin finite element method is a very attractive method for partial differential equations because of its flexibility and efficiency in terms of mesh and shape functions. And the higher order of convergence can be achieved without over many iterations. Such a method was first proposed and analyzed in the early 1970s as a technique to seek numerical solutions of partial differential equations. The discontinuous Galerkin finite element method becomes a very attractive tool for the initial problems of the ODEs and the initial-boundary problems of PDEs; see [22–27].The rest of this paper is constructed as follows. The fractional derivative space is introduced in Section2. In Section 3, the discontinuous Galerkin finite element scheme is introduced. The existence and uniqueness of numerical solution are proved in Section 3. And the error estimate of the discontinuous Galerkin finite element approximation is studied in Section 4. Finally in Section 5, numerical examples are also taken to show the efficiency of the theoretical results.


## 2. Fractional Derivative Space

In this section, we firstly introduce the fractional integral, fractional Caputo derivative, and their properties.Definition 2.1.
Theαth order left and right Riemann-Liouville integrals of function u(x) are defined as follows:
(2.1)Da,x-αu(x)=1Γ(α)∫axu(s)(x-s)1-αds,Dx,b-αu(x)=1Γ(α)∫xbu(s)(s-x)1-αds,
where α>0 and a<x<b.Definition 2.2.
Theαth order left and right Caputo derivatives of function u(x) are defined in
(2.2)DCa,xαu(x)=Da,xα-nDnu(x),DCx,bαu(x)=Dx,bα-n(-D)nu(x),
where n-1<α<n∈Z+ and a<x<b.Definition 2.3.
Letα>0; define a fractional derivative space Jsα,p(Ω)(Ω=(a,b)) as
(2.3)Jsα,p(Ω)={u(x)∈Hn(Ω):DCa,xαu(x),DCx,bαu(x)∈Lp(Ω),n-1≤α<n},
endowed with the seminorm
(2.4)|u|Jsα,p=|(DCa,xαu(x),DCx,bαu(x))1/2|Lp(Ω)
and the norm
(2.5)‖u(x)‖Jsα,p=(∑k≤[α]‖Dβu(x)‖p+|u(x)|Jsα,pp)1/p,β∈Z,0≤β≤[α]-1.With the help of Fourier Transform, we can conclude that the following three expressions:∫Ω|DCa,xαu(x)·CDx,bαu(x)|dΩ, ∫Ω(DCa,xαu(x))2dΩ, and ∫Ω(DCx,bαu(x))2dΩ, are equivalent [20]. So (∫Ω(DCa,xαu(x))2dΩ)1/2 and (∫Ω(DCx,bαu(x))2dΩ)1/2 can also be recognized as the seminorms of the fractional space Jsα,p, and when we use the seminorm of fractional space Jsα,p, there is no difference among them. Now we restrict our discussion in case p=2, and the following notations are used, Jsα,2(Ω) being rewritten as Jsα(Ω) with norm ∥·∥Jsα or ∥·∥α and seminorm |·|Jsα or |·|α. We denote Js,0α(Ω) as the closure of C0∞(Ω) under its norm.The following lemmas are useful for our discussions later on.Lemma 2.4 (see [19]).
Letn-1≤α<n; if u∈Js,0α([a,b]), then
(2.6)(Da,x-αDCa,xαu)(x)=u(x),(Dx,b-αDCx,bαu)(x)=u(x).Lemma 2.5 (see [20]).
Foru∈Js,0α([a,b]), 0<β<α, then
(2.7)DCa,xαu(x)=DCa,xα-βDCa,xβu(x).Lemma 2.6.
Foru∈C02α([a,b]),v∈C0α([a,b]),n-1≤α<n, then
(2.8)(DCa,x2αu,v)=(-1)n(DCa,xαu,DCx,bαv).Lemma 2.7 (see [21]).
Foru∈Js,0α(Ω), 0<β<α, then
(2.9)|u|β≤|u|α,‖u‖β≤‖u‖α.
For u∈H0α(Ω), then
(2.10)‖u‖L2(Ω)≤c|u|H0α.
For 0<β<α, then
(2.11)‖u‖H0β(Ω)≤C|u|H0α.Lemma 2.8.
Letμ>0. The following mapping properties hold:
(2.12)DCa,xα:Jsμ(Ω)⟶L2(Ω)is a bounded linear operator,DCx,bα:Jsμ(Ω)⟶L2(Ω)is a bounded linear operator.
The proof is similar to that in [19].


## 3. The Space-Time Discontinuous Galerkin Finite Element Approximation

In this section, we formulate a fully discrete discontinuous Galerkin finite element method for a type of nonlinear spacial fractional Fokker-Planck equation.Problem 1 (nonlinear fractional Fokker-Planck equation).
For0<β<1,
(3.1)ut=[-DC0,xβb(x,t,u)+DC0,x2βa(x,t,u)]u,(x,t)∈Ω×(0,T],u(x,t)=ϕ(x,t),x∈∂Ω×(0,T],u(x,0)=g(x),x∈Ω¯,
where the Ω is a bounded domain. For positive constants m and M, the coefficients a and b satisfy
(3.2)0<m<a(u)<M,0<m<b(u)<M.Throughout the paper, we always assume that the following mild Lipschitz continuity conditions ona,b are satisfied: there exists a positive constant L such that for x∈Ω,t∈(0,T], and s,r∈R, there have
(3.3)|a(x,t,s)-a(x,t,r)|≤L|s-r|,|b(x,t,s)-b(x,t,r)|≤L|s-r|.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.4)Sh={v∈H0β(Ω):v∣K∈Pr-1(K),K∈Kh},
where Pr-1(K) is a set of polynomials with degree r-1 on a given domain K. And the functions in Sh are continuous on Ω.Let us consider a partition of time domainI=(0,T], 0=t0<t1<⋯<tN=T, and In=(tn,tn+1], kn=tn+1-tn,n=0,1,…N-1, k=maxkn. On each time slab In, we define a discrete function space ℑk as (3.5)Ik={v:In⟶R;v∣t∈In∈Pq(In)}.
The functions in the space ℑk can be discontinuous at the time node tn, but is at least left continuous and right continuous. And the functions in the space ℑk are polynomials, whose degree is no more than q-1.Define a discrete function spaceTkh on Ω×I as
(3.6)Tkh={v:Ω×I⟶R∣∀t∈I,v(⋅,t)∈Ik,∀x∈Ω,v(x,⋅)∈Sh}.
Moreover, the space Tkh can be verified as follows:
(3.7)Tkh={v∣v(x,t)=∑j=0q-1tjχj(x),χj∈Sh},
where q is a positive integer. In other words, for each t∈In the functions in Tkh are the elements of Sh, and for each x∈Ω piecewise polynomial functions in t of degree q-1 with possible discontinuities at the nodes tn,n=0,1,…,N-1. Set Tkhn={v∣Ω×In:v∈Tkh}.We introduce some norms on different spaces which will be used later on. AndL2(I,Jsβ(Ω)) is equipped with the norms (3.8)‖v‖0,β=(∫I‖v‖β2dt)1/2,‖v‖∞,β=max⁡I⁡‖v‖β=max⁡0≤t≤tN⁡‖v‖β.AndL2(In,Jsβ(Ω)) is equipped with the norms (3.9)‖v‖0,β,n=(∫In‖v‖β2dt)1/2,‖v‖∞,β,n=max⁡In⁡‖v‖β=max⁡tn<t≤tn+1⁡‖v‖β.In order to derive a variational form of Problem1, we assume that u is a sufficiently smooth solution of Problem 1, then multiply an arbitrary v∈Tkh to obtain the integration formulation (3.10)∫I(∂u∂t,v)dt=∫I(DC0,x2βa(u)u,v)dt-∫I(DC0,xβb(u)u,v)dt,
where (·,·) denotes the inner product on L2(Ω). Integrating by parts in the right, noting that for all t∈I, v∈J0,sβ(Ω), and the discontinuous property at the time node tn,n=0,…,N-1, one has (3.11)∑n=0N-1{∫In(∂u∂t+DC0,xβb(u)u,v)+(DC0,xβa(u)u,DCx,bβv)dt+([u]n,vn+)}=0.
The notation [u]n denotes the jump of the function u at the time node tn; that is,(3.12)[u]n=lim⁡ε→0⁡[u(⋅,tn+|ε|)-u(⋅,tn-|ε|)].Using the superscripts “–” and “+” for left and right limits, respectively, the jump is described by (3.13)[u]n=u(⋅,tn+)-u(⋅,tn-).
It shows the discontinuity of the scheme. The inner product ([u]n,vn+) also denotes the transport process during different time-space slabs.Thus, we define(3.14)A(u,v)=∫In{(∂u∂t,v)+(DC0,xβa(u)u,CDx,bβv)+(DC0,xβb(u)u,v)}dt+([u]n,vn+).Definition 3.1.
For allv∈Tkh, the function u∈Tkh is a variational solution of Problem 1 provided that
(3.15)∑n=0N-1A(u,v)=0,∀v∈Tkh.Now we are ready to describe a fully discrete space-time finite element method to solve the nonlinear problem1, where the Galerkin finite element method is used in spacial domain and discontinuous finite element approximation is used in time domain.Applying the same notations with Definition3.1, the space-time discontinuous Galerkin scheme for Problem 1 can be now formulated as follows. Find U∈Tkhn, satisfying, for all v∈Tkhn,n=0,1,…,N-1,
(3.16)A(U,v)=0,∀v∈Tkh.The value ofU on t0- is replaced by the initial condition g(x). The term including U0 can be moved to the right-hand side of (3.16). Once it is computed, the value of U on t0+ is different from that on u0-, and it is the error introduced by the discretization in the numerical scheme.Noting the discontinuity at each nodetn,n=0,1,…,N-1 in Thk, the computation of uh can be decoupled in each time slab. Once uh is known on tn-, this value is taken as an initial condition for the time slab In and the following equation needs to be solved:(3.17)(Un+1-,vn+1-)+∫In{-(U,∂v∂t)+(DC0,xβa(U)U,DCx,bβv)+(DC0,xβb(U)U,v)}dt=(Un-,vn+).Next, we investigate the uniqueness and existence about the numerical scheme. First, we give the scheme in detail. For a fixed integerq≥1, let li(τ)i=1q be the Lagrange polynomials associated with the abscissa 0<τ1<⋯<τq=1; that is, (3.18)li(τ)=∏j=1,…,qi≠jτ-τjτi-τj.For the quadrature in time slab(tn,tn+1], we use the Radau quadrature rule. For a given function g(τ),τ∈[0,1], the following approximation holds: (3.19)∫01g(τ)dτ≈∑j=1qwjg(τj),
where wj=∫01lj(τ)dτ. This quadrature rule is exact for all polynomials of degree ≤2q-2. Using the linear transformation t=tn+τkn, which maps [0,1] into I¯n, one gets
(3.20)tn,j=tn+τjkn,j=1,…,q,tn,q=tn+1,ln,j(t)=lj(τ),t=tn+τkn,wn,j=∫tntn+1ln,jdt=kn∫01lj(τ)dτ=knwj,j=1,…,q.Then the Radau rule inI¯n can be got by(3.21)∫Ing(t)dt=∑j=1q∫01lj(τ)dτg(tn+τjkn)=∑j=1qwn,jg(tn,j).We choose{ln,i(t)}j=1q as the basis functions for the piecewise polynomial function space Pq-1(In). Then U∣In is uniquely determined by the functions Un,j=Un,j(x),U(x)∈Sh, such that (3.22)U(x,t)=∑j=1qln,j(t)Un,j(x),(x,t)∈Ω×In.Takingv=ln,i(t)ψ(x)∈Vhkn into (3.17), where ψ∈Shn, we have (3.23)(ln,q(tn,q)Un,q,ln,i(tn,q)ψ(x))-∫In(∑j=1qln,jUn,j(x),ln,i′(t)ψ(x))dt+∫In(DC0,xβa(∑k=1qln,kUn,k(x))∑j=1qln,jUn,j(x),DCx,bβ(ln,i(t)ψ(x)))dt+∫In(DC0,xβb(∑k=1qln,kUn,k(x))∑j=1qln,jUn,j(x),ln,i(t)ψ(x))dt=(ln,1(tn,1)Un,1,ln,i(tn,1)ψ(x)).On time slabIn=(tn,tn+1], we define a Lagrange interpolation operator În:C(In)→Pq-1(In), such that (3.24)Îny(tn,j)=y(tn,j),j=1,2,…,q,
where the interpolation points are Radau points. It is easy to see that for given x∈Ω, Înu(x,t)∈Pq-1(Tn) and Înu(x,t)=u(x,tn+1) are available. On time slab In=(tn,tn+1], tn,1=tn,tn,q=tn+1,  with the definition of interpolation operator, the discontinuous Galerkin scheme can be rewritten in detail as follows: (3.25)δqi(Un,q,ψ(x))-∫In(∑j=1qln,jUn,j(x),ln,i′(t)ψ(x))dt+∫In(DC0,xβa(∑k=1qln,kUn,k(x))∑j=1qln,jUn,j(x),CDx,bβ(ln,i(t)ψ(x)))dt+∫In(DC0,xβb(∑k=1qln,kUn,k(x))∑j=1qln,jUn,j(x),ln,i(t)ψ(x))dt=ln,i(tn)(Un,ψ(x)).Next we introduce a lemma which is useful to prove the existence and estimate the error.Consider aq×q matrix W, (3.26)Wij=wn,jln,i′(tn,j)=wjli′(τj),M=eqeqT-N,
where eq=(0,0,…,1).It is clear thatW,M are independent of Kn. And if YT=(yn,1,…,yn,q)∈Rq, then (3.27)YTMY=∑j=1qδqiyn,qyn,i-∑j=1qwn,jln,i′(tn,j)yn,jyn,i.Lemma 3.2 (see [27]).
LetM̃,D be the matrixes
(3.28)M̃=D-1/2MD1/2,D=diag⁡{τ1,…,τq}.
If λ=(α0/2)min{w1/τ1,w2/τ2,…,wq-1/τq-1,1+wq}>0, then there holds
(3.29)XTM̃X≥λ|X|2,∀X∈Rq,
where |x|=∑j=1qXi2.In order to prove the existence and uniqueness, we need to define a new exhibition forU by
(3.30)U=∑j=1qτj1/2ln,j(T)Un,j(x),(x,t)∈Ω×In,
where Ũ=τj-1/2Un,j∈Shn. Then U∣In is uniquely determined by the function Ũ=Ũ(x)∈Shn.We chooseψ=τi-1/2ϕ, where ϕ∈Shn, and use the new expression U to obtain the following results: (3.31)δqi(Ũ,ϕ)+∫In(-∑j=1qτj1/2ln,jŨn,j,∑k=1qτi-1/2ln,kln,i′(tn,k)ϕ)dt+∫In(∑j=1qτj1/2DC0,xβa(∑m=1qτm1/2ln,mŨn,m)ln,jŨn,j,DCx,bβτi-1/2ln,iϕ)dt+∫In(∑j=1qτj1/2DC0,xβb(∑m=1qτm1/2ln,mŨn,m)ln,jŨn,j,τi-1/2ln,iϕ)dt=ln,i(tn+)τ-1/2(Ũn,ϕ).Set thata(Ũ)=a(∑m=1qτm1/2ln,mŨn,m) and Wn,j=∫Inln,j2(t)dt, the above equation can be rewritten as (3.32)-∑j=1qWn,jτj1/2τi-1/2ln,i′(tn,j)(Ũn,j,ϕ)+∑j=1qWn,jτj1/2τi-1/2(DC0,xβa(Ũ)Ũn,j,DCx,bβϕ)+∑j=1qWn,jτj1/2τi-1/2(DC0,xβb(Ũ)Ũn,j,ϕ)+δqi(Ũ,ϕ)=ln,i(tn+)τ-1/2(Ũn,ϕ).Theorem 3.3.
LetUn be given in Shn, then for the sufficiently small kn, there exists {Ũn,j}j=1q∈Shn satisfying (3.32); hence, (3.17) has a unique solution U∈Thk.Proof of Theorem3.3.
The vector space(Shn)q is a Hilbert space with finite dimension. For all x=(x1,x2,…,xq)T, ϕ=(ϕ1,ϕ2,…,ϕq)T∈Shn, equipped with the norm(3.33)‖|x|‖=(∑j=1q‖xi2‖)1/2,‖|x|‖β=(∑j=1q‖xi2‖β)1/2.
Define a mapF from (Shn)q to itself by
(3.34)(F(v)i,ϕ)=(vi,ϕ)+δqi(vq,ϕ)-∑j=1qWn,jτj1/2τi-1/2ln,i′(vi,ϕ)+Wn,i(DC0,xβ(a(v)vi),DCx,bβϕ)+Wn,i(DC0,xβb(v)vi,ϕ)-(Un,τi-1/2ln,i(tn)ϕ).
Since (3.35)(vi,ϕ)+δqi(vq,ϕ)-∑j=1qWn,jτj1/2τi-1/2ln,i′(vi,ϕ)-(Un,τi-1/2ln,i(tn)ϕ)
is a linear system, a(x,t,u),b(x,t,u) are all continuous maps, and the map F is a continuous map from (Shn)q into itself. So according to the Brouwer fixed point theorem, there exists at least one fixed point vi such that (F(v)i,ϕ)=(vi,ϕ). So (3.32) has at least one solution, denoted by v. Next we investigate the uniqueness of the solution.
Letv and v* be two solutions of (3.32), setting ϕ=vi-vi* and summing i from 1 to q, then we have
(3.36)∑i=1q(F(v-v*)i,vi-vi*)=δqi(vq-vq*,vi-vi*)-∑i=1q∑j=1qWn,jln,i′(tn,j)τj1/2τi-1/2(vj-vj*,vi-vi*)+∑i=1q∑j=1qWn,jτj1/2τi-1/2(DC0,xβa(v-v*)(vj-vj*),DCx,bβ(vi-vi*])+∑i=1q∑j=1qWn,jτj1/2τi-1/2(DC0,xβb(v-v*)(vj-vj*),(vi-vi*))+∑i=1q(vi-vi*,vi-vi*).
According to Lemma 3.2, we can see that for the first two terms of the right-hand side of (3.36), there exists a constant α, such that (3.37)δqi(vq-vq*,vi-vi*)-∑i=1q∑i=1qWn,jln,i′(tn,j)τj1/2τi-1/2(vj-vj*,vi-vi*)=(M̃(v-v*),v-v*)≥α‖|v-v*|‖2≥Cα‖|v-v*|‖β2.
With the help of bounded assumption of a and b, we have (3.38)∑i=1q∑j=1qWn,jτj1/2τi-1/2(DC0,xβa(v-v*)(vj-vj*),DCx,bβ[(vi-vi*)])≤∑i=1q∫In|(DC0,xβa(v-v*)∑j=1qτj1/2τi-1/2ln,j(t)(vj-vj*),ln,i(t)DCx,bβ(vi-vi*))|dt≤∑i,j=1qM∫In|τj1/2τi-1/2ln,i(t)ln,j(t)(DC0,xβ(vj-vj*),DCx,bβ(vi-vi*))|dt≤Ckn|‖vk-vk*‖|β2,∑i=1q∑j=1qWn,jτj1/2τi-1/2(DC0,xβb(v-v*)(vj-vj*),(vi-vi*))≤∑i=1q∫In|(∑j=1qτj1/2τi-1/2ln,j(t)DC0,xβb(v-v*)(vj-vj*),ln,i(t)(vi-vi*))|dt≤∑i,j=1qM∫In|τj1/2τi-1/2ln,i(t)ln,j(t)DC0,xβ(vj-vj*),(vi-vi*)|dt.
According to the boundedness of fractional operator DC0,xβ, inequality ∥DC0,xβ(vi-vi*)∥≤C∥vi-vi*∥ holds. Furthermore, we have (3.39)∑i,j=1qM∫In|τj1/2τi-1/2ln,i(t)ln,j(t)(DC0,xβ(vj-vj*),(vi-vi*))|dt≤∑j=1qτj1/2τj-1/2‖Wn,k(vk-vk*,vk-vk*)‖=M∑k=1qWn,k‖vk-vk*‖2≤Ckn|||vk-vk*|||β2.
The fifth term of the right-hand side of (3.36) is estimated in
(3.40)∑j=1q(vi-vi*,vi-vi*)=|||v-v*|||2≤|||v-v*|||β2.
And by the Brouwer fixed point theorem, one gets(3.41)(F(v-v*)i,vi-vi*)=(vi-vi*,vi-vi*),
furthermore,
(3.42)∑i=1q(F(v-v*)i,vi-vi*)=∑i=1q(vi-vi*,vi-vi*)=|||vi-vi*|||2.
By using inequalities (3.37)–(3.40) and (3.42), one has
(3.43)C(α-kn)|||vi-vi*|||β2≤0.
Choosingα≥kn, we have
(3.44)|||vi-vi*|||β2≤0.
That is impossible. So the uniqueness is proved.


## 4. Error Estimation

Now we turn to analyze the error estimate of the D-G scheme.Letuhn be the approximate solution on time slab In. Denote Phun as a Ritz-Galerkin projection operator defined as follows: (4.1)(DC0,xβa(u)(un-Phun),DCx,bβv)=0,(DC0,xβa(u0)(un-Phun),DCx,bβv)=0.Andρn=Phun-un,θn=uhn-Wn, then one has (4.2)uhn-un=uhn-Wn+Wn-un=θn+Wn-un,
where W:I→H0β(Ω) is defined by (4.3)W(x,t)=ÎnqPhu(x,t),(x,t)∈Ω×In.
It can be seen that W is an element of Thkn.Lemma 4.1.
Leta(u), b(u) be smooth functions on Ω, and 0<m≤a(u), b(u)≤M,0<β<1, and Phun is defined as above, then
(4.4)‖Phun-un‖≤c2hn1+β‖u‖1+β.Proof of Lemma4.1.
For allϕ∈L2(Ω), w is the solution of the following equation:
(4.5)-a(w)DC0,x2βw=ϕ,w∈Ω,w=0,w∈∂Ω.
So the next equation holds(4.6)‖w‖2β≤c2‖ϕ‖.
And for allχ∈Sh, by the aid of approximation properties of Sh and the weak form, we can derive
(4.7)(Phun-un,ϕ)=-(Phun-un,a(w)DC0,x2βw)=(DCx,bβa(w)(Phun-un),DC0,xβw)=(DCx,bβa(w)(Phun-un),DC0,xβ(w-χ))≤M‖Phun-un‖β‖w-χ‖β≤M‖Phun-un‖βinf⁡χ∈Sh‖w-χ‖β≤c2hr-β‖u‖rhβ‖w‖β=c2hr‖u‖r‖ϕ‖.
So
(4.8)‖Phun-un‖=sup⁡0≠ϕ∈L2(Ω)⁡(Phun-un,ϕ)∥ϕ∥≤c2hr∥u∥r.Lemma 4.2.
LetPhun be defined as (4.1), then ∥D(Phun)∥ and ∥DC0,xβPhun∥ are bounded, where 0≤β≤1.Proof of Lemma4.2.
Using the properties of the norm, then the following inequality is valid:(4.9)‖DC0,xβPhun‖L∞≤‖DC0,xβ(Phun-Ihun)‖L∞+‖DC0,xβ(Ihun)‖L∞.
According to the estimate of interpolation, we get
(4.10)‖DC0,xβ(Phun-Ihun)‖L∞≤c(‖DC0,xβρn‖L∞+‖DC0,xβ(Ihun-un)‖L∞).
From Lemmas 2.8 and 4.1, and the inverse estimate, we know that ∥DC0,xβρn∥L∞≤c∥un∥. Therefore, it is easy to see that
(4.11)‖DC0,xβ(Ihun)‖L∞≤c‖DC0,xβun‖L∞,‖DC0,xβ(Ihun-un)‖L∞≤c‖DC0,xβun‖L∞.
So
(4.12)‖DC0,xβ(Phun)‖L∞≤c‖DC0,xβun‖L∞+c‖un‖.
The boundedness ∥DC0,xβPhun∥ can be derived from above. The proof of ∥DPhun∥≤c(u) can be similarly given.LetÎnq:C(In)→Pq(In) be the usual Lagrange interpolation operator at the Radau points on In=(tn,tn+1]; that is, (4.13)(Înqy)(tn,j)=y(tn,j),j=1,2,…,q,
where tn,j,j=1,2,…,q are the Radau points. Therefore, we can see that (4.14)(Înqu)(x,⋅)∈Pq-1(In),(Înqu)(x,tn+1)=u(x,tn+1),x∈Ω.Lemma 4.3.
The interpolation operationÎnq on Ph(In) has the following property:
(4.15)‖Înqu-u‖≤chs‖u‖s,∀u∈Hs∩H0β.Lemma 4.4.
ForW defined above, we have the following error estimates:
(4.16)max⁡In⁡‖u-W‖β≤Cknqmax⁡In‖u(q)‖β+chs-βmax⁡In‖u‖s,2≤s≤r+1,‖u-W‖0,β,n≤Cknq‖u(q)‖0,0,n+Chns-β‖u‖0,s,n,2≤s≤r.Proof of Lemma4.4.
Since(4.17)‖u-W‖β=‖u-Phu+Phu-W‖β≤‖u-Phu‖β+‖Phu-W‖β,
according to Lemma 4.1, we have
(4.18)‖u-Phu‖β+‖Phu-W‖β≤hs-β‖u‖s+‖Phu-W‖β≤hs-β‖u‖s+Cknq‖(Phu)(q)‖β≤hs-β‖u‖s+Cknq‖u(q)‖β≤hs-β‖u‖s+Cknq‖u(q)‖β.
The first inequality can be similarly proved.We are now ready to prove convergence result. PuttingW into D-G finite element formula (3.17), and setting θ=U-W, we can get the basic error equation as follows: (4.19)(θn+1-,vn+1-)-(θn-,vn+)-∫In(θ,vt)dt+∫In(DC0,xβa(U)θ-b(U)(U-W),DCx,bβv)dt=-(Wn+1-,vn+1-)+(Wn-,vn+)+∫In(W,vt)dt-∫In(DC0,xβa(U)W+b(U)W,DCx,bβv)dt,
where W0=Phu0 and θ0=U0-W0=u0-Phu0.Theorem 4.5.
Letu be the solution of Problem 1. If U is a solution of the discontinuous Galerkin scheme (3.17), then the following error estimation holds:
(4.20)‖u-U‖∞,β=Cmax⁡0≤n≤N-1⁡hsmax⁡In⁡(‖ut‖s+‖u‖s+h-β‖u‖s)+Cmax⁡0≤n≤N-1{knqmax⁡In(‖u(q+1)‖+‖u(q)‖+‖Dx2βu(q)‖+‖Dxβu(q)‖)},
where 2≤s≤r.Proof of Theorem4.5.
It is necessary to rewriteW and discrete error θ as follows:
(4.21)W=∑j=1qln,j(t)Phun,j,θ=U-W=∑j=1qln,j(t)(Un,j-Phun,j),
where un,j=u(x,tn,j). We also denote a(∑j=1qln,j(t)Un,j),b(∑j=1qln,j(t)Un,j)  as a(U) and b(U). Noting that tn,q=tn+1, tn,1=tn, so ln,q(tn,q)=1 and ln,1(tn,1)=1. Setting Wn,i=∫Inln,i2(t)dt and v=ln,i(t)ϕ(x),ϕ∈Shn, then the basic error equation can be rearranged as follows:
(4.22)-δqi(Phun,q,ϕ)+ln,i(tn)(Phun,ϕ)+(∑j=1qWn,jln,i′(tn,j)Phun,j,ϕ)-(un-,ln,i(tn)ϕ)-(∑j=1qWn,jDC0,xβa(U)Phun,j,DCx,bβϕ)+∑j=1qWn,j(b(U)(t)Phun,j,DCx,bβϕ)-∫I-n(u,ln,i′(t)ϕ)dt+∫In(DC0,xβ(a(u)u)+b(u)u,ln,i(t)DCx,bβϕ)dt+(un+1-,ln,i(tn+1)ϕ)=(R1,ϕ)+(R2,ϕ)+(R3,DCx,bβϕ)+(R4,DCx,bβϕ)+(R5,ϕ),
where
(4.23)R1=∑j=1qWn,jln,i′(tn,j)(Phun,j-un,j)-δqi(Phun,q-un,q)-ln,i(tn)(un+-Phun+),R2=∑j=1qWn,jln,i′(tn,j)un,j-∫I-nln,i′(t)udt,R3=-∑j=1qWn,jDC0,xβa(U)Phun,j+∫Inln,i(t)DCx,bβa(u)udt,R4=-∑j=1qWn,jb(U)(Phun,j)+∫Inln,i(t)b(u)udt,R5=ln,i[un-Phun]=ln,i(tn){(un+-Phun+)-(un--Phun-)}.
Letϕ=θn,i and sum i from 1 to q, then the error equation can be derived as
(4.24)(θn+1-,θn+1-)-(θn-,θn+)-∫In(θ,θt)dt+∫In(DC0,xβa(U)θ,DCx,bβθ)dt+∫In(b(U)θ,DCx,bβθ)dt=∑i=1q(R1,θn,i)+(R2,θn,i)+(R3,DCx,bβθn,i)+(R4,DCx,bβθn,i)+(R5,θn,i).
Note that(4.25)-∫In(θ,θt)dt=-12∫Inddt‖θ‖2=-12‖θn+1-‖2+12‖θn+‖2.
Consider the bounded assumption ofa(u) and b(u); the following inequality can be got:
(4.26)∫In(DC0,xβa(U)θ,DCx,bβθ)dt≥m∫In(DC0,xβθ,DCx,bβθ)dt=m∫In‖DC0,xβθ‖2dt.
Sinceb(u) is bounded, and from Lemma 2.8, where DCx,bβ is a linear bounded map, the following inequality is established:
(4.27)∫In(b(U)θ,DCx,bβθ)dt≥m∫In(θ,DCx,bβθ)dt≥Cm∫In‖θ‖‖DCx,bβθ‖dt≥Cm∫In‖θ‖2dt.
The first and second terms of the left side of (4.24) can be seen as
(4.28)(θn+1-,θn+1-)=‖θn+1-‖2,(θn-,θn+)≤‖θn-‖⋅‖θn+‖≤‖θn-‖2+‖θn+‖24.
Taking (4.24)–(4.27) into the left side of (4.22) yields
(4.29)12‖θn+1-‖2+12‖θn+‖2+m∫In‖DC0,xβθ‖2dt+Cm∫In‖θ‖2dt≤‖θn-‖2+‖θn+‖24+∑i=1q{(R1,θn,i)+(R2,θn,i)+(R3,DCx,bβθn,i)+(R4,DCx,bβθn,i)+(R5,θn,i)}.
As toR1, for i=1,2,…,q, there exists
(4.30)δqi-∑j=1qWn,jln,i′(tn,i)-ln,i(tn+)=δqi-∫Inln,i′(t)dt-ln,i(tn+)=0.
This means that there exists constantsCij and C, which are independent of n, such that
(4.31)∥R1∥=‖∑j=1qCij(ρn,j-ρn,j-1)‖≤‖∫InCρt(τ)dτ‖≤C∫Inhs‖ut(τ)‖sdτ≤Chskn1/2‖ut‖0,s,n.
In order to estimateR2, we apply the interpolation operator. Let Îqn:C(În)→Pq(În) be the interpolation operator on the time slab In, whose order is less than q. The interpolation points involve not only Radau points, but also tn, and these points satisfy that Înqy(tn,j)=y(tn,j),Înqy(tn+)=y(tn+),j=1,2,…,q. Then, for every x∈Ω, such that ln,i′Înqu is a polynomial whose order is 2q-2,
(4.32)‖R2‖=‖∑j=1qWn,jln,i′(tn,j)un,j-∫Inln,i′(t)udt‖=‖∫Inln,i′(t)(Înqu-u)dt‖≤Cknq+1/2‖u(q+1)‖0,0,n.
As toR3, when kn is sufficiently small, with the aid of Hölder inequality, we obtain
(4.33)(R3,DCx,bβθn,i)=-(DC0,xβR3,θn,i)≤‖∑j=1qWn,jDC0,x2β(a(U)-a(u))un,j‖⋅‖θn,i‖+‖∑j=1qWn,jDC0,x2βa(u)un,j-∫Inln,i(t)DC0,x2βa(u)udt‖⋅‖θn,i‖=(‖R31‖+‖R32‖)‖θn,i‖,
where
(4.34)‖R31‖=‖∑j=1qWn,jDC0,x2β(a(U)-a(u))un,j‖≤∫Inln,i2(t)‖DC0,x2βθ‖2dt+LC∫Inln,i2(t)‖DC0,x2β(W-u)‖dt≤LC∫Inln,i2(t)‖DC0,x2βθ‖2dt+Cknq+1/2‖DC0,x2βu(q)‖0,0,n≤LC∫Inln,i2(t)‖θ‖2dt+Cknq+1/2‖DC0,x2βu(q)‖0,0,n,‖R32‖=‖∑j=1qWn,jDC0,x2βa(u)un,j-∫Inln,i(t)DC0,x2βa(u)udt‖≤M‖∑j=1qWn,jDC0,x2βÎnqun,j-∫Inln,i(t)DC0,x2βudt‖≤M‖∫Inln,l2(t)(Înq-I)DC0,x2βudt‖≤CMknq+1/2‖DC0,x2βu(q)‖0,0,n.
Almost similar to the estimation ofR3, the analysis of R4 can be got by the help of boundedness of b(u),u, and Phun.
(4.35)(R4,DCx,bβθn,i)≤‖∑j=1qWn,jDC0,xβ(b(U)-b(u))Phun,j‖‖θn,i‖+|∑j=1qWn,jDC0,xβb(u)Phun,j-∫Inln,i(t)DC0,xβb(u)udt|‖θn,i‖≤C{∫Inln,j‖DC0,xβθ‖dt+knq+1/2‖DC0,xβu(q)‖0,0,n+(M+L)hs-βkn1/2‖u‖0,s,n}‖θn,i‖.
The last term containsR5 which is estimated as follows:
(4.36)‖(R5,θn+)‖≤C‖(un+-Phun+)‖2+C‖un--Phun-‖2+14‖θn+‖2≤Chn2s‖u‖s2+Chn-12s‖u‖s2+14‖θn+‖2,
where hn is the space partition in the time slab In.
Taking the estimations fromR1 to R5 into (4.29), and setting Wn,j=knWj, since ∥θ∥0,0,n=(∑j=1qWn,j∥θn,j∥2)1/2=kn1/2(∑j=1qWj∥θn,j∥2)1/2, the following inequality holds:
(4.37)12‖θn+1-‖2+m∫In‖DC0,xβθ‖2dt+Cm∫In‖θ‖2dt‖θ‖0,2,n≤‖θn-‖2+∑i=1q{‖R1‖2+‖R2‖2+‖DC0,xβR3‖2+‖DC0,xβR4‖2}+C∑i=1q‖θn,i‖2+Chn2s‖u‖s2+Chn-12s‖u‖s2.
Taking the interpolation representation ofθ on time slab In into the error equation (4.22), and multiplying τi-1/2 in both side of the error equation, we have
(4.38)-∑j=1qWn,jln,i′(tn,j)τj1/2τi-1/2(θ̃n,j,ϕ)+Wn,jτj1/2τi-1/2(DC0,xβa(U)θ̃n,j,DCx,bβϕ)+Wn,jτj1/2τi-1/2(b(U)θ̃n,j,DCx,bβϕ)+δqi(θ̃n,j,ϕ)=τi-1/2[(θn-,ln,i(tn)ϕ)+(R1,ϕ)+(R2,ϕ)+(R3,DCx,bβϕ)+(R4,DCx,bβϕ)+(R5,ϕ)].
Letϕ=θ̃n,i, and sum from i=1 to q, with Lemma 3.2, we know that
(4.39)∑i=1qδqi(θ̃n,q,θ̃n,i)-∑i,j=1qτi-1/2τj1/2(θ̃n,j,θ̃n,i)=(M̃θ̃n,θ̃n)≥α∑j=1q‖θ̃n,j‖2,
where θ̃n=(θ̃n,1,θ̃n,2,…,θ̃n,q).
The second term of the left side of (4.42) is estimated as
(4.40)∑i=1qWn,i(DC0,xβa(U)θ̃n,i,DCx,bβθ̃n,i)≥m∑i=1qWn,j(DC0,xβθ̃n,i,DCx,bβθ̃n,i)≥mC∑i=1q‖DC0,xβθ̃n,i‖2,∑i=1qWn,i(b(U)θ̃n,i,DCx,bβθ̃n,i)≥m∑i=1qWn,i(θ̃n,i,DCx,bβθ̃n,i)≥mC∑i=1qWn,i‖θ̃n,i‖⋅‖DCx,bβθ̃n,i‖≥mCkn1/2∑i=1q‖DC0,xβθ̃n,i‖2.
From (4.39)-(4.40), the following estimation can be derived:
(4.41)α∑i=1q‖θ̃n,i‖2+mc∑i=1q‖DC0,xβθ̃n,i‖2+mCkn1/2∑i=1q‖DC0,xβθ̃n,i‖2≤C(∑i=1q‖θ̃n,i‖2)1/2{‖θn-‖2+(∑i=1q(‖R1‖2+‖R2‖2+‖DC0,xβR3‖2+‖DC0,xβR4‖2)∑i=1q)1/2+‖[un-Phun]‖}.
Introduce the estimations fromR1 to R4 and the inequality containing ∥[un-Phun]∥ into the above expression, by the definition and characteristics of the norm of the fractional derivative space and the relation of ∥θ∥0,β,n, there holds
(4.42)∥θ∥0,β,n≤Ckn{‖θn-‖β+hns‖ut‖0,s,n+hns-β‖u‖0,s,n+knq(‖u(q+1)‖0,0,n2)+‖u(q)‖0,s,n+‖CD0,x2βu(q)‖0,0,n+‖CD0,xβu(q)‖0,0,n+(hns+hn-1s)‖u‖shns-β‖u‖0,s,n}.
Introducing the estimations fromR1 to R5 into (4.37), based on the boundedness of the fractional operation, one has
(4.43)‖θn+1-‖2≤Ckn2q(‖u(q+1)‖0,0,n2+‖u(q)‖0,0,n2+‖CD0,x2βu(q)‖0,0,n2+‖DC0,xβu(q)‖0,0,n2)+(hn2s+hn-12s)‖u‖s2+hn2s(‖ut‖0,s,n2+hn-2s‖ut‖0,s,n2)+C‖θ‖0,β,n2+‖θn-‖2.
From inequalities (4.42) and (4.43), there holds
(4.44)‖θn+1-‖2≤Ckn2q(1+kn)(‖u(q+1)‖0,0,n2+‖DC0,x2βu(q)‖0,0,n2+‖DC0,xβu(q)‖0,0,n2)+(1+kn)(hn2s+hn-12s)‖u‖s2+(1+kn)hn2s‖ut‖0,s,n2+(1+kn)‖θn-‖β2.
According to the above inequality, we can estimateθn-. With the θj,0≤j≤n,
(4.45)‖θn+1-‖β2≤∏j=0nC(1+kj)‖θ0‖β2+C∑m=0n∏j=m+1n(1+kj){Ckm2q(‖u(q+1)‖0,0,m2+‖DC0,x2βu(q)‖0,0,m2+‖DC0,xβu(q)‖0,0,m2)+(hm2s+hm-12s)‖u‖s2+h2s‖ut‖0,s,n‖DC0,x2βu(q)‖0,0,m2+(1+kn)hn2s-2β‖u‖0,s,n+kn2q(1+kn)‖u(q)‖0,0,n2}.
For a fixedn, if kj→0,j=1,2,…,n, we choose k=max0≤j≤nkj, then there holds
(4.46)∏j=0n(1+kj)≤∏j=0n(1+k)=(1+k)n≤eCtn+1.
SettingC1=eCtn+1, since θ0-=u0-Phu0, we get
(4.47)‖θn+1-‖2≤C‖u0-Phu0‖β2+C1∑m=0n{km2q(‖u(q+1)‖0,0,m2+‖DC0,x2βu(q)‖0,0,m2+‖DC0,xβu(q)‖0,0,m2)+(hm2s+hm-12s)‖u‖s2+hm2s‖ut‖0,s,m+h2s-2β‖u‖s2‖u(q+1)‖0,0,m2+(1+kn)hn2s-2β‖u‖0,s,n+kn2q(1+kn)‖u(q)‖0,0,n2}.
Taking it into (4.43), we have
(4.48)‖θ‖0,β,n≤Ckn1/2h0‖u0‖s+C1kn1/2∑m=0n{km2q(‖uq+1‖0,0,m2+‖DC0,x2βu(q)‖0,0,m2+‖DC0,xβu(q)‖0,0,m2)+(hm2s+hm-12s)‖u‖s2+hm2s‖ut‖0,s,m+(1+kn)hn2s-2β‖u‖0,s,n+kn2q(1+kn)‖u(q)‖0,0,n2}.
Note that ∥U-u∥∞,β,n≤∥θ∥∞,β,n+∥W-u∥∞,β,n. By Lemma 2.8, the estimate of ∥W-u∥ is proved. Therefore, the proof of Theorem 4.5 is ended.


## 5. Numerical Examples

In this section, we present numerical results for the Galerkin approximations which support the theoretical analysis derived in Section4.LetSh denote a uniform partition on spacial domain [0,a]  and Xh the space of continuous piecewise linear functions on Sh. In order to implement the discontinuous Galerkin finite element approximation, we adapt the finite elements scheme in space domain and the discontinuous finite element scheme along time domain. We associate shape function of space Xh with the standard basis of hat functions on the uniform grid of size h=1/n and adapt the same shape functions along time axis. In our scheme, the finite element trial and test spaces for Problem 1 are chosen to be the same.Example 5.1.
Consider the following equation:(5.1)∂u∂t=(-DC0,x1.6x1.6Γ(4.6)-DC0,x0.8x0.8Γ(3.8))u,0≤x≤1,0≤t≤1,u(x,0)=x2,0≤x≤1,u(0,t)=0,u(1,t)=et,0≤t≤1.
The exact solution of the equation is u(x,t)=e-tx2.
If we selectkn=ch and note u(x,0) is smooth enough, then we have the following numerical results presented in Table 1.
Table1 includes numerical calculations over a regular partition of [0,1]. We can see that, when the size of grid becomes smaller, the finite element approximation becomes better. We can also observe that the experimental rates of convergence agree well with the theoretical rates for the numerical solution.Table 1
Numerical error result for Example5.1.
h∥u-uh∥∞,0Cvge. rate1/100.1220e-003—1/200.7433e-0040.71471/400.3244e-0041.19621/800.1642e-0040.98281/1600.7880e-0051.05881/3200.3952e-0050.9955Example 5.2.
We consider another space fractional nonlinear Fokker-Planck differential equation. Theu(x,t)=e-tx3 is the exact solution of the following equation:
(5.2)∂u∂t=(-DC0,x1.812etuΓ(5.8)x1.2+DC0,x0.96etuΓ(4.9)x2.1)u(x,t),x∈(0,2),t∈[0,1),u(x,0)=x3,0≤x≤2,u(0,t)=0,u(1,t)=et,0≤t≤1.
We can also selectkn=ch. Table 2 shows the error results at different sizes of space grid. We can still observe that the experimental rates of convergence are inline with the theoretical convergence rates.Table 2
Numerical error result for Example5.2.
h∥u-uh∥∞,0Cvge. rate1/100.0231—1/200.01270.86101/400.00670.93151/800.00340.96601/1600.00170.98301/3200.00090.9915Example 5.3.
Equationu(x,t)=e-t(x2-x1.9) is the exact solution of the following problem:
(5.3)∂u∂t=-DC0,x1.6(x-0.4u+x1.5)uΓ(4.6)+x0.9Γ(3.5)Γ(4.6)Γ(1.9)-x22e-t,0≤x≤1,0≤t≤1,u(x,0)=x2-x1.9,0≤x≤1,u(0,t)=0,u(1,t)=0,0≤t≤1.
Table 3 includes numerical calculations over a regular partition of [0,1]. We can see that the experimental numerical is valid and inline with the theoretical result. We can also see that this method is valid for more generalized fractional diffusion-type equation, and the rate of convergence depends on the highest order of the equation and the order of the shape function.Table 3
Numerical error result for Example5.3.
h∥u-uh∥∞,0Cvge. rate1/100.6125e-003—1/200.4199e-0030.54471/400.1765e-0031.25001/800.0985e-0030.84251/1600.0527e-0030.9003



---

*Source: 279038-2010-12-15.xml*