Finite difference/local discontinuous Galerkin method for solving the fractional diffusion-wave equation

# Finite difference/local discontinuous Galerkin method for solving the fractional diffusion-wave equation

Leilei Wei111Corresponding author. E-mail addresses: leileiwei09@gmail.com.
College of Science, Henan University of Technology, Zhengzhou, Henan 450001, P.R. China
\CJKtilde

Abstract: In this paper a finite difference/local discontinuous Galerkin method for the fractional diffusion-wave equation is presented and analyzed. We first propose a new finite difference method to approximate the time fractional derivatives, and give a semidiscrete scheme in time with the truncation error , where is the time step size. Further we develop a fully discrete scheme for the fractional diffusion-wave equation, and prove that the method is unconditionally stable and convergent with order , where is the degree of piecewise polynomial. Extensive numerical examples are carried out to confirm the theoretical convergence rates.

Key words: Fractional diffusion-wave equation; Time fractional derivative; Local discontinuous Galerkin method; Stability.

Mathematics Subject Classi?cation: 65M12; 65M06; 35S10

## 1 Introduction

Fractional calculus, which might be considered as an extension of classical calculus, attracts much attention in recent decades. Fractional order partial differential equations (FPDEs) have been frequently used to solve many scientific problems in various fields, such as quantitative finance, engineering, biology, chemistry, hydrology, and so on [1, 15, 18, 19, 31, 42, 43].

However, analytical solutions for the majority of fractional partial differential equations, which are too complex and cannot expressed explicitly, are very difficult to be applied in the science and engineering, so it is a good choice to use numerical methods to finding numerical solutions for fractional partial differential equations, and has very important theoretical and practical significance. The existed methods solving the FPDEs include finite difference methods [2, 6, 9, 10, 13, 22, 25, 26, 27, 37, 38, 46, 47, 49], finite element methods[8, 11, 12, 16, 17, 33, 50], spectral methods[3, 23, 24], discontinuous Gakerkin methods [51, 52], homotopy perturbation method and the variational method [14, 29, 32, 36, 44, 34, 48].

In this paper we consider the following fractional diffusion-wave equation

 ∂αu(x,t)∂tα−∂2u(x,t)∂x2=f(x,t),      (x,t)∈[a,b]×[0,T],u(x,0)=u0(x),    ∂u(x,0)∂t=u1(x),  x∈[a,b], (1.1)

where is a parameter describing the order of the fractional time, are given smooth functions. We do not pay attention to boundary condition in this paper; hence the solution is considered to be either periodic or compactly supported.

The time fractional derivative in the equation (1.1), uses the Caputo fractional partial derivative of order , defined as [10]

 ∂αu(x,t)∂tα=1Γ(2−α)∫t0∂2u(x,s)∂s2ds(t−s)α−1, t>0, 1<α<2, (1.2)

where is the Gamma function.

The fractional diffusion-wave equation is obtained by replacing the first- or second-order time derivative of the classical diffusion or wave equation with a fractional derivative of order , and can be used to interpolate the diffusion equation and wave equation and model many of the mechanical responses and acoustics accurately.

The rest of this paper is constructed as follows. In the section 2 some basic notations and theoretic results are introduced. Then in section 3 we construct our finite difference/discontinuous Galerkin method for the fractional diffusion-wave equation, and stability and error analysis are given. Numerical results are presented in section 4, and the concluding remarks is included in the final section.

## 2 Notations and auxiliary results

In this section we introduce some notations and definitions that will be used later in the following sections.

Let be a finite domain, and a partition is given by

 a=x12

we denote the cell by for and the cell lengths .

We denote by and the values of at , from the right cell and from the left cell ,respectively.

The piecewise-polynomial space is defined as the space of polynomials of the degree up to in each cell , i.e.

 Vkh={v:v∈Pk(Ij),x∈Ij,j=1,2,⋯N}.

For error estimates, we will be using two projections in one dimension , denoted by , i.e., for each ,

 ∫Ij(Pω(x)−ω(x))v(x)=0,∀v∈Pk(Ij), (2.1)

and special projection , i.e., for each ,

 ∫Ij(P+ω(x)−ω(x))v(x)=0,∀v∈Pk−1(Ij), and P+ω(x+j−12)=ω(xj−12) ∫Ij(P−ω(x)−ω(x))v(x)=0,∀v∈Pk−1(Ij), and P−ω(x−j+12)=ω(xj+12). (2.2)

