Local discontinuous Galerkin method for a nonlinear time-fractional fourth-order partial differential equation

Accepted Manuscript Local discontinuous Galerkin method for a nonlinear time-fractional fourth-order partial differential equation

Yanwei Du, Yang Liu, Hong Li, Zhichao Fang, Siriguleng He

PII: DOI: Reference:

S0021-9991(17)30366-2 http://dx.doi.org/10.1016/j.jcp.2017.04.078 YJCPH 7344

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

12 October 2015 17 April 2016 30 April 2017

Please cite this article in press as: Y.W. Du et al., Local discontinuous Galerkin method for a nonlinear time-fractional fourth-order partial differential equation, J. Comput. Phys. (2017), http://dx.doi.org/10.1016/j.jcp.2017.04.078

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Local discontinuous Galerkin method for a nonlinear time-fractional fourth-order partial differential equation

∗

Yanwei Du, Yang Liu† , Hong Li, Zhichao Fang, Siriguleng He (School of Mathematical Sciences, Inner Mongolia University, Hohhot, 010021, China)

Abstract: In this article, a fully discrete local discontinuous Galerkin (LDG) method with high-order temporal convergence rate is presented and developed to look for the numerical solution of nonlinear time-fractional fourth-order partial differential equation (PDE). In the temporal direction, for approximating the fractional derivative with order α ∈ (0, 1), the weighted and shifted Gr¨ unwald difference (WSGD) scheme with second-order convergence rate is introduced and for approximating the integer time derivative, two step backward Euler method with second-order convergence rate is used. For the spatial direction, the LDG method is used. For the numerical theories, the stability is derived and a priori error results are proved. Further, some error results and convergence rates are calculated by numerical procedure to illustrate the effectiveness of proposed method. Keywords: Nonlinear time-fractional fourth-order problem; WSGD scheme; LDG method; High-order scheme; Caputo fractional derivative 2000 Mathematics Subject Classification: 65M60, 65N30

1

Introduction

Fractional differential equations, which mainly cover temporal fractional derivatives [3, 7, 8, 9, 4, 11, 12, 13, 33, 27, 36, 38, 30, 28], spatial fractional derivatives [5, 6, 14, 29, 15, 17, 19, 20, 21, 26, 35, 37] and space-time fractional derivatives [16, 22, 10, 31], have been found in a lot of problems in many scientific fields. So, ones need to look for the analytical or numerical solutions of fractional differential equations and to study their properties. The numerical methods of differential equations with fractional derivatives of Caputo type and Riemann-Liouville type were studied and developed by many authors. In these numerical methods, we can find that many approximate formulas were used to discretize Caputo and Riemann-Liouville fractional derivatives. In [23], Meerschaert and Tadjeran considered some finite difference schemes for advection-dispersion flow equations with Riemann fractional derivative and the first-order approximate result with O(Δt) was obtained. In 2006, Sun and Wu [54] gave some approximations covering (2 − α)order for Caputo fractional derivatives. In 2007, Lin and Xu [24] showed spectral method with finite difference approximation for the time Caputo fractional diffusion equation and gave the detailed proof of L1 approximation with order (2 − α). Up to now, the L1-approximation has been used to get the numerical solutions for fractional PDEs. Lin et al. [25] carried out spectral approximations with L1-formula for some fractional PDEs. In [36], finite difference methods combined with L1-approximations were studied ∗ Foundation item: Supported by the National Natural Science Fund (11301258, 11361035, 11501311), the Natural Science Fund of Inner Mongolia Autonomous Region (2014BS0101), the Scientific Research Projection of Higher Schools of Inner Mongolia(NJZZ12011, NJZY14013). † Corresponding E-mail: [email protected](Y.Liu)

1

and analyzed for solving PDEs with some fractional derivatives. In [34, 16, 28], L1-formulas combined with finite element analysis were shown to look for the numerical solutions for some temporal and spacetime fractional PDEs. Li et al. [30] considered finite element method with L1-approximation for time fractional maxwell’s equations. Finite volume (element) methods were also provided for fractional PDEs with Riemann-Liouville fractional derivative in [18, 19] and Caputo fractional derivative in [14]. In the literatures [56, 57, 49], some (local) discontinuous Galerkin method were discussed and developed by some authors. Recently, some high-order approximations for fractional derivatives of Caputo type or Riemann-Liouville type were proposed and discussed in [1, 2, 39, 40, 41, 50, 52, 33, 32]. By assembling the Gr¨ unwald-Letnikov difference operators with different weights and shifts, Tian et al. [1] presented some high-order approximations to discretize the Riemann-Liouville space-fractional derivatives. Following this idea in [1], Wang and Vong [2] provided a second-order accuracy formula to approximate the time-fractional derivative and established a compact finite difference scheme for solving the modified anomalous fractional diffusion equations. But, we find that these high-order schemes [1, 2] have not been used in the LDG method. LDG method proposed by Cockburn and Shu [58] is a kind of important numerical methods and concerned by a lot of scholars. So far, many integer order PDEs [59, 60, 61, 62, 63, 64, 65, 66] have been solved by LDG methods. Recently, some fractional partial differential equations [55, 56, 57, 49] have been studied by LDG methods. In [55], Deng and Hesthaven applied LDG methods to solve fractional ordinary differential equations. In [57], authors discussed LDG approximation for time-fractional Schr¨ odinger equation. In [56] and [49], authors gave some detailed analysis on LDG methods for time-fractional fourth-order problems with L1-approximation, respectively. Here, we will use the LDG methods to solve the following nonlinear fourth-order problem with time Caputo fractional derivative ∂u(x, t) ∂ α u(x, t) ∂ 4 u(x, t) + + + f (u(x, t)) = g(x, t), ∂t ∂tα ∂x4 u(x, 0) = u0 (x), x ∈ [a, b], α

(x, t) ∈ [a, b] × (0, T ],

(1.1)

u(x,t) with α ∈ (0, 1) is the fractional derivative of Caputo type defined later in the next section. where ∂ ∂t α The source term g and the initial value u0 are two given smooth functions. The boundary condition is not given since the solution is considered to be periodic. The nonlinear reaction term f (u) satisfies the assumed conditions: |f (u)| ≤ C|u|; The first-order derivative of f (u) with respect to u is bounded, that is to say, there exists a positive constant C such that |f  (u)| ≤ C. The important applications of the fourth-order model equations [44, 45] can be found in many scientific and engineering fields, such as modeling formation of grooves on a flat surface [42, 50], wave propagation in beams [43], mathematical modeling of fourth-order subdiffusion systems [46, 50, 51, 52] and so forth. Because it is hard to find out the exact solutions of fourth-order fractional PDEs by the analytic methods, more and more fourth-order fractional problems have been solved by some numerical methods such as finite difference methods [46, 47, 48, 50], finite element methods [51, 52, 53], LDG methods [56, 49] and so on. In [56], Wei and He considered a fully discrete LDG method for solving linear fourth-order problems with time-fractional derivative, and obtained the temporal convergence order with O(Δt2−α ). Very recently, Guo, Wang and Vong [49] also discussed the analysis of fully discrete LDG methods for some linear time-fractional fourth-order problems and obtained the temporal convergence results with O(Δt2−α ). Here our purpose is to discuss a higher-order approximation (called WSGD method [1, 2]) for temporal fractional derivative combined with LDG method for nonlinear fractional fourth-order partial differential equation. Compared with the results obtained in [56, 49], the current result in the temporal direction is O(Δt2 ). We derive and analyze some stable results and a priori estimates for LDG method, then provide some numerical results to verify our analysis for our convergence results covering second-order convergence rate. In this article, we use C to denote a positive constant which may have a different value in each occurrence.

2

The organizational structure of this article is as follows. We introduce the definitions of fractional derivatives, approximate formula and some lemmas in section 2. In section 3, we give the LDG scheme for the Eq.(1.1) and prove the stability. In section 4, we derive the high-order error results of LGD method. In section 5, we show the detailed numerical algorithm and calculate the errors and convergent rates to confirm our theoretical analysis. In section 6, we make some concluding remarks.

2

Fractional derivatives and some lemmas

2.1

Fractional derivatives

In many literatures, we can get the following definitions on fractional derivatives of Caputo type and Riemann-Liouville type. Definition 2.1 :The α-order (0 < α < 1) fractional derivative of Caputo type for the function u(t) is defined by ∂ α u(t) . = ∂tα

C α 0 ∂t u(t)

1 Γ(1 − α)

=



t 0

u (τ ) dτ, (t − τ )α

(2.1)

