A high-order stabilizer-free weak Galerkin ﬁnite element method on nonuniform time meshes for subdi ﬀ usion problems

: We present a stabilizer-free weak Galerkin ﬁnite element method (SFWG-FEM) with polynomial reduction on a quasi-uniform mesh in space and Alikhanov’s higher order L2-1 σ scheme for discretization of the Caputo fractional derivative in time on suitable graded meshes for solving time-fractional subdi ﬀ usion equations. Typical solutions of such problems have a singularity at the starting point since the integer-order temporal derivatives of the solution blow up at the initial point. Optimal error bounds in H 1 norm and L 2 norm are proven for the semi-discrete numerical scheme. Furthermore, we have obtained the values of user-chosen mesh grading constant r , which gives the optimal convergence rate in time for the fully discrete scheme. The optimal rate of convergence of order O ( h k + 1 + M − 2 ) in the L ∞ ( L 2 )-norm has been established. We give several numerical examples to conﬁrm the theory presented in this work.


Introduction
In this work, we consider the following time fractional diffusion equation (TFDE) with zero Dirichlet boundary value: where Ω ⊂ R 2 be a bounded polygonal domain with its boundary ∂Ω, J = (0, T ] is the time interval with the final time T , A(x) = (a i j (x)) 2×2 is a symmetric positive definite matrix-valued function in Ω, g and f are given sufficiently smooth and C 0 D α T denotes the Caputo fractional derivative of order α [1], ∂s ds, α ∈ (0, 1). (1.2) In particular, we assume that the function A(x) satisfies where k 1 and k 2 are positive constants.
We introduce the left-hand sided and right-hand sided Riemann-Lioville fractional integrals and the left-hand sided and right-hand sided Riemann-Lioville fractional derivatives, respectively, defined by [1] 0 D −α t u(x, t) = 1 Γ(α) t 0 (t − s) α−1 u(x, s)ds, α ≥ 0, (1.4) (1.7) We will use the following identity from [1] in our analysis In recent decades, many physical phenomenon or processes in science and engineering have been modeled using fractional partial differential equations since they describe better memory effect and hereditary properties (see, [2][3][4][5] and references therein).The TFDE considered in this paper is obtained from the classical diffusion problem by taking a fractional derivative of order α in the place of the firstorder time derivative.TFDEs are obtained from a fractional Fick law describing transport equations with long memory [6].The TFDE is formulated by a linear integro-differential equation.Several analytical and numerical techniques have been considered for solving TFDE in the literature.The Fourier transform method, the Laplace transform method, the Mellin transform method and the Green function technique are some examples of the analytical approach (see [1,7] and references therein).In general, the analytical solutions of fractional differential or partial differential equations are not available; thus one must construct a numerical method for solving this equation.Many papers have investigated robust numerical approximations to the time fractional diffusion equations or sub-diffusion problems [8][9][10][11][12][13][14][15][16][17][18].The TFDEs with the Caputo derivatives have been solved by the widely used L1 method in [8-11, 19, 20].The compact difference methods [12,14,16,21] have been considered to improve the order of time approximation, and spectral methods [9,17,18] have been proposed for improving the accuracy in the space discretization.High order numerical methods [22][23][24][25][26][27] have been presented for solving fractional partial differential equations.However, typical solutions of (1.1) have an initial layer when t → 0 and the first order time derivative of the solution blows up at t = 0. Since the solution is singular at t = 0, one should construct a reliable and efficient numerical scheme in time [1].Furthermore, the forcing function f in our model problem may blow up at the starting point t = 0 even though u(x, t) is a continuous function in time.Thus, one needs to pay attention to the case f (x, 0) when proposing a numerical method for solving the problem (1.1).As a remedy, we use the Alikhanov's L2 − 1 σ method on graded meshes to approximate the time derivative in (1.1) for higher order convergence.
Our aim of this paper is to combine the L2-1 σ method [28] due to Alikhanov for time discretization of the Caputo fractional derivative on graded temporal meshes with the SFWG-FEM in space on quasi-uniform meshes to develop a robust scheme to solve the TFDE (1.1).The weak Galerkin finite element method first appeared in [29] for solving the second order elliptic problems.The novelty of this method lies on the definitions of weak function space and weak derivatives such as weak gradient operator on this weak function space.This renders the method more applicable in defining finite element space and mesh constructions; for example, one uses completely discontinuous function spaces in the approximation and polygonal meshes.In the last decade, this method has been applied to solve a variety of differential equations, e.g., parabolic problems [30,31], second order elliptic problems [32][33][34][35][36][37][38][39], Stokes equations [31], Maxwell equations [40] and fractional differential equations [41,42].The stabilized-free weak Galerkin method is introduced in [43] for solving second-order elliptic partial differential equations.This new method does not have a stabilization component, which is essential in discontinuous approximations to enforce the jump across the element boundaries.Also, the stabilization term is needed to ensure the coercivity of the WG-FEM.The absence of stabilization terms in numerical methods makes the formulation much simpler and provides much flexibility in implementation.The main idea of the SFWG-FEM is to use a higher degree of polynomials in computing weak differential operators for the strong connectivity of weak functions on element boundaries.In [44] and [45], the authors applied the SFWG-FEM to second-order elliptic equations and proved that the method has supercloseness convergence in an energy norm and L 2 -norm on triangular meshes.For more discussions of the method, we refer the interested readers to [43,46,47].Further, in [48] the authors derived the semi-discrete SFWG-FEM formulation in space and presented the stability results and error estimates.Then, they established the fully discrete SFWG-FEM by discretization of the time using the L2-1 σ formula to approximate the fractional time Caputo derivative in uniform meshes.However, Alikhanov's higher order L2-1 σ formula has a rate of convergence 3 − α on a uniform mesh if the solution is smooth enough, while this rate of convergence reduces to first order for the case where the typical solutions with a weak singularity at the starting point.As pointed out in [49, Theorem 2.1], the smooth solutions exist only under acceptably restrictive data.These observations and facts lead to using the L2-1 σ formula on graded meshes in time and SFWG-FEM in space in this paper.We show the order of convergence result that gives second-order convergence in time and an optimal rate of convergence in space.
The remaining parts of the paper are organized as follows.Section 2 introduces a semi-discrete SFWG-FEM.Stability and error analysis of the semi-discrete method are presented.The fully-discrete SFWG-FEM and its stability and error analyses are given in Section 3. We present several numerical examples to support the theory in Section 4. Finally, we close the paper with some concluding remarks in Section 5.
Throughout the paper, we use the following notations.u represents the L 2 -norm of a function u in Ω and (u, v) := Ω uv dx denotes the L 2 inner product.For D ⊂ Ω, the classical Sobolev spaces are denoted by H s (Ω) with s ≥ 0 an integer and the corresponding norm u s,D and semi-norms |u| s .Further, the space L ∞ (J; H k (Ω)) with the norm is used.When D = Ω, we do not write D in the subscript.

