Convergence and superconvergence analysis of finite element methods for the time fractional diffusion equation

Journal Pre-proof Convergence and superconvergence analysis of finite element methods for the time fractional diffusion equation

Meng Li, Dongyang Shi, Lifang Pei

PII:

S0168-9274(19)30374-5

DOI:

https://doi.org/10.1016/j.apnum.2019.12.023

Reference:

APNUM 3737

To appear in:

Applied Numerical Mathematics

Received date:

29 June 2019

Revised date:

22 December 2019

Accepted date:

23 December 2019

Please cite this article as: M. Li et al., Convergence and superconvergence analysis of finite element methods for the time fractional diffusion equation, Appl. Numer. Math. (2020), doi: https://doi.org/10.1016/j.apnum.2019.12.023.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2020 Published by Elsevier.

Convergence and superconvergence analysis of finite element methods for the time fractional diffusion equation Meng Lia , Dongyang Shia,∗ , Lifang Peia a

School of Mathematics and Statistics, Zhengzhou University, Zhengzhou 450001, P. R. China

Abstract In this paper, two classes of finite element methods (FEMs) for the two-dimensional time fractional diffusion equation (TFDE) with non-smooth solution are proposed and analyzed. In the temporal direction, we adopt nonuniform L2-1σ method, and in the spatial direction, the conforming bilinear element and the nonconforming modified quasi-Wilson element are utilized. For the proposed conforming and nonconforming FEMs, by using new fractional discrete Gr¨onwall inequalities, the theoretical analysis including the L2 -norm error estimates and H 1 -norm superclose results are given in details. Furthermore, by virtue of the interpolated postprocessing techniques, the global H 1 -norm superconvergence results are presented. Finally, some numerical results illustrate the correctness of theoretical analysis. Keywords: Two dimensional TFDE, nonuniform L2-1σ formula, Bilinear FEM, Modified quasi-Wilson FEM, Superconvergence.

1. Introduction In this paper, we consider the following two-dimensional TFDE, C α 0 Dt u(x, t)

− Δu(x, t) = f (x, t), (x, t) ∈ Ω × (0, T ], u(x, t) = 0, (x, t) ∈ ∂Ω × (0, T ], u(x, 0) = u0 (x), x ∈ Ω,

(1.1) (1.2) (1.3)

where the domain Ω in R2 is a bounded polygonal region with boundary ∂Ω, u0 (x) and f (x, t) are given functions, and C0 Dαt is the αth-order left-sided Caputo derivative operator defined by  t ∂u(x, s) C α ds, t > 0, (1.4) D u(x, t) = ω1−α (t − s) 0 t ∂s 0 

This work was supported by NSF of China (No. 11701522, 11771163, 11671160, 11801527) and China Postdoctoral Science Foundation (No. 2018M632791). ∗ Corresponding author. Email addresses: [email protected] (M. Li), shi [email protected] (D. Shi∗ ), Preprint submitted to Elsevier

with ωβ (t) = tβ−1 /Γ(β). The model (1.1)-(1.3) has been proved to be one of the very valuable tools in modeling complex systems, such as glassy and disordered media [1]; see [2, 3] for further in physical and biological background and applications. The solutions for the fractional PDEs are typically less regular than ones for classical parabolic PDEs, which arise as the limiting case when α verges to 1. Sakamoto and Yamamoto [4] presented the lack of regularity of u for more general initial function u0 and a linear source term f . Stynes et al. [5] applied the theory of sectorial operators to show that u can be a smooth function of t only if the initial function u0 and the source term f satisfy some restrictive compatibility conditions. For the model (1.1) − (1.3) , we define λi and χi , i = 1, 2, . . . be the eigenvalues and eigenfunctions for the following problem Δχi = λi Xi on Ω, with boundary condition Xi = 0 on ∂Ω, where the eigenfunctions are normalised by Xi L2 . By using the theory of sectorial operators [4], for each ζ ∈ R, define the fractional power Δζ of the operator Δ with the domain ∞    2 λ2ζ |(g, X )| < ∞ , D(Δζ ) := g ∈ L2 (Ω) : i i i=1

and set the norm gΔζ :=

∞ 

2 λ2ζ i |(g, Xi )|

1/2

i=1

for j ∈ N and t ∈ (0, T ] if u0 ∈ D(Δ

j+2

) and ∂lt f (·, t) ∈ D(Δ j ), such that

u0 Δ j+2 + δkt f (·, t)Δ j+1 ≤ Cu0 , f ,

k = 1, 2,

(1.5)

where Cu0 , f > 0 is a constant independent of t. Then, from [6], the system (1.1)-(1.3) has a unique solution u, such that u(·, t) j + C0 Dαt u(·, t) j ≤ L,

∂lt u(·, t) j ≤ L(1 + tα−l ),

l = 0, 1, 2, 3,

(1.6)

for all t ∈ (0, T ], where L > 0 is a constant independent of t and υ j :=

 |β|≤ j

 12 |D υ| dx . β

Ω

2

For the homogeneous subdiffusion problem, Jin et al. [7] adopted the Laplace transform and the theory of sectorial operators drawing a conclusion that the classical uniform L1 formula has accuracy of order O(tnα−1 τ) and O(tn−1 τ) for smooth and nonsmooth initial data, respectively. To deal with the initial singularity for the Caputo-type fractional subdiffusion problems, Stynes et al. [8] studied the L1 formula based on the uniform and graded time grids, respectively. In their works, the discrete maximum principle and the direct analysis of local truncation error were utilized to build the error estimates under some regularity assumptions. For a class of timefractional reaction-diffusion problems, Huang and Stynes [9] also considered the discrete scheme 2

based on the well-known L1 discretisation in time on a graded mesh and a direct discontinuous Galerkin FEM in space on a uniform mesh. Recently, Liao etc. [10] constructed a discrete fractional Gr¨onwall-type inequality, and illustrated its use by showing the stability and convergence of the numerical schemes (permitting the use of nonuniform time steps) for fractional reactionsubdiffusion problems. Furthermore, many works about the numerical methods of the time fractional PDEs based on nonuniform meshes have been existing (see [9, 11, 12, 13, 14, 15, 16]). Very recently, Huang and Stynes [17] studied the superconvergence of the FEM of the multi-term time fractional diffusion equation. In this paper, the authors constructed the discrete scheme based on the L1 scheme on a temporal graded mesh and the bilinear FEM in spatial direction. This work is the first job to study the superconvergence of the FEM of the time fractional partial differential equation with the non-smooth solution. Motivated by above results, the main work of this paper is to study the nonuniform L2-1σ method (considered as the extension of the uniform L2-1σ method [18]), and the conforming and nonconforming FEMs for the system (1.1)-(1.3) with non-smooth solution. By means of the new discrete fractional Gr¨onwall-type inequality [10], we focus on the optimal L2 -norm optimal error estimate and H 1 -norm superconvergence analysis of the discrete schemes. The theoretical analysis shows that for the conforming bilinear FEM, the L2 -norm error estimate, and the H 1 norm superclose and superconvergence results are all sharp in the temporal direction. For the nonconforming modified quasi-Wilson FEM, we show the L2 -norm optimal error estimate is sharp in the temporal direction. Although we cannot obtain the optimal-in-time H 1 -norm error estimate, it is still meaningful to study the superclose and superconvergence properties of the nonconforming FEM. During the proof, we find that the H 1 -norm error estimate is optimal for tn not close to t = 0. A brief outline of the paper follows. In Section 2, we mainly introduce the nonuniform L2-1σ method, and the discrete fractional Gr¨onwall inequality and global consistency errors. In Section 3, the conforming bilinear FEM is considered. We obtain the optimal L2 -norm error estimate, and the H 1 -norm superclose and superconvergence results of the constructed discrete scheme, which are shown to be sharp in the temporal direction. In Section 4, we construct the nonconforming modified quasi-Wilson FEM of the system (1.1)-(1.3). The theoretical analysis of the constructed nonconforming FEM, including the the optimal L2 -norm error estimate, and the H 1 -norm superclose and superconvergence results, are obtained in details. Finally, Section 5 presents several numerical tests to verify the correctness of the theoretical analysis. In Section 6, some conclusions are drawn. In this paper, Cν , Cu , CΩ and CΛ denote positive constants which may depend on u, f , Ω and T , but independent of the mesh; they can take different values in different places. 2. Nonuniform L2-1σ method 2.1. The discrete formula of Caputo derivative Let us denote the time levels 0 = t0 < t1 < t2 < . . . < tN = T with the time-step τi = ti − ti−1 , 1 ≤ i ≤ N. Define ρi = τi /τi+1 , 1 ≤ i ≤ N − 1 be the step size ratios, and τ = max1≤i≤N τi be the maximum time-step size. We also denote the weighted time by tn−ϑ = ϑtn−1 + (1 − ϑ)tn , ϑ ∈ [0, 1). Then, let Π1,i v be the linear Lagrange interpolation function with respect to the notes ti−1 and ti , and Π2,i v be the quadratic Lagrange interpolation function with respect to the notes ti−1 , ti and 3

ti+1 . Denote Π p,i v(t) = v(t) − (Π p,i v)(t), p = 1, 2 be the corresponding interpolant error functions. Obviously, we have (Π1,i v) (t) =

∇τ vi , τi

(Π2,i v) (t) =

∇τ vi 2(t − ti−1/2 ) ρi ∇τ vi+1 − ∇τ vi , + τi τi (τi + τi+1 )