where Γ(·) is Gamma function. Definition 2.2 :The α-order (0 < α < 1) fractional derivative of Riemann-Liouville type for the function u(t) is defined as α 0 ∂t u(t) =

d 1 Γ(1 − α) dt



t 0

u(τ ) dτ. (t − τ )α

(2.2)

Lemma 2.3 [39] The relationship between Caputo fractional derivative and Riemann-Liouville fractional derivative can be given by α 0 ∂t u(t)

α =C 0 ∂t u(t) +

u(0)t−α . Γ(1 − α)

(2.3)

α From the relationship (2.3), we easily see that 0 ∂tα u(t) =C 0 ∂t u(t) when u(0) = 0.

2.2

Time-fractional derivative’s approximation

To give the fully discrete analysis of Eq.(1.1), we should approximate both integer and fractional derivatives. The grid points in the time interval [0, T ] are labeled as ti = iτ , i = 0, 1, 2, . . . , M , where τ = T /M is the time step parameter. We define z n = z(tn ) for a smooth function on [0, T ]. For 0 < α < 1, when the initial value u0 (x) = 0, we use lemma 2.3 and do a simple calculation at t = tn+1 to get the approximate formula [1, 2] with second-order accuracy n+1 ∂ α z(x, tn+1 )  b(i) = z(x, tn+1−i ) + O(τ 2 ), α ∂tα τ i=0

(2.4)

where b(i) =

⎧ ⎨ α+2 g α , 2 2

if i = 0,

0

⎩ α+2 g α + i

−α α 2 gi−1 ,

(2.5)

if i > 0,

(α)

and gi = Γ(i − α)/(Γ(−α)Γ(i + 1)). Then, we present two lemmas to be used later in our article.

3

Lemma 2.4 [1, 2] Let {b(i)}∞ i=1 be defined as in (2.5). Then for any positive integer L and real vector (v 0 , v 1 , · · · , v L ) ∈ RL+1 , it holds that N  n  n=0

 b(i)v n−i v n ≥ 0.

(2.6)

i=0

Lemma 2.5 Let {b(i)}∞ i=1 be defined as in (2.5). Then for any positive integer i, b(i) is bounded. Proof. Using the definition of giα and natures of Gamma function, we can easily get giα =

⎧ ⎨1,

if i = 0,

⎩ i−1−α g α , i−1 i

if i > 0.

(2.7)

Since | i−1−α | < 1, we have |giα | ≤ 1. According to the definition of b(i), we know that b(i) is bounded. i Remark 2.6 In the current literatures on Caputo fractional derivatives covering order α, 0 < α < 1, most of them are approximated by L1-formula with (2 − α)-order convergence rate. So the studies of the high-order approximate methods for fractional derivatives will be important and significant. Based on the discussions in Refs. [1, 2], we can find that a second-order convergence result for fractional derivative was obtained and the related numerical theories and results were discussed by using finite difference methods. However, the related numerical theories based on finite element methods and LDG methods have not been considered. For finite difference methods, the results can be arrived at by considering the positive definite matrix (also see the results in lemma 2.4). However, for finite element methods and LDG methods, derivations of the results are more difficult.

3

The analysis of stability for fully discrete scheme

We assume the considered domain Ω = [a, b] is constituted by cells Ij = [xj− 12 , xj+ 12 ], for j = 1, . . . N , where a = x 12 < x 32 < . . . < xN + 12 = b, and denote the cell lengths Δxj = xj+ 12 − xj− 12 , 1 ≤ j ≤ N , h = max Δxj . At xj+ 12 , we denote the values of u by u+ and u− from the right cell Ij+1 and the j+ 1 j+ 1 1≤j≤N

2

2

− u− by [u]j+ 12 . left cell Ij , respectively, and the value of u+ j+ 12 j+ 12 In this section, we introduce the numerical scheme for solving Eq.(1.1). By the introduced auxiliary variables, the problem (1.1) can be rewritten as the first-order system: p =ux , q =px , (3.1)

s =qx , ∂ α u(x, t) ∂u + sx + f (u(x, t)) =g(x, t). + ∂tα ∂t

Based on the system (3.1), by considering initial value u0 (x) = 0, we can obtain the following weak form:

4

Case I: n = 0  1  b(i) τα

i=0

N 

+

j=1

Ω

(s(x, t1 )v − )j+ 12 



E0 vdx +

+ Ω



Ω

 Ω

E1 vdx +



s(x, t1 )wdx +

Ω

 u(x, t1 ) − u(x, t0 ) vdx − s(x, t1 )vx dx τ Ω Ω   − (s(x, t1 )v + )j− 12 = (−f (u(x, t0 )))vdx + g(x, t1 )vdx 

u(x, t1−i )vdx +



Ω

E2 vdx,

q(x, t1 )wx dx −

Ω

q(x, t1 )ρdx +



Ω

p(x, t1 )ρx dx −

Ω

Ω

N  j=1

 p(x, t1 )ξdx +

N  j=1

 Ω

Ω

u(x, t1 )ξx dx −

((q(x, t1 )w− )j+ 12 − (q(x, t1 )w+ )j− 12 ) = 0,

(3.2)

((p(x, t1 )ρ− )j+ 12 − (p(x, tn+1 )ρ+ )j− 12 ) = 0,

N  ((u(x, t1 )ξ − )j+ 12 − (u(x, t1 )ξ + )j− 12 ) = 0. j=1

Case II: n ≥ 1 n+1  i=0

+

b(i) τα

N  j=1

 Ω

(s(x, tn+1 )v − )j+ 12 

 +

Ω

 Ω

 3u(x, tn+1 ) − 4u(x, tn ) + u(x, tn−1 ) vdx − s(x, tn+1 )vx dx 2τ Ω Ω  − (s(x, tn+1 )v + )j− 12 = (−2f (u(x, tn )) + f (u(x, tn−1 )))vdx 

u(x, tn+1−i )vdx +

g(x, tn+1 )vdx +

Ω



s(x, tn+1 )wdx +





Ω

E0 vdx +

 Ω

E3 vdx +

q(x, tn+1 )wx dx −

j=1

 Ω

q(x, tn+1 )ρdx +



Ω

 Ω

p(x, tn+1 )ξdx +

Ω

N 

p(x, tn+1 )ρx dx −

N  j=1

Ω

Ω

E4 vdx,

((q(x, tn+1 )w− )j+ 12 − (q(x, tn+1 )w+ )j− 12 ) = 0,

(3.3)

((p(x, tn+1 )ρ− )j+ 12 − (p(x, tn+1 )ρ+ )j− 12 ) = 0,

N  u(x, tn+1 )ξx dx − ((u(x, tn+1 )ξ − )j+ 12 − (u(x, tn+1 )ξ + )j− 12 ) = 0, j=1

where E0 =

n+1  i=0

b(i) ∂ α u(x, tn+1 ) u(x, t ) − = O(τ 2 ), n+1−i τα ∂tα

u(x, t1 ) − u(x, t0 ) − ut (x, t1 ) = O(τ ), τ E2 = f (u1 ) − f (u0 ) = f  (uξ )(u(x, t1 ) − u(x, t0 )) = f  (uξ )(ut (x, tπ )τ ) = O(τ ), E1 =

(3.4)

3u(x, tn+1 ) − 4u(x, tn ) + u(x, tn−1 ) − ut (tn+1 ) = O(τ 2 ), 2τ E4 = f (u(x, tn+1 )) − (2f (u(x, tn )) − f (u(x, tn−1 ))) = O(τ 2 ).

E3 =

For formulating LDG scheme, we define the following piecewise-polynomial space Vhk = {v : v ∈ P k (Ij ), x ∈ Ij , j = 1, 2, . . . N },

(3.5)

where P k (Ij ) denotes the set of polynomials of degree up to k defined on the cell Ij . Let un+1 , pn+1 , qhn+1 , sn+1 ∈ Vhk be the approximation of u(., tn+1 ), p(., tn+1 ), q(., tn+1 ), s(., tn+1 ), h h h 5

respectively. We can define a fully discrete LDG scheme as follow: Case I: n = 0  1  b(i) i=0

+

τα

N  j=1



Ω

 Ω

 Ω

Ω

 u1h − u0h vdx − s1h vx dx τ Ω Ω   + 0 1 − (sh v )j− 12 = (−f (uh ))vdx + g 1 vdx,

uh1−i vdx +

(s 1h v − )j+ 12

s1h wdx + qh1 ρdx + p1h ξdx +

 Ω

 Ω

 Ω



Ω

qh1 wx dx −