For the above projections and , we have [5, 35, 40, 41]

 ∥ωe∥+h∥ωe∥∞+h12∥ωe∥τh≤Chk+1, (2.3)

where or .

The notations are used: the scalar inner product on be denoted by , and the associated norm by . If , we drop . In the present paper we use to denote a positive constant which may have a different value in each occurrence.

## 3 The schemes

In this section, we first present a finite difference method to approximate the time fractional derivatives, and then give the implicit fully discrete scheme with space discretized by the local discontinuous Galerkin method. Stability and convergence are detailed analysis.

### 3.1 Time fractional derivative discretization

We divide the interval uniformly with a time step size , , be the mesh points.

Let , and from the fact

 v(x,ti)=∂u(x,ti)∂t=3u(x,ti)−4u(x,ti−1)+u(x,ti−2)2Δt+rn1,

where the truncation error , we can obtain

 ∂αu(x,tn)∂tα =1Γ(2−α)∫tn0∂v(x,s)∂sds(tn−s)α−1 (3.1) =1Γ(2−α)n−1∑i=0∫ti+1ti∂v(x,s)∂sds(tn−s)α−1 =1Γ(2−α)n−1∑i=0∫ti+1tiv(x,ti+1)−v(x,ti)Δtds(tn−s)α−1+rn2 =(Δt)2−αΓ(3−α)n−1∑i=0bn−i−1v(x,ti+1)−v(x,ti)Δt+rn2 =(Δt)1−αΓ(3−α)[v(x,tn)+n−1∑i=1(bn−i−bn−i−1)v(x,ti)−bn−1v(x,t0)]+rn2 =(Δt)1−αΓ(3−α)[3u(x,tn)−4u(x,tn−1)+u(x,tn−2)2Δt +n−1∑i=1(bn−i−bn−i−1)3u(x,ti)−4u(x,ti−1)+u(x,ti−2)2Δt −bn−1v(x,t0)]+rn3,

where

 b0=1,    bi=(i+1)2−α−i2−α,i=1,2,3,⋯

when , we take by Taylor expansion.

Similar to the proof in [24],the truncation error , so satisfied

 |rn3|≤C(Δt)3−α.

It is easy to check that

 bi>0,i=1,2⋯,n.1=b0>b1>b2>⋯>bn,bn→0(n→∞). (3.2)

Substituting (3.1) into (1.1), we have

 3u(x,tn)−β∂2(x,tn)∂x2= n−1∑i=1(bn−i−1−bn−i)(3(x,ti)−4u(x,ti−1)+u(x,ti−2)) +2Δtbn−1v(x,t0)+βf(x,tn)+4u(x,tn−1) −u(x,tn−2)+βrn3,

where .

Let be the numerical approximation to , , the problem (1.1) can be discretized by the following scheme

 3un−β∂2un∂x2= n−1∑i=1(bn−i−1−bn−i)(3ui−4ui−1+ui−2) (3.3) +2Δtbn−1v0+βfn+4un−1−un−2,

where We know

 |βrn3|≤C(Δt)3,

however, by Taylor expansion we have

 |u(x,t−1)−u(x,0)+Δtu1(x)|≤C(Δt)2,

therefore the truncation error is in scheme (3.3).

### 3.2 Fully discrete schemes

In this subsection we present the fully discrete LDG scheme for the problem (1.1) based on the semidiscrete scheme (3.3).

We rewrite Eq. (1.1) as a first-order system:

 p=ux,      ∂αu(x,t)∂tα−px=f(x,t). (3.4)

Let be the approximations of , respectively, . We define a fully discrete local discontinuous Galerkin scheme as follows: find such that for all test functions

 3∫Ωunhϕdx+β(∫Ωpnhϕxdx−N∑j=1((ˆpnhϕ−)j+12−(ˆpnhϕ+)j−12))=n−1∑i=1(bn−i−1−bn−i)∫Ω(3uih−4ui−1h+ui−2h)ϕdx+2Δtbn−1∫Ωv0hϕdx+4∫Ωun−1hϕdx−∫Ωun−2hϕdx+β∫Ωfnϕdx,∫Ωpnhwdx+∫Ωunhwxdx−N∑j=1((ˆunhw−)j+12−(ˆunhw+)j−12)=0, (3.5)

The initial conditions are taken as the projections of , respectively,

 (3.6)

The “hat” terms in (3.5) in the cell boundary terms from integration by parts are the so-called “numerical fluxes”, which are single valued functions defined on the edges and should be designed based on different guiding principles for different PDEs to ensure stability. It turns out that we can take the simple choices such that

 ˆunh=(unh)−,  ˆpnh=(pnh)+. (3.7)