(2.1)

where ∇τ vi = vi − vi−1 . For convenience, we denote ωβ (t) = tβ−1 /Γ(β) and Qn (t) = −ω2−α (tn−ϑ − t). It is obvious that, for 0 ≤ t ≤ tn−ϑ , Qn (t) < 0 and for 0 ≤ t < tn−ϑ , Qn (t) = ω1−α (tn−ϑ − t) > 0,

Qn (t) = −ω−α (tn−ϑ − t) > 0,

Q n (t) = ω−α−1 (tn−ϑ − t) > 0.

Then, the nonuniform L2-1σ approximation to the Caputo derivative C0 Dαtn+σ u is defined by [13]  tn−ϑ n−1  ti  α n−ϑ   = Qn (s)(Π1,n v) (s)ds + Qn (s)(Π2,i v) (s)ds (Dτ v) tn−1

n = a(n) 0 ∇τ v +

i=1

n−1  

ti−1

 (n) (n) i i+1 i ∇ v + ρ b ∇ v − b ∇ v a(n) , i n−i τ n−i τ n−i τ

(2.2)

i=1

where



tn−ϑ

a(n) 0

1 = τn

b(n) n−i

2 = τi (τi + τi+1 )

Qn (s)ds,

tn−1



ti

ti−1

a(n) n−i

1 = τi



ti

Qn (s)ds,

1 ≤ i ≤ n − 1;

ti−1

(s − ti− 12 )Qn (s)ds,

1 ≤ i ≤ n − 1.

Then, we obtain the following compact form (Dτα v)n−ϑ =

n 

i A(n) n−i ∇τ v ,

(2.3)

i=1

A(n) n−i

is defined as where the discrete convolution kernel n ≥ 2, ⎧ (n) ⎪ ⎪ a0 + ρn−1 b(n) ⎪ 1 , ⎪ ⎪ ⎨ (n) (n) (n) An−i = ⎪ an−i + ρi−1 b(n) ⎪ n−i+1 − bn−i , ⎪ ⎪ ⎪ ⎩a(n) − b(n) , n−1 n−1

(1) follows: A(1) 0 = a0 for n = 1, and for

i = n, 2 ≤ i ≤ n − 1, i = 1.

2.2. Discrete fractional Gr¨onwall inequality and global consistency errors In the following discrete fractional Gr¨onwall inequality, a complementary discrete kernel P(n) n−k is involved and satisfies [10] n  ( j) P(n) 1 ≤ i ≤ n ≤ N. (2.4) n− j A j−i ≡ 1, j=i

From (2.4), we have P(n) 0 =

1

; A(n) 0

P(n) n− j =

n 1 

A(0j)

(i) (n) (A(i) i− j−1 − Ai− j )Pn−i ,

i= j+1

4

1 ≤ j ≤ n − 1.

(2.5)

Lemma 2.1. (See [13].) Assume that the discrete convolution kernels A(n) n−i hold the following three properties: (n) A1. The discrete kernel is monotone, which is to say, A(n) i−2 ≥ Ai−1 > 0 for 2 ≤ i ≤ n ≤ N.

A2. There exists a constant πA > 0, holding that  ti 1 ω1−α (tn − s) (n) ds, An−i ≥ πA ti−1 τi

1 ≤ i ≤ n ≤ N.

(2.6)

A3. There exists a positive constant ρ, such that the step ratios ρi ≤ ρ, 1 ≤ i ≤ N − 1. N−1 are non-negative constants independent of the Suppose the offset parameter ϑ ∈ [0, 1), and {λl }l=0 time steps. Meanwhile, Assume further that there exists a constant Λ (independent of the time-step √α N−1 sizes) such that Λ ≥ l=0 λl , and the maximum step size satisfies τ ≤ 1/ 2πA Γ(2 − α)Λ. Then, if N N the non-negative sequences {ξi }i=1 and {vi }i=0 satisfy that n 

i 2 A(n) n−i ∇τ (v ) ≤

i=1

we have

n 

λn−k ( vn−ϑ )2 + vn−ϑ gn ,

1 ≤ n ≤ N,

(2.7)

k=1

i    j P(i) g , vn ≤ 2Eα (2 max(1, ρ)πA Λtnα ) v0 + max i− j 1≤i≤n

1 ≤ n ≤ N.

(2.8)

j=1

Remark 2.1. About Lemma 2.1, we here state the whole result as the same as Ref. [13]. However, we in this paper only consider the case of λk = 0, k = 0, 1, . . . , n − 1. According to the detailed proof of Lemma 2.1, if λk = 0, the discrete inequality (2.8) holds in a simpler form, which only requires A1 and A2 but no restrictions on time steps, 

v ≤ 2 v + max n

0

1≤i≤n

i 

 j P(i) g , i− j

1 ≤ n ≤ N.

(2.9)

j=1

From [10], we know that if A1 holds, P(n) n− j is non-negative and well-defined. Moreover, if A2 holds, then one obtains n  P(n) 1 ≤ n ≤ N. (2.10) n− j ≤ πA ω1+α (tn ), j=1

We here show two assumptions, which are going to be used repeatedly in the next theoretical analysis: M1. The parameter ϑ = α/2 = 1 − σ, and the maximum time-step ratio ρ = 7/4. M2. There is a positive constant Cγ , such that τi ≤ Cγ τ min{1, ti1−1/γ } for 1 ≤ i ≤ N, ti ≤ Cγ ti−1 and τi /ti ≤ Cγ τi−1 /ti−1 for 2 ≤ i ≤ N. 5

Remark 2.2. (1) In order to achieve second-order temporal accuracy, we have to ask for the assumption M1. Once we choose ϑ  α/2, the consistency error would be limited to an order of O(τ2−α n ), even for the model with smooth solutions. (2) In M2, the parameter γ ≥ 1 controls the concentrated extent of the time levels near t = 0. It is obvious that M2 holds with γ = 1, in which case, the mesh is quasi-uniform. As γ increases, the initial step sizes will become smaller than the later ones. Here, the graded mesh tk = (k/N)γ T is a natural choice to satisfy the assumption M2. From [13], if M1 holds, we have the following useful properties of the discrete kernel A(n) n−i (0 ≤ i ≤ n): P1. A(n) n−i are bounded, such that

 tn 24 ω1−α (tn − s)ds, 11τn tn−1  ti 4 ≥ ω1−α (tn − s)ds, 11τi ti−1

A(n) 0 ≤ A(n) n−i P2. A(n) n−i are monotone, such that A(n) n−i−1

−

A(n) n−i

≥ (1 +

ρi )b(n) n−i

1 + 5τi



ti

(2.11) 1 ≤ i ≤ n;

(ti − s)Qn (s)ds,

1 ≤ i ≤ n − 1;

(2.12)

(2.13)

ti−1

(n) (n) (n) P3. For n ≥ 2, we have A(n) 0 − A1 ≥ ϑ(2A0 − A1 ).

The property P1 implies A2 holds with πA = 11/4. Meanwhile, the property P2 ensures A1 valid, and has been used to establish a stronger estimate in the error analysis. For a fixed function u(t), denote the local consistency error by Υn−ϑ = C0 Dαtn−ϑ u − (Dτα u)n−ϑ =

n 

Υn−ϑ , i

1 ≤ n ≤ N,

(2.14)

i=1

where

 Υn−ϑ = i = Υn−ϑ n

ti



Qn (s)Π 2,i v (s)ds,

t  i−1tn−ϑ



 Qn (s)Π 1,n v (s)ds,

1 ≤ i ≤ n − 1 ≤ N − 1, 1 ≤ n ≤ N.

tn−1

Lemma 2.2. ([13]). Assume (1.5) and the mesh condition M1 hold. Then the global consistency error satisfies n  τα   1 1 α−3 3 α α−3 3−α j−ϑ + t max P(n) |Υ | ≤ C τ + (t − t ) t τ , 1 ≤ n ≤ N. (2.15) v i 1 1 2 i−1 i n− j α 1 − α 2≤i≤n j=1 6

Lemma 2.3. ([13]). Denote R n−ϑ = un−ϑ − un−ϑ , 1 ≤ n ≤ N be the local truncation error. Assume (1.5) and the mesh condition M1 hold. Then the global consistency error satisfies n 

j−ϑ P(n) | ≤ Cv n− j |R

 τ2α

j=1

1

α

 α−2 2 + tnα max ti−1 τi , 2≤i≤n

1 ≤ n ≤ N.

(2.16)

Remark 2.3. From Lemmas 2.2 and 2.3, following similar analysis as [13], we also obtain that n 

j−ϑ P(n) q ≤ Cv n− j Υ

j=1

 τα 1

α

+ t1α−3 τ32 +

 1 α−3 3−α max(ti − t1 )α ti−1 τi , 1 − α 2≤i≤n

1 ≤ n ≤ N and q = 0, 1, 2, (2.17)

and n 

j−ϑ P(n) q ≤ Cv n− j R

 τ2α 1

α

j=1

 α−2 2 + tnα max ti−1 τi ,

1 ≤ n ≤ N and q = 0, 1, 2, 3, 4.

2≤i≤n

(2.18)

Lemma 2.4. Under the condition M1, we have

vn−ϑ (Dτα v)n−ϑ ,

h

1  (n) A ∇τ (vi 20 ). 2 i=1 n−i n

≥

(2.19)