N  j=1

p1h ρx dx −

N  j=1

u1h ξx dx −

Ω

((q h1 w− )j+ 12 − (q h1 w+ )j− 12 ) = 0,

(3.6)

1 ρ− ) 1 + ((p

j+ 12 − (ph ρ )j− 12 ) = 0, h

N 

1 ξ − ) 1 − (u

1 ξ + ) 1 ) = 0. ((u j+ 2 j− 2 h h j=1

Case II: n ≥ 1 n+1  i=0

+ 

b(i) τα

N  j=1

Ω

 Ω

 Ω

 Ω

 3un+1 − 4unh + un−1 h h vdx − sn+1 vx dx h 2τ Ω Ω   n+1 + n−1 n 1 = − (s v ) (−2f (u ) + f (u ))vdx + g n+1 vdx, j− 2 h h h

un+1−i vdx + h

n+1 − (s v )j+ 12 h

sn+1 wdx + h qhn+1 ρdx + pn+1 ξdx + h

 Ω

 Ω

 Ω



Ω

qhn+1 wx dx − pn+1 ρx dx − h

un+1 ξx dx − h

N  j=1

N  j=1 N  j=1

Ω

n+1 − n+1 + ((q w )j+ 12 − (q w )j− 12 ) = 0, h h

(3.7)

n+1 − n+1 + ((p ρ )j+ 12 − (p ρ )j− 12 ) = 0, h h

n+1 − n+1 + ((u ξ )j+ 12 − (u ξ )j− 12 ) = 0. h h

Based on the system (3.6)-(3.7), we can obtain the following theorem of stability. Theorem 3.1 For the LDG scheme (3.6)-(3.7), the following stable inequality for uL h ∈ Vh holds, for sufficiently small τ  L−1 1  n+1 2 2 uL  ≤ C τ g  , h

(3.8)

n=1

where C is a positive constant which does not depend on time step length Δt and space mesh h. , ξ = −sn+1 , w = pn+1 , ρ = qhn+1 in the scheme (3.7), we can Proof. Taking the test functions v = un+1 h h h

6

obtain n+1 



1 [(un+1 2 + 2un+1 − unh 2 ) − (unh 2 + 2unh − un−1 2 ) h h h 4τ Ω i=0  N  n−1 2 n+1 n+1 n+1 n+1 − n+1 n+1 + n + un+1 − 2u + u  ] − s (u ) dx + (s (uh ) )j+ 12 − (s (uh ) )j− 12 x h h h h h h h b(i) τα

Ω

 = Ω

 Ω

un+1−i un+1 dx + h h

(−2f (unh ) + f (un−1 ))un+1 dx + h h

sn+1 pn+1 dx + h h

qhn+1 2 +  −

Ω

 Ω

 Ω

qhn+1 (pn+1 )x dx − h

pn+1 (qhn+1 )x dx − h

pn+1 sn+1 dx − h h

 Ω

N  j=1

j=1



g n+1 un+1 dx, h

Ω N  j=1

n+1 n+1 − n+1 n+1 + ((q (ph ) )j+ 12 − (q (ph ) )j− 12 ) = 0, h h

(3.9)

n+1 n+1 − n+1 n+1 + ((p (qh ) )j+ 12 − (p (qh ) )j− 12 ) = 0, h h

un+1 (sn+1 )x dx + h h

N  j=1

n+1 n+1 − n+1 n+1 + ((u (sh ) )j+ 12 − (u (sh ) )j− 12 ) = 0. h h

Combining the above four equations in (3.9), we have n+1  i=0

− −

b(i) τα

(unh 2



1 [(un+1 2 + 2un+1 − unh 2 ) h h 4τ

+ 2unh − un−1 2 ) + un+1 − 2unh + un−1 2 ] h h h

N     n+1 − n+1 n+1 + (un+1 + (qhn+1 pn+1 s ) − (u s ) )− − (qhn+1 pn+1 )+ h h h h h h j+ 1 j− 1 j+ 1 j− 1 2

2

2

j=1

2

N  

N     n+1 n+1 − n+1 n+1 + n+1 n+1 − n+1 n+1 +    1 − (s 1 1 − (u 1 (s + ( u (u ) ) (u ) ) (s ) ) (s ) ) j+ j− j+ j− h h h h h h h h 2 2 2 2

j=1

−

Ω

un+1−i un+1 dx + qhn+1 2 + h h

N   j=1

+



j=1

N   j=1

= Ω



n+1 n+1 − n+1 n+1 + (q (ph ) )j+ 12 − (q (ph ) )j− 12 − h h

(−2f (unh ) + f (un−1 ))un+1 dx + h h

 Ω

N   j=1

n+1 n+1 − n+1 n+1 + (p (qh ) )j+ 12 − (p (qh ) )j− 12 h h



g n+1 un+1 dx. h (3.10)

The “hat” terms in (3.2) in the cell boundary terms from integration by parts are the so-called “numerical fluxes”, in order to ensure stability, we can make the following choices n+1 n+1 n+1 n+1 u = (un+1 )− , s = (sn+1 )+ , q = (qhn+1 )− , p = (pn+1 )+ . h h h h h h h

7

(3.11)

Substituting the fluxes (3.11) into (3.10), we have n+1  i=0

− −

b(i) τα



+ 2unh − un−1 2 ) + un+1 − 2unh + un−1 2 ] h h h 2

2

2

j=1

N  

