On the Ulam stability and existence of $L^p$-solutions for fractional differential and integro-differential equations with Caputo-Hadamard derivative

1.   Introduction

We will present a weak Galerkin finite element method for the following parabolic singularly perturbed convection-reaction-diffusion problem:

where $\varOmega = (0, 1)^2$, $0 < \varepsilon\ll 1$, and $T > 0$ is a constant. Assume ${\bf b} = {\bf b}(x, y), \; c = c(x, y)$, and $u^0 = u^0(x, y)$ are sufficiently smooth functions on ${\varOmega}$, and

for some constants $\beta_1, \beta_2$, and $c_{0}.$ The parabolic convection-dominated problem (1.1) has been utilized in a broad range of applied mathematics and engineering including fluid dynamics, electrical engineering, and the transport problem ^[1,2].

In general, the solution of the problem (1.1) will have abrupt changes along the boundary. In other words, the solution exhibits boundary/interior layers near the boundary of $\varOmega$. We are only interested in the boundary layers by excluding the interior layers which can be accomplished by assuming some extra compatible conditions on the data; see, e.g., ^[1,3]. The standard numerical schemes including the finite element method give unsatisfactory numerical results due to the boundary layers. Some nonphysical oscillations in the numerical solution can occur even on adapted meshes, and it is not easy to solve efficiently the resulting discrete system ^[4]. There are many numerical schemes for solving convection-dominated problems efficiently and accurately in the literature. These methods include Galerkin finite element methods ^[5,6,7], weak Galerkin finite element methods (WG-FEMs) ^[8,9,10], the streamline upwind Petrov-Galerkin (SUPG) ^[11,12], and the discontinuous Galerkin (DG) methods ^[13,14,15]. Among these numerical methods, the standard WG-FEM introduced in ^[16] is also an effective and flexible numerical algorithm for solving PDEs. Recently, the WG methods demonstrate robust and stable discretizations for singularly perturbed problems. In fact, while the WG-FEM and the hybridizable discontinuous Galerkin share something in common, the WG-FEM seems more appropriate for solving the time dependent singularly perturbed problems since the inclusion of the convective term in the context of hybridized methods is not straightforward and makes the analysis more subtle. Errors estimates of arbitrary-order methods, including the virtual element method (VEM), are typically limited by the regularity of the exact solution. A distinctive feature of the WG-FEM lies in its use of weak function spaces. Moreover, hybrid high-order (HHO) methods have similar features with WG-FEMs. In fact, the reconstruction operator in the HHO method and the weak gradient operator in WG methods are closely related, and that the main difference between HHO and WG methods lies in the choice of the discrete unknowns and the design of the stabilization operator ^[17]. Notably, in weak Galerkin methods, weak derivatives are used instead of strong derivatives in variational form for underlying PDEs and adding parameter free stabilization terms. Considering the application of the WG method, various PDEs arising from the mathematical modeling of practical problems in science are solved numerically via WG-FEMs using the concept of weak derivatives. There exist many papers on such PDEs including elliptic equations in ^[16,18,19], parabolic equations ^[20,21,22], hyperbolic equations ^[23,24], etc.

However, to the best of the author's knowledge, there is no work regarding the uniform convergence results of the fully-discrete WG-FEM for singularly perturbed parabolic problems on layer-adapted meshes. This paper uses three layer-adapted meshes defined through mesh generating functions, namely, Shishkin-type meshes, Bakhvalov-Shishkin type meshes and Bakhvalov-type meshes given in ^[25]. The error estimates in this work show that one has optimal order of convergence for Bakhvalov-Shishkin type meshes and Bakhvalov-type meshes while almost optimal convergence for Shishkin-type meshes. The main ingredient of the error analysis is the use of the vertices-edges interpolation of Lin ^[26]. The main advantage of this interpolation operator is that we have sharper error bounds compared with the classical interpolation operators. For the sake of simplicity, the Crank–Nicolson method is used for time discretization. This scheme yields optimal order estimates for fully-discrete WG-FEM. As an alternative, one can apply a discontinuous Galerkin method in time and present optimal order estimates for the fully-discrete scheme ^[27].

The rest of the paper is organized as follows. In Section 2, we introduce some notations and recall some definitions. The regularity of the solution is also summarized and three layer-adapted meshes have been introduced in Section 2. Also, we define the weak gradient and weak convection operators, and using these concepts we define our bilinear forms. In Section 3, the semi-discrete WG-FEM and its stability results have been presented. Section 4 introduces a special interpolation operator and analyses interpolation error estimates. Section 5 presents error analysis of the semi-discrete WG-FEM for the problem (1.1) on the layer-adapted meshes. In Section 6, we apply the Crank-Nicolson scheme on uniform time mesh in time to obtain the fully-discrete WG-FEM, and prove uniform error estimates on the layer-adapted meshes. In Section 7, we conduct some numerical examples to validate the robustness of the WG-FEM for the problem (1.1). Summary on the contributions of this work are presented in Section 8.

