Discontinuous Galerkin method for fractional convection-diffusion equations

# Discontinuous Galerkin method for fractional convection-diffusion equations

## Abstract

We propose a discontinuous Galerkin method for convection-subdiffusion equations with a fractional operator of order defined through the fractional Laplacian. The fractional operator of order is expressed as a composite of first order derivatives and fractional integrals of order , and the fractional convection-diffusion problem is expressed as a system of low order differential/integral equations and a local discontinuous Galerkin method scheme is derived for the equations. We prove stability and optimal order of convergence for subdiffusion, and an order of convergence of is established for the general fractional convection-diffusion problem. The analysis is confirmed by numerical examples.

f

ractional convection-diffusion equation, fractional Laplacian, fractional Burgers equation, discontinuous Galerkin method, stability, optimal convergence

{AMS}

26A33, 35R11, 65M60, 65M12.

## 1 Introduction

In this paper we consider the following convection diffusion equation

 {∂u(x,t)∂t+∂∂xf(u)=ε(−(−Δ)α/2)u(x,t),x∈R,t∈(0,T],u(x,0)=u0(x),x∈R (1)

where is assumed to be Lipschitz continuous and the subdiffusion is defined through the fractional Laplacian , which can be defined using Fourier analysis as[22, 33, 37, 24],

 −ˆ(−Δ)α/2u(ξ)=(2π)α|ξ|α^u(ξ) (2)

Equivalently, the fractional Laplacian can also be defined by a singular integral[3, 33],

 −(−Δ)α2u(x)=cα∫|z|>0u(x+z)−u(x)|z|1+αdz (3)

Equation (1), also known as a fractional conservation law, can be viewed as a generalization of the classical convection diffusion equation. During the last decade, it has arisen as a suitable model in many application areas, such as geomorphology[25, 26, 2], overdriven detonations in gases[14, 1], signal processing[5], anomalous diffusion in semiconductor growth [36] etc. With the special choice of , it is recognized as a fractional version of the viscous Burgers’ equation. Fractional conservation laws, especially Burgers’ equation, has been studied by many authors from a theoretical perspective, mainly involving questions of well-posedness and regularity of the solutions[27, 3, 31, 30, 20]. In the case of , the solution is in general not smooth and shocks may appear even with smooth initial datum. Similar to the classical scalar conservation laws, an appropriate entropy formulation is needed to guarantee well-posedness[3, 4]. For the case , the existence and uniqueness of a regular solution has been proved in [3, 30]. In this case, the nonlocal term serves as subdiffusion term, and smooth out discontinuities from initial datum.

Numerical studies of partial differential equations with nonlocal operator have attracted a lot of interest in recent years. Liu, Deng et al. have worked on numerical methods for fractional diffusion problems with fractional Laplacian operators or Riesz fractional derivatives [17, 40, 29, 19]. Briani, Cont, Matache, et al. have considered numerical methods for financial model with fractional Laplacian operators[9, 16, 32]. However, for conservation laws with fractional Laplacian operators, the development of accurate and robust numerical methods remains limited. Droniou is the first to analyze a general class of difference methods for fractional conservation laws[21]. Azeras proposed a class of finite difference schemes for solving a fractional antidiffusive equation[6]. Bouharguane proposed a splitting methods for the nonlocal Fowler equation in [8]. For the case , Cifani et al. [12, 13] applied the discontinuous Galerkin method to the fractional conservation law and degenerate convection diffusion equations and developed some error estimation. Unfortunately, their numerical results failed to confirm their analysis.

The discontinuous Galerkin method is a well established method for classical conservation laws [28]. However, for equations containing higher order spatial derivatives, discontinuous Galerkin methods cannot be directly applied[15, 39]. A careless application of the discontinuous Galerkin method to a problem with high order derivatives could yield an inconsistent method. The idea of local discontinuous Galerkin methods [15] for time dependent partial differential equations with higher derivatives is to rewrite the equation into a first order system and then apply the discontinuous Galerkin method to the system. A key ingredient for the success of this methods is the correct design of interface numerical fluxes. These fluxes must be designed to guarantee stability and local solvability of all the auxiliary variables introduced to approximate the derivatives of the solution.

In this paper, we will consider fractional convection diffusion equations with a fractional Laplacian operator of order . Especially for , it is conceptionally similar to a fractional derivative with an order between 1 and 2. To obtain a consistent and high accuracy method for this problem, we shall rewrite the fractional operator as a composite of first order derivatives and a fractional integral, and convert the fractional convection diffusion equation into a system of low order equations. This allows us to apply the local discontinuous Galerkin method.

This paper is organized as follows. In Section 2, we introduce some basic definitions and state a few central lemmas. In Section 3, we will derive the discontinuous Galerkin formulation for the fractional convection-diffusion law. In section 4, we develop stability and convergence analysis for fractional diffusion and convection-diffusion equations. In Section 5 some numerical examples are carried out to support our analysis. Conclusions are offered in Section 6.

## 2 Background and definitions

Apart from the definitions of the fractional Laplacian based on the Fourier and integral form above, it can also be defined using ideas of fractional calculus[22, 33, 40], as

 −(−Δ)α/2u(x)=∂α∂|x|αu(x)=−−∞Dαxu(x)+xDα∞u(x)2cos(απ2) (4)

where and refer to the left and right Riemann-Liouville fractional derivatives, respectively, of ’th order. This definition is also known as a Riesz derivative. In this paper, we will based our developments and analysis mainly on this definition and will, in preparation, introduce a few definitions and properties of fractional integral/derivative to set the stage.

The left and right fractional integrals of order are defined as

 −∞Iαxu(x)=1Γ(α)∫x−∞(x−t)α−1u(t)dt (5)
 xIα∞u(x)=1Γ(α)∫∞x(t−x)α−1u(t)dt (6)

This allows for the definition of the left and right Riemann-Liouville fractional derivative of order as,

 −∞Dαxu(x)=1Γ(n−α)(ddx)n∫x−∞(x−t)n−1−αu(t)dt=Dn(−∞In−αxu(x)) (7)
 xDα∞u(x)=1Γ(n−α)(−ddx)n∫∞x(t−x)n−1−αu(t)dt=(−D)n(xIn−α∞u(x)) (8)

Mirroring this, Caputo’s left and right fractional derivatives of order are defined as

 C−∞Dαxu(x)=1Γ(n−α)∫x−∞(x−t)n−1−α(ddt)nu(t)dt=−∞In−αx(Dnu(x)) (9)
 CxDα∞u(x)=1Γ(n−α)∫∞x(t−x)n−1−α(−ddt)nu(t)dt=xIn−α∞((−D)nu(x)) (10)

The two definitions, the Riemann-Liouville fractional derivative and Caputo’s fractional derivatives, are naturally related and they are equivalent provided all derivatives less than order of disappear when .

Fractional integrals and derivatives satisfy the following properties,

{lemma}

(Linearity[34]).

 −∞Iαx(λf(x)+μg(x))=λ−∞Iαxf(x)+μ−∞Iαxg(x) (11) −∞Dαx(λf(x)+μg(x))=λ−∞Dαxf(x)+μ−∞Dαxg(x) (12)
{lemma}

(Semigroup property[34]). For

 −∞Iα+βxf(x)=−∞Iαx(−∞Iβxf(x))=−∞Iβx(−∞Iαxf(x)) (13)
{lemma}

[34] Suppose ,when ( ), then

 Extra open brace or missing close brace (14) Extra open brace or missing close brace (15)

From Lemma 2, we get, for and a continuous function with

 −(−Δ)α/2u(x) = −12cos(απ/2)d2dx2(−∞I2−αxu(x)+xI2−α∞u(x)) (16) = −12cos(απ/2)ddx(−∞I2−αxdu(x)dx+xI2−α∞du(x)dx) = −12cos(απ/2)(−∞I2−αxd2u(x)dx2+xI2−α∞d2u(x)dx2)

For , we define

 Δ−s/2u(x)=−−∞D−sxu(x)+xD−s∞u(x)2cos((2−s)π2)=−∞Isxu(x)+xIs∞u(x)2cos(sπ2)

When , we have

 −(−Δ)α/2u(x)=d2dx2(Δ(α−2)/2u)=Δ(α−2)/2(d2udx2)=ddx(Δ(α−2)/2dudx) (17)

To carry out the analysis, we introduce appropriate fractional spaces.
{definition}(Left fractional space[23] ). We define the following seminorm

 |u|JαL(R)=∥−∞Dαxu∥L2(R) (18)

and the norm

 ∥u∥JαL(R)=(|u|2JαL(R)+∥u∥2L2(R))12 (19)

and let denote the closure of with respect to .

{definition}

(Right fractional space[23] ). We define the following seminorm

 |u|JαR(R)=∥xDα∞u∥L2(R) (20)

and the norm

 ∥u∥JαR(R)=(|u|2JαR(R)+∥u∥2L2(R))12 (21)

and let denote the closure of with respect to .

{definition}

(Symmetric fractional space[23]) . We define the following seminorm

 |u|JαS(R)=∣∣(−∞Dαxu,xDα∞u)L2(R)∣∣12 (22)

and the norm

 ∥u∥JαS(R)=(|u|2JαS(R)+∥u∥2L2(R))12 (23)

and let denote the closure of with respect to .

Using these definitions, we recover the following result,

 (−∞Iαxu,u)R=(u,xIα∞u)R (24)
{lemma}

[23]

 (−∞Iαxu,xIα∞u)R=cos(απ)|u|2J−αL(R)=cos(απ)|u|2J−αR(R) (25)

From Lemma 2-2, we recover
{lemma}For all ,

 (Δ−su,u)R=|u|2J−sL(R)=|u|2J−sR(R) (26)

Generally, we consider the problem in a bounded domain instead of . Hence, we restrict the definition to the domain .
{definition} Define the spaces as the closures of under their respective norms.

For these fractional spaces, we have the following Theorem[18],
{theorem} If , then is embedded into , and is embedded into both of them.

{lemma}

(Fractional Poincaré-Friedrichs)[23] For , we have

 ∥u∥L2(Ω)⩽C|u|JμL,0,

and for , we have

 ∥u∥L2(Ω)⩽C|u|JμR,0.

From the definition of the left and right fractional integral, we obtain the following result,
{lemma} Suppose the fractional integral is defined in . Let , then

 xIαbf(x)y=b−x=0Iαyg(y) (27)
{lemma}

Suppose is a smooth function . is a discretization of the domain with interval width , is an approximation of in . i, is a polynomial of degree up to order , and . is the degree of the polynomial. Then for , we have

 ∥Δα/2u(x)−Δα/2uh(x)∥L2⩽Chk+1

and for we have

 ∥−(−Δ)α2u(x)+(−Δ)α2uh(x)∥L2⩽Chk+1−n (28)

where is a constant independent of .

{proof}

Case ,

First, we consider the approximation error for a fractional integral .

First recall that from classical approximation theory. Suppose . We only need to consider .

 ri = 1Γ(−α)∫xa(x−s)−α−1O(hk+1)ds ⩽ (b−a)−αO(hk+1)Γ(1−α)

Recalling Lemma 2, the case is proved.

The case , recall that and use the result for case , we obtain,

 ∥−(−Δ)α2u(x)+(−Δ)α2uh(x)∥L2⩽Chk+1−n (29)

This proves the Lemma.
{lemma}(Inverse properties[11]) Suppose is a finite element space spanned by polynomials up to degree . For any , there exists a positive constant in dependent of and , such that

 ∥∂xuh∥⩽Ch−1∥uh∥,∥uh∥Γh⩽Ch−12∥uh∥

## 3 LDG scheme for the fractional convection-diffusion equation

Let us consider the fractional convection-diffusion equation with . From (4) and Lemma 2, we know that a direct approximation of the fractional laplacian operator is of order . To obtain a high order discontinuous Galerkin scheme for the fractional derivative, we rewrite the fractional derivative as a composite of first order derivatives and fractional integral and convert the equation to a low order system.

Following (17), we introduce two variables , and set

 q = Δ(α−2)/2p p = √ε∂∂xu

Then, the fractional convection-diffusion problem can be rewritten as

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂u∂t+∂∂xf(u)=√ε∂∂xqq=Δ(α−2)/2pp=√ε∂∂xu

Consider the equation in . Given a partition , we denote the mesh by , .

We consider the solution in a polynomial space , which is certainly embedded into the fractional space according to theorem 2. The piece-wise polynomial space on the mesh is defined as,

 Vh={v:v∈Pk(Ij),x∈Ij}

We seek an approximation to such that, for any , we have

 Unknown environment '% (30)

To complete the LDG scheme, we introduce some notations and the numerical flux. Define

 u±(xj)=limx→x±ju(x),¯ui=u++u−2,\llbracketu\rrbracket=u+−u−

and the numerical flux as,

 ^u=hu(u−,u+),^q=hq(q−,q+),^fh=^f(u−h,u+h)

For the high order derivative part, a good choice is the alternating direction flux[15, 39], defined as,

 ^ui+12=u−i+12,^qi+12=q+i+12,0⩽i⩽K−1

or

 ^ui=u+i+12,^qi=q−i+12,1⩽i⩽K

For the nonlinear part, , any monotone flux can be used [28].

Applying integration by parts to (30), and replacing the fluxes at the interfaces by the corresponding numerical fluxes, we obtain,

 ((uh)t,v)Ii+(^fhv−√ε^qhv)|x−i+12x+i−12−(f(uh)−√εqh,vx)Ii = 0 (31) (qh,w(x))Ii−(Δ(α−2)/2ph,w(x))Ii = 0 (32) (ph,z(x))Ii−√ε^uhz|x−i+12x+i−12+√ε(uh,zx)Ii = 0 (33) (uh(x,0),v(x))Ii−(u0(x),v(x))Ii = 0 (34)
###### Remark \thetheorem

Originally, the problem is defined in . However, for numerical purposes, we assume there existing a domain in which  has compact support and restrict the problem to this domain . As a consequence, we impose homogeneous Dirichlet boundary conditions for all . to recover

For the flux at the boundary, we use the flux introduced in [10], defined as

 ^uK+12=u(b,t),^qK+12=q−K+12+βh\llbracketuK+12\rrbracket

or

 ^u12=u(a,t),^q12=q−12+βh\llbracketu12\rrbracket

where is a positive constant.

## 4 Stability and error estimates

In the following we shall discuss stability and accuracy of the proposed scheme, both for the fractional diffusion problem and the more general convection-diffusion problem.

### 4.1 Stability

In order to carry on analysis of the LDG scheme, we define

 B(u,p,q;v,w,z) = ∫T0K∑i=1(ut,v)Iidt+∫T0K∑i=1(^fv−√ε^qv)|x−i+12x+i−12dt (36) −∫T0K∑i=1(f(u)−√εq,vx)Iidt+∫T0K∑i=1(q,w(x))Iidt −∫T0K∑i=1(Δ(α−2)/2p,w(x))Iidt−∫T0K∑i=1√ε^uz|x−i+12x+i−12dt +∫T0K∑i=1(p,z(x))Iidt+∫T0K∑i=1√ε(u,zx)Iidt−∫T0L(v,w,z)dt

where contains the boundary term, defined as

 L(v,w,z)=√εu(a,t)z+12−√εβhu(b,t)v−K+12dt+√εu(b,t)z−K+12=0 (37)

If is a solution, then for any . Considering the fluxes and the flux at the boundaries we recover,

 B(u,p,q;v,w,z) = ∫T0K∑i=1(ut,v)Iidt−∫T0K∑i=1(f(u),vx)Iidt (38) +∫T0K∑i=1(√εq,vx)Iidt+∫T0K∑i=1√ε(u,zx)Iidt+∫T0K∑i=1(q,w(x))Iidt −∫T0K∑i=1(Δ(α−2)/2p,w(x))Iidt+∫T0K∑i=1(p,z(x))Iidt −∫T0K−1∑i=1^fi+12\llbracketv\rrbracketi+12dt+∫T0K−1∑i=1√εq+i+12\llbracketv\rrbracketi+12dt +∫T0K−1∑i=1√εu−i+12\llbracketz\rrbracketi+12dt−∫T0(^f12v+12−^fK+12v−K+12)dt +∫T0√ε(q+12v+12−q−K+12v−K+12)dt+∫T0√εβhu−K+12v−K+12dt
{lemma}

Set in (38), and define . Then the following result holds,

 B(u,p,q;u,−p,q) = ∥u(x,T)∥2−∥u0∥2+∫T0(Δ(α−2)/2p,p)dt+∫T0√εβh(u−K+12)2dt +∫T0Φ(u)12−Φ(u)K+12−(^fu)12+(^fu)K+12dt +∫T0K−1∑j=1(\llbracketΦ(u)\rrbracketj+12−^f\llbracketu\rrbracketj+12)dt
{proof}

Set in (38), and consider the integration by parts formula , to recover at an interface

 K∑i=1(√εq,vx)Ii+K∑i=1(√εu,zx)Ii+K−1∑i=1√εq+i+12\llbracketv\rrbracketi+12+K−1∑i=1√εu−i+12\llbracketz\rrbracketi+12 = K∑i=1√ε(uq)|x−i+12x+i−12+K−1∑i=1√εq+i+12\llbracketu\rrbracketi+12+K−1∑i=1√εu−i+12\llbracketq\rrbracketi+12 = −√εq+12u+12+√εq−K+12u−K+12

Then

 B(u,p,q;u,−p,q) = ∫T0K∑i=1(ut,u)Iidt−∫T0K∑i=1(f(u),ux)Iidt (39) −∫T0K∑i=1(Δ(α−2)/2p,p)Iidt−∫T0K−1∑i=1^fi+12\llbracketv\rrbracketi+12dt −∫T0(^f12u+12−^fK+12u−K+12)dt+∫T0√εβh(u−K+12)2dt

Define , then

 K∑i=1(f(u),ux)Ii=K∑i=1Φ(x)|x−i+12x+i−12=−K−1∑i=1\llbracketΦ(u)\rrbracketi+