N  

)+ (qhn+1 )− )j+ 12 − ((pn+1 )+ (qhn+1 )+ )j− 12 ((pn+1 h h

N  

)− )j+ 12 − ((qhn+1 )− (pn+1 )+ )j− 12 ((qhn+1 )− (pn+1 h h

j=1

= Ω



)− (sn+1 )− )j+ 12 − ((un+1 )− (sn+1 )+ )j− 12 ((un+1 h h h h

N  

(−2f (unh ) + f (un−1 ))un+1 dx + h h

 Ω

2



((sn+1 )+ (un+1 )− )j+ 12 − ((sn+1 )+ (un+1 )+ )j− 12 h h h h

j=1

−

1 [(un+1 2 + 2un+1 − unh 2 ) h h 4τ

N     n+1 − n+1 n+1 + s ) − (u s ) )− − (qhn+1 pn+1 )+ (un+1 + (qhn+1 pn+1 h h h h h h j+ 1 j− 1 j+ 1 j− 1

j=1

−

un+1−i un+1 dx + qhn+1 2 + h h

N  

j=1

+

Ω

(unh 2

j=1

+



(3.12)

 

g n+1 un+1 dx. h

Simplifying the equation above, we get n+1  i=0

b(i) τα

 Ω

un+1−i un+1 dx + qhn+1 2 + h h

1 [(un+1 2 + 2un+1 − unh 2 ) h h 4τ (3.13)

2 ) + un+1 − 2unh + un−1 2 ] − (unh 2 + 2unh − un−1 h h h   ))un+1 dx + g n+1 un+1 dx. = (−2f (unh ) + f (un−1 h h h Ω

Ω

Multiplying (3.13) by 4τ , we arrive at n+1 

4τ

1−α

 b(i) Ω

i=0

un+1−i un+1 dx + 4τ qhn+1 2 + [(un+1 2 h h h

+ 2un+1 − unh 2 ) − (unh 2 + 2unh − un−1 2 ) + un+1 − 2unh + un−1 2 ] h h h h   ))un+1 dx + 4τ g n+1 un+1 dx. =4τ (−2f (unh ) + f (un−1 h h h Ω

(3.14)

Ω

Using Cauchy-Schwarz inequality and Young inequality, we have n+1 

4τ 1−α b(i)

 Ω

i=0

≤(unh 2

+

2unh

−

un+1−i un+1 dx + (un+1 2 + 2un+1 − unh 2 ) h h h h un−1 2 ) h

+

Cτ (unh 2

+

un−1 2 h

+

un+1 2 h

+ g

(3.15) n+1 2

 ).

Sum with respect to n from 1 to L − 1 to get L−1 2 2 L (uL  ) + 4τ 1−α h  + 2uh − uh

≤(u1h 2 + 2u1h − u0h 2 ) + Cτ

L−1 

L−1  n+1 

Ω

n=1 i=0

(unh 2 +

n=1

 b(i)

un+1−i un+1 dx h h

un−1 2 h

8

(3.16) +

un+1 2 h

+ g n+1 2 ).

According to Lemma 2.4, we obtain L−1  n+1 

 b(i)

n=1 i=0

Ω

un+1−i un+1 dx ≥ 0. h h

(3.17)

Therefore, we have L−1 2 2 L (uL  ) h  + 2uh − uh

≤(u1h 2 + 2u1h − u0h 2 ) + Cτ

L−1  n=1

(3.18)

(unh 2 + un−1 2 + un+1 2 + g n+1 2 ). h h

The equation above can be written as L−1 2 2 L (1 − Cτ )(uL  ) ≤(u1h 2 + 2u1h − u0h 2 ) + Cτ u0h 2 h  + 2uh − uh

+ Cτ

L−1 

L−1 

n=1

n=1

(unh 2 + 2unh − un−1 2 ) + Cτ h

g n+1 2 .

(3.19)

Taking sufficiently small τ so that 1 − Cτ > 0 and using Gronwall lemma, we get L−1 2 2 L (uL  ) h  + 2uh − uh

≤(u1h 2 + 2u1h − u0h 2 ) + Cτ u0h 2 + Cτ

L−1 

g n+1 2 .

(3.20)

n=1

For n = 0, setting v = u1h , ξ = −s1h , w = p1h , ρ = qh1 in the scheme (3.6), using a similar process to the case n > 0, we obtain u1h 2 ≤ u0h 2 + Cτ 2 g 1 2 .

(3.21)

Based on (3.20) and (3.21), noting u0 (x) = 0, we can easily get the results of Theorem 3.1. Remark 3.2 We here consider the initial value u(x, 0) = 0. Otherwise, we only need to make a transform w(x, t) = u(x, t) − u0 (x), then use some triangle inequalities to get  L−1 12  L (g n+1 2 + uh0 2 ) . uL h  = wh + uh0  ≤ C τ

(3.22)

n=1

For (3.22), we use triangle inequality to get 1  L−1  L n+1 2 2  ≤ w  + u  ≤ Cu  + C τ g  . uL h0 h0 h h

(3.23)

n=1

4

The estimation of error

For the need of error estimates, we will use the standard L2 −projection of a function w(x) with k + 1 continuous derivatives into space Vhk , denoted by P, i.e., for each j  Ij

(Pw(x) − w(x))v(x) = 0, ∀v ∈ P k (Ij )

(4.1)

9

+

and special projection P − into Vhk , i.e., for each j  

Ij

Ij

(P + w(x) − w(x))v(x) = 0, ∀v ∈ P k−1 (Ij ), P + w(x+ ) = w(xj− 12 ), j− 1 2

(4.2) (P − w(x) − w(x))v(x) = 0, ∀v ∈ P k−1 (Ij ), P − w(x+ ) = w(xj+ 12 ). j+ 1 2

Based on the projection P, ones can get the following approximation results 1

we  + hwe ∞ + h 2 we τh ≤ C0 hk+1 ,

(4.3)

+

where we = Pw − w or we = P − w − w, τh denotes the set of boundary points of all elements Ij . we τh denotes the L2 −norm of we on τh , which is we τh =

N  j=0

((we )+ )2 + (we )− )2 j+ 1 j+ 1 2

 12

2

.

(4.4)

For simplifying the process of estimates, we denote = u(x, tn+1 ) − un+1 = P − en+1 − (P − u(x, tn+1 ) − u(x, tn+1 )), en+1 u u h = s(x, tn+1 ) − sn+1 = P + en+1 − (P + s(x, tn+1 ) − s(x, tn+1 )), en+1 s s h en+1 = q(x, tn+1 ) − qhn+1 = P − en+1 − (P − q(x, tn+1 ) − q(x, tn+1 )), q q

(4.5)

= p(x, tn+1 ) − pn+1 = P + en+1 − (P + p(x, tn+1 ) − p(x, tn+1 )). en+1 p p h Subtracting (3.7) from (3.3), and using the fluxes (3.11), we can obtain the error equations Case: n = 0    1 1 1  eu − e0u 1−i vdx − b(i) e vdx + e1s vx dx u τ α i=0 τ Ω Ω Ω  N  + ((e1s )+ v − )j+ 12 − ((e1s )+ v + )j− 12 + (f (u0 ) − f (u0h ))vdx Ω

j=1

 +

Ω

 + Ω

 + 

Ω

+ Ω

e1s wdx + e1q ρdx + e1p ξdx +

 Ω

 Ω

 Ω

e1q wx dx −

e1p ρx dx −

N  j=1

N  j=1

e1u ξx dx −

N  j=1

(((e1q )− w− )j+ 12 − ((e1q )− w+ )j− 12 ) (4.6)

(((e1p )+ ρ− )j+ 12 − ((e1p )+ ρ+ )j− 12 )

(((e1u )− ξ − )j+ 12 − ((e1u )− ξ + )j− 12 )

(E0 + E1 + E2 )vdx = 0.

10

Case: n ≥ 1    n+1 1  3en+1 − 4enu + en−1 u u n+1−i vdx − b(i) e vdx + en+1 vx dx u s τ α i=0 2τ Ω Ω Ω  N  + − n+1 + + 1 − ((e 1 + + ((en+1 ) v ) ) v ) (2f (un ) − 2f (unh ) − f (un−1 ) + f (un−1 ))vdx j+ 2 j− 2 s s h Ω

j=1

 +

Ω

 + Ω

 + 

Ω

+ Ω

en+1 wdx + s en+1 ρdx + q en+1 ξdx + p



en+1 wx dx − q

Ω

 Ω

 Ω

N  (((en+1 )− w− )j+ 12 − ((en+1 )− w+ )j− 12 ) q q j=1

en+1 ρx dx − p

j=1

en+1 ξx dx − u

(4.7)

N 

N  j=1

(((en+1 )+ ρ− )j+ 12 − ((en+1 )+ ρ+ )j− 12 ) p p

(((en+1 )− ξ − )j+ 12 − ((en+1 )− ξ + )j− 12 ) u u

(E0 + E3 + E4 )vdx = 0.

in L2 -norm. In the next theorem, we will mainly discuss the estimate of P − en+1 u Theorem 4.1 For the error system (4.6)-(4.7) with u(t) ∈ C 2 [0, T ], the following error inequality for P − en+1 holds, for sufficiently small τ u 1

 ≤ C(τ 2 + hk+1 + τ − 2 −α hk+1 ), P − en+1 u

(4.8)

where C is a positive constant which does not depend on time step length Δt and space mesh h. Proof. For n ≥ 1, Eq.(4.7) can be written as   n+1 1  3P − en+1 − 4P − enu + P − en−1 u u − n+1−i vdx b(i) P e vdx + u α τ i=0 2τ Ω Ω  N  − P + en+1 vx dx + ((P + en+1 )+ v − )j+ 12 − ((P + en+1 )+ v + )j− 12 s s s Ω

 + Ω

 + Ω

 + Ω

j=1

P + en+1 wdx + s P − en+1 ρdx + q P + en+1 ξdx + p



Ω

 Ω

 Ω

P − en+1 wx dx − q

N  (((P − en+1 )− w− )j+ 12 − ((P − en+1 )− w+ )j− 12 ) q q j=1

P + en+1 ρx dx − p

N  j=1

P − en+1 ξx dx − u

N  j=1

(((P + en+1 )+ ρ− )j+ 12 − ((P + en+1 )+ ρ+ )j− 12 ) p p

(((P − en+1 )− ξ − )j+ 12 − ((P − en+1 )− ξ + )j− 12 ) u u

11

 = Ω

+

+

(−2f (un ) + 2f (unh ) + f (un−1 ) − f (un−1 ))vdx + h

N 



Ω N  j=1



Ω N  j=1

(P p(x, tn+1 ) − p(x, tn+1 ))ξdx +

j=1

 Ω



+

Ω N 

Ω

(P − q(x, tn+1 ) − q(x, tn+1 ))wx dx

(4.9)

(P + p(x, tn+1 ) − p(x, tn+1 ))ρx dx

((((P + p(x, tn+1 ) − p(x, tn+1 )))+ ρ− )j+ 12 − (((P + p(x, tn+1 ) − p(x, tn+1 )))+ ρ+ )j− 12 )

 +



((((P − q(x, tn+1 ) − q(x, tn+1 )))− w− )j+ 12 − (((P − q(x, tn+1 ) − q(x, tn+1 )))− w+ )j− 12 )

(P − q(x, tn+1 ) − q(x, tn+1 ))ρdx +

+

−

Dt P − u(x, tn+1 )vdx

(((P + s(x, tn+1 ) − s(x, tn+1 )))+ v − )j+ 12 − (((P + s(x, tn+1 ) − s(x, tn+1 )))+ v + )j− 12

(P + s(x, tn+1 ) − s(x, tn+1 ))wdx +

+

−

Ω

  n+1 1  − b(i) (P u(x, t ) − u(x, t ))vdx − (P + s(x, tn+1 ) − s(x, tn+1 ))vx dx n+1−i n+1−i τ α i=0 Ω Ω

j=1

−



Ω

(P − u(x, tn+1 ) − u(x, tn+1 ))ξx dx

((((P − u(x, tn+1 ) − u(x, tn+1 )))− ξ − )j+ 12 − (((P − u(x, tn+1 ) − u(x, tn+1 )))− ξ + )j− 12 )

 −

Ω

(E0 + E3 + E4 )vdx,

where Dt P − u(x, tn+1 ) =

3(P − u(x, tn+1 ) − u(x, tn+1 )) − 4(P − u(x, tn ) − u(x, tn )) + (P − u(x, tn−1 ) − u(x, tn−1 )) . 2τ

, ξ = −P − en+1 , w = P − en+1 , ρ = P − en+1 , then we get Take the test functions v = P − en+1 u s p q   n+1 3P − en+1 − 4P − enu + P − en−1 1  u u − n+1−i − n+1 P − en+1 b(i) P e P e dx + dx + P − en+1 2 u u u q α τ i=0 2τ Ω Ω   − n+1 = (−2f (un ) + 2f (unh ) + f (un−1 ) − f (un−1 ))P e dx + Dt P − u(x, tn+1 )P − en+1 dx u u h Ω

n+1 

Ω



1 b(i) (P − u(x, tn+1−i ) − u(x, tn+1−i ))P − en+1 dx u τ α i=0 Ω   dx − (E0 + E3 + E4 )P − en+1 dx. + (P − q(x, tn+1 ) − q(x, tn+1 ))P − en+1 q u

+

Ω

Ω

(4.10)

12

The equation above can be written as  n+1 1  1 [(P − en+1 b(i) P − en+1−i P − en+1 dx + 2 + 2P − en+1 − P − enu 2 ) u u u u τ α i=0 4τ Ω − (P − enu 2 + 2P − enu − P − en−1 2 ) + P − en+1 − 2P − enu + P − en−1 2 ] + P − en+1 2 u u u q   = (−2f (un ) + 2f (unh ) + f (un−1 ) − f (un−1 ))P − en+1 dx + Dt P − u(x, tn+1 )P − en+1 dx u u h Ω

Ω



n+1 

(4.11)

1 b(i) (P − u(x, tn+1−i ) − u(x, tn+1−i ))P − en+1 dx u τ α i=0 Ω   dx − (E0 + E3 + E4 )P − en+1 dx. + (P − q(x, tn+1 ) − q(x, tn+1 ))P − en+1 q u

+

Ω

Ω

Multiplying (4.11) by 4τ , we get 4τ 1−α

n+1 

 b(i) Ω

i=0

P − en+1−i P − en+1 dx + [(P − en+1 2 + 2P − en+1 − P − enu 2 ) u u u u

2 ) + P − en+1 − 2P − enu + P − en−1 2 ] + 4τ P − en+1 2 − (P − enu 2 + 2P − enu − P − en−1 u u u q   =4τ (−2f (un ) + 2f (unh ) + f (un−1 ) − f (un−1 ))P − en+1 dx + 4τ Dt P − u(x, tn+1 )P − en+1 dx (4.12) u u h Ω Ω  n+1  + 4τ 1−α b(i) (P − u(x, tn+1−i ) − u(x, tn+1−i ))P − en+1 dx u Ω

i=0

 + 4τ Ω

(P − q(x, tn+1 ) − q(x, tn+1 ))P − en+1 dx − 4τ q

 Ω

(E0 + E3 + E4 )P − en+1 dx. u

According to the Cauchy-Schwarz inequality and Young inequality, we have 4τ

1−α

n+1 

 b(i) Ω

i=0

P − en+1−i P − en+1 dx + [(P − en+1 2 + 2P − en+1 − P − enu 2 ) u u u u

− (P − enu 2 + 2P − enu − P − en−1 2 ) + P − en+1 − 2P − enu + P − en−1 2 ] + 4τ P − en+1 2 u u u q ≤Cτ [ − 2f (un ) + 2f (unh ) + f (un−1 ) − f (un−1 )2 + P − en+1 2 + Dt P − u(x, tn+1 )2 u h 2

2

2

−

2

+ E0  + E3  + E4  ] + 2τ P q(x, tn+1 ) − q(x, tn+1 ) + + C[τ 1−2α

n+1 

(4.13)

2τ P − en+1 2 q

|b(i)|2 P − u(x, tn+1−i ) − u(x, tn+1−i )2 + τ P − en+1 2 ]. u

i=0

Eq.(4.13) can be simplified as 4τ 1−α

n+1 

 b(i) Ω

i=0

P − en+1−i P − en+1 dx + [(P − en+1 2 + 2P − en+1 − P − enu 2 ) u u u u

2 )] − (P − enu 2 + 2P − enu − P − en−1 u )2 + P − en+1 2 + Dt P − u(x, tn+1 )2 ≤Cτ [ − 2f (un ) + 2f (unh ) + f (un−1 ) − f (un−1 u h 2

2

2

−

2

+ E0  + E3  + E4  + P q(x, tn+1 ) − q(x, tn+1 ) ] + Cτ 1−2α

n+1 

|b(i)|2 P − u(x, tn+1−i ) − u(x, tn+1−i )2 .

i=0

13

(4.14)

Use Cauchy-Schwarz inequality with triangle inequality to get  − 2f (un ) + 2f (unh ) + f (un−1 ) − f (un−1 )2 h 



= − 2f (uθn )(un − unh ) + f (uθn−1 )(un−1 − un−1 ) h ≤(P − enu 2

+

P − en−1 2 u

−

(4.15) 2

−

2

+ P u(x, tn ) − u(x, tn ) + P u(x, tn−1 ) − u(x, tn−1 ) ).

Utilizing the inequality above and summing with respect to n from 1 to L − 1, we have 4τ 1−α

L−1  n+1 

 b(i)

n=1 i=0

Ω

2 − L − L−1 2 P − en+1−i P − en+1 dx + (P − eL  ) u u u  + 2P eu − P eu

≤(P − e1u 2 + 2P − e1u − P − e0u 2 ) + Cτ

L−1 

[P − enu 2 + P − en−1 2 u

n=1

+ P − u(x, tn ) − u(x, tn )2 + P − u(x, tn−1 ) − u(x, tn−1 )2 + P − en+1 2 + P − q(x, tn+1 ) − q(x, tn+1 )2 + E0 2 u  tL + E3 2 + E4 2 ] + C (P − u(x, t) − u(x, t))t 2 ds

(4.16)

t0

+ Cτ 1−2α

L−1  n+1 

|b(i)|2 P − u(x, tn+1−i ) − u(x, tn+1−i )2 .

n=1 i=0

According to Lemma 2.4, we obtain L−1  n+1 

 b(i)

n=1 i=0

Ω

P − en+1−i P − en+1 dx ≥ 0. u u

(4.17)

So, we get 2 − L − L−1 2  ) (P − eL u  + 2P eu − P eu

≤(P − e1u 2 + 2P − e1u − P − e0u 2 ) + Cτ

L−1 

[P − enu 2 + P − en−1 2 u

n=1

+ P − u(x, tn ) − u(x, tn )2 + P − u(x, tn−1 ) − u(x, tn−1 )2 + P − en+1 2 + P − q(x, tn+1 ) − q(x, tn+1 )2 + E0 2 u  tL + E3 2 + E4 2 ] + C (P − u(x, t) − u(x, t))t 2 ds

(4.18)

t0

+ Cτ 1−2α

L−1  n+1 

|b(i)|2 P − u(x, tn+1−i ) − u(x, tn+1−i )2 .

n=1 i=0

We rewrite (4.18) as 2 − L − L−1 2  ) (1 − Cτ )(P − eL u  + 2P eu − P eu

≤(P − e1u 2 + 2P − e1u − P − e0u 2 ) + Cτ

L−1 

[P − enu 2

n=1

+ 2P − enu − P − en−1 2 + P − u(x, tn ) − u(x, tn )2 u + P − u(x, tn−1 ) − u(x, tn−1 )2 + P − q(x, tn+1 ) − q(x, tn+1 )2  tL 2 2 2 + E0  + E3  + E4  ] + C (P − u(x, t) − u(x, t))t 2 ds t0

+ Cτ P − e0u 2 + Cτ 1−2α

L−1  n+1 

|b(i)|2 P − u(x, tn+1−i ) − u(x, tn+1−i )2 .

n=1 i=0

14

(4.19)

Using the Gronwall inequality, we arrive at 2 − L − L−1 2 (P − eL  ) u  + 2P eu − P eu

≤(P − e1u 2 + 2P − e1u − P − e0u 2 ) + Cτ

L−1 

[P − u(x, tn ) − u(x, tn )2

n=1

+ P − u(x, tn−1 ) − u(x, tn−1 )2 + P − q(x, tn+1 ) − q(x, tn+1 )2  tL + E0 2 + E3 2 + E4 2 ] + C (P − u(x, t) − u(x, t))t 2 ds

(4.20)

t0

+ Cτ P − e0u 2 + Cτ 1−2α

L−1  n+1 

|b(i)|2 P − u(x, tn+1−i ) − u(x, tn+1−i )2 .

n=1 i=0

For n = 0, we have 1 2τ 1−α b(0)P − e1u 2 + P − e1u 2 2 ≤Cτ 2 [P − u(x, t0 ) − u(x, t0 )2 + E0 2 + E1 2 + E2 2 ]  t1 (P − u(x, t) − u(x, t))t 2 ds + 2τ P − q(x, t1 ) − q(x, t1 )2 + Cτ t0

+ C[τ 2−2α

1 

(4.21) |b(i)|2 P − u(x, t1−i ) − u(x, t1−i )2 ]

i=0

≤C(τ 2 h2k+2 + τ 6 + τ 4 ) + C(τ 2 h2k+2 + τ h2k+2 ) + Cτ 2−2α h2k+2 ≤C(τ 4 + h2k+2 ). For n ≥ 1, we can easily get 2 − L − L−1 2 (P − eL  ) + 2τ u  + 2P eu − P eu

l−1 

P − en+1 2 q

(4.22)

n=1 4

2k+2

≤C(τ + h

+τ

−1−2α 2k+2

h

).

So, we finish the proof of Theorem 4.1. Remark 4.2 From the conclusion of theorem, by a simple calculation with triangle inequality and (4.3), we easily get the L2 -error un − unh . Based on the results of theorem, ones can easily see that the convergence order in the temporal direction is 2, which is higher than the results with O(τ 2−α ) obtained by L1-approximation in Refs. [56, 49]. In what follows, we will implement the numerical procedure by using the current algorithm in this article to verify theoretical results in Theorem 4.1.

5

Numerical tests

In this section, we will show the detailed process of numerical calculations based on LDG method with second-order WSGD approximation, then give some numerical results to illustrate the feasibility and effectiveness of the proposed method.

5.1

Numerical algorithm

Choosing Legendre polynomials Pl as local basis functions, we can express our approximate solution uh as follows: uh (x, t) =

N  k 

uj,l ϕj,l (x),

(5.1)

j=1 l=0

15

where ⎧ ⎨P (2(x − x )/Δx ), l j j ϕj,l (x) = ⎩0,

if x ∈ Ij ,

(5.2)

else.

To carry out the numerical test, we approximate auxiliary variables sh , qh and ph in the same way and choose the parameter k = 1 in (5.1). In (3.6) and (3.7), we choose v = w = ρ = ξ = ϕj,l and substitute (5.1) into the two equations to get the numerical procedure. The computational process mainly contains two steps: Step I: The case n = 0 ⎡ ⎢ ⎢ ⎢ ⎣

( τ1 +

b(0) τ α )A

B A

⎤ ⎡ 1⎤ 1  U R ⎥ ⎥ ⎢S ⎥ ⎢ ⎥ ⎢  1 ⎥ ⎢ 0 ⎥ ⎥ ⎢  1⎥ = ⎢  ⎥ −B⎦ ⎣Q ⎦ ⎣ 0 ⎦ 0 P 1 A ⎤⎡

C A

C

(5.3)

 1 can be calculated by the initial vectors U  0 and the vectors G  1 on the source In (5.3), the right term R term g. Since the coefficient matrix is invertible, we can easily solve the equation and find the solution  1, S  1, Q  1 , P 1 ]T . vector [U Step II: The case n ≥ 1 ⎡ ⎢ ⎢ ⎢ ⎣

3 + ( 2τ

b(0) τ α )A

C

B A

⎤⎡

C A

⎤ ⎡ n+1 ⎤  n+1  U R ⎥ ⎢S ⎢ 0 ⎥ n+1 ⎥  ⎥⎢ ⎥ ⎥ ⎢ ⎥ ⎢  n+1 ⎥ = ⎢  ⎥ ⎦ ⎣ 0 ⎦ −B⎦ ⎣Q 0 P n+1 A

(5.4)

 1, Q  1 , P 1 ]T obtained in step I, the value R  2 can  1, S Based on the initial vectors and the solution vector [U 2 2 2 2 T  ,S  ,Q  , P ] can be uniquely solved in (5.4). Similarly, we can have be gotten and further the vector [U n+1  n+1  n+1  n+1 T   n+1 arrived at by the the solution vectors [U ,S ,Q ,P ] on the basis of the vector value R j j  j  j T  solution vectors [U , S , Q , P ] , j = 0, 1, · · · n. In (5.3) and (5.4), the used vectors are as follows: 1 = R  n+1 R

b(1)   0  0  1 AU − F + G , τ τα n+1 2 n 1  n−1  b(i)  n+1−i  n+1 , AU = AU − − AU − 2F n + F n−1 + G α τ 2τ τ i=1

1

+

 n = (sn , sn , sn , sn , · · · , sn , sn )T ,  n = (un , un , un , un , · · · , un , un )T , S U 1,0 1,1 2,0 2,1 N,0 N,1 1,0 1,1 2,0 2,1 N,0 N,1 n n n n n n n T n n n n n n  Q = (q1,0 , q1,1 , q2,0 , q2,1 , · · · , qN,0 , qN,1 ) , P = (p1,0 , p1,1 , p2,0 , p2,1 , · · · , pN,0 , pnN,1 )T ,  n = ((g n , ϕ1,0 ), (g n , ϕ1,1 ), (g n , ϕ2,0 ), (g n , ϕ2,1 ), · · · , (g n , ϕN,0 ), (g n , ϕN,1 ))T , G F n = ((f (unh ), ϕ1,0 ), (f (unh ), ϕ1,1 ), (f (unh ), ϕ2,0 ), (f (unh ), ϕ2,1 ), · · · , (f (unh ), ϕN,0 ), (f (unh ), ϕN,1 ))T .

16

(5.5)

where ⎛

⎞

Δx1

⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ A=⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛

−1 ⎜ ⎜−1 ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ B=⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜1 ⎝ 1

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

Δx1 /3 ..

. ..

. ..

. ΔxN

1 −1

1 1 ..

.

..

.

−1 −1 .. . .. .

ΔxN /3 ⎞

..

.

..

.

..

.

..

.

−1 −1

1 −1

−1 −1

2N ×2N

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ 1 −1⎟ ⎟ 1 −1⎟ ⎟ −1 1 ⎟ ⎠ −1 −1

⎛ −1 −1 ⎜ ⎜ 1 −1 ⎜ ⎜1 1 −1 ⎜ ⎜ −1 −1 1 ⎜ ⎜ C=⎜ .. ⎜ . ⎜ ⎜ . .. ⎜ ⎜ ⎜ ⎝

2N ×2N

⎞ 1 ⎟ −1⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ .. ⎟ . ⎟ ⎟ .. ⎟ . ⎟ 1 −1 −1⎟ ⎠ −1 1 −1 1 −1

−1 −1 .. . .. .

..

.

..

.

1 −1

2N ×2N

(5.6)

5.2

Numerical example

To implement the numerical calculations, we consider the following nonlinear fourth-order PDE with fractional derivative ∂u(x, t) ∂ α u(x, t) ∂ 4 u(x, t) + + + u3 (x, t) − u(x, t) = g(x, t), ∂t ∂tα ∂x4 u(x, 0) = 0, x ∈ [0, 2π], 2−α

(x, t) ∈ [0, 2π] × [0, 1],

(5.7)

2t ) sin(x) + t6 sin3 (x). Based on the known source where the known source term is g(x, t) = (2t + Γ(3−α) term, we can easily find that the exact solution is u(x, t) = t2 sin(x). In Table 1, we firstly take changed spatial meshes h = 2π/5, 2π/10, 2π/20, 2π/40 and fixed time step τ = 1/100, then check that the spatial convergence rate is close to 2, which is not changed as the values 0.01, 0.5, 0.99 of α and is consistent with the theoretical result covering the index k = 1. From the calculated results in Table 1, ones can easily know that the errors in L2 and L∞ -norms attain optimal order of accuracy for α = 0.01, 0.5, 0.99. In Table 2, for checking the temporal convergence rate, we fix the spatial step parameter h = 2π/200 and take changed time meshes τ = 1/20, 1/40, 1/60, 1/80, then we can see that the convergence order in L2 and L∞ -norms with α = 0.01 in temporal direction is close to 2. The similar convergence results with α = 0.5, 0.99 are also arrived at in Table 2. The approximate results in Table 2 illustrate that we can get second-order accuracy in the temporal direction by using our method, which also tells ones the current approximation can keep the second-order convergence rate for changed α. For seeing well the convergence rate, we show two figures on convergence rate based on the calculated data. With parameter α = 0.7 and time mesh τ = 1/160, we describe the spatial convergence order in Fig.1. Again, Fig.2 displays the time convergence order with α = 0.4 and h = 2π/200. Based on the discussions for the numerical results in Tables 1-2 and Figs. 1-2, we can see that the current method studied in this paper can well numerically solve the nonlinear fractional fourth-order

17

Table 1: The errors and convergence orders for u with time mesh τ = 1/100 α 0.01

0.5

0.99

h

u − uh L2

Rate

u − uh L∞

Rate

2π/5 2π/10 2π/20 2π/40

1.3753e-001 3.0072e-002 6.9448e-003 1.4106e-003

2.1932 2.1144 2.2996

7.5835e-002 1.6658e-002 3.9294e-003 8.1003e-004

2.1867 2.0838 2.2783

2π/5 2π/10 2π/20 2π/40

1.3750e-001 3.0112e-002 6.9822e-003 1.4461e-003

2.1910 2.1086 2.2715

7.5920e-002 1.6680e-002 3.9505e-003 8.3020e-004

2.1863 2.0780 2.2505

2π/5 2π/10 2π/20 2π/40

1.3728e-001 3.0151e-002 7.0297e-003 1.4935e-003

2.1868 2.1007 2.2347

7.6016e-002 1.6720e-002 3.9958e-003 8.7431e-004

2.1847 2.0651 2.1923

Table 2: The errors and convergence orders for u with space mesh h = 2π/200 α 0.01

0.5

0.99

Δt

u − uh L2

Rate

u − uh L∞

Rate

1/20 1/40 1/60 1/80

8.4185e-003 2.3397e-003 1.0485e-003 5.7271e-004

1.8472 1.9797 2.1020

5.3592e-003 1.4771e-003 6.6194e-004 3.6276e-004

1.8592 1.9796 2.0906

1/20 1/40 1/60 1/80

7.9140e-003 2.1627e-003 9.6108e-004 5.2094e-004

1.8716 2.0003 2.1288

5.0969e-003 1.3805e-003 6.1345e-004 3.3374e-004

1.8844 2.0004 2.1160

1/20 1/40 1/60 1/80

7.5132e-003 2.0096e-003 8.8596e-004 4.7762e-004

1.9025 2.0199 2.1477

5.1028e-003 1.3597e-003 6.0079e-004 3.2568e-004

1.9080 2.0144 2.1285

PDE. Furthermore, the results also show that the WSGD operator can be well combined with LDG method.

6

Conclusion

In this paper, based on WSGD operator [1, 2] for fractional derivative, we present a LDG method to solve the nonlinear time-fractional fourth-order partial differential problem. Compared with the L1approximation with order (2 − α) used in many references, the current approximate results, which are not affected by the changed parameters α, can arrive at order 2 by considering the idea of WSGD operator. The stability is proved in detail and the detailed process of error estimates is derived. To confirm the accuracy and capability of our method, some numerical results are calculated and analyzed. Future work will include the analysis of LDG method for two-dimensional fractional problems.

18

Figure 1: Convergence rate for u with time mesh τ = 1/160

α=0.7 −3 L2 L −4

∞ 2

h

log(error)

−5

−6

−7

−8

−9 −4.5

−4

−3.5

−3

−2.5

−2

log(h)

Figure 2: Convergence rate for u with space mesh h = 1/200

α=0.4 −3 L2 L∞ −4

2

τ

log(error)

−5

−6

−7

−8

−9 −4.5

−4

−3.5

−3 log(τ)

19

−2.5

−2

Acknowledgements Authors thank the reviewers and editor very much for their valuable comments and constructive suggestions to improve our work. This work is supported by the National Natural Science Fund (11301258, 11361035, 11501311), the Natural Science Fund of Inner Mongolia Autonomous Region (2014BS0101), the Scientific Research Projection of Higher Schools of Inner Mongolia(NJZZ12011, NJZY14013).

References [1] W.Y. Tian, H. Zhou, W. Deng, A class of second order difference approximations for solving space fractional diffusion equations. Math. Comput. 84 (2015) 1703-1727. [2] Z.B. Wang, S.W. Vong, Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation, J. Comput. Phys. 277 (2014) 1-15. [3] A. Atangana, On the stability and convergence of the time-fractional variable order telegraph equation. J. Comput. Phys. 293 (2015) 104-114. [4] A. Atangana, D. Baleanu, Numerical solution of a kind of fractional parabolic equations via two difference schemes, Abst. Appl. Anal. 2013 (2013), Article ID 828764, 8 pages. [5] P.D. Wang, C.M. Huang, An energy conservative difference scheme for the nonlinear fractional Schr¨ odinger equations, J. Comput. Phys. 293 (2015) 238-251. [6] P.D. Wang, C.M. Huang, A conservative linearized difference scheme for the nonlinear fractional Schr¨ odinger equation, Numer. Algor. 69(3) (2015) 625-641. [7] A.H. Bhrawy, M.A. Zaky, Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation, Nonlinear Dynamics, 80(1-2) (2015) 101-116. [8] A.H. Bhrawy, D. Baleanu, F. Mallawi, A new numerical technique for solving fractional sub-diffusion and reaction sub-diffusion equations with a non-linear source term. Thermal Science, 19(suppl.1) (2015) 25-34. [9] S.B. Yuste, J. Quintana-Murillo, A finite difference method with non-uniform timesteps for fractional diffusion equations, Comput. Phys. Commun. 183(12) (2012) 2594-2600. [10] R.H. Nochetto, E. Ot´ arola, A.J. Salgado, A PDE approach to space-time fractional parabolic problems, arXiv:1404.0068 [math.NA]. [11] K. Mustapha, M. McLean, Piecewise-linear, discontinous Galerkin method for a fractional diffusion equation, Numer. Algorithm 56 (2011) 159-184. [12] B. Jin, R. Lazarov, Y.K. Liu, Z. Zhou, The Galerkin finite element method for a multi-term time-fractional diffusion equation, J. Comput. Phys. 281 (2015) 825-843. [13] B. Jin, R. Lazarov, Z. Zhou, Error estimates for a semidiscrete finite element method for fractional order parabolic equations, SIAM J. Numer. Anal. 51 (2013) 445-466. [14] A.J. Cheng, H. Wang, K.X. Wang, A Eulerian-Lagrangian control volume method for solute transport with anomalous diffusion, Numer. Methods Partial Differ. Equ. 31(1) (2015) 253-267. [15] W.P. Bu, Y.F. Tang, J.Y. Yang, Galerkin finite element method for two-dimensional Riesz space fractional diffusion equations. J. Comput. Phys. 276 (2014) 26-38. [16] W. Bu, Y. Tang, Y. Wu, et al., Finite difference/finite element method for two-dimensional space and time fractional Bloch-Torrey equations, J. Comput. Phys. 293 (2015) 264-279. [17] F.H. Zeng, F.W. Liu, C.P. Li, K. Burrage, I. Turner, V. Anh, A Crank-Nicolson ADI spectral method for a two-dimensional Riesz space fractional nonlinear reaction-diffusion equation, SIAM J. Numer. Anal. 52(6) (2014) 2599-2622. [18] H. Hejazi, T. Moroney, F. Liu, Stability and convergence of a finite volume method for the space fractional advection-dispersion equation, J. Comput. Appl. Math. 255 (2014) 684-697. [19] P. Zhuang, F. Liu, I. Turner, Y.T. Gu, Finite volume and finite element methods for solving a one-dimensional space-fractional Boussinesq equation, Appl. Math. Model. 38 (2014) 3860-3870. [20] C. Tadjeran, M. Meerschaert, A second-order accurate numerical method for the two-dimensional fractional diffusion equation, J. Comput. Phys. 220 (2007) 813-823. [21] L.B. Feng, P. Zhuang, F. Liu, I. Turner, Stability and convergence of a new finite volume method for a two-sided space-fractional diffusion equation, Appl. Math. Comput. 257 (2015) 52-65.

20

[22] Q. Yang, I. Turner, F. Liu, et al., Novel numerical methods for solving the time-space fractional diffusion equation in two dimensions, SIAM J. Sci. Comput. 33(3) (2011) 1159-1180. [23] M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection-dispersion flow equations, J. Comput. Appl. Math. 172 (2004) 65-77. [24] Y.M. Lin, C.J. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys. 225(2) (2007) 1533-1552. [25] Y.M. Lin, X.J. Li, C.J. Xu, Finite difference/spectral approximations for the fractional Cable equation, Math. Comput. 80 (2011) 1369-1396. [26] F. Liu, P. Zhuang, I. Turner, K. Burrage, V. Anh, A new fractional finite volume method for solving the fractional diffusion equation, Appl. Math. Model. 38 (2014) 3871-3878. [27] C.M. Chen, F. Liu, V. Anh, I. Turner, Numerical methods for solving a two-dimensional variable-order anomalous subdiffusion equation, Math. Comput. 81 (2012) 345-366. [28] Y.J. Jiang, J.T. Ma, High-order finite element methods for time-fractional partial differential equations, J. Comput. Appl. Math. 235(11) (2011) 3285-3290. [29] J.T. Ma, J.Q. Liu, Z.Q. Zhou, Convergence analysis of moving finite element methods for space fractional differential equations, J. Comput. Appl. Math. 255 (2014) 661-670. [30] J.C. Li, Y.Q. Huang, Y.P. Lin, Developing finite element methods for maxwell’s equations in a cole-cole dispersive medium, SIAM J. Sci. Comput. 33(6) (2011) 3153-3174. [31] C.P. Li, Z.G. Zhao, Y.Q. Chen, Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion, Comput. Math. Appl. 62(3) (2011) 855-875. [32] C.P. Li, H.F. Ding, Higher order finite difference method for the reaction and anomalous-diffusion equation, Appl. Math. Model. 38 (15-16) (2014) 3802-3821. [33] H.F. Ding, C.P. Li, High-order compact difference schemes for the modified anomalous subdiffusion equation. arXiv preprint arXiv 1408.5591 (2014). [34] N. Zhang, W. Deng, Y. Wu, Finite difference/element method for a two-dimensional modified fractional diffusion equation, Adv. Appl. Math. Mech. 4 (2012) 496-518. [35] H. Wang, K. Wang, T. Sircar, A direct O(N log 2 N ) finite difference method for fractional diffusion equations, J. Comput. Phys. 229 (2010) 8095-8104. [36] Y.N. Zhang, Z.Z. Sun, H.L. Liao, Finite difference methods for the time fractional diffusion equation on non-uniform meshes, J. Comput. Phys. 265 (2014) 195-210. [37] E. Sousa, C. Li, A weighted finite difference method for the fractional diffusion equation based on the Riemann-Liouville derivative, Appl. Numer. Math. 90 (2015) 22-37. [38] M.R. Cui. Compact finite difference method for the fractional diffusion equation, J. Comput. Phys. 228 (2009) 7792-7804. [39] C.C. Ji, Z.Z. Sun, A high-order compact finite difference scheme for the fractional sub-diffusion equation, J. Sci. Comput. in press. [40] G.H. Gao, H.W. Sun, Z.Z. Sun, Stability and convergence of finite difference schemes for a class of timefractional sub-diffusion equations based on certain superconvergence, J. Comput. Phys. 280 (2015) 510-528. [41] W.L. Li, D. Xu, Finite central difference/finite element approximations for parabolic integro-differential equations, Computing, 90 (2010) 89-111. [42] K.B. Oldham, J. Spainer, The Fractional Calculus, Academic Press, New York, 1974. [43] O.P. Agrawal, A general solution for a fourth-order fractional diffusion-wave equation defined in a bounded domain, Comput. Struct. 79 (2001) 1497-1501. [44] H. Jafari, M. Dehghan, K. Sayevand, Solving a fourth-order fractional diffusion wave equation in a bounded domain by decomposition method, Numer. Methods Partial Diff. Equ. 24 (2008) 1115-1126. [45] A. Golbabai, K. Sayevand, Fractional calculus-a new approach to the analysis of generalized fourth-order diffusion-wave equations, Comput. Math. Appl. 67 (2011) 2227-2231. [46] X.L. Hu, L.M. Zhang, On finite difference methods for fourth-order fractional diffusion-wave and subdiffusion systems, Appl. Math. Comput. 218 (2012) 5019-5034. [47] X.L. Hu, L.M. Zhang, A new implicit compact difference scheme for the fourth-order fractional diffusion-wave system, Int. J. Comput. Math. 91(10) (2014) 2215-2231.

21

[48] C.C. Ji, Z.Z. Sun, Z.P. Hao, Numerical algorithms with high spatial accuracy for the fourth-order fractional sub-diffusion equations with the first dirichlet boundary conditions, J. Sci. Comput. 2015, in press. [49] L. Guo, Z. Wang, S. Vong, Fully discrete local discontinuous Galerkin methods for some time-fractional fourth-order problems, Int. J. Comput. Math. 2015, in press. [50] S. Vong, Z. Wang, Compact finite difference scheme for the fourth-order fractional subdiffusion system, Adv. Appl. Math. Mech. 6(4) (2014) 419-435. [51] Y. Liu, Y.W. Du, H. Li, J.C. Li, S. He. A two-grid mixed finite element method for a nonlinear fourth-order reaction diffusion problem with time-fractional derivative, Comput. Math. Appl. 2015, in press. [52] Y. Liu, Y.W. Du, H. Li, S. He, W. Gao, Finite difference/finite element method for a nonlinear time-fractional fourth-order reaction-diffusion problem, Comput. Math. Appl. 70 (2015) 573-591. [53] Y. Liu, Z.C. Fang, H. Li, S. He, A mixed finite element method for a time-fractional fourth-order partial differential equation, Appl. Math. Comput. 243 (2014) 703-717. [54] Z.Z. Sun, X.N. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193-209. [55] W.H. Deng, J.S. Hesthaven, Local discontinuous Galerkin methods for fractional ordinary differential equations, arXiv:1403.5759 [math.NA]. [56] L. Wei, Y. He, Analysis of a fully discrete local discontinuous Galerkin method for time-fractional fourthorder problems, Appl. Math. Model. 38 (2014) 1511-1522 [57] L. Wei, Y. He, X. Zhang, et al. Analysis of an implicit fully discrete local discontinuous Galerkin method for the time-fractional Schr¨ odinger equation, Finite Elem. Anal. Des. 59 (2012) 28-34. [58] B. Cockburn, C.W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion systems, SIAM J. Numer. Anal. 35(6) (1998) 2440-2463. [59] X. Meng, C.W. Shu, B. Wu, Superconvergence of the local discontinuous Galerkin method for linear fourth order time dependent problems in one space dimension, IMA J. Numer. Anal. 32(4) (2012) 1294-1328. [60] M. Baccouch, The local discontinuous Galerkin method for the fourth-order Euler-Bernoulli partial differential equation in one space dimension. Part I: Superconvergence error analysis, J. Sci. Comput. 59(3) (2014) 795-840. [61] Y. Xu, C.W. Shu, Optimal error estimates of the semidiscrete local discontinuous Galerkin methods for high order wave equations, SIAM J. Numer. Anal. 50(1) (2012) 79-104. [62] Y. Xia, Y. Xu, C.W. Shu, Local discontinuous Galerkin methods for the Cahn-Hilliard type equations, J. Comput. Phys. 227(1) (2007) 472-491. [63] J. Yan, S. Osher, A local discontinuous Galerkin method for directly solving Hamilton-Jacobi equations, J. Comput. Phys. 230(1) (2011) 232-244. [64] F. Gao , J. Qiu, Q. Zhang, Local discontinuous Galerkin finite element method and error estimates for one class of Sobolev equation, J. Sci. Comput. 41(3) (2009) 436-460. [65] J. Zhu, J. Qiu, Local DG method using WENO type limiters for convection-diffusion problems, J. Comput. Phys. 230(11) (2011) 4353-4375. [66] Q. Zhang, F. Gao, A fully-discrete local discontinuous Galerkin method for convection-dominated Sobolev equation, J. Sci. Comput. 51(1) (2012) 107-134.

22