Semi-discrete scheme
We start with discretization of (1.1)only in space using a SFWG finite element method.Let T h be a partition of the domain Ω consisting of polygons with a set of shape-regular requirements given in [29].Let E h be the set of all edges in T h , and E 0 h = E \ ∂Ω be the set of all interior edges.For any element K∈ T h , h K is the diameter of K and the mesh size h = max K∈T h h K .
We first formulate the weak form of (1.1):For any t ∈ (0, T ], we seek a solution u(x, t) ∈ H 1 0 (Ω) satisfying where π h g is the L 2 projection of the function g on the weak Galerkin finite element space defined below.
For given integers k ≥ 1, the weak Galerkin finite element space S h associated with T h is defined by and its subspace S 0 h is given by where P k (K) is the space of polynomials with degree of at most k on an element K and P k−1 (e) is the space of polynomials with degree of at most k − 1 on an edge e ⊂ ∂K.Observe that ω b of v = {ω 0 , ω b } is allowed to have different value from the trace of ω 0 on the boundary of each element.
The configuration of the SFWG-FEM scheme consists of P k (K), P k−1 (e), P j (K) Using the definition of the classical definition of weak gradient given by (2.4), the configuration of the SFWG-FEM P k (K), P k−1 (e), P j (K) d delivers only sub-optimal spatial accuracy in both energy norm and the L 2 norm, presented in Table 1 ( [50]).
Table 1.Weak gradient calculated by (2.4), The standard definition for a weak gradient ∇ w u ∈ P j (K) 2 of a weak function u = {u 0 , u b } is a unique polynomial on each K ∈ T h satisfying [29,43,51] (2.4) Following the ideas in [52], we shall use a different definition of weak gradient given by (2.5) for the SFWG-FEM using the configuration P k (K), P k−1 (e), P j (K)

