A high-order stabilizer-free weak Galerkin finite element method on nonuniform time meshes for subdiffusion problems

1.   Introduction

In this work, we consider the following time fractional diffusion equation (TFDE) with zero Dirichlet boundary value:

where $\varOmega\subset \mathbb{R}^2$ be a bounded polygonal domain with its boundary $\partial \varOmega$, $J = (0, T]$ is the time interval with the final time $T$, ${A(x) = (a}_{ij}(x))_{2\times 2}$ is a symmetric positive definite matrix-valued function in $\varOmega$, $g$ and $f$ are given sufficiently smooth and ${ }_0^C D_{T}^\alpha$ denotes the Caputo fractional derivative of order $\alpha$ ^[1],

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]

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 $\alpha$ in the place of the first-order 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,9,10,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\to 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_\sigma$ 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$_\sigma$ 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_\sigma$ formula to approximate the fractional time Caputo derivative in uniform meshes. However, Alikhanov's higher order L2-$1_\sigma$ formula has a rate of convergence $3-\alpha$ 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_\sigma$ 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. $\Vert u \Vert$ represents the $L^2$-norm of a function $u$ in $\Omega$ and $(u, v): = \int_{\Omega} u v\, d x$ denotes the $L^2$ inner product. For $D\subset \varOmega$, the classical Sobolev spaces are denoted by $H^s(\Omega)$ with $s\geq0$ an integer and the corresponding norm $\Vert u\Vert _{s, D}$ and semi-norms $\vert u\vert_s$. Further, the space $L^\infty(J; H^k(\Omega))$ with the norm

is used. When $D = \varOmega$, we do not write $D$ in the subscript.

2.   Semi-discrete scheme

We start with discretization of (1.1) only in space using a SFWG finite element method. Let $\mathcal{T}_h$ be a partition of the domain $\Omega$ consisting of polygons with a set of shape-regular requirements given in ^[29]. Let $\mathcal{E}_h$ be the set of all edges in $\mathcal{T}_h$, and $\mathcal{E}_h^0 = \mathcal{E}\setminus \partial \Omega$ be the set of all interior edges. For any element $K$$\in \mathcal{T}_h$, $h_K$ is the diameter of $K$ and the mesh size $h = \max_{{K}\in \mathcal{T}_h}{h_{{K}}}$.

We first formulate the weak form of (1.1): For any $t\in (0, T]$, we seek a solution $u(x, t)\in H_0^1(\varOmega)$ satisfying

where $\pi_h g$ is the $L^2$ projection of the function $g$ on the weak Galerkin finite element space defined below.

For given integers $k \geq 1$, the weak Galerkin finite element space $S_h$ associated with $\mathcal{T}_h$ is defined by

and its subspace $S_h^0$ is given by

where $\mathbb{P}_k(K)$ is the space of polynomials with degree of at most $k$ on an element $K$ and $\mathbb{P}_{k-1}(e)$ is the space of polynomials with degree of at most $k-1$ on an edge $e\subset \partial {K}$. Observe that $\omega_b$ of $v = \{\omega_0, \omega_b\}$ is allowed to have different value from the trace of $\omega_0$ on the boundary of each element. The configuration of the SFWG-FEM scheme consists of $\left(\mathbb{P}_k({K}), \mathbb{P}_{k-1}(e), \mathbb{P}_{j}({K})^{d}\right)$ with $j = d+k-1$. Using the definition of the classical definition of weak gradient given by (2.4), the configuration of the SFWG-FEM $\left(\mathbb{P}_k({K}), \mathbb{P}_{k-1}(e), \mathbb{P}_{j}({K})^{d}\right)$ delivers only sub-optimal spatial accuracy in both energy norm and the $L^{2}$ norm, presented in Table 1 (^[50]).

The standard definition for a weak gradient $\nabla_{w} u\in\left[\mathbb{P}_{j}({K})\right]^{2}$ of a weak function $u = \left\{u_{0}, u_{b}\right\}$ is a unique polynomial on each ${K} \in \mathcal{T}_{h}$ satisfying ^[29,43,51]

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 $\left(P_{k}({K}), P_{k-1}(e), \left[P_{j}({K})\right]^{2}\right)$. 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.

For each $\omega = \{\omega_0, \omega_b\}\in S_h\cup H^1(\varOmega),$ we define the weak gradient operator $\nabla_{w} \omega \in [\mathbb{P}_{j}({K})]^2$ of the weak function $\omega$ on ${K}$ as a unique polynomial on ${K}$ satisfying the following equation