2.   Preliminaries and weak Galerkin finite element method

Let $S$ be a measurable subset of $\Omega$. We shall use the classical Sobolev spaces $W^{r, q}(S), H^r(S) = W^{r, 2}(S), H_0^r(S), L^q(S) = W^{0, q}(S)$ for negative integers $r > 0$ and $1 \leq q \leq \infty$, and $(\cdot, \cdot)_S$ for the $L^2$ inner product on $S$. The semi-norm and norm on $H^r(S)$ are denoted by $\vert \cdot \vert_{{r}, S}$ and $\|\cdot\|_{{r}, S}$, respectively. If $S = \Omega$, we do not write $S$ in the subscript. Throughout the study, we shall use $C$ as a positive generic constant, which is independent of the mesh parameters $h$ and $\varepsilon$.

2.1.   Decomposition of the solution

This section deals with the introduction of a decomposition of the solution which provides a priori information on the exact solution and its derivatives. Based on this solution decomposition, we construct layer-adapted meshes. As we noted in the introduction, the solution of (1.1) exhibits typically two exponential boundary layers at $x = 1$ and $y = 1$. The following lemma gives some information on the solution decomposition and bounds on the solution of (1.1) and its derivatives.

Lemma 2.1. Let $k$ be positive integer and $l\in (0, 1)$. Assume that the solution $u$ of the problem (1.1) belongs to the space $C^{k+l}(Q_T)$ where $Q_T: = \varOmega\times (0, T]$. Assume further that the solution $u$ can be decomposed into a smooth part $u_R$ and layer components $u_{L_0}, u_{L_1}$, and $u_{L_1}$ with

where the smooth and layer parts satisfy

for any $(x, y)\in \overline{\varOmega}, \; t\in [0, T]$, and positive integers $i, \; j, \; r$ with $i+j+2r\leq k$, and $C$ only depends on ${\bf b}, c$, and $f$ and is independent of $\varepsilon$. Here, $u_R$ is the regular part of $u$, $u_{L_0}$ is an exponential boundary layer along the side $x = 1$ of $\Omega$, $u_{L_1}$ is an exponential boundary layer along the side $y = 1$, while $u_{L_2}$ is an exponential corner layer at $(1, 1)$.

Proof. Under some smoothness conditions and strong imposed compatibility requirements, Shishkin proved this solution decomposition; see, ^[1]. □

2.2.   Layer-adapted meshes

Let $N_x$ and $N_y$ be positive integers. For the sake of simplicity, we assume that $N_x = N_y = N$ is an even integer number. We shall construct the tensor product mesh $\mathcal{T}_N = \{T_{ij}\}_{i{, j = 1, 2\dots, N}}$ in $\overline{\varOmega}$ consisting of elements $T_{ij} = I_i\times K_j$ with the intervals $I_i = (x_{i-1}, x_i)$ and $K_j = (y_{j-1}, y_j)$, where the mesh points are defined by

Since the construction of the meshes in both directions is similar, the mesh construction in $x$-variable is given here.

We define the transition parameter as

where $\sigma\geqslant p+1$ is a positive constant. Here, $p$ is the degree of the polynomials used in the approximation space. The function $\varphi$ obeys the conditions

Assumption 1. Throughout this article, we assume that $\varepsilon \leq C N^{-1}$ such that $\tau = \frac{(p+1) \varepsilon}{\beta_{1}} \varphi(1/2)$, since otherwise the analysis can be carried out using uniform mesh.

Let the mesh points $x_i$ be equally distributed in $[0, x_{N/2}]$ with $N/2$ intervals and distributed $[x_{N/2}, 1]$ with $N/2$ intervals using the mesh generating function defined by

For example, as in ^[25], the Shishkin-type (S-type) meshes can be deduced by $\varphi(1/2) = \ln N$ while Bakhvalov-type meshes (B-type) can be recovered by taking $\varphi(1/2) = \ln (1/\varepsilon)$.

We will use the mesh characterizing function $\psi$ defined by $\psi = e^{-\varphi}$, which is an essential tool in our analysis below.

Following ^[2], we list some famous adaptive meshes including S-type, Bakhvalov-Shishkin meshes (BS-mesh), and B-type in Table 1.

Similarly, we define the transition point in the $y$-direction as

We first split the domain $\varOmega$ into four subdomains as in Figure 1:

Clearly, the mesh is uniform in $\varOmega_r$ with a mesh size of $\mathcal{O}(N^{-1})$, highly anisotropic in $\varOmega_x$ and $\varOmega_y$, while it is very fine in $\varOmega_{xy}$.