2
. With this new definition, our proposed SFWG-FEM has optimal order spatial convergence rates in the energy and L 2 norms, presented in Table 2.
Table 2. Weak gradient calculated by (2.5), For each ω = {ω 0 , ω b } ∈ S h ∪ H 1 (Ω), we define the weak gradient operator ∇ w ω ∈ [P j (K)] 2 of the weak function ω on K as a unique polynomial on K satisfying the following equation where n is the unit outward normal vector to ∂K.In (2.5), we let ω 0 = ω and Here, π b is the L 2 -projection onto P k−1 (e) for any e ∈ E h .Let π 0 be the L 2 -projection onto the space P k (K) for any K ∈ T h and define π h u = {π 0 u, π b u} ∈ S h for any u ∈ H 1 (K).For any K ∈ T h , we denote by Π h : [L 2 (K)] 2 → [P j (K)] 2 the L 2 -orthogonal projection defined by K (Π h τ − τ)σ dx = 0, ∀σ ∈ [P j (K)] 2 .From now on we take j = k + 1.
We have the relation between the weak gradient and the projection Π h stated in the following lemma.
Lemma 2.1.For any v ∈ H 1 (K) and K ∈ T h , one has Proof.From (2.5) and integration by parts, one has for r where we used the definitions of projections π 0 π b and Π h .The proof is completed.
For simplicity, we will use the following notations.
(z, w) zw dx, and z, w For u = {u 0 , u b } ∈ S h and ω = {ω 0 , ω b } ∈ S h , we define the bilinear form as Based on the weak form (2.1), the semi-discrete SFWG-FEM for the problem (1.1) is to find a numerical solution for all ω = {ω 0 , ω b } ∈ S 0 h .We define two energy norms ) (2.9) The following lemma shows the equivalence of the two norms defined above.