where $\textbf{n}$ is the unit outward normal vector to $\partial {K}$. In (2.5), we let $\omega_0 = \omega$ and $\omega_b = \omega$ if $\omega\in H^1(\varOmega)$. Here, $\pi_b$ is the $L^2$-projection onto $\mathbb{P}_{k-1}(e)$ for any $e\in \mathcal{E}_h$. Let $\pi_0$ be the $L^2$-projection onto the space $\mathbb{P}_{k}({K})$ for any ${K}\in \mathcal{T}_h$ and define $\pi_h u = \{ \pi_0 u, \pi_b u\}\in S_h$ for any $u\in H^1({K})$. For any ${K}\in \mathcal{T}_h$, we denote by $\Pi_h: [L^2({K})]^2\to [\mathbb{P}_{j}({K})]^2$ the $L^2$-orthogonal projection defined by $\int_{{K}} (\Pi_h {\bf \tau}-{\bf \tau}) {\bf \sigma} \; dx = 0,$ $\forall {\bf \sigma}\in [\mathbb{P}_{j} ({K})]^2$. From now on we take $j = k+1$.

We have the relation between the weak gradient and the projection $\Pi_h$ stated in the following lemma.

Lemma 2.1. For any $v\in H^1({K})$ and ${K}\in \mathcal{T}_h$, one has

Proof. From (2.5) and integration by parts, one has for ${\bf r}\in [\mathbb{P}_{k+1}({K})]^2$

where we used the definitions of projections $\pi_0\; \pi_b$ and $\Pi_h$. The proof is completed. □

For simplicity, we will use the following notations.

For $u = \{u_0, u_b\}\in S_h$ and $\omega = \{\omega_0, \omega_b\}\in 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 $u_h(t) = \{u_{0, h}(t), u_{b, h}(t)\}: (0, T]\to S_h^0$ such that

for all $\omega = \left\{\omega_0, \omega_b\right\} \in S_h^0$.

We define two energy norms $\Vert \cdot\Vert_E$ and $\Vert \cdot \Vert_{1, h}$, respectively: For $\omega = \{\omega_0, \omega_b\}\in S_h^0$,

The following lemma shows the equivalence of the two norms defined above.

Lemma 2.2. For any $\omega = \left\{\omega_0, \omega_b\right\} \in S_h^0$, there are two constants $a_1, a_2 > 0$ such that

Proof. To avoid repetition, the interested reader is refereed to [52, Lemma 3.2] for the proof. □

Lemma 2.3. There is a constant $C > 0$ such that, for any $\omega\in S_h^0$