Let $h_i^x: = x_i-x_{i-1}, \; {i = 1, \dots, N}\quad h_{{j}}^y: = y_j-y_{j-1}, \; {j = 1, \dots, N}$ be the mesh sizes of the subintervals. For the sake of simplicity, we assume that $\beta_1 = \beta_2 = \beta$. Then, one has $h_i^x = h_{{j}}^y$, and we simply write $h_i, \; i = 1, \dots, N$ for simplicity. The following technical lemmas show the smallness of the boundary layer-functions and the basic properties of the mesh sizes of the layer-adapted meshes.

Lemma 2.2. ^[28] Denote by $\Theta_i = \min \left\{h_i / \varepsilon, 1\right\} e^{-\alpha\left(1-x_i\right) / \sigma \varepsilon}$ for $i = N / 2+1, \ldots, N$. There exists a constant $C > 0$ independent of $\varepsilon$ and $N$ such that

Lemma 2.3. ^[28] For the layer-adapted meshes we considered here, we have

Moreover, for $i = N / 2+1, \cdots, N,$

and, for the B-type mesh,

where $C > 0$ is a constant independent of $\varepsilon$ and $N$.

2.3.   The numerical method

A weak formulation of the problems (1.1) and (1.2) is to look for $u\in H_0^1(\varOmega)$ such that

where the bilinear from ${A}(\cdot, \cdot)$ is defined by

Based on the weak formulation (2.8), we define the WG-FEM on the layer-adapted mesh. Let $p$ be a positive integer. We define a local WG-FE space $V(p, K)$ on each $K\in \mathcal{T}_N$ given by

where $\mathbb{Q}_p(K)$ is the polynomials of degree $p$ on $K$ in both variables, and $\mathbb{P}_p(e)$ denotes the polynomials of degree $p$ on the edge $e$.

Defining the WG finite element space $V_N$ globally on $\mathcal{T}_N$ as

and its subspace

The weak gradient operator $\nabla_w u_N \in [\mathbb{Q}_{p-1}(K)]^2$ can be defined on $K$ as

where ${\bf n}$ represents the outward unit normal $\partial K$, $(w, v)_K$ denotes the inner product on $K$ for functions $w$ and $v$, and $\langle w, v\rangle_{\partial K}$ is the $L^{2}-$inner product on $\partial K$.

The weak convection operator ${\bf b}\cdot \nabla_w u_N \in \mathbb{Q}_{p}(K)$ can be defined on $K$ as

For simplicity, we adapt

For $u_N = \{u_0, u_b\}\in V_N$ and $v_N = \{v_0, v_b\}\in V_N$, the bilinear form ${\cal A}_{w}(\cdot, \cdot)$ is given by

where $s_d(\cdot, \cdot)$ and $s_c(\cdot, \cdot)$ are bilinear forms defined by

with $\partial K^{+} = \{ z \in \partial K: {\bf b}(z) \cdot {\bf n}(z) \geq 0\}$ and ${\rho_K}$ is the penalty term given by

Given a mesh rectangle $K$, its dimensions parallel to the $x$ and $y$-axes are written as $h_{x, K}$ and $h_{y, K}$, respectively.

Lemma 2.4. ^[29] For all $K \in {\cal T}_{N}$ with $h_K = \min \{h_{x, K}, h_{y, K}\}$, there exists a constant $C$ depending only on $p$ such that

We next formulate our semi-discrete WG scheme as follows (Algorithm 1).

3.   The semi-discrete scheme and stability analysis

This section is devoted to establishing the stability results of the WG-FEM defined by (2.15). Define the energy norm $\Vert \cdot\Vert_E$ on the weak function space $V_N$ for $v_N = \{v_0, v_b\}\in V_N$,

Define also an $H^1$ equivalent norm on $V_N$ by

The equivalence of these two norms on $V_N^0$ is given in the next lemma.

Lemma 3.1. For $v_N\in V_N^0$, one has

Proof. For $v_N = \{v_0, v_b\}\in V_N^0$, it follows from the weak gradient operator (2.10) and integration by parts that;

Choosing $w = \nabla_w v_N$ in (3.3) and using the Cauchy-Schwartz inequality and the trace inequality (4), we arrive at

Hence, we get

Summing over all $K\in \mathcal{T}_N$ yields

From the penalty term (2.13), we get

and

As a result, for $v_N\in V_N^0$, we have

Taking $w = \nabla v_0$ in (3.3) and using the Cauchy-Schwartz inequality, we get

where we have again used the trace inequality (4).

Consequently,

Summing over all $K\in \mathcal{T}_N$ yields

Therefore, we have

which implies

From (3.4) and (3.5), we obtain the desired conclusion, which completes the proof. □

We shall show the coercivity of the bilinear form ${\cal A}_{w}(\cdot, \cdot)$ in $\Vert \cdot \Vert_E$ norm on $V_N^0$.

Lemma 3.2. For any $v_N\in V_N^0$, the following inequality holds:

Proof. For $v_N = \{v_0, v_b\}, w_N = \{w_0, w_b\}\in V_N^0$, using the weak convection derivative (2.11) and integration by parts gives

and

where we use the facts that $v_N, w_N\in V_N^0$, and $\big \langle\boldsymbol{b} \cdot \boldsymbol{n}v_b, w_b\big \rangle = 0$ in the last inequality.

Combining (3.7) and (3.8), and taking $v_N = w_N$, we obtain

From (3.9), we have

which yields (3.6) with $C = \min \left\{c_{0}, \frac{1}{2}\right\}$. The proof is completed. □

Therefore, the existence of a unique solution to (2.15) follows from the coercivity property of the bilinear form (3.6). As a result of the two lemmas above, the bilinear form is also coercive in the $\Vert \cdot \Vert _\varepsilon$-norm in the sense that for any $v_N\in V_N^0$, there holds

Lemma 3.3. If $f \in {L}^{2}(\Omega)$ for each $t\in (0, T)$, then there is a constant $C > 0$ independent of $t$ and mesh size $h$ such that the solution $u_{N}(t) = \{u_0(t), u_b(t)\}$ defined in (2.15) satisfies

Proof. Choosing $v = u_{N}(t)$ in (2.15) gives that

Using the Cauchy-Schwarz inequality and the coercivity (3.6) of the bilinear form ${\cal A}_{w}(\cdot, \cdot)$,

Integrating the above inequality with respect to the time variable $t$, we arrive at

Using the Gronwall's inequality gives the desired conclusion. We are done. □

4.   Interpolation estimates

First, we define "vertices-edges' interpolation $\mathcal{P}_K v$ of a function $v$ on an element $K$ as follows. Let $\hat{K}: = (-1, 1)\times (-1, 1)$ be the reference element with the vertices $\hat{a_i}$ and the edges $\hat{e}_i$ for $i = 1, \dots, 4$. For $\hat{v}(\cdot, \cdot) \in C(\hat{K})$, the approximation $\hat{\mathcal{P}}: C(\hat{K})\to \mathbb{Q}_p(\hat{K})$ is given by

The approximation operator $\hat{\mathcal{P}}$ is well-defined ^[30]. Thus, we can define a continuous global interpolation operator $\mathcal{P}_N: C(\overline{\varOmega})\to V_N$ by writing

where the bijective mapping $G_K:\hat{K}\to K$ is given by $G_K(\xi, \eta) = (x, y) = (x_K+h_K^x \xi/2, y_K+h_K^y \eta/2)$.

This interpolation operator has the following stability estimate ^[30]

Since our approximation operator $(\mathcal{P}_N v)|_{K}$ is continuous on $K$, we have $\{ (\mathcal{P}_N v)|_{K}, (\mathcal{P}_N v)|_{e} \} \in V_N$ for $e\in \partial K$. By the trace theorem, we will denote this by $\mathcal{P}_N v$.

Lemma 4.1. ^[31] For any $\phi\in H^1(K)$,

We recall some technical results from ^[30].

Lemma 4.2. For any $K \in \mathcal{T}_N$ and $q \in[1, \infty]$, there exists a constant $C$ such that the vertices-edges-element approximant $\mathcal{P}_{N}\phi$ satisfies

for all $\phi \in W^{p+1, q}(K)$.

Lemma 4.3. Let $s \in [1, p]$. The following estimates hold for any $K \in \mathcal{T}_{N}$:

for all $\psi \in H^{r+1}(K)$.

A careful inspection of the proof of Lemma 4.3 in ^[30] reveals that the following results also hold true.

Lemma 4.4. For $K \in \mathcal{T}_{N}$ and $s \in [1, p+1-\ell]$ with $\ell = 1, 2, \dots, p+1$, there exists a constant $C$ such that the vertices-edges-element approximant $\mathcal{P}_{N}\phi$ satisfies

for all $\phi \in H^{r+1}(K)$.

Lemma 4.5. Let the assumption of Lemma 2.1 hold such that $u = u_R+u_L, \; u_L = u_{L_0}+u_{L_1}+u_{L_2}.$ Then there is a constant $C > 0$ such that the following interpolation error estimates are satisfied:

Proof. The proof of (4.6) follows from Lemma 4.2 and the solution decomposition (2.2) of Lemma 2.1.

Using the decay bound (2.3) of $u_{L_0}$ and the fact that $\varphi(1/2)\geq \ln N$, we have

which shows (4.7) for $u_{L_0}$. Similar arguments can be applied to the layer functions $u_{L_1}$ and $u_{L_2}$. Thus, we complete the proof of (4.7). For (4.8) and (4.9), we will prove for $u_{L_0}$ since the other two parts follow similarly. One can use the decay bound (2.3) of $u_{L_0}$ to obtain

We now shall prove (4.9). Appealing (2.3), one gets