Lemma 2.3.
There is a constant C > 0 such that, for any ω ∈ S 0 Proof.The proof follows from from the definition the bilinear form A h (•, •) given by (2.6) and (1.3).
In order to analyze convergence and stability properties of the time-fractional partial differential equation, we introduce some function spaces.
Let H α (J) denote the classical fractional Sobolev space with the norm • H α (J) .We denote the spaces of infinitely differentiable functions and compactly supported infinitely differentiable functions on J by C ∞ (J) and C ∞ 0 (J), respectively.Let 0 C ∞ (J) be the the space of infinitely differentiable functions on the interior of J with compact support in J. Now, the closure of the space 0 C ∞ (J) in the norm • H α (J) for α ∈ (0, 1) defines the Sobolev space 0 H α (J).The following Sobolev spaces can be found in [53] and [54].The set A denotes the closure of a set A.
For a nonnegative real number α and the Sobolev space Y with the norm • Y , define the space .
We also adapt the following notations: The following theorem ensures the well-posedness of the semi-discrete SFWG-FEM (2.7) and the stability estimate for t > 0.
Theorem 2.1.For any α ∈ (0, 1) and f ∈ L q (J; L 2 (Ω)), the solution u h of the semi-discrete problem (2.7) satisfies the following stability estimate with q = 2 (2.12) ), then using (1.8) and the coercivity of the bilinear form where we have used that u(x, 0) = g(x) and A h (u h , u h ) ≥ C u h 2 E .Employing Lemma 2.5 to the first term on LHS of (2.13), knowing that u h 2 E ≥ 0, and using Cauchy-Schwarz inequality on the right side of (2.13), we obtain Integrating the above expression with respect to t, we get By Lemma 2.4, we obtain Employing Hölder inequality on the left side of the above inequality gives where q = 2 1 + α and p = 2 1 − α .With the aid of the facts that (see, e.g., [53]) and the embedding theorem from [56] from which the desired result follows after dividing both sides of the above inequality by u 0,h H α/2 (J;L 2 (Ω)) .Thus, the proof is completed.
Lemma 2.7.Let u be the solution of the problem (1.1).Then, for any ω = {ω 0 , ω b } ∈ S 0 h , we have where Z(u, ω) ) Proof.Using integration by parts, we have for any where we have used that ω b | ∂Ω = 0 and A∇u is single value on E h .Using the definition of projection Π h , one has On the other hand, it follows from the definition of weak gradient (2.5) that From (2.20), (2.21) and using the definitions of Π h and π b , we have Plugging the above equation (2.22) into (2.19)yields which gives (2.15).Thus, we complete the proof.
Next, we present an error equation for the discretization error e h (t) = {e 0,h (t), e b,h (t)} := {π 0 u(t) − u 0,h (t), π b u(t) − u b,h (t)}.This error equation plays an important role in our error analysis.Lemma 2.8.Assume that u(t) is the solution of the problem (1.1) and u h (t) is the solution of the semidiscrete problem (2.7) for t ∈ (0, T ].Then, we have for any Proof.Multiplying the first equation in (1.1) by a test function From (2.15) in Lemma 2.7, we obtain where we have used that Subtracting the first equation in (2.7) from (2.24), we get which is the desired result (2.23).Thus, we complete the proof.
For ω ∈ H 1 (T ), we have the trace inequality (see, e.g., [29]) (2.25) The following approximation results will be used in the sequel.
The optimal error rate in the L 2 -norm can be derived by introducing the elliptic (Ritz) projection In fact, E h w is the WG-FEM solution of the corresponding elliptic equation that has exact solution w.
Note that E h w = {E 0 w, E b w} for any w ∈ H 1 0 (Ω), where E 0 w is the WG-FEM solution in the inside of elements and E b w is the trace on E h .
The following error bounds for the elliptic projection E h follow from [52, Theorem 4.4 ] and [52, Theorem 4.5].
Lemma 2.11.[52] Suppose that w ∈ H 1 0 (Ω) ∩ H k+1 (T h ).Then, the elliptic projection E h defined by (2.36) has the following error estimates We now state and prove the error estimate of a semi-discrete scheme in L 2 -norm.
Theorem 2.3.Let u ∈ H k+1 (Ω) and u h be the solution of the problem (1.1) and the semi-discrete problem (2.7), respectively.Then for any fixed t ∈ (0, T ], there holds (2.37) Proof.We first split the error e h (t) := ξ h (t) + θ h (t), where (2.39) It follows from Lemma 2.11 that We shall bound ξ h using the semi-discrete problem (2.7), the definitions of the elliptic projection (2.36) and the projection π 0 as follows.For any φ = {φ 0 , φ b } ∈ S 0 h , we have Choosing φ = ξ h in the above equation, we get we get Applying the arithmetic mean inequality yields , where we have used the equivalent of the norms (2.10) in the last inequality.From (2.40), we obtain which proves the desired result.

Temporal discretization
In order to present semi-discrete numerical scheme, we discretize (1.1) in time direction.We investigate time discretization using the well known L2 − 1 σ formula on graded meshes to deal with the singularity of solution at t = 0.
Let M > 0 be a natural number.We define temporal graded meshes by setting t m = T ( m M ) r for m = 0, 1, . . ., M, where r ≥ 1 is a user-chosen grading constant.Similar temporal graded meshes have been investigated in the literature, e.g, [59,60].The graded constant r will influence the convergence rate; thus we take it such a way that the rate is optimal as presented below.