Proof. The proof follows from from the definition the bilinear form $A_h(\cdot, \cdot)$ 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^\alpha(J)$ denote the classical fractional Sobolev space with the norm $\Vert \cdot \Vert_{H^\alpha(J)}$. We denote the spaces of infinitely differentiable functions and compactly supported infinitely differentiable functions on $J$ by $C^\infty(J)$ and $C^\infty_0(J)$, respectively. Let ${}_0C^\infty(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 ${}_0C^\infty(J)$ in the norm $\Vert \cdot \Vert_{H^\alpha(J)}$ for $\alpha\in (0, 1)$ defines the Sobolev space ${}_0H^\alpha(J)$. The following Sobolev spaces can be found in ^[53] and ^[54]. The set $\overline{A}$ denotes the closure of a set $A$.

Definition 2.1. ^[53] Define $H_r^\alpha(0, T) : = \overline{C_0^{\infty}(0, T)}$ in the norm $\|\cdot\|_{H_r^\alpha(0, T)}$ defined by

where $|v|_{H_r^\alpha(0, T)} = \| { }_{t}^R D_{T}^\alpha v \|_{L^2(0, T)}$.

Definition 2.2. ^[53] Define $H_c^\alpha(0, T): = \overline{C_0^{\infty}(0, T)}$ in the norm $\|\cdot\|_{H_c^\alpha(0, T)}$ defined by

where $|v|_{H_c^\alpha(0, T)} = \Big\vert \big({ }_0^R D_{T}^\alpha v, { }_{t}^R D_{T}^\alpha v\big)_{L^2(0, T)}\Big\vert^{1/2}$.

Lemma 2.4. ^[53] For $\alpha\in (0, 1)$, the spaces $H^\alpha_r(J)$, $H^\alpha_c(J)$, $H^\alpha_0(J)$ are equal and their norms $\Vert \cdot \Vert_{H^\alpha_r(J)}$, $\Vert \cdot \Vert_{H^\alpha_c(J)}$ and $\Vert \cdot \Vert_{H^\alpha_0(J)}$ are equivalent.

Lemma 2.5. ^[55] If $0 < \alpha < 2,$ $\alpha\neq 1$, $u, v\in H^{\alpha/2}(J)$, then we have

Lemma 2.6. ^[54] For $\alpha > 0$ and $u\in H^\alpha(J)$, we have

where $\hat{u}$ is the extension of $u$ by zero outside of $J$.

For a nonnegative real number $\alpha$ and the Sobolev space $Y$ with the norm $\Vert \cdot \Vert_Y$, define the space

with respect to the norm

We also adapt the following notations:

The broken ${\text { Sobolev space }} H^s\left(\mathcal{T}_h\right)$ is denoted by $H^s\left(\mathcal{T}_h\right): = \left\{v \in L^2(\Omega):\left.v\right|_K \in H^s(K), \quad \forall K \in \mathcal{T}_h\right\}$.

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 $\alpha\in (0, 1)$ and $f\in L^{q}(J; L^2(\varOmega))$, the solution $u_h$ of the semi-discrete problem (2.7) satisfies the following stability estimate with $q = \cfrac{2}{1+\alpha}$.

Proof. Choosing $v = u_h = {\{ u_{0, h}, u_{b, h}\}}$ in (2.7), then using (1.8) and the coercivity of the bilinear form $A_h(\cdot, \cdot)$, we get

where we have used that $u(x, 0) = g(x)$ and $A_h(u_h, u_h)\geq C \Vert u_h\Vert^2_E.$ Employing Lemma 2.5 to the first term on LHS of (2.13), knowing that $\Vert u_h \Vert_E^2\geq 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 = \cfrac{2}{1+\alpha}$ and $p = \cfrac{2}{1-\alpha}$.

With the aid of the facts that (see, e.g., ^[53])

and the embedding theorem from ^[56]

one can conclude that

from which the desired result follows after dividing both sides of the above inequality by $\Vert u_{0, h}\Vert_{H^{\alpha/2}(J; L^2(\varOmega))}$. Thus, the proof is completed. □

Lemma 2.7. Let $u$ be the solution of the problem (1.1). Then, for any $\omega = \{\omega_0, \omega_b\}\in S_h^0$, we have

where $Z(u, \omega): = Z_1(u, \omega)+Z_2(u, \omega)+Z_3(u, \omega)$ with

Proof. Using integration by parts, we have for any $\omega = \{\omega_0, \omega_b\}\in S_h^0$

where we have used that $\omega_b\vert_{\partial \varOmega} = 0$ and ${A}\nabla u$ is single value on $\mathcal{E}_h$.

Using the definition of projection $\Pi_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 $\Pi_h$ and $\pi_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)\}: = \{\pi_0 u(t)-u_{0, h}(t), \pi_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 semi-discrete problem (2.7) for $t\in (0, T]$. Then, we have for any $\omega = \{\omega_0, \omega_b\}\in S_h^0$

Proof. Multiplying the first equation in (1.1) by a test function $\omega_0$ of $\omega = \{\omega_0, \omega_b\}\in S_h^0$ yields

From (2.15) in Lemma 2.7, we obtain

where we have used that $\big({ }_0^C D_{T}^\alpha u, \omega_0\big)_{ \mathcal{T}_{h}} = \big({ }_0^C D_{T}^\alpha \pi_0 u, \omega_0\big)_{ \mathcal{T}_{h}}$ which follows from the definition of projection $\pi_0$.

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 $\omega\in H^1(T)$, we have the trace inequality (see, e.g., ^[29])

The following approximation results will be used in the sequel.

Lemma 2.9. ^[29] Let $u$ be the solution of the problem (1.1) and assume that $\mathcal{T}_h$ is shape regular. Then, the $L^2$ projections $\pi_0$ and $\Pi_h$ satisfy

Lemma 2.10. Assume that $w\in H^{k+1}(\varOmega)$. Then, for any $\omega = \{\omega_0, \omega_b\}\in V_h^0$, one has

Proof. Using Cauchy-Schwarz inequality, (1.3), Lemmas 2.1 and Lemma 2.9 we obtain

It follows from the Cauchy-Schwarz inequality, the trace inequality (2.25), (1.3), (2.10) and Lemma 2.9 that

where we have used that

From Cauchy-Schwarz inequality, the trace inequality (2.25), (1.3), (2.10) and Lemma 2.9 that

where $\Pi_{k-1}$ is the $L^ 2$ projection operator onto the space $[\mathbb{P}_{k-1}({K})]^2$ on each ${K} \in \mathcal{T}_h$ and we have used the fact $\Vert \pi_b \omega_0-\omega_0\Vert _{L^2(\partial {K})}\leq \Vert \pi_0 \omega_0-\omega_0\Vert _{L^2(\partial {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 $\Vert \cdot \Vert_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, \frac{\partial u}{\partial t}\in L^2(0, t; H^{k+1}(\mathcal{T}_h))$ for any fixed $t\in J$. Then, we have the error estimate

where $C(g, u, \frac{\partial u}{\partial t})$ depends on the norms of $g$, $u$ and $\frac{\partial u}{\partial t}$.

Proof. For a fixed $t\in J$, choosing $\omega = e_h$ in (2.23) yields

The coercivity property (2.11) of the bilinear form $A_h(\cdot, \cdot)$ and Lemma 2.10 imply that

Hence, we have

Since $e_h(0) = 0$, one has ${ }_0^C D_{T}^\alpha e_{0, h} = { }_0^R D_{T}^\alpha e_{0, h}$, thus integrating (2.33) with respect to $t$ and using Lemma 2.5, we have

Now, taking $v = \frac{\partial e_h}{\partial t}$ in (2.23) gives

or, equivalently

where we have used the fact that ${ }_0^R D_{T}^\alpha e_{0, h} = { }_0^R D_{T}^{\alpha-1}\frac{\partial e_{0, h}}{\partial t} = { }_0^R D_{T}^{-\beta}\frac{\partial e_{0, h}}{\partial t}$ with $\beta: = 1-\alpha$ (see, e.g., ^[1,57]).

Hence, integrating on the time domain, using the adjoint property of the Riemann-Lioville fractional integrals and then Lemma 2.6, we get

where $\widehat{ e_{0, h}}$ is the zero extension of ${ e_{0, h}}$ in $(0, T]$. Using the fact that $\cos(\pi \beta/2)$ is non-negative for $\beta\in (0, 1)$ and Lemma 2.10, we obtain

Appealing (2.34), one has

Note that $u(t) = u(0)+\int_0^t \frac{\partial u}{\partial t}(s) = g+\int_0^t \frac{\partial u}{\partial t}(s)$. Thus, we get

which combined with (2.35) gives

We complete the proof. □

The optimal error rate in the $L^2$-norm can be derived by introducing the elliptic (Ritz) projection $E_h:H^1_0(\varOmega)\to S_h^0$ defined as follows. For any $w\in H_0^1(\varOmega)$ ^[58]

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\in H^1_0(\varOmega)$, where $E_0 w$ is the WG-FEM solution in the inside of elements and $E_b w$ is the trace on $\mathcal{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\in H_0^1(\varOmega)\cap H^{k+1}(\mathcal{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\in H^{k+1}(\varOmega)$ and $u_h$ be the solution of the problem (1.1) and the semi-discrete problem (2.7), respectively. Then for any fixed $t\in (0, T]$, there holds

Proof. We first split the error $e_h(t): = \xi_h(t)+\theta_h(t)$, where

It follows from Lemma 2.11 that

We shall bound $\xi_h$ using the semi-discrete problem (2.7), the definitions of the elliptic projection (2.36) and the projection $\pi_0$ as follows. For any $\phi = \{\phi_0, \phi_b\}\in S_h^0$, we have

Choosing $\phi = \xi_h$ in the above equation, we get

Since $\xi_0(0) = 0$, one has ${ }_0^C D_{T}^\alpha \xi_0 = { }_0^R D_{T}^\alpha \xi_0$ and Lemma 2.5 implies that

Recalling that

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. □

3.   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_\sigma$ 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(\cfrac{m}{M})^r$ for $m = 0, 1, \dots, M$, where $r\geq1$ 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.

One can easily show that

and

For $m = 0, 1, \dots, M-1$ and $\sigma\in[0, 1]$, let $t_{m+\sigma} = t_m+\sigma \tau_{m+1}$. Following the ides of ^[28], the Caputo derivative ${ }_0^C D_{T}^\alpha u$ in (1.1) at $t = t_{m+\sigma}$ is approximated using the $L2-1_\sigma$ formula.

Here, $c_{0, 0} = \tau_1^{-1} d_{0, 0}, \; c_{m, -1} = 0$ and for $m \geq 1$ one has

where

The following technical lemmas will be useful in the later analysis.

Lemma 3.1. [28, Corollary 1] Assume that $0 \leq \sigma \leq 1$. we have the following coercivity property of a function $Z(\cdot, t)$ on the graded temporal mesh $\left\{t_m\right\}_{m = 0}^M$,

Lemma 3.2. [61, Lemma 7] Let $w \in C[0, T] \cap C^{3}(0, T]$. Suppose that $\left|w^{(d)}(t)\right| \leq C(1+t^{\alpha-d})$ for $d = 0, 1, 2, 3$ and $t \in(0, T]$. Then

where

and

Lemma 3.3. [61, Lemma 5] Assume that $1 \geq \sigma \geq 1-\alpha / 2$. Then for any function $\left\{V^{j}\right\}_{j = 0}^{M}$, one has

Lemma 3.4. [61, Lemma 1] Suppose $\sigma = 1-\alpha / 2$. Assume that $\tau_{j+1} \leq C t_{j}$ for $j \geq 2$ and $\tau_{1} / \tau_{2} \leq \rho$, where $\rho$ is any fixed positive constant. Then for any function $v(t) \in C^{3}(0, T]$, one has

where $\psi_{v}^{j+\sigma}$ and $\psi_{v}^{j, s}$ are given by (3.4) and (3.5), respectively.

4.   Fully discrete numerical scheme

We now formulate the fully discrete L2-$1_\sigma$ SFWG-FEM for the problem (2.1) as follows: find $u_h^{m+1} = \{u_{0, h}^{m+1}, u_{b, h}^{m+1} \}\in S_h^0$ such that

where $u_h^{m, \sigma } = \sigma u_h^{m+1 }+(1-\sigma)u_h^{m}$ and $f^{m+\sigma } = f(\cdot, t_{m+\sigma})$ for $m = 0, 1, \dots, M-1.$

We will prove that the $L^2$ stability of the fully discrete L2-$1_\sigma$ SFWG-FEM (4.1). First, we give a Poincare-type inequality in the WG finite element space $S_h^0$.

Lemma 4.1. ^[62] There is a positive constant $C$ independent of the mesh such that

Lemma 4.2. Let $\{u_h^j\}_{j = 0}^M$ be the solution of the fully discrete problem (4.1). Then

Proof. For $m = 0, \dots, M-1$, choosing $v_h = u_h^{m, \sigma }$ in (4.1) yields

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_\sigma$ SFWG-FEM (4.1). For the error analysis, similar to (2.38), we split the error

where $\xi_h^m = u^m-E_hu^m$ and $\theta_h^m = E_h u^m-u_h^m$. The estimation of $\xi_h^m$ follows from Lemma 2.9, and thus we shall estimate $\theta_h^m$ as follows. From (1.1) and (4.1), for $m = 0, \dots, M-1$ we have $\forall v_h = \{v_{0, h}, v_{b, h}\}\in S_h^0,$

where $\Phi^{m+\sigma } = { }_0^C D_{T}^{\alpha} u^{m+\sigma }-\delta _{t_{m+\sigma }}^\alpha u^m$ and $\mathcal{R}^{m+\sigma } = \nabla\cdot ({A} \nabla u^{m+\sigma })-\nabla\cdot ({A} \nabla u^{m, \sigma })$.

The following theorem shows convergence of the L2-$1_\sigma$ SFWG-FEM in the norm $L^\infty(L^2)$.

Theorem 4.1. Assume that $\sigma = (2-\alpha)/2$. Let $u^m$ and $u_h^m$ be the solutions of (1.1), (4.1), respectively. Assume that $u \in L^{\infty}\left(J; H_0^1(\Omega) \cap H^{k+1}(\Omega)\right), { }_0^C D_{T}^\alpha u \in L^{\infty}\left(J; H_0^1(\Omega) \cap H^{k+1}(\Omega)\right)$ and $\left\|\cfrac{\partial^lu(\cdot, t)}{\partial t^l}\right\| \leq C(1+t^{\alpha-l})$ for $l = 0, 1, 2, 3$. Then there exists a constant $C$ such that

Proof. Let $m \in \{0, 1, \dots, M\}$ be a fixed number. Taking $v_h = \theta_h^{m, \sigma }$ in (4.4) yields

From Lemma 3.1 and Cauchy-Schwarz inequality, one has

where $\beta: = \delta _{t_{m+\sigma }}^\alpha\xi_0+ \mathcal{R}^{m+\sigma } + \Phi^{m+\sigma }$. Now using the fact that the stability of the $L^2$ projection in $L^2$-norm ^[63] and the triangle inequality, we have

Thus, Lemma 3.3 and $u_{0, h}^0-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 $\Vert E_0 \Phi^{i+\sigma}\Vert \leq \Vert \Phi^{i+\sigma}\Vert$ and Lemma 2.9 give that

Now, from Lemma 3.4 and Lemma 3.2 we get

We next consider the second term $\max _{0 \leq i \leq M-1}\left\{t_{i+\sigma}^\alpha\left\|\mathcal{R}^{i+\sigma}\right\|\right\}$ in (4.6). When $i = 0$, the assumption $\Vert u(t)\Vert \leq C$ implies that

where we have used that $t_1 = M^{-r}$. Taylor's theorem [61, Lemma 9] and the assumption that $\left\|\cfrac{\partial^2 u(\cdot, t)}{\partial t^2}\right\| \leq C(1+t^{\alpha -2})\leq C t^{\alpha -2}$, for each $i\geq 1$ yield

where we have used (3.1) and (3.2). Therefore,

This result, together with (4.9), gives

From (4.6)–(4.10), one can conclude that

Combining (4.11), Lemmas (2.9) and (4.3) gives the desired result (4.5). The proof is now completed. □

Corollary 4.1. Assume that $r \alpha \geq 2$. Then the error between the numerical solution $u_{h}^{\mathrm{m}}$ computed by $L2-1_\sigma$ SFWG-FEM scheme (4.1) satisfies

Proof. From Theorem 4.1 and $r \geq 2 / \alpha$, one has $\left\| u^m-u_{0, h}^m\right\| \leq C\left(h^{k+1}+M^{-2}\right)$ for each $m$. The result now follows. The proof is now completed. □

The convergence order $O\left(M^{-2}\right)$ in time for the L2-$1_\sigma$ SFWG-FEM proved in Corollary 4.1 is higher-order than the $O\left(M^{-(2-\alpha)}\right)$ temporal rate obtained by the well-known L1-WG-FEM. The $O\left(h^{k+1}\right)$ spatial convergence rate of the numerical methods is optimal in $L^2(\Omega)$.

5.   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\times (0, 1]$. In the following examples, we have first divided the domain $\varOmega$ 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_\sigma$ SFWG-FEM scheme (4.1) on graded meshes. In the numerical solutions, the L2-$1_\sigma$ SFWG-FEM solution $u_h^m = \{ u_{0, h}^m, u_{b, h}^m\}$ is computed by the scheme (4.1). We compute the temporal errors $e_h = \{e_0, e_b\} = \{u^m-u_{0, h}^m, u^m-u_{b, h}^m \}$ in the $L^{\infty}\left(L^2\right)$-norm using the following formula

where we have used $5$-point Gaussian quadrature rules on time mesh interval for each $m$ to approximate $\left\|u\left(x, t_m\right)-u_{0, h}^m\right\|$. The order of convergence (OC) is computed using the following formula

Example 5.1. Let ${A} = 1$ in (1.1). We take $f(x, y, t)$ so that the problem (1.1) has the following solution

where $\theta$ is a fractional number.

If $\theta \in (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 = \frac{1}{2^m}, \quad 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 $\alpha = 0.4$ and $\theta = 2$. The orders of convergence listed in this table show that the proposed method has the optimal convergence order of $\mathcal{O}(h^k)$ in the energy norm and $\mathcal{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 $\mathcal{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 $\theta\in (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$.

Example 5.2. Let $K = 0.1$ in (1.1). The function $f(x, y, t)$ is taken so that the exact solution of the problem (1.1) is

where $\theta$ 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 $\Vert u^m-u_{0, h}^m\Vert_{L^\infty(L^2), M}$ of L2-$1_\sigma$ SFWG-FEM and their rates of convergence in Table 8 and in Table 9 for various values of $\alpha$. 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 $\mathcal{O}(M^{-2})$ is only obtained by the values of the grading constant $r\geq 2/\theta$.

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.

Example 5.3. Let $\varOmega = (0, 1)^2, \; A = 1$, $T = 1$ in (1.1) and $f(x, y, t)$ be given such that the exact solution

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.

6.   Conclusions

In this paper, we studied the SFWG-FEM in space and L2-$1_\sigma$ 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.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The authors would like to express their deep gratitude to the editor and anonymous referees for their valuable comments and suggestions that improve the paper.

Conflict of interest

The authors declare that they have no conflicts of interest.