The proof of (4.10) is a little longer. Using an inverse estimate yields

Hence, we shall estimate $\| \mathcal{P}_{N} u_{L}\|_{L^2(\Omega_{r})}$. With the help of the stability estimate (4.5) and the decay bound (2.3) of $u_{L_0}$, we have

If $i < N / 2$, then the sum can be small as a function of $\varepsilon$ but not small if $i = N / 2$. For $i = 1, \ldots, N / 2-1$ and $x \in\left[x_{i-1}, x_i\right]$, we have

and when $i = N / 2$, again using the fact $\varphi(1/2)\geq \ln N$,

Thus,

which proves (4.10). To prove (4.11), we use (4.7) and (4.10) to obtain

On the set $\varOmega_r\cup \varOmega_y$, from the triangle inequality, one obtains

For $I_1$, we have

For $I_2$, using (4.5) and the decay property (2.3) of $u_{L_0}$,

where we have used that $\varphi(1/2)\geq \ln N.$ Applying Lemma 4.2 and the decay property (2.3), we have for any $K\subset\Omega_{x}\cup \Omega_{xy}$

where we have used $\cfrac{h_{x, K}}{\varepsilon}\geq CN^{-1}\geq Ch_{y, K}$ and Lemma 2.2. Similarly, one can show that (4.11) holds for $u_{L_1}$ as well.

For $K\subset \varOmega_r\cup \varOmega_x\cup \varOmega_y$, one can prove as above $\sum_{K\subset \varOmega_r\cup\varOmega_x\cup \varOmega_y}\left\|u_{L_2}- \mathcal{P}_N u_{L_2}\right\|_{L^2(K)}\leq C N^{-(p+1)}$. For $K\subset \varOmega_{xy}$, we obtain

Thus, we complete the proof of (4.11). The proof of the lemma is now completed. □

Lemma 4.6. Let $u\in H^{p+1}(\varOmega)$. Assume that the conditions of Lemma 4.5 hold. Then, we have

for $l = 1, 2$, where $u_L$ denotes $u_{L_0}, \; u_{L_1}$, or $u_{L_2}$.

Proof. The first and second estimates follow from Lemma 4.3, Lemma 4.4, and the fact that $\max\{h_{x, K}, h_{y, k}\}\leq C N^{-1}$. From the triangle inequality and (4.8) and (4.10) of Lemma 4.5, we have

where we have used that $\varepsilon\leq CN^{-1}$. This completes the proof of the third and fifth inequalities for $l = 1$. An inverse inequality and (4.9) and (4.10) of Lemma 4.5 lead to

where again we have used $\varepsilon N\leq C$. This proves the third and fifth inequalities for $l = 2$. Using Lemma 4.4 with $r = p+1-\ell$ for $\ell = 1, 2$ and the decay bound (2.3) of $u_{L_0}$, one can show that for any $K\subset \varOmega_x\cup \varOmega_{xy}$,

Similarly, one can prove that the result holds for $u_{L_1}$, too.

For $K\subset \varOmega_r\cup \varOmega_x\cup \varOmega_y$, we get $\sum_{K\subset \varOmega_r\cup\varOmega_x\cup \varOmega_y}\left\|\cfrac{\partial^\ell (u_{L_2}- \mathcal{P}_N u_{L_2})}{\partial x^\ell}\right\|_{L^2(K)}\leq C \varepsilon^{1/2-\ell}N^{-(p+1)}$. For $K\subset \varOmega_{xy}$, we obtain

which completes the proof of the fourth inequality. Similarly, one can prove that the last inequality holds true. The proof is now finished.

□

5.   Error analysis of the semi-discrete WG-FEM

Unlike the classical numerical methods such as FEM and the SUPG, the proposed WG-FEM does not have Galerkin orthogonality property. This results in some inconsistency errors in the error bounds. We first give a useful error equation in the following lemma.

Lemma 5.1. ^[31] Let $u$ solve the problem (1.1). For $v_N = \{v_0, v_b\}\in V_N^0$,

where

The following error equation $e_N = \{e_0, e_b\}: = \{ \mathcal{P}_N u-u_0, \mathcal{P}_N u-u_b\}$ will be needed in the error analysis.

Lemma 5.2. Let $u$ and $u_N$ be the solutions of (1.1) and (2.15), respectively. For $v_N = \{v_0, v_b\}\in V_N^0$, one has

where $E(u, v_N): = E_d(u, v_N)+E_c(u, v_N)+E_r(u, v_N)$, which are defined by (5.4), (5.5), and (5.6), respectively.

Proof. Multiplying (1.1) by a test function $v_N = \{v_0, v_b\}\in V_N^0$, we arrive at

With the help of (5.1), (5.2), and (5.3), the above equation becomes

Since $\mathcal{P}_N u$ is continuous in $\varOmega$, we get

Therefore, we have