Fully discrete numerical scheme
We now formulate the fully discrete L2-1 σ SFWG-FEM for the problem (2.1) as follows: find where u m,σ h We will prove that the L 2 stability of the fully discrete L2-1 σ SFWG-FEM (4.1).First, we give a Poincare-type inequality in the WG finite element space S 0 h .Lemma 4.1.[62] There is a positive constant C independent of the mesh such that Lemma 4.2.Let {u j h } M j=0 be the solution of the fully discrete problem (4.1).Then Invoking Lemma 3.1, we obtain Using Cauchy-Schwarz inequality, the Young's inequality and Lemma 4.1, we obtain Using this inequality in (4.2) gives Thus, from Lemma 3.3, one has which completes the proof.
We now prove an error estimate for the fully discrete L2-1 σ SFWG-FEM (4.1).For the error analysis, similar to (2.38), we split the error where The estimation of ξ m h follows from Lemma 2.9, and thus we shall estimate θ m h as follows.From (1.1) and (4.1), for m = 0, . . ., M − 1 we have where ).The following theorem shows convergence of the L2-1 σ SFWG-FEM in the norm L ∞ (L 2 ).Theorem 4.1.Assume that σ = (2 − α)/2.Let u m and u m h be the solutions of ( Assume that u ∈ L ∞ J; Then there exists a constant C such that Proof.Let m ∈ {0, 1, . . ., M} be a fixed number.Taking From Lemma 3.1 and Cauchy-Schwarz inequality, one has where β := δ α t m+σ ξ 0 + R m+σ + Φ m+σ .Now using the fact that the stability of the L 2 projection in L 2norm [63] and the triangle inequality, we have Thus, Lemma 3.3 and u 0 0,h − E 0 u 0 = 0 imply that Therefore, one has We shall estimate each term in (4.6) individually as follows.We first note that E 0 Φ i+σ ≤ Φ i+σ and Lemma 2.9 give that We next consider the second term max 0≤i≤M−1 t α i+σ R i+σ in (4.6).When i = 0, the assumption u(t) ≤ C implies that where we have used that t 1 = M −r .Taylor's theorem [61,Lemma 9] and the assumption that ∂ 2 u(•, t) where we have used (3.1) and (3.2).Therefore, This result, together with (4.9), gives The convergence order O M −2 in time for the L2-1 σ SFWG-FEM proved in Corollary 4.1 is higherorder than the O M −(2−α) temporal rate obtained by the well-known L1-WG-FEM.The O h k+1 spatial convergence rate of the numerical methods is optimal in L 2 (Ω).

Numerical results
In this section, we present some numerical experiments to show that our theoretical results are valid.For easy implementation and to avoid much algebra, we take T = 1 and Q T = [0, 1] 2 × (0, 1].In the following examples, we have first divided the domain Ω into squares and then further divided these squares into triangles by connecting the diagonal lines (as illustrated in Figure 1).We use the space of piecewise linear polynomials on the triangles and constant polynomials on the edges.We aim to verify the spacial and temporal accuracy of the L2-1 σ SFWG-FEM scheme (4.1) on graded meshes.In the numerical solutions, the L2-1 σ SFWG-FEM solution u m h = {u m 0,h , u m b,h } is computed by the scheme (4.1).We compute the temporal errors e h = {e 0 , where we have used 5-point Gaussian quadrature rules on time mesh interval for each m to approximate u (x, t m ) − u m 0,h .The order of convergence (OC) is computed using the following formula . where θ is a fractional number.
If θ ∈ (0, 1), then the solution u will have a weak singularity at t = 0.That is, the integer-order derivatives of u with respect to the time will blow up at the starting point.
We first fix M = 2000 which is small enough to avoid the time discretization errors and take the mesh size h = 1 2 m , m = 1, 2, 3, 4, 5, 6.We present the spatial errors at t = 1 in the energy and L 2 norms and their orders of convergence is shown in Table 3 when α = 0.4 and θ = 2.The orders of convergence listed in this table show that the proposed method has the optimal convergence order of O(h k ) in the energy norm and O(h k+1 ) in the L 2 norm, as stated by Corollary 4.1.Then, we fix the mesh size h in space so that the errors in time dominate the errors in space.The results in Table 4, Table 5 and Table 6 show that the rate of convergence is of order O(M −2 ), which is in perfect agreement with Corollary 4.1.In these tables, when r = r opt , we observe that the computed OC is slightly bigger than the expected ones.The reason for this is that the optimum value of the grading constant yields a better approximation in L 2 norm for the problem having a weak singularity at the initial point.Observe that if the regularity parameter θ ∈ (0, 1), then one cannot achieve the optimal convergence rate using the uniform mesh or non-optimal grading constant r due to the singularity of the solution at t = 0.
where θ is a regularity parameter.
In this example, we take M = 2000 so that the spatial error dominates the error in time.Table 7 lists the computed errors at t = 1 in the energy and L 2 norms and their orders of convergence.These rates of convergence displayed are in good agreement with the theory predicted by Corollary 4.1.These rates show that the optimal orders of convergence are obtained.
Next, we fix the spatial step size h = 1/300 to ensure that the temporal errors dominate the errors in space.We compute u m − u m 0,h L ∞ (L 2 ),M of L2-1 σ SFWG-FEM and their rates of convergence in Table 8 and in Table 9 for various values of α.These rates of convergence displayed are in good agreement with the theory predicted by Corollary 4.1.Further, Table 8 displays that the optimal rate of convergence O(M −2 ) is only obtained by the values of the grading constant r ≥ 2/θ.In order to show the advantage of the SFWG method in dealing with the triangular meshes with hanging nodes, we report the measured errors in the energy and L 2 norms in Table 10.The results show the flexibility of the SFWG-FEM in treating hanging nodes in meshes, while the standard finite element method cannot be easily applied without any hp refinement, which makes the computations formidably expensive.

Conclusions
In this paper, we studied the SFWG-FEM in space and L2-1 σ method in time for the fractional diffusion problems on graded meshes in time.We derived optimal error estimates of semi-discrete in the L 2 and H 1 norms and fully discrete numerical schemes in the L 2 norm.Because of the singularity at the initial time, graded meshes in time were used, and the optimal values of the grading constant gave the second order convergence in time.Various examples were carried out to verify the theory presented in this work.We will investigate the method of this paper to fractional diffusion problems with time-dependent diffusion coefficient in future work.

Figure 1 .
Figure 1.Uniform triangulation of the domain.

Figure 2 .
Figure 2. Triangulation of the domain with hanging nodes.
For a fixed t ∈ J, choosing ω = e h in (2.23) yields ∂K) and above estimate(2.30).Combining the estimates (2.28), (2.29) and (2.31) finishes the proof.We now study the estimate for the error e h (t) of the numerical scheme (2.7) in the • E -norm.Theorem 2.2.Let u be the solution of the problem (1.1) and u h be the solution of the semi-discrete problem (2.7).Suppose that u, ∂u ∂t ∈ L 2 (0, t; H k+1 (T h )) for any fixed t ∈ J.Then, we have the C 0 D α T e 0,h , e 0,h + A h (e h , e h ) = Z(u, e h ).The coercivity property (2.11) of the bilinear form A h (•, •) and Lemma 2.10 imply that

Table 3 .
Spatial error histories and the rates of convergence at t = 1 of the SFWG-FEM using P 2 element for Example 5.1 on triangular meshes.

Table 4 .
Time error histories and the rates of convergence of the SFWG-FEM for Example 5.1 on triangular meshes.

Table 5 .
Time error histories and the rate convergence of the SFWG-FEM for Example 5.1 on triangular meshes.

Table 6 .
Time error histories and the rate convergence of the SFWG-FEM for Example 5.1 on triangular meshes.

Table 8 .
Time error histories and the rate convergence of the SFWG-FEM for Example 5.2 on triangular meshes.

Table 9 .
Time error histories and the rate convergence of the SFWG-FEM for Example 5.2 on triangular meshes.In Example 5.3 we use triangular meshes with hanging nodes to illustrate the advantage of the SFWG method in treating irregular meshes shown in Figure 2.

Table 10 .
Spatial error histories and the rates of convergence at t = 1 of the SFWG-FEM using P 2 element for Example 5.3 on triangular meshes with hanging nodes.