3. Convergence and superconvergence analysis of the bilinear FEM In this section, the convergence and superconvergence analysis for the bilinear FEM will be given in details. 3.1. Bilinear FEM Denote Th1 be a family of rectangular meshes of Ω with Ω = ∪K∈Th1 K. For each K ∈ Th1 , denote OK = (xK , yK ) be the center of K, and hKx and hyK be the perpendicular distances between the center OK and two sides of K paralleling to the two coordinate planes, respectively. A1K = (xK − hKx , yK − hyK ), A2K = (xK + hKx , yK − hyK ), A3K = (xK + hKx , yK + hyK ) and A4K = (xK + hKx , yK + hyK ) represent the four vertices of K. Denote hK = max{hKx , hyK } and h = maxK∈Th1 hK . Then, we define the following finite element spaces by Wh = {v; v|K ∈ Q11 (K)},

W0h = {v ∈ Wh ; v|∂Ω = 0},

(3.1)

where Qi j (K) = span{xr y s ; (x, y) ∈ K, 0 ≤ r ≤ i, 0 ≤ s ≤ j}. Next, we define the L2 -projector Ph : L2 (Ω) → W0h by (Ph w, vh ) = (w, vh ),

∀vh ∈ W0h .

(3.2)

As the result of the quasi-uniform of the meshes, the operator Ph holds that [19] ∇Ph v0 ≤ CΩ ∇v0 , 7

∀v ∈ H01 (Ω).

(3.3)

Define the Ritz projection Rh1 : H01 (Ω) → W0h by (∇Rh1 w, ∇vh ) = (∇w, ∇vh ),

∀vh ∈ W0h ,

(3.4)

∀v ∈ H 2 (Ω) ∩ H01 (Ω).

(3.5)

The operator Rh1 holds that v − Rh1 v0 + hv − Rh1 vh ≤ CΩ h2 |v|2 ,

Define Ih : H 2 (Ω) → W0h be the associated interpolation operator such that Ih | K = I K ,

IK v(di ) = v(di ),

i = 1, 2, 3, 4.

Following lemmas are necessary for the stability and convergence analysis [20]. Lemma 3.1. Let u ∈ H 2 (Ω) ∩ H01 (Ω). Under anisotropic meshes, there holds u − Ih u0 + hu − Ih u1 ≤ CΩ h2 u2 .

(3.6)

Furthermore, if u ∈ H 3 (Ω) ∩ H01 (Ω), we have (∇(u − Ih u), ∇vh ) ≤ CΩ h2 |u|3 vh 1 , ∀vh ∈ Vh .

(3.7)

Lemma 3.2. For u ∈ H 3 (Ω) ∩ H01 (Ω), we have Rh1 u − Ih u1 ≤ CΩ h2 u3 .

(3.8)

We first introduce the following weak formulation of (1.1)-(1.3), which is to find u ∈ H01 (Ω), such that (C0 Dαt u, v) + (Δu, v) = ( f, v), ∀v ∈ H01 (Ω), (3.9) n−ϑ = ϑU n + (1 − ϑ)U n−1 , 1 ≤ n ≤ N. Then, with the initial value u(x, 0) = u0 (x), x ∈ Ω. Denote U we arrive at the following fully discrete scheme, which is to find U n ∈ W0h , 1 ≤ n ≤ N, such that