Subtracting (2.15) from (5.8) gives (5.7), which completes completed. □

Lemma 5.3. Let $\mathcal{P}_N u$ be the vertex-edge-cell interpolation of the solution $u$ of the problem (1.1). Then, there holds

Proof. The solution decomposition (2.1) implies that

Using Lemma 4.2 with $q = \infty$ and Lemma 2.1,

Next, we examine the layer parts one by one. Let $(x, y)\in K\subset \varOmega_r\cup\varOmega_y$. From the $L^\infty$ stability property (4.5) of the interpolation operator, one has

Let $(x, y)\in K \subset \Omega_{x} \cup \Omega_{xy}$. The stability property (4.5) and Lemma 4.2 with $q = \infty$ yield

Similarly, we can derive the estimates on the other layer components $u_{L_1}$ and $u_{L_2}$. Combining the above estimates gives the desired conclusion.

Thus, we complete the proof. □

We recall the following trace inequality. For any $u \in H^1(K)$, one has

Lemma 5.4. Let $u\in H^{p+1}(\varOmega)$ and $\rho_K$ be given by (2.13). Assume that the conditions of Lemma 4.5 hold. Then, one has

Proof. For the sake of simplicity, we use the following notations. Let $\zeta_R: = u_R- \mathcal{P}_N u_R$ and $\zeta_L: = u_L- \mathcal{P}_N u_L$ represent the interpolation errors of the regular and layer components of the solution. Hence, by the triangle inequality, we have

With the help of the trace inequality (5.9), we have

Now, appealing the definition (2.13) of stabilization parameter $\rho_K$ and Lemma 4.6 gives

where we have used that $\varepsilon N < 1$.

Using once again the trace inequality (5.9), we have

Now, appealing the definition (2.13) of stabilization parameter $\rho_K$ and Lemma 4.6 again reveals that

Plugging (5.12) and (5.11) into (5.10) yields

Consequently, we have

which completes the proof.

□

Now, we shall prove the error bounds for the consistency errors.

Lemma 5.5. (A priori bounds) Assume that ${\cal T}_{N}$ is the tensor product mesh as defined in Section 2. Then, for $u\in H^{k+1}(\Omega)$ and $v_N\in V_{N}^0,$ we have

Proof. It follows from the Cauchy-Schwarz and Holder inequalities that

Now, it then follows from Lemma 4.6 that

Next, we consider the term ${\cal S}_{2}$. From the Cauchy-Schwarz inequality and Lemma 5.4, we have

Combining (5.16) and (5.17), we get

From (5.5) and (5.6) and using $\langle u-I_N u, {\bf b}\cdot {\bf n}v_b\rangle = 0$, we arrive at

Now, the Hölder inequality and Lemma 5.3 lead us to write

The Cauchy Schwartz and inverse inequalities give

Appealing the Cauchy Schwartz inequality on $\varOmega\setminus\varOmega_r$, we have

Using the error bounds (5.20) and (5.21) in (5.19), we obtain

where we have used the fact that $N^{-1} (\ln N)^{1/2} \max \left|\psi^{\prime}\right| = N^{-1}(\ln N)^{1/2} \left|\psi^{\prime}(0)\right| < C$.

Since $u$ and $\mathcal{P}_N u$ are continuous, we conclude that $\Vert u- \mathcal{P}_N u\Vert_{L^\infty(e)}\leq \Vert u- \mathcal{P}_N u\Vert_{L^\infty(K)}$ for any $e\subset \partial K \in \mathcal{T}_N$. Then, the Hölder inequality and Lemma 5.3 imply that

where we have used that $\varepsilon \ln N < C$ and $(N^{-1} \max \left|\psi^{\prime}\right|)^{2(p+1)}\ln N \leq N^{-(2p+1)} (\max \left|\psi^{\prime}\right|)^{2(p+1)}$.

Using the Cauchy Schwartz inequality and (4.6) and (4.11) of Lemma 4.5, we obtain

The proof is completed. □

By letting $v_N = e_{N}$ in (5.7), we obtain

It then follows from the estimates (5.13) and (5.14), together with Young's inequality and 3.10, that

As a result,

Then, by integrating from $0$ to $t$, we have

This result is collected in the following theorem.

Theorem 5.1. (Semi-discrete estimate) Let $u\in H^{k+1}(\Omega)$ be the solution of (1.1)-(1.2) and $u_{N}\in V_{N}^0$ be the solution of (2.15). Then, we have

6.   Fully-discrete WG finite element scheme

In this section, we shall use the Crank-Nicolson scheme on uniform time mesh in time to derive the fully discrete approximation of the problem (1.1) and (1.2). For a given partition $0 = t_0 < t_1 < \cdots < t_M = T$ of the time interval $J = [0, T]$ for some positive integer $M$ and step length $\tau = \frac{T}{M}$, we define