We remark that the choice for the fluxes (3.7) is not unique. In fact the crucial part is taking and from opposite sides [40, 5].

### 3.3 Stability and Convergence

In order to simplify the notations and without lose of generality, we consider the case in its numerical analysis.

###### Theorem 3.1.

For periodic or compactly supported boundary conditions, the fully-discrete LDG scheme (3.5) is unconditionally stable, and there exists a positive constant depending on , such that

 ∥unh∥≤C(∥u0h∥+Δt∥u1(x)∥),     n=1,2⋯,M. (3.8)
###### Proof.

Taking in scheme (3.5), we obtain

 3∥unh∥2+β∥pnh∥2+βN∑j=1(Ψ(unh,pnh)j+12−Ψ(unh,pnh)j−12+Θ(unh,pnh)j−12)
 =n−1∑i=1(bn−i−1−bn−i)∫Ω(3uih−4ui−1h+ui−2h)unhdx+2Δtbn−1∫Ωv0hunhdx+4∫Ωun−1hunhdx−∫Ωun−2hunhdx, (3.9)

where

 Ψ(unh,pnh)=(pnh)−(unh)−−ˆpnh(unh)−−ˆunh(pnh)−,Θ(unh,pnh)=(pnh)−(unh)−−(pnh)+(unh)+−ˆpnh(unh)−+ˆpnh(unh)+−ˆunh(pnh)−+ˆunh(pnh)+.

If we take the fluxes (3.7), after some manual calculation, we can easily obtain

Then based on the equation (3.9), we can get

 3∥unh∥2+β∥pnh∥2=n−1∑i=1(bn−i−1−bn−i)∫Ω(3uih−4ui−1h+ui−2h)unhdx+2Δtbn−1∫Ωv0hunhdx+4∫Ωun−1hunhdx−∫Ωun−2hunhdx≤n−1∑i=1(bn−i−1−bn−i)(3∥uih∥+4∥ui−1h∥+∥ui−2h∥)∥unh∥+2Δtbn−1∥v0h∥∥unh∥+4∥un−1h∥∥unh∥+∥un−2h∥∥unh∥,

that is

 3∥unh∥≤n−1∑i=1(bn−i−1−bn−i)(3∥uih∥+4∥ui−1h∥+∥ui−2h∥)+2Δtbn−1∥v0h∥+4∥un−1h∥+∥un−2h∥. (3.10)

We will prove the Theorem 3.1 by mathematical induction. When , we can obtain

 3∥u1h∥≤2Δt∥v0h∥+4∥u0h∥+∥u−1h∥ (3.11)

Notice that

 ∫Iju−1hvdx=∫IjP(u(x,0)−Δtu1(x))vdx=∫Iju0hvdx−Δt∫Iju1(x)vdx,

for any . Taking , we can obtain

 ∥u−1h∥2Ij=∫Iju0hu−1hdx−Δt∫Iju1(x)u−1hdx≤∥u0h∥2Ij+14∥u−1h∥2Ij+(Δt)2∥u1(x)∥2Ij+14∥u−1h∥2Ij,

summing over from to , we can get

 ∥u−1h∥≤C(∥u0h∥+Δt∥u1(x)∥). (3.12)

Similar to the proof of (3.12), we can easily obtain

 ∥v0h∥≤∥u1(x)∥. (3.13)

By using (3.11),(3.12) and (3.13), it is easily to know that there exists a positive constant , such that

 ∥u1h∥≤C(∥u0h∥+Δt∥u1(x)∥). (3.14)

Now suppose the following inequality holds

 ∥umh∥≤C(∥u0h∥+Δt∥u1(x)∥),m=2,3⋯K, (3.15)

we need to prove

Let in the inequality (3.10), we can obtain

 3∥uK+1h∥≤K∑i=1(bK−i−bK+1−i)(3∥uih∥+4∥ui−1h∥+∥ui−2h∥)+2ΔtbK∥v0h∥+4∥uKh∥+∥uK−1h∥.

Using (3.12), (3.13) and (3.15), we can obtain the following inequality easily

 ∥uK+1h∥≤C(∥u0h∥+Δt∥u1(x)∥).

This finishes the proof of the stability result. ∎

###### Theorem 3.2.

Let be the exact solution of problem (1.1), which is sufficiently smooth such that with . Let be the numerical solution of the fully discrete LDG scheme (3.5), then there holds the following error estimate:

 ∥u(x,tn)−unh∥≤C(hk+1+(Δt)2),n=1,⋯,M, (3.16)