n−ϑ , ∇vh ) = (Ph f n−ϑ , vh ), (Dτα U)n−ϑ , vh + (∇U



∀vh ∈ Wh ,

(3.10)

with the initial value U 0 = Ih u0 . Define the discrete Laplacian Δh : W0h → W0h by [21] (Δh v, w) = −(∇v, ∇w),

∀v, w ∈ W0h .

(3.11)

There holds that [21] Δh Rh1 v = Ph Δv,

∀v ∈ H 2 (Ω).

(3.12)

n−ϑ = Ph f n−ϑ . (Dτα U)n−ϑ − Δh U

(3.13)

The fully discrete scheme can be written as

8

3.2. Sharp L2 -norm optimal error estimate In this subsection, we intend to give the sharp L2 -norm optimal error estimate of the full discrete scheme. Theorem 3.1. Assume (1.5) holds for j = 2. Then, if M1 holds, the fully discrete solution U n is convergent with respect to the following discrete L2 -norm,  τα  1 1 n n α α−3 3−α α α−2 2 2 + max(ti − t1 ) ti−1 τi + tn max ti−1 τi + h . U − u 0 ≤ CΛ (3.14) 2≤i≤n α 1 − α 2≤i≤n In particular, if M2 also holds, we have U n − un 0 ≤ CΛ (τmin{γα,2} + h2 ).

(3.15)

Proof. For convenience, we denote un − U n = (un − Rh1 un ) + (Rh1 un − U n ) = ηn + ξn .

(3.16)

Then, from (1.1), (3.12) and (3.13), we obtain ξn−ϑ = (Rh1 − Ph )(Dτα u)n−ϑ + Ph (Dτα u)n−ϑ − Ph Δ (Dτα ξ)n−ϑ − Δh un−ϑ − Ph f n−ϑ = −Ph (I − Rh1 )(Dτα u)n−ϑ + Ph (C0 Dαtn−ϑ u − Υn−ϑ ) − Ph (Δun−ϑ − R n−ϑ ) − Ph f n−ϑ = −Ph (Dτα η)n−ϑ + Ph (R n−ϑ − Υn−ϑ ).

(3.17)

Multiplying (3.17) by  ξn−ϑ and integrating over Ω, and then by using Green formula, we obtain that ((Dτα ξ)n−ϑ ,  ξn−ϑ ) + ∇ ξn−ϑ 20 = −((Dτα η)n−ϑ ,  ξn−ϑ ) + (R n−ϑ − Υn−ϑ ,  ξn−ϑ ). (3.18) By virtue of Cauchy-Schwarz inequality and (3.5), and removing the positive term ∇ ξn−ϑ 20 , we derive ((Dτα ξ)n−ϑ ,  ξn−ϑ ) ≤ Ch2 |(Dτα u)n−ϑ |2  ξn−ϑ 0 + (R n−ϑ 0 + Υn−ϑ 0 ) ξn−ϑ 0 . (3.19) From Lemma 2.4, (3.19) and (1.6), it holds that n 

i 2 2 C α n−ϑ A(n) |2 ) ξn−ϑ 0 + 2(R n−ϑ 0 + Υn−ϑ 0 ) ξn−ϑ 0 n−i ∇τ (ξ 0 ) ≤ 2C Ω h (|0 Dtn−ϑ u|2 + |Υ

i=1

≤



 2CΩ Lh2 + 2CΩ h2 |Υn−ϑ |2 + 2R n−ϑ 0 + 2Υn−ϑ 0  ξn−ϑ 0 .

(3.20)

By Lemma 2.1, we get i     2 2 j−ϑ j−ϑ j−ϑ P(i) Lh + 2C h |Υ | + 2R  + 2Υ  ξn 0 ≤ 2 ξ0 0 + max 2C . Ω Ω 2 0 0 i− j 1≤i≤n

(3.21)

j=1

Then, by using Lemma 3.2, (2.17) and (2.18), we obtain   τα 1 1 n α−3 3 α α−3 3−α α α−2 2 2 + t 1 τ2 + max(ti − t1 ) ti−1 τi + tn max ti−1 τi + h . ξ 0 ≤ CΛ 2≤i≤n α 1 − α 2≤i≤n 9

(3.22)

Therefore, by virtue of the triangular inequality and (3.5), we conclude (3.14). For convenience, we denote β := min{2, γα}. The mesh assumption M2 implies that τ1 ≤ Cγ τγ . In addition, we obtain

β t1α−3 τ32 = t1α−3 τ3−α−β (t2 − t1 )α τβ2 ≤ t1α−3 τ3−α−β (t2 − t1 )α Cγ τ min{1, t21−1/γ } 2 2

β ≤ (Cγ )3−α+β τ3−α−β t2α−3+α (τ min{1, t21−1/γ } ≤ (Cγ )3−α+β t2α−β/γ (τ2 /t2 )3−α−β τβ 2 ≤ (Cγ )3−α+β t2α−β/γ τβ .

(3.23)

Moreover, it holds that β 1−1/γ β α−3 3−α β α− γ β (ti − t1 )α ti−1 τi ≤ tiα+α−3 τ3−α−β τ min{1, t } ≤ (C ) t τ, C γ γ i i i

which implies that

α−3 3−α (ti − t1 )α ti−1 τi ≤ (Cγ )β tiα−min{(3−α)/γ,α} τmin{2,γα} .

(3.24)

Furthermore, we have β β

β α−2 2 α−2 2−β β 2+β−α τi ≤ ti−1 τi Cγ τ min{1, ti1−1/γ } ≤ (Cγ )2+β−α (ti )α− γ +β−2 τ2−β (ti )α− γ τβ . ti−1 i τ ≤ (C γ )

(3.25)

From (3.22), (3.23), (3.24) and (3.25), if tn ≤ T , we obtain ξn 0 ≤ CΛ (τmin{γα,2} + h2 ).

(3.26)

Then, the triangular inequality and (3.5) imply that (3.15) holds.  3.3. Sharp H 1 -norm superconvergence analysis In order to obtain the H 1 -norm superconvergence result, we first introduce the following lemma. Lemma 3.3. (See [12].) Assume that the discrete convolution kernels A(n) n−i hold A1-A3. Suppose N−1 the offset parameter ϑ ∈ [0, 1), and {λl }l=0 are non-negative constants independent of the time steps. Meanwhile, Assume further that there exists a constant Λ (independent of the time-step √ N−1 sizes) such that Λ ≥ l=0 λl , and the maximum step size satisfies τ ≤ 1/ α 2πA Γ(2 − α)Λ. Then, if N N N the non-negative sequences {gi }i=1 , {si }i=1 and {vi }i=0 satisfy that n 

i 2 A(n) n−i ∇τ (v )

i=1

≤

n 

λn−k (vn−ϑ )2 + vn−ϑ gn + (sn )2 ,

1 ≤ n ≤ N,

(3.27)

k=1

or (vn )2 ≤ (v0 )2 +

n  j=1

P(n) n− j

j 

λ j−k (vn−ϑ )2 +

n 

j−ϑ j P(n) g + n− j v

j=1

k=1

n 

j 2 P(n) n− j (s ) ,

1 ≤ n ≤ N,

(3.28)

j=1

we have i     α/2 i j P(i) g + π Γ(1 − α) max {t s } , vn ≤ 2Eα (2 max(1, ρ)πA Λtnα ) v0 +max A i i− j 1≤i≤n

1≤i≤n

j=1

10

1 ≤ n ≤ N. (3.29)

Theorem 3.2. Assume (1.5) holds for j = 3. Then, if M1 holds, the fully discrete solution U n is convergent with respect to the following H 1 -norm,  τα  1 α−3 3−α α−2 2 max(ti − t1 )α ti−1 τi + tnα max ti−1 τi + h 2 . (3.30) |U n − Rh1 un |1 ≤ CΛ 1 + t1α−3 τ32 + 2≤i≤n α 1 − α 2≤i≤n In particular, if M2 also holds, we have |U n − Rh1 un |1 ≤ CΛ (τmin{γα,2} + h2 ).

(3.31)

ξn−ϑ and integrating over Ω, and then by using Green formula Proof. Multiplying (3.17) by −Δh and (3.11), we obtain that ξn−ϑ 20 = (Ph (Dτα η)n−ϑ , Δh ξn−ϑ ) − (Ph (R n−ϑ − Υn−ϑ ), Δh ξn−ϑ ) ((Dτα ∇ξ)n−ϑ , ∇ ξn−ϑ ) + Δh

(3.32)

Obviously, we have (Dτα η)n−ϑ = (Dτα η)n−ϑ − C0 Dαtn−ϑ η + C0 Dαtn−ϑ η = −Υn−ϑ + Rh1 Υn−ϑ + C0 Dαtn−ϑ η.

(3.33)

Therefore, by using Green formula, Cauchy-Schwarz inequality and Young’s inequality, and utilizing the definition of Rh1 (3.4), we obtain ξn−ϑ ) = −(Υn−ϑ , Δh ξn−ϑ ) + (Rh1 Υn−ϑ , Δh ξn−ϑ ) + (C0 Dαtn−ϑ η, Δh ξn−ϑ ) (Ph (Dτα η)n−ϑ , Δh C0 Dαtn−ϑ η20 n−ϑ n−ϑ ξn−ϑ 20 ≤ |Υ |1 |ξ |1 + + Δh 4 (CΩ )2 L2 4 ξn−ϑ 20 . h + Δh ≤ |Υn−ϑ |1 | ξn−ϑ |1 + (3.34) 4 Similarly, we have (Ph (R n−ϑ − Υn−ϑ ), Δh ξn−ϑ ) ≤ (|R n−ϑ |1 + |Υn−ϑ |1 )| ξn−ϑ |1 .

(3.35)

Substituting (3.34) and (3.35) into (3.32), and by virtue of Lemma 2.4, we get n 

i2 n−ϑ A(n) |1 + 2|Υn−ϑ |1 )| ξn−ϑ |1 + n−i ∇τ (|ξ |1 ) ≤ (|R

i=1

(CΩ )2 L2 4 h. 4

(3.36)

Then, by Lemma 3.3, we have 

|ξ |1 ≤ 2 |ξ |1 + max n

0

1≤i≤n

i  j=1

j−ϑ P(i) |1 i− j (|R

+ 2|Υ

j−ϑ

√  CΩ L 11Γ(1 − α) max{tiα/2 }h2 , |1 ) + 1≤i≤n 4

1 ≤ n ≤ N. (3.37)

Therefore, by using Lemma 3.2, (2.17) and (2.18), we conclude  τα  1 1 n α−3 3 α α−3 3−α α α−2 2 2 + t 1 τ2 + max(ti − t1 ) ti−1 τi + tn max ti−1 τi + h . |ξ |1 ≤ CΛ 2≤i≤n α 1 − α 2≤i≤n

(3.38)

Hence, (3.30) holds. If the mesh assumption M2 holds, from (3.23)-(3.25), we get (3.31).  By using Theorem 3.2, we obtain the following superclose results. 11

Theorem 3.3. Assume (1.5) holds for j = 3. Then, if M1 holds, the fully discrete solution U n is convergent with respect to the following H 1 -norm, U n − Ih un 1 ≤ CΛ

 τα 1

α

 1 α−3 3−α α−2 2 max(ti − t1 )α ti−1 τi + tnα max ti−1 τi + h 2 . 2≤i≤n 1 − α 2≤i≤n

+ t1α−3 τ32 +

(3.39)

In particular, if M2 also holds, we have U n − Ih un 1 ≤ CΛ (τmin{γα,2} + h2 ).

(3.40)

U n − Ih un 1 ≤ U n − Rh1 un 1 + Rh1 un − Ih un 1 .

(3.41)

Proof. Obviously, we have

Then, by virtue of (3.22) and (3.30), we get U n − Rh1 un 1 ≤ CΛ

 τα 1

α

+ t1α−3 τ32 +

 1 α−3 3−α α−2 2 τi + tnα max ti−1 τi + h 2 . max(ti − t1 )α ti−1 2≤i≤n 1 − α 2≤i≤n

(3.42)

By using Lemma 3.2, (3.41) and (3.42), it follows that U n − Ih un 1 ≤ CΛ

 τα 1

α

+ t1α−3 τ32 +

 1 α−3 3−α α−2 2 max(ti − t1 )α ti−1 τi + tnα max ti−1 τi + h 2 . 2≤i≤n 1 − α 2≤i≤n

(3.43)

Thus, (3.39) is obtained. If the mesh assumption M2 also holds, we derive (3.40).  In order to gain the global superconvergence result, a new family of meshes Γ2h is constructed,  consists of four neighboring elements belonging to Γh (see Fig. 1). Define and the new element K

 Figure 1: New element K.

the interpolation postprocessing operator I2h satisfying  I2h u|K ∈ Q22 (K), I2h Ih u = I2h u, ∀u ∈ H 2 (Ω), I2h u − u1 ≤ CΩ h2 u3 , ∀u ∈ H 3 (Ω) ∩ H01 (Ω), I2h vh 1 ≤ CΩ vh 1 , ∀vh ∈ V0h ,  denotes the space of biquadratic polynomials on K.  where Q22 (K) 12

(3.44) (3.45) (3.46) (3.47)

Theorem 3.4. Assume (1.5) holds for j = 3. Then, if M1 holds, the fully discrete solution U n is convergent with respect to the following H 1 -norm, I2h U n − un 1 ≤ CΛ

 τα 1

α

+ t1α−3 τ32 +

 1 α−3 3−α α−2 2 max(ti − t1 )α ti−1 τi + tnα max ti−1 τi + h 2 . 2≤i≤n 1 − α 2≤i≤n

(3.48)

In particular, if M2 also holds, we have I2h U n − un 1 ≤ CΛ (τmin{γα,2} + h2 ).

(3.49)

Proof. Utilizing (3.44)-(3.47) results in I2h U n − un 1 ≤ I2h U n − I2h Ih un 1 + I2h Ih un − un 1 ≤ CΩ U n − Ih un 1 + CΩ h2 u3 . Therefore, (3.48) and (3.49) are obtained by using (3.39) and (3.40), respectively.  Remark 3.1. This method is motivated by the work given in [22]. We notice that this method can be used to deal with the H 1 -norm error estimate of all the conforming FEMs. However, it still cannot handle the nonconforming FEMs as the result of the existence of the consistency term. In the following section, we will construct the nonconforming FEMs of the system (1.1)-(1.3). For the constructed scheme, we will get the sharp L2 -norm optimal error estimate. However, for the H 1 -norm error estimate, although we can get the corresponding superclose and superconvergence results, it is still suboptimal in the temporal direction. 4. Convergence and superconvergence analysis of the modified quasi-Wilson FEM In this section, we will give the detailed convergence analysis of the modified quasi-Wilson FEM for the system (1.1)-(1.3). 4.1. Modified quasi-Wilson element Let Kˆ = [−1, 1] × [−1, 1] be a reference element on ξ − η plane. Denote Aˆ 1 = (−1, −1), ˆ Define lˆ1 = Aˆ 1 Aˆ 2 , lˆ2 = Aˆ 2 Aˆ 3 , Aˆ 2 = (1, −1), Aˆ 3 = (1, 1), Aˆ 1 = (−1, 1) be the four vertices of K. ˆ Then, we define the modified quasi-Wilson FE lˆ3 = Aˆ 3 Aˆ 4 and lˆ4 = Aˆ 4 Aˆ 1 be four edges of K. ˆ ˆ ˆ {K, P, Σ} as follows   ˆ ˆ ψ(η) , Σˆ = {ˆv1 , vˆ 2 , vˆ 3 , vˆ 4 , vˆ 5 , vˆ 6 }, Pˆ = span Ni (ξ, η), i = 1, 2, 3, 4, ψ(ξ), where Ni (ξ, η) = (ξ4 , η4 ) = (−1, 1),

1 (1 4

+ ξi ξ)(1 + ηi η), (ξ1 , η1 ) = (−1, −1), (ξ2 , η2 ) = (1, −1), (ξ3 , η3 ) = (1, 1), ˆ = − 1 + 1 s 2 − 7 s 4 + 7 s6 , ψ(s) 30 2 6 10

and vˆ i = vˆ (Aˆ i ),

i = 1, 2, 3, 4;

1 vˆ 5 = ˆ |K| 13

 Kˆ

∂2 vˆ dξdη, ∂ξ2

1 ˆ |K|

 Kˆ

∂2 vˆ dξdη. ∂η2

ˆ h vˆ can be expressed uniquely as ˆ the associated interpolation Π For vˆ ∈ H 2 (K), ˆ h vˆ (ξ, η) = Π

4 

ˆ v5 + ψ(η)ˆ ˆ v6 . Ni (ξ, η)ˆvi + ψ(ξ)ˆ

(4.1)

i=1

 Denote Th2 be a convex quadrilateral decomposition of Ω, such that Ω = K∈Th2 K satisfying the regular condition [23]. For K ∈ Th2 , we denote Ai (xi , yi ), i = 1, 2, 3, 4 be its four vertices, and let hK =diam(K), h = maxK∈Th2 hK . Obviously, there exists a invertible mapping F K∗ : Kˆ → K, such that 4 4   x= Ni (ξ, η)xi , y = Ni (ξ, η)yi . i=1

i=1

ˆ = K and F ∗ (Aˆ i ) = Ai , i = 1, 2, 3, 4. Then, we define the following finite Thus, we have F K∗ (K) K element space and associated interpolation operator   Vh := v : vˆ |Kˆ = v|K ◦ F K∗ , ∀K ∈ Th2 , v(A) = 0, ∀node A ∈ ∂Ω , (4.2) and Πh : v ∈ H 2 (Ω) → vh ∈ Vh ,

Πh |K = ΠK ,

ˆ v ◦ (F K∗ )−1 . ΠK v = Πˆ

(4.3)

For any vh ∈ Vh , from [24], we can split vh = v¯ h + v˜ h , in which v¯ h and v˜ h denote the conforming  and nonconforming parts, respectively. It is obvious that  · h = ( K| · |21,K )1/2 is a norm over Vh .  Denote ∇h be the gradient operator piecewisely and (γ1 , γ2 )h = K K γ1 · γ2 dx. From [25], we have the following features on quadrilateral meshes. Lemma 4.1. For any vh ∈ Vh , K ∈ Th2 , there hold  q˜vh dX = 0, ∀q ∈ P1 (K), K  ∂˜vh ∂˜vh dX = dX = 0, ∀q ∈ P2 (K), q q ∂x ∂y K K vh 2h = ¯vh 2h + ˜vh 2h , ˜vh h ≤ CΩ hvh h ,   ∂u vh ds ≤ CΩ h2 |u|3 vh h , ∂n ∂K K∈T

(4.4) (4.5) (4.6) (4.7)

h2

where n is unit normal vector on ∂K. Lemma 4.2. Assume that u ∈ H 4 (Ω) ∩ H01 (Ω). Then for any vh ∈ Vh and K ∈ Th2 , we have     ∂u   ≤ C h3 |u| v  . ds (4.8) v h Ω 4 h h   ∂n ∂K K Lemma 4.3. Let u ∈ H 3 (Ω) ∩ H01 (Ω). Then, for any vh ∈ Vh , there holds (∇h (u − Πh u), ∇h vh )h ≤ CΩ h2 |u|3 vh h . 14

(4.9)

Let Rh2 : H 2 (Ω) → Vh be a Ritz projection operator, such that (∇h (u − Rh2 u), ∇h vh )h = 0,

∀vh ∈ Vh .

(4.10)

u − Rh2 uh ≤ CΩ hu2 .

(4.11)

Then, from [23], it follows that u − Rh2 u0 ≤ CΩ h2 u2 ,

n−ϑ = ϑU n + (1 − ϑ)U n−1 , 1 ≤ n ≤ N. Then, we arrive at the following fully discrete Denote U scheme, which is to find U n ∈ Vh , 1 ≤ n ≤ N, such that

n−ϑ , ∇h vh )h = ( f n−ϑ , vh )h , (Dτα U)n−ϑ , vh h + (∇h U

∀vh ∈ Vh ,

(4.12)

with the initial value U 0 = Rh2 u0 . 4.2. L2 -norm optimal error estimate From (1.1)-(1.3) and (4.12), we get the following error equation

(Dτα ξ)n−ϑ , vh h + (∇h ξn−ϑ , ∇h vh )h = − ((Dτα η)n−ϑ , vh )h − (∇h ηn−ϑ , ∇h vh )h − (∇h R n−ϑ , ∇h vh )h   ∂ un−ϑ vh ds. − (Υn−ϑ , vh )h + ∂K ∂n K (4.13)

Theorem 4.1. Assume (1.5) holds for j = 4. Then, if M1 holds, the fully discrete solution U n is convergent with respect to the following discrete L2 -norm, U n − un 0 ≤ Ch2 |u|2 + CΛ

 τα 1

α

+

 1 α−3 3−α α−2 2 max(ti − t1 )α ti−1 τi + tnα max ti−1 τi + tnα h2 . 2≤i≤n 1 − α 2≤i≤n

(4.14)

In particular, if M2 also holds, we have U n − un 0 ≤ CΛ (τmin{γα,2} + h2 ). ξn−ϑ in (4.13), we obtain Proof. Substituting vh =  α n−ϑ n−ϑ  n−ϑ 2  0 (Dτ ξ) , ξ h + ∇h ξ ξn−ϑ )h − (∇h R n−ϑ , ∇h ξn−ϑ )h − (Υn−ϑ ,  ηn−ϑ , ∇h = −((Dτα η)n−ϑ ,  ξn−ϑ )h − (∇h ξn−ϑ )h   ∂ un−ϑ n−ϑ ξ ds. + ∂n ∂K K

(4.15)

(4.16)

By using Lemma 2.4, we have

1  (n) ξn−ϑ 2h ≥ ξn−ϑ 2h . ξn−ϑ h + ∇h A ∇τ (ξi 20 ) + ∇h (Dτα ξ)n−ϑ ,  2 i=1 n−i n



15

(4.17)

It is obvious that (Dτα η)n−ϑ = C0 Dαtn−ϑ η −



C α 0 Dtn−ϑ η

 − (Dτα η)n−ϑ = C0 Dαtn−ϑ η − Υn−ϑ + Rh2 Υn−ϑ .

(4.18)

By Cauchy-Schwarz inequality and (4.18), we obtain   ((Dτα η)n−ϑ ,  ξn−ϑ )h ≤ C0 Dαtn−ϑ η0 + Υn−ϑ 0 + Rh2 Υn−ϑ 0  ξn−ϑ 0 .

(4.19)

Assume Rh2 v0 ≤ CΩ v0 . Then, by using (1.6), (4.11) and (4.19), we get   ((Dτα η)n−ϑ ,  ξn−ϑ )h ≤ CΩ Lh2 + (1 + CΩ )Υn−ϑ 0  ξn−ϑ 0 .

(4.20)

By virtue of (4.10), we have

(∇h ξn−ϑ )h = 0. ηn−ϑ , ∇h

(4.21)

Applying Green’s formula obtains (ΔR

n−ϑ

, ξn−ϑ )h =

 K

∂K

∂R n−ϑ n−ϑ ξn−ϑ )h . ξ ds − (∇h R n−ϑ , ∇h ∂n

Hence, from Lemma 4.2 and inverse inequality, we obtain     ∂R n−ϑ n−ϑ  n−ϑ n−ϑ n−ϑ n−ϑ    ξ ds − (ΔR , ξ )h  (∇h R , ∇h ξ )h =  ∂n ∂K K ξn−ϑ h + |R n−ϑ |2  ξn−ϑ 0 ≤ CΩ h3 |R n−ϑ |4  ≤ CΩ h2 |R n−ϑ |4  ξn−ϑ 0 + |R n−ϑ |2  ξn−ϑ 0 .

(4.22)

(4.23)

Due to the Cauchy-Schwarz inequality, we have (Υn−ϑ ,  ξn−ϑ )h ≤ Υn−ϑ 0  ξn−ϑ 0 . Then by using (1.6), Lemma 4.2 and inverse inequality again, the last term on the RHS of (4.16) can be bounded by     ∂ un−ϑ n−ϑ    ≤ CΩ h3 | ds un−ϑ |4  ξn−ϑ h ≤ CΩ h2 | un−ϑ |4  ξn−ϑ 0 ≤ CΩ Lh2  ξn−ϑ 0 . (4.24) ξ  ∂n ∂K K Substituting (4.17), (4.20), (4.21) and (4.23) into (4.16) yields   1  (n) ξn−ϑ 2h ≤ (CΩ + CΩ L)h2 +Υn−ϑ 0 +CΩ h2 |R n−ϑ |4 +|R n−ϑ |2  An−i ∇τ (ξi 20 )+∇h ξn−ϑ 0 . (4.25) 2 i=1 n

ξn−ϑ 2h , (4.25) can be written as Removing the positive term ∇h n 

  i 2 2 n−ϑ 2 n−ϑ n−ϑ n−ϑ 0 . A(n) ∇ (ξ  ) ≤ 2(C + C L)h + 2Υ  + 2C h |R | + 2|R | Ω Ω 0 Ω 4 2 ξ 0 n−i τ

i=1

16

(4.26)

By using Lemma 2.1, we get i     n 0 2 j−ϑ 2 j−ϑ j−ϑ ξ 0 ≤ 2 ξ 0 + max P(i) + C L)h + 2Υ  + 2C h |R | + 2|R | 2(C . (4.27) Ω Ω 0 Ω 4 2 i− j 1≤i≤n

j=1

By virtue of (2.17) and (2.18), we have  τα  1 α−3 3−α α−2 2 max(ti − t1 )α ti−1 τi + tnα max ti−1 τi + tnα h2 . ξn 0 ≤ CΛ 1 + t1α−3 τ32 + 2≤i≤n α 1 − α 2≤i≤n If the mesh assumption M2 holds, from (3.23)-(3.25), it follows that ξn 0 ≤ CΛ (τmin{γα,2} + h2 ).

(4.28)

(4.29)

Then, by virtue of (4.29) and (4.11), we have completed the proof of Theorem 4.1.  4.3. H 1 -norm superconvergence analysis In this subsection, we give in details the H 1 -norm superclose and superconvergence results of the discrete scheme (4.12). To this end, we first show the following lemmas. Lemma 4.4. ([12].) Under the same conditions of Lemma 2.2, we have n  Cv Cv j−ϑ 2 2α−6 6 2α P(n) | ≤ 2 τα1 + t1α−6 τ62 + max tiα ti−1 τi /τi−1 , 1 ≤ n ≤ N. n− j |Υ 2≤i≤n α 1 − α j=1 Lemma 4.5. Under the same conditions of Lemma 2.2, we have n  τ3α   1 j−ϑ 2 α 2α−4 4 P(n) |R | ≤ C + t max t τ 1 ≤ n ≤ N. v n i , i−1 n− j 2≤i≤n α2 j=1 Proof. By using the Taylor formula, it is easy to verify that  t j−ϑ  j−ϑ  (s − t j−1 )v (s)ds − (1 − ϑ) R = −ϑ

tj

t j−1

t j−ϑ

(4.30)

(4.31)

(t j − s)v (s)ds.

(4.32)

2 ≤ j ≤ N.

(4.33)

Under the regularity assumption, we have |R 1−ϑ |2 ≤ Cv

τ2α 1 , α2

4 |R j−ϑ |2 ≤ Cv t2α−4 j−1 τ j ,

From [13], we have P(n) n−1 ≤

1 A(1) 0

≤ 3Γ(2 − α)τα1 .

(4.34)

Combining (4.33), (4.34) and (2.10) yields n n   j−ϑ 2 α 1−ϑ 2 i−ϑ 2 P(n) |R | ≤ 3Γ(2 − α)τ |R | + max |R | P(n) 1 n− j n− j 2≤i≤n

j=1

≤ Cv

 τ3α 1

α2



2α−4 4 + tnα max ti−1 τi , 2≤i≤n

which completes the proof.  17

j=2

(4.35)

Remark 4.1. From Lemmas 4.4 and 4.5, following similar analysis as [12], we also obtain that n 

j−ϑ 2 P(n) 0 ≤ n− j Υ

j=1

and

n 

Cv α α−6 6 Cv max tα t2α−6 τ6 /τ2α , τ1 + t 1 τ 2 + 2 α 1 − α 2≤i≤n i i−1 i i−1

j−ϑ 2 P(n) |q ≤ C v n− j |R

j=1

 τ3α 1 α2

 2α−4 4 + tnα max ti−1 τi , 2≤i≤n

1 ≤ n ≤ N,

1 ≤ n ≤ N and q = 2, 4.

(4.36)

(4.37)

Theorem 4.2. Assume (1.5) holds for j = 4. Then, if M1 holds, the fully discrete solution U n is convergent with respect to the following H 1 -norm, α/2−3 3 α−3 3 α α−2 2

τ2 + tnα/2 h2 + max tiα/2 ti−1 τi /τi−1 + tnα/2 max ti−1 τi . (4.38) ξn h ≤ CΛ τα/2 1 + t1 2≤i≤n

2≤i≤n

In particular, if M2 also holds, we have ξn h ≤ CΛ (h2 + τmin{γα/2,2} ).

(4.39)

Proof. Substituting vh = (Dτα ξ)n−ϑ in (4.13), we obtain  α n−ϑ 2

(Dτ ξ) 0 + ∇h ξn−ϑ , ∇h (Dτα ξ)n−ϑ ηn−ϑ , ∇h (Dτα ξ)n−ϑ )h − (∇h R n−ϑ , ∇h (Dτα ξ)n−ϑ )h = −((Dτα η)n−ϑ , (Dτα ξ)n−ϑ )h − (∇h   ∂ un−ϑ α n−ϑ n−ϑ α n−ϑ (Dτ ξ) ds. − (Υ , (Dτ ξ) )h + ∂K ∂n K

(4.40)

Similar as Lemma 2.4, we have n−ϑ

1  (n) ∇h ξ , ∇h (Dτα ξ)n−ϑ ≥ A ∇τ (∇h ξi 20 ). 2 i=1 n−i n

(4.41)

Similar as (4.18)-(4.20), we have ((Dτα η)n−ϑ , (Dτα ξ)n−ϑ )h By virtue of (4.10), we have

  2 n−ϑ ≤ CΩ Lh + (1 + CΩ )Υ 0 (Dτα ξ)n−ϑ 0 .

(∇h ηn−ϑ , ∇h (Dτα ξ)n−ϑ )h = 0.

(4.42)

(4.43)

Applying Green’s formula, Lemma 4.2 and inverse inequality, and by virtue of the CauchySchwarz inequality and Young’s inequality, we obtain       ∂R n−ϑ α n−ϑ (∇h R n−ϑ , ∇h (Dτα ξ)n−ϑ )h  =  (Dτ ξ) ds − (ΔR n−ϑ , (Dτα ξ)n−ϑ )h  ∂n ∂K K (4.44) ≤ CΩ h3 |R n−ϑ |4 (Dτα ξ)n−ϑ h + ΔR n−ϑ 0 (Dτα ξ)n−ϑ 0   ≤ CΩ h2 |R n−ϑ |4 + |R n−ϑ |2 (Dτα ξ)n−ϑ 0 . 18

By using Cauchy-Schwarz inequality, we have (Υn−ϑ , (Dτα ξ)n−ϑ )h ≤ Υn−ϑ 0 (Dτα ξ)n−ϑ 0 .

(4.45)

Thanks to Lemma 4.2 and inverse inequality, there yields     ∂ un−ϑ α n−ϑ   (Dτ ξ) ds ≤ CΩ Lh2 (Dτα ξ)n−ϑ 0 .  ∂n ∂K K

(4.46)

Substituting (4.41)-(4.46) into (4.40), we obtain that n 

i 2 α n−ϑ 2 A(n) 0 n−i ∇τ (∇h ξ 0 ) + (Dτ ξ)

i=1



≤ 4CΩ Lh + 2CΩ Υ 2

n−ϑ

0 + 4Υ

n−ϑ

0 + 2CΩ h |R 2

n−ϑ

|4 + 2|R

n−ϑ



(4.47)

|2 (Dτα ξ)n−ϑ 0 .

By using Young’s inequality, we obtain that n 

  i 2 4 n−ϑ 2 4 n−ϑ 2 n−ϑ 2 A(n) ∇ (∇ ξ  ) ≤ C + Υ  + h |R | + |R | h h Λ 0 0 4 2 . n−i τ

(4.48)

i=1

From (4.48), we have n 

P(n) n− j

j=1

j 

i 2 A(n) n−i ∇τ (∇h ξ 0 )

≤ CΛ

i=1

n 

4

j−ϑ 2 P(n) 0 + h4 |R j−ϑ |24 + |R n−ϑ |22 , n− j h + Υ

(4.49)

j=1

Exchanging the summation order arrives at n  j=1

P(n) n− j

j 

i 2 A(n) n−i ∇τ (∇h ξ 0 )

i=1

=

n 

∇τ (∇h ξi 20 )

n 

( j) n 2 0 2 P(n) n− j A j−i = ∇h ξ 0 − ∇h ξ 0 .

(4.50)

j=i

i=1

It is obvious from (2.10) that n 

4 4 α 4 P(n) n− j h ≤ h πA ω1+α (tn ) ≤ C v tn h .

(4.51)

j=1

Substituting (4.50) and (4.51), we get C Cv v 3α α−6 6 α 4 max tα t2α−6 τ6 /τ2α ∇h ξn 20 ≤CΛ 2 (τα1 + h4 τ3α 1 + τ 1 ) + t 1 τ2 + C v tn h + α 1 − α 2≤i≤n i i−1 i i−1  2α−4 4 2α−4 4 + Cv h4 tnα max ti−1 τi + Cv tnα max ti−1 τi . 2≤i≤n

(4.52)

2≤i≤n

Thus, if there exists a positive constant h0 , when h < h0 , it holds that α/2−3 3 α−3 3 α α−2 2

τ2 + tnα/2 h2 + max tiα/2 ti−1 τi /τi−1 + tnα/2 max ti−1 τi . ∇h ξn 0 ≤ CΛ τα/2 1 + t1 2≤i≤n

19

2≤i≤n

(4.53)

In particular, if the mesh assumption M2 also holds, we have τ1 ≤ Cγ τγ and τ2 ≤ Cγ τγ . Thus, we have α/2−3 3 τα/2 τ2 ≤ Cγ τγα/2 , (4.54) 1 + t1 and α−3 3 α −3 3−α −3 3−α tiα/2 ti−1 τi /τi−1 = ti3α/2 ti−1 τi (τi /ti )α /(τi−1 /ti−1 )α ≤ Cγ ti3α/2 ti−1 τi

≤ Cγ ti3α/2−3 τ3−α ≤ Cγ ti3α/2−3 τ3−α−β (τ min{1, ti1−1/γ })β i i ≤ Cγ tiα/2−β/γ (τi /ti )3−α−β τβ ≤ Cγ timax{0,α/2−(3−α)/γ} τβ ,

2 ≤ i ≤ n.

(4.55)

Moreover, there holds 1−1/γ β α−2 2 ti−1 τi ≤ Cγ tiα−2 τ2−β }) ≤ Cγ tiα−β/γ (τi /ti )2−β τβ ≤ Cγ timax{0,α−2/γ} τβ . i (τ min{1, ti

(4.56)

Therefore, it follows from (4.53)-(4.56) that ∇h ξn 0 ≤ CΛ (h2 + τmin{γα/2,2} ).

(4.57)

By using (4.28) and (4.29), the proof is completed.  In what follows, we focus on the superclose analysis of u in H 1 -norm by using a new technique. Theorem 4.3. Under the same conditions of Theorem 4.2, there holds U n − Πh un h = O(h2 + τmin{γα/2,2} ).

(4.58)

Proof. Using the triangle inequality yields U n − Πh un h ≤ U n − Rh un h + Rh un − Πh un h = ξn h + Rh un − Πh un h .

(4.59)

By virtue of the known high accuracy result of bilinear element, and thanks to Lemma 4.3 and (4.39), one arrives at Rh un − Πh un 2h = (∇h (Rh un − Πh un ), ∇h (Rh un − Πh un ))h = (∇h (un − Πh un ), ∇h (Rh un − Πh un ))h ≤ CΩ h2 |u|3 Rh un − Πh un h .

(4.60)

From (4.59) and (4.60), we can obtain (4.58). Thus, the proof is complete.  Next, we show the global superconvergence result of u in H 1 -norm. Define the postprocessing operator Π2h , such that (Π2h u)|K ∈ Q2 (K ) and Π2h u(Ai ) = u(Ai ), i = 1, 2, . . . , 9, where Q2 (K ) is a biquadratic polynomial space, K ∈ T2h consists of the four small elements K j in Th , j = 1, 2, 3, 4, and u(Ai ) is the value of u on the nodes Ai and Ai are nodes of K j , j = 1, 2, 3, 4. From [26], we have Π2h Πun = Π2h un ,

Π2h un − un h ≤ Ch2 |un |3 ,

Π2h vh h ≤ Cvh h ,

∀vh ∈ Vh .

(4.61)

un − Π2h U n h ≤ un − Π2h Πh un h + Π2h Πh un − Π2h U n h ≤ un − Π2h un h + CΠh un − U n h ≤ CΩ h2 |un |3 + C(h2 + τmin{γα/2,2} ) ≤ C(h2 + τmin{γα/2,2} ).

(4.62)

Therefore, we get

We state the result (4.62) in the following theorem. 20

Theorem 4.4. Under the same conditions of Theorem 4.2, there holds un − Π2h U n h = O(h2 + τmin{γα/2,2} ).

(4.63)

Remark 4.2. From the proof process of Theorem 4.2, it is noticed that the suboptimal conclusion of the H 1 -norm bound is resulting from (4.54). But the H 1 -norm error estimate is optimal for tn not close to t = 0. The result is similar as Ref. [27]. Meanwhile, it may be difficult to overcome this deficiency if we use traditional analytical techniques. We can use the local grid refinement technique, not a best method, to ensure τ1 ≤ Cγ τ2γ and τ2 ≤ Cγ τ2γ . Recently, there are two different works to obtain the H 1 -norm error estimates. The first one is Huang and Stynes [22]. We have used this method in Section 3. However, for the nonconforming FEMs, as the result that   ∂ α n−ϑ n−ϑ



ξn−ϑ α n−ϑ ξn−ϑ h , (Dτ ξ) , Δh ξ = (Dτ ξ) ds − (Dτα ∇h ξ)n−ϑ , ∇h h ∂K ∂n K Huang- Stynes’ method still cannot obtain the temporal optimal error estimate as we analyze the H 1 -norm superconvergence. Another one is Ren etc. [12]. This work adopts the time-space split method to derive the sharp H 1 -norm error estimate result. This method has been used in [14]. However, we should estimate the boundedness of |(Dτα Uτ )n−ϑ |2 when we use FEMs, to obtain which we must multiply the temporal error equation by Δ2 en−ϑ τ , integrate over Ω and then use the fourth-order Green n formula, where Uτ is the temporal semi-discrete solution and enτ = un − Uτn . For the Dirichlet ∂en problem, it is inappropriate to use the fourth-order Green formula, as the result that enτ and ∂nτ may not equal to zero on ∂Ω. To sum up, it is still an open problem to obtain the sharp H 1 -norm error estimate for the nonconforming FEMs. But, we do think it is still meaningful to study the sharp L2 -norm error estimate of the nonconforming FEM, which has a great difference with one of the conforming FEM. Moreover, although it is suboptimal in temporal direction for the H 1 -norm error estimate result, we still can obtain the superclose and superconvergence properties of the constructed discrete scheme in the spatial direction. Remark 4.3. The analytical methods and corresponding results are also suitable for the semilinear reaction diffusion equations. 5. Numerical experiments In this section, some numerical results are presented to verify our theoretical analysis. We here only show the error estimates based on the discrete scheme (4.12). Example 1. Consider the problem (1.1)-(1.3) with Ω = [0, 1] × [0, 1] and T = 1. The right-hand term f is chosen such that the system (1.1)-(1.3) has an exact solution u = ω1+α (t) sin πx sin πy. 21

In order to confirm our error estimates and convergence orders on spatial direction, including the L2 -norm and H 1 -norm, we choose the graded time steps N = 64, and uniform rectangular partition with different nodes (M + 1) in each direction is utilized in our computation. The corresponding results with different α are listed in Tables 1-3, which are all in agreement with the theoretical analysis. From Tables 1-3, we can easily find the superconvergence property of the numerical method proposed in this work. Table 4 presents the error estimates uN −U N 1 and corresponding temporal convergence orders with the different fractional parameters α. In order to obtain the desired convergence orders √ as Theorem 4.1, the numerical tests in Table 4 assume the harsh time-space ratio conditions M = N, which in the end take more computing time during the process (see the CPU times in Table 4). The check processes of L2 -norm convergence order don’t need this harsh assumption, as the result of uN − U N 0 = O(τ2 + h2 ). Moreover, the H 1 -norm superconvergence property of the numerical scheme from Theorem 4.4 shows similar approach as L2 -norm. Therefore, it is noticed that we only need the assumption M = N in Tables 5-7, which results in the less CPU times in every computing processes. In the end, we check the convergence and superconvergence properties of the modified quasiWilson FEM under the quadrilateral meshes. In the tests, we assume M = N, α = σ and the graded parameter γ is chosen satisfying γα = 2. The results of L2 -norm error estimate and H 1 norm superconvergence results with different α are shown in Figs. 3 amd 4, respectively. Table 1: The errors and spatial convergence orders with α = 0.2, N = 64, α = 0.8 and γopt = 2/α.

M uN − U N 0 Order uN − U N 1 Order N uN − I2h U N 1 Order

2 1.2772e-01 1.0709e+00 8.3323e-01 -

4 3.1572e-02 2.7359 5.3837e-01 1.3461 2.4938e-01 2.3615

6 1.4010e-02 2.4148 3.5974e-01 1.1982 1.1209e-01 2.3766

8 7.8763e-03 2.2916 2.7005e-01 1.1411 6.3288e-02 2.2745

10 5.0396e-03 2.2252 2.1613e-01 1.1099 4.0574e-02 2.2154

12 3.4991e-03 2.1838 1.8015e-01 1.0900 2.8202e-02 2.1773

6. Conclusion In this paper, we focus on the numerical schemes for the two-dimensional time fractional diffusion equations with non-smooth solution. In the spatial direction, the nonuniform L2-1σ method is adopted to discrete the Caputo-type fractional derivative. In the spatial direction, two kinds of FEMs are proposed and analyzed, including the bilinear element and the modified quasi-Wilson element. Detailed theoretical analysis, including L2 -norm, H 1 -norm superclose and H 1 -norm superconvergence error estimates, is derived by using new analytical technique and new discrete fractional Gr¨onwall-type inequalities. Finally, some numerical results illustrate the correctness of theoretical analysis. 22

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

Figure 2: Quadrilateral meshes.

=0.2 =0.5 =0.8 slope=-2

10 -1

10 -2

L2 -norm error

0

10 -3

10 -4

10 -5

10 -6 10 0

10 1

10 2

M(N)

Figure 3: L2 -norm error estimate for the nonconforming FEM under quadrilateral meshes.

23

1

Table 2: The errors and spatial convergence orders with α = 0.5, N = 64, α = 0.8 and γopt = 2/α.

M uN − U N 0 Order uN − U N 1 Order N uN − I2h U N 1 Order

2 1.2771e-01 1.0709e+00 8.3336e-01 -

4 3.1565e-02 2.7362 5.3837e-01 1.3463 2.4940e-01 2.3617

6 1.4006e-02 2.4149 3.5974e-01 1.1982 1.1210e-01 2.3766

8 7.8736e-03 2.2918 2.7005e-01 1.1411 6.3296e-02 2.2743

10 5.0375e-03 2.2256 2.1613e-01 1.1099 4.0579e-02 2.2154

12 3.4974e-03 2.1843 1.8015e-01 1.0900 2.8207e-02 2.1770

Table 3: The errors and spatial convergence orders with α = 0.8, N = 64, α = 0.8 and γopt = 2/α.

M uN − U N 0 Order uN − U N 1 Order N uN − I2h U N 1 Order

2 1.2794e-01 1.0708e+00 8.3121e-01 -

4 3.1647e-02 2.7346 5.3837e-01 1.3461 2.4916e-01 2.3585

6 1.4044e-02 2.4146 3.5974e-01 1.1982 1.1200e-01 2.3764

8 7.8955e-03 2.2916 2.7005e-01 1.1411 6.3236e-02 2.2745

10 5.0515e-03 2.2256 2.1613e-01 1.1099 4.0541e-02 2.2154

12 3.5071e-03 2.1843 1.8015e-01 1.0900 2.8181e-02 2.1769

Table 4: The errors uN − U N 1 , temporal convergence orders and CPU times with α = 0.8, γopt = 2/α and M =

(M, N) α = 0.2 Order CPU time α = 0.5 Order CPU time α = 0.8 Order CPU time

(2, 4) 5.3838e-01 ≈ 0.24s 5.3839e-01 ≈ 0.23s 5.4036e-01 ≈ 0.28s

(4, 16) 1.3515e-01 1.9941 ≈ 0.54s 1.3520e-01 1.9936 ≈ 0.46s 1.3538e-01 1.9969 ≈ 0.56s

24

(8, 64) 3.3798e-02 1.9996 ≈ 9.07s 3.3801e-02 2.0000 ≈ 8.46s 3.3803e-02 2.0018 ≈ 8.80s

(16, 256) 8.4495e-03 2.0000 ≈ 2534.47 8.4503e-03 2.0000 ≈ 2551.67s 8.4507e-03 2.0000 ≈ 2548.93s

√

N.

N N Table 5: The errors uN −U N 0 and uN − I2h U 1 , temporal convergence orders and CPU times with α = 0.2, α = 0.8, γopt = 2/α and M = N.

(M, N) uN − U N 0 Order N uN − I2h U N 1 Order CPU time

(4, 4) 3.1424e-02 2.4982e-01 ≈ 0.29s

(8, 8) 7.8419e-03 2.0026 6.3384e-02 1.9787 ≈ 0.29s

(16, 16) 1.9605e-03 2.0000 1.5899e-02 1.9952 ≈ 0.84s

(32, 32) 4.9008e-04 2.0001 3.9780e-03 1.9988 ≈ 5.16s

N N Table 6: The errors uN −U N 0 and uN − I2h U 1 , temporal convergence orders and CPU times with α = 0.5, α = 0.8, γopt = 2/α and M = N.

(M, N) uN − U N 0 Order N uN − I2h U N 1 Order CPU time

(4, 4) 3.0949e-02 2.5126e-01 ≈ 0.25s

(8, 8) 7.7818e-03 1.9917 6.3552e-02 1.9832 ≈ 0.34s

(16, 16) 1.9454e-03 2.0000 1.5940e-02 1.9953 ≈ 0.70s

(32, 32) 4.8620e-04 2.0004 3.9886e-03 1.9987 ≈ 4.43s

N N Table 7: The errors uN −U N 0 and uN − I2h U 1 , temporal convergence orders and CPU times with α = 0.8, α = 0.8, γopt = 2/α and M = N.

(M, N) uN − U N 0 Order N uN − I2h U N 1 Order CPU time

(4, 4) 3.0315e-02 2.5323e-01 ≈ 0.22s

(8, 8) 7.7790e-03 1.9624 6.3560e-02 1.9943 ≈ 0.26s

(16, 16) 1.9457e-03 1.9993 1.5939e-02 1.9956 ≈ 0.59s

(32, 32) 4.8644e-04 2.0000 3.9879e-03 1.9989 ≈ 4.59s

References [1] J. P. Bouchaud, A. Georges, Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications, Phys. Rep. 195 (4) (1990) 127–293. [2] R. L. Magin, Fractional calculus in bioengineering, Begell House Redding, 2006. [3] R. Hilfer, Applications of fractional calculus in physics, World Scientific, Singapore, 2000.

25

10 0 =0.2 =0.5 =0.8 slope=-2

H1 -norm superconvergence

10 -1

10 -2

10 -3

10 -4 10 0

10 1

10 2

M(N)

Figure 4: H 1 -norm superconvergence for the nonconforming FEM under quadrilateral meshes. [4] K. Sakamoto, M. Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems, J. Math. Anal. Appl. 382 (1) (2011) 426–447. [5] M. Stynes, Too much regularity may force too much uniqueness, Frac. Calc. Appl. Anal. 19 (6) (2016) 1554– 1562. [6] N. Kopteva, Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions, Math. Comput. 88 (319) (2019) 2135–2155. [7] B. Jin, R. Lazarov, Z. Zhou, An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data, IMA J. Numer. Anal. 33 (1) (2016) 691–698. [8] M. Stynes, E. O’Riordan, J. L. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal. 55 (2) (2017) 1057–1079. [9] C. Huang, M. Stynes, N. An, Optimal L∞ (L2 ) error analysis of a direct discontinuous Galerkin method for a time-fractional reaction-diffusion problem, BIT Numer. Math. 58 (3) (2018) 661–690. [10] H.-l. Liao, W. McLean, J. Zhang, A discrete Gr¨onwalll inequality with applications to numerical schemes for subdiffusion problems, SIAM J. Numer. Anal. 57 (1) (2019) 218–237. [11] P. Lyu, S. Vong, A high-order method with a temporal nonuniform mesh for a time-fractional Benjamin-BonaMahony equation, J. Sci. Comput. 80 (3) (2019) 1607–1628. [12] J. Ren, H.-L. Liao, J. Zhang, Z. Zhang, Sharp H 1 –norm error estimates of two time–stepping schemes for reaction–subdiffusion problems, arXiv preprint arXiv:1811.08059. [13] H.-l. Liao, W. McLean, J. Zhang, A second-order scheme with nonuniform time steps for a linear reactionsudiffusion problem, arXiv preprint arXiv:1803.09873. [14] X. Li, L. Zhang, H.-l. Liao, Sharp H 1 -norm error estimate of a cosine pseudo-spectral scheme for 2D reactionsubdiffusion equations, Numer. Algor.doi:10.1007/s11075-019-00722-w. [15] D. Li, C. Wu, Z. Zhang, Linearized Galerkin FEMs for nonlinear time fractional parabolic problems with nonsmooth solutions in time direction, J. Sci. Comput. 80 (1) (2019) 403–419. [16] M. Li, J. Zhao, C. Huang, S. Chen, Nonconforming virtual element method for the time fractional reactionsubdiffusion equation with non-smooth, J. Sci. Comput. 81 (3) (2019) 1823–1859. [17] C. Huang, M. Stynes, Superconvergence of a finite element method for the multi-term time-fractional diffusion problem, 2019, https://www.researchgate.net/publication/336104219.

26

[18] A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput. Phys. 280 (2015) 424–438. [19] J. H. Bramble, J. E. Pasciak, O. Steinbach, On the stability of the l2 projection in h1 (ω), Math. Comput. 71 (237) (2002) 147–156. [20] D. Shi, H. Zhu, The superconvergence analysis of an anisotropic finite element, J. Syst. Sci. Complex 18 (4) (2005) 478–487. [21] V. Thomee, Galerkin Finite Element Methods for Parabolic Problems, 2nd Edition, Springer Series in Computational Mathematics 25, Springer Berlin Heidelberg, 1997. [22] C. Huang, M. Stynes, Optimal spatial H 1 -norm analysis of a finite element method for a time-fractional diffusion equation, J. Comput. Appl. Math. 367 (2020) 112435. [23] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam,, 1978. [24] D. Shi, S. Chen, A kind of improved Wilson arbitrary quadrilateral elements, Numer. Math. J. Chin. Univ. 16 (2) (1994) 161–167. [25] D. Shi, L. Pei, Nonconforming quadrilateral finite element method for a class of nonlinear Sine-Gordon equations, Appl. Math. Comput. 219 (17) (2013) 9447–9460. [26] Q. Lin, J. Lin, Finite Element Methods: Accuracy and Improvement, Science Press, Beijing,, 2007. [27] C. Huang, M. Stynes, A direct discontinuous Galerkin method for a time-fractional diffusion equation with a Robin boundary condition, Appl. Numer. Math. 135 (2019) 15–29.

27