where the sequence $\{\omega^n\}_{n = 0}^M \subset L^2({\cal D})$. For simplicity, we denote $\xi(\cdot, t_n)$ by $\xi^n$ for a function $\xi:[0, T]\rightarrow L^2(\Omega).$ We now state our fully discrete weak Galerkin finite element approximation. Find $U^n_{N} = \{U_0^n, U_b^n\}\in V_N^0$ such that

with $U^0_{N} = { \mathcal{P}}_Nu(0)$ and $f^{n+\frac{1}{2}} = \frac{1}{2}(f^{n+1}+f^{n})$.

The following lemma shows that the Crank-Nicolson scheme is unconditionally stable in the $L^2$ norm.

Lemma 6.1. Let $f\in C(0, T; L^2(\Omega))$. Then, we have the following stability estimate for the fully-discrete scheme (6.1):

Proof. Choosing $v = U_N^{n+1/2}$ in (6.1), and using the Cauchy-Schwarz inequality, we get

where we have used that ${\cal A}_{w}(U_N^{n+\frac{1}{2}}, U_N^{n+\frac{1}{2}})\geq C \Vert U_N^{n+\frac{1}{2}}\Vert_\varepsilon ^{2}$. Using the fact that $a^2-b^2 = (a-b)(a+b)$ and the coercivity of the bilinear form, we have

Let $1\leq j \leq M$ be an integer. We sum the above inequality from $n = 1$ to $n = j$:

Recalling that $M\tau = T$, the result follows. The proof is now completed.

□

Next, we shall present the convergence analysis. To begin, we prove the error estimate of the discretization error ${\cal P}_N u(t_n)-U^n_{N}$. To this end, we need to derive an error equation involving the error $e^n_{N}: = U^n_{N}-{\cal P}_{N}u^n$.

We formulate the error equation for $e_{N}^n$ in the following lemma.

Lemma 6.2. For $v_{N} = \{v_0, v_b\} \in V_{N}^0,$ we have

where $\xi^n = \frac{1}{\tau}(\mathcal{P}_N u^{n+1}- \mathcal{P}_N u^n)-\partial_{t}u^{n+\frac{1}{2}}$; and $E(u, v_N) = E_d(u, v_N)+E_c(u, v_N)+E_r(u, v_N).$

Proof. From (1.1), one obtains the following equation:

On each element $T\in {\cal T}_{N},$ for $v_N = \{v_0, v_b\}\in V_{N}^0,$ we test equation (6.4) against $v_0$ to arrive at

Using a similar argument in deriving (5.8), one can show that

where we have used that ${\cal S}_{d}({\cal P}_{N}u^{n+\frac{1}{2}}, v_N) = 0\; \; \text{ and }\; \; {\cal S}_{c}({\cal P}_{N}u^{n+\frac{1}{2}}, v_N) = 0$ since ${\cal P}_{N}u$ is continuous in $\Omega$. Thus, we get

Subtracting (6.1) from (6.6) gives the conclusion. We complete the proof. □

Lemma 6.3. Let $u\in H^{k+1}(\Omega)$. Assume that $u$ and $U^n_{N}$ are the solutions (1.1), (1.2), and (6.1), respectively. One has for $n = 1, 2, \ldots, M$,

Proof. Choosing $v_N = e_{N}^{n+\frac{1}{2}}$ in (6.3) and by the coercivity property (3.6), we find

or, equivalently,

We can express the term $\xi^n = (\partial_\tau \mathcal{P}_N u^n-\partial_\tau u^n)+ (\partial_\tau u^n-\partial_t u^{n+\frac{1}{2}}) = : T_1+T_2$. We write

and

From (6.9) and (6.10), we obtain

Hence, with the aid of the Cauchy-Schwarz and the Poincare inequality, $\mathcal{W}_1$ in (6.8) can be estimated as follows.

where we have used the Young's inequities in the second inequality, and the estimate (6.11) and Lemma 4.5 in the second estimates of the righthand side. Applying Lemma 5.5 and Young's inequality, we obtain the estimate of the term $I_2$ in the righthand side of (6.8) as follows:

Combining (6.8)–(6.13) yields

Let $1 \leq j \leq M$. Using the fact that $e_N^0 = 0$, we sum the above expression from $n = 0$ to $n = j-1$ to obtain

We complete the proof. □

Theorem 6.1. Let $u\in H^{k+1}(\Omega)$. Assume that $u$ and $U^n_{N}$ are the solutions (1.1), (1.2), and (6.1), respectively. One has for $n = 1, 2, \ldots, M$,

Proof. Choosing $v = \partial_{\tau}e_N^{n}$ in (6.3) and by coercivity (3.6), we find

or, equivalently,

Thus, we have

Because $e_0^0 = 0$, we sum up the above term from $n = 0$ to $n = m-1$ for any fixed $m$ to get

From (6.11), we have

Observe that

Similar to (6.13), one has

It follows from Lemma 5.5, the Cauchy-Schwarz inequality, and Young's inequality that

From (6.16), (6.17), and (6.18) together with $\tau M = T$, we have

Combining (6.14), (6.15), and (6.19) yields that

Finally, using (6.7), we obtain

which completes the proof.

□

7.   Numerical Experiments

This section presents various numerical examples for the fully-discrete Crank-Nicolson WG finite element method. We used MATLAB R2020A in our the calculations. We also used the 5-point Gauss-Legendre quadrature rule for evaluating of all integrals. All the calculations were calculated using MATLAB R2016a. The systems of linear equations resulting from the discrete problems were solved by lower-upper (LU) decomposition.

We apply the fully-discrete WG-FEM on the adaptive meshes shown in Table 1. We choose $\sigma = p+1$ and calculate the energy-norm $\Vert e_N^n\Vert_E$ and the $L^2$-norm error $\Vert e_0^n\Vert$, where $e_N = \{e_0^n, e_b^n\} = \{u-U_0^{n}, u-U_b^n\}$ is the error using $N$ intervals in each direction. The order of convergence (OC) is computed by the formula

The numerical errors and the order of convergences in space are also tested. In order for the space error to dominate the errors, we take $\tau = N^{-2}$ for $N$ element in each direction. We list the errors in the energy norm and $L^2$-norm and the order of convergence in Tables 2 and 3, respectively. These numerical results show that the order of convergence is of order $p$ and of order $p+1$ in the energy and $L^2$ norms, respectively, which support the stated error estimates in Theorem 6.1.

Example 7.1. Let ${\bf b} = (1, 1)$ and $T = 1$ in the problem (1.1). We choose $f$ and $u^0$ such that the exact solution is

where $\kappa(z) = 1-e^{-(1-z)/\varepsilon}$.

In Figure 2, we plot the numerical solutions of the WG-FEM using the $\mathbb{P}_1$ element on the three layer-adapted meshes given in Figure 1 for $\varepsilon = 10^{-5}$ and $N = 32$.

We next present the temporal convergence rate for Example 7.1. In order for the temporal error to dominate the error, we take $N = 256$ and $\varepsilon = 10^{-5}$, and use the $\mathbb{P}_3$ element. We report the results in the $L^2$-norm and the energy norm in Tables 4 and 5, respectively. We see that the order of convergence in time is of order $\mathcal{O}(\tau^2)$, which verifies the theoretical estimate claimed in Theorem 6.1.

Lastly, we test the robustness of the WG-FEM method with respect to the small parameter $\varepsilon$ for Example 7.1. We take $N = 256$ and use the $\mathbb{P}_1$ element for the values of $\varepsilon = 10^{-r}, \; r = 5, 6, \dots, 10$. The results are reported in Table 6. These results show that the WG-FEM is robust with respect to the perturbation parameter $\varepsilon$.

The order of convergence via loglog plot in the energy norm and $L^2$ norm are plotted in Figures 3 and 4, respectively, for Example 7.2. We observe that the order of convergence of order $p$ and of order $p+1$ in the energy and $L^2$ norms, respectively, which support the stated error estimates in Theorem 6.1 as in Example 7.1. To test the temporal error, we choose $N = 256$ and $\varepsilon = 10^{-5}$, and use the $\mathbb{P}_1$ element. We present the results in the $L^2$-norm and the energy norm in Tables 7 and 8, respectively. We see that the order of convergence in time is of order $\mathcal{O}(\tau^2)$ as claimed in Theorem 6.1.

Example 7.2. Let ${\bf b} = (1+x, 2-y)$, $c = (1+x^2+y^2)$, and $T = 1$ in the problem (1.1). We take $f$ and $u^0$ such that the exact solution is

We also test the WG-FEM for Example 7.2 for the robustness against $\varepsilon$. The results are presented in Table 9 for $N = 256$ and the $\mathbb{P}_1$ element for the values of $\varepsilon = 10^{-r}, \; r = 5, 6, \dots, 10$. Again, one sees that the WG-FEM is the parameter-uniform method.

8.   Conclusions

In this paper, we present the Crack-Nicolson- WG-FEM applied to the singularly perturbed parabolic convection-dominated problems in 2D. We use the Crack-Nicolson scheme in time on uniform mesh and the WG-FEM in space on three layer-adapted meshes: Shishkin, Bakhvalov-Shishkin, and Bakhvalov meshes. We prove (almost) uniform error estimates of order $p$ in the energy norm and second order estimate in time. With the use of a special interpolation operator, the error analysis of the semi-discrete WG-FEM and the fully discrete WG-FEM have been carried out. Various numerical examples are conducted to validate the convergence rate of the proposed method.

Use of AI tools declaration

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

Conflict of interest

The authors declare there is no conflicts of interest.