where C is a constant depending on .

###### Proof.

By Taylor expansion we know

 |u(x,t−1)−u(x,0)+Δtu1(x)|≤C(Δt)2,

here is a positive constant depending on . Then by using the property (2.3), we can obtain the following estimate which will be used later,

 ∥u(x,t−1)−u−1h∥≤C((Δt)2+hk+1). (3.17)

It is easy to verify that the exact solution of PDE (1.1) satisfies

 3∫Ωu(x,tn)ϕdx+β(∫Ωp(x,tn)ϕxdx−N∑j=1((p(x,tn)ϕ−)j+12−(p(x,tn)ϕ+)j−12))=n−1∑i=1(bn−i−1−bn−i)∫Ω(3u(x,ti)−4u(x,ti−1)+u(x,ti−2))ϕdx+2Δtbn−1∫Ωv(x,t0)ϕdx+4∫Ωu(x,tn−1)ϕdx−∫Ωu(x,tn−2)ϕdx+β∫Ωf(x,tn)ϕdx+β∫Ωrn3ϕdx,∫Ωp(x,tn)wdx+∫Ωu(x,tn)wxdx−N∑j=1((u(x,tn)w−)j+12−(u(x,tn)w+)j−12)=0, (3.18)

, for .

Denote

 enu=u(x,tn)−unh=P−enu−(P−u(x,tn)−u(x,tn)),enp=p(x,tn)−pnh=P+enp−(P+p(x,tn)−p(x,tn)). (3.19)

Subtracting (3.5) from (3.18), and with the fluxes (3.7) we can obtain the error equation:

 3∫Ωenuϕdx+β(∫Ωenpϕxdx−N∑j=1(((enp)+ϕ−)j+12−((enp)+ϕ+)j−12))−n−1∑i=1(bn−i−1−bn−i)∫Ω(3eiu−4ei−1u+ei−2u)ϕdx−4∫Ωen−1uϕdx+∫Ωen−2uϕdx+β∫Ωrn3ϕdx+∫Ωenpwdx+∫Ωenuwxdx−N∑j=1(((enu)−w−)j+12−((enu)−w+)j−12)=0. (3.20)

Using (3.19), the error equation (3.20) can be written as follows:

 3∫ΩP−enuϕdx+β(∫ΩP+enpϕxdx−N∑j=1(((P+enp)+ϕ−)j+12−((P+enp)+ϕ+)j−12))+∫ΩP+enpwdx+∫ΩP−enuwxdx−N∑j=1(((P−enu)−w−)j+12−((P−enu)−w+)j−12)=n−1∑i=1(bn−i−1−bn−i)∫Ω(3P−eiu−4P−ei−1u+P−ei−2u)ϕdx+4∫ΩP−en−1uϕdx−∫ΩP−en−2uϕdx−β∫Ωrn3ϕdx
 −n−1∑i=1(bn−i−1−bn−i)∫Ω(3(P−u(x,ti)−u(x,ti))−4(P−u(x,ti−1)−u(x,ti−1))+(P−u(x,ti−2)−u(x,ti−2)))ϕdx−4∫Ω(P−u(x,tn−1)−u(x,tn−1))ϕdx+∫Ω(P−u(x,tn−2)−u(x,tn−2))ϕdx+3∫Ω(P−u(x,tn)−u(x,tn))ϕdx+β(∫Ω(P+p(x,tn)−p(x,tn))ϕxdx−N∑j=1((((P+p(x,tn)−p(x,tn)))+ϕ−)j+12−(((P+p(x,tn)−p(x,tn)))+ϕ+)j−12))+∫Ω(P+p(x,tn)−p(x,tn))wdx+∫Ω(P−u(x,tn)−u(x,tn))wxdx−N∑j=1(((P−u(x,tn)−u(x,tn))−w−)j+12−((P−u(x,tn)−u(x,tn))−w+)j−12). (3.21)

Taking the test functions in (3.21), using the properties (2.1)-(2), then the following equality holds,

Therefore, we obtain

 ∥P−enu)∥≤n−1∑i=1(bn−i−1−bn−i)(3∥P−eiu∥+4∥P−ei−1u∥+∥P−ei−2u∥)+4∥P−en−1u∥+∥P−en−2u∥+β∥rn3∥+n−1∑i=1(bn−i−1−bn−i)(3∥P−u(x,ti)−u(x,ti)∥+4∥P−u(x,