Error estimations of a weak Galerkin finite element method for a linear system of $\ell\geq 2$ coupled singularly perturbed reaction-diffusion equations in the energy and balanced norms

1.   Introduction

We consider the following system of $\ell\geq 2$ coupled singularly perturbed reaction-diffusion boundary value problems. Find $\textbf{u} \in (C^2(0, 1)\cap C[0, 1])^\ell$ such that

where $\mathcal{E} = \mbox{diag}(\varepsilon_1^2, \dots, \varepsilon_\ell^2)$ with small perturbation parameters $0 < \varepsilon_i < < 1, \; i = 1, 2\dots, \ell$, the vector-valued function $\textbf{g} = (g_1, g_2, \dots, g_\ell)^T$ and the reaction coefficients matrix $\textbf{A} = (a_{kl})_{k, l = 1}^\ell$ are twice continuously differentiable on $[0, 1]$ and the constants $\bf{u_0}$ and $\bf{u_1}$ are given. The exact solution of (1.1) is the vector $\textbf{u} = (u_1, u_2, \dots, u_\ell)^T$. We assume that $\textbf{A}:[0, 1]\to \mathbb{R}^{(\ell, \ell)}$ and the vector valued function $\textbf{g}:[0, 1]\to \mathbb{R}^{\ell}$ are independent of the perturbation parameters $\mathcal{E}$ and the reaction matrix $\textbf{A}$ is strongly diagonally dominant with

Then the condition (1.2) implies that $\textbf{A}$ is an $M$-matrix and its inverse is positive definite and bounded in the maximum norm (see e.g., ^[1]). Under these assumptions, the problem (1.1) has a unique solution $\textbf{u} = (u_1, \dots, u_\ell)^T\in (C^2(0, 1)\cap C[0, 1])^\ell$. In this paper, without loss of generality we will assume the most general case

which always can be done, if necessary, by renumbering of the equations in the system.

It has been well-known that standard numerical methods including finite difference (FD) methods and finite element methods (FEM) are inefficient and inaccurate when applied to singularly perturbed problems (SPPs) on uniform meshes. The solutions of SPPs have boundary or/and interior layers which are very thin regions in which the solution or its derivative change abruptly. The width of the layers depends on the perturbation parameter and these layers are not resolved unless a large number of mesh points are used which is computationally expensive. As a remedy, fitted mesh methods based on layer-adapted meshes have been proposed and studied during recent years for solving boundary layer problems. The construction of these meshes require a priori knowledge on the bounds for the solution and its derivative. These meshes are finer in the part of boundary layers and coarser in the outside of region of the boundary layers. The well-known layer-adapted meshes are piecewise uniform Shishkin meshes ^[2] and Bakhvalov-type meshes ^[3]. We refer the readers to the books ^[1,2,4] and references therein for more details.

Unlike the non-coupled SPPs, the boundary layer behaviour of the solution to each equation in the system can be dramatically different and complicated. Each solution in the coupled system may have a sublayer corresponding to each of the perturbation parameter in the domain if the perturbation parameter in each equation has a different magnitude. This renders the construction of numerical methods very subtle. Shiskin ^[5] considered coupled system of two reaction diffusion equations on an infinite strip and he proved that the finite difference method is a robust method and has the rate of convergence $\mathcal{O}(N^{-1/4})$ on piecewise uniform meshes when the perturbation parameters are small and different each other. Later, it is shown that the method has higher order convergence in ^[6,7,8,9] using piecewise uniform Shishkin mesh. The finite element method has been developed and analyzed for SPPs in the papers ^[7,10,11]. Recently, the numerical solution of system of reaction diffusion problems have been presented in ^[12,13,14] and reference therein.

Although there has been increasing interest in the numerical solution of coupled system of two singularly perturbed differential equations, few articles discuss the numerical solution of coupled system of more than two singularly perturbed differential equations. Kellogg et al. ^[15] considered a system of {$\ell\geq 2$} reaction-diffusion equations in two dimension with the same perturbation parameter for each equation in the system. A system of {$\ell\geq 2$} reaction-diffusion equations each of which has a different perturbation parameter has been studied and analyzed in one dimension in ^[9]. In general, the reaction coefficient matrix $\textbf{A}$ is assumed to be diagonally dominant with positive diagonal and nonpositive off-diagonal elements in most of the papers. However, this condition is weakened in ^[9]. Usually, error analyses in FEMs or DG methods have been analyzed in the energy norms derived from corresponding variational formulations. Unfortunately, these norms are too weak to capture the boundary layers of SPPs of reaction-diffusion type, see, e.g., ^{[14,16,17,18]}. Up to now, the only work on the balanced error estimates of finite element method for a system of $\ell \geq 2$ coupled singularly perturbed reaction-diffusion two-point boundary value problems is the paper of Lin and Stynes ^[14]. They have proved that the classical FEM using quadratic $C^1$ splines is order of $\mathcal{O}(N^{-1} \ln N)$ in the balanced norm provided that each perturbation parameter is equal to the same small number. {A new FEM is presented for SPPs of reaction-diffusion type in the weighted and balanced norm in ^[19]. The convergence analysis of the classical FEM in a balanced norm on Bakhvalov-type rectangular meshes has been studied in ^[20].} The analysis is much more involved and complicated when each parameter in the system is different. The error analysis in the balanced norm, to the best of the authors' knowledge, is not studied in the literature when the perturbation parameters are different. In this paper, to fill this gap, we derive the error estimates of a weak Galerkin method for a system of $\ell \geq 2$ coupled SPPs of reaction-diffusion type in the energy and balanced norms when each equation in the system has a different parameter. The classical FEM on the balanced norm using $C^0$-elements is open ^[21].

Wang and Ye ^[22] first introduced the weak Galerkin finite element method (WG-FEM) and presented for the second order elliptic equation. The key in the WG finite element scheme is to use the weak functions and weak derivatives on the completely discontinuous piecewise polynomials spaces. Since then, many papers have been devoted to WG finite element methods including the implementation results in ^[23], parabolic problems in ^[24], the Maxwell equations ^[25], the Stokes equations ^[26], the Helmholtz equations with high wave numbers in ^[27] and the multi-term time fractional diffusion equations in ^[28]. In ^[29], a discrete gradient and divergence operators have been introduced for convection-dominated problems. A uniformly convergent weak Galerkin finite element method on Shishkin mesh for convection-diffusion problem in one dimension has been presented in ^[30]. Uniform convergence of the WG-FEM on Shishkin mesh for SPPs of convection-dominated type has been studied in 2D in ^[31] and in 1D in ^[32] and singularly perturbed reaction-convection-diffusion problems with two parameters has been analyzed in ^[33]. Uniform convergence of a weak Galerkin method on Bakhvalov-type mesh for singularly perturbed convection-diffusion problem has been analyzed in ^[34,35] and nonlinear singularly perturbed reaction-diffusion problems in ^[36]. Supercloseness in an energy norm of a WG-FEM on a Bakhvalov-type mesh for a singularly perturbed two-point boundary value problem has been demonstrated in ^[37] and superconvergence results in ^[38]. The WG-FEM for two coupled system of SPPs of reaction-diffusion type has been presented in the energy norm in ^[39]. We wish to study a robust WG-FEM for the coupled systems of SPPs of reaction-diffusion type. Thus, the main aim of this paper is to construct a uniformly convergent WG-FEM for the problem (1.1).

The paper is organized as follows. In Section 2, we present and study a decomposition of the exact solution and a uniform Shishkin mesh. We introduce the WG-FEM in Section 3. Stability properties of the proposed method have been demonstrated in Section 4. Error analysis in the energy and balanced norm is presented in Sections 5 and 6, respectively. In Section 7, the numerical results are conducted to confirm the theory in the previous sections. Finally. conclusion is given in Section 8.

In this work, by $C$ we mean a generic constant independent of $N$ and the perturbation parameters $\varepsilon_i, \; i = 1, \dots, \ell$ which may not be the same at each occurrence. Constants with subscript such as $C_c$ are fixed numbers and also do not depend on $\varepsilon_i, \; i = 1, \dots, \ell$, and the mesh parameter $N$.

2.   Preliminaries

In this section, we first give a decomposition of the analytical solution of the linear system (1.1). Then we will derive the bounds for the solution and its derivatives. Next, a piecewise-uniform Shishkin mesh is constructed. Sobolev spaces with the related norms and some basic notations are introduced at the end of this section.

2.1.   Properties of the solution

The solution of the system (1.1) can be decomposed as $\textbf{u} = \textbf{R}+ \textbf{L}$ where $\textbf{R}$ is the regular solution part and $\textbf{L}$ is the layer parts. In light of (1.2), there is a constant $\rho\in (0, 1)$ such that

Define $\alpha = \alpha(\rho)$ by

For the future reference, we set

and define $\ell \times \ell$ matrix $\Gamma = (c_{ij})$ by

We assume that the matrix $\Gamma$ is inverse monotone, that is, $\Gamma^{-1}$ exists and

We first provide the stability of the solution of (1.1) from ^[9].

Lemma 2.1. Assume that $\boldsymbol{u}$ is the solution of (1.1) and the reaction coefficient matrix $\boldsymbol{A}$ has strictly positive diagonal elements $a_{ii} > 0$ for $i = 1, 2, \dots, \ell$. Let the matrix $\Gamma$ be inverse monotone. Then the solution to each equation in the system has the following bounds

Proof. We refer the reader to ^[9] for the detailed proof. □

The following theorem states that the coercivity of $\bf{A}$ and inverse monotone property (2.3) of the matrix $\Gamma$ are related.

Theorem 2.2. ^[40] Assume that the reaction coefficient matrix $\boldsymbol{A}$ has strictly positive diagonal elements $a_{ii} > 0$ for $i = 1, 2, \dots, \ell$ and the matrix $\Gamma$ is inverse monotone. Then, there is a constant diagonal matrix $\bf{D}$ with positive elements and a positive constant $\beta$ such that

Remark 2.1.

(1) If the matrix $\boldsymbol{A}$ has the property (1.2), then the matrix $\Gamma$ is a strongly diagonally dominant $L_0$ matrix which implies that the matix $\Gamma$ is inverse monotone.

(2) If $\boldsymbol{A}$ and $\boldsymbol{g}$ are twice continuously differentiable, then the above stability result guarantees the existence of a unique solution $\boldsymbol{u}\in C^4[0, 1]^\ell$.

(3) The reaction matrix is assumed to be strongly diagonally dominant with positive diagonal elements and nonpositive off-diagonal elements in most of existence papers on coupled system of SPPs with the exception ^[9,41]. This assumption implies that the operator $\mathcal{L}$ is inverse monotone and satisfies the maximum principle which is a useful tool in finite difference method. In this paper, the assumptions on $\boldsymbol{A}$ are weakened and we consider problems in a more general setting.

(4) Since the form of system (1.1) and the matrix $\Gamma$ do not change when a constant positive diagonal matrix is applied on the left, Theorem 2.2 implies that we can assume, without loss of generality, the reaction matrix $\bf{A}$ is coercive if it has positive diagonal elements. That means that there exists $\eta > 0$ such that

We have to consider the solution decomposition consisting of smooth and layer components because of the boundary layers. Thus, we will use the following decomposition of $\textbf{u}$ in the forthcoming analysis.

where $\textbf{R}$ is the smooth part, $\textbf{L}_L$ and $\textbf{L}_R$ are the boundary layer parts, and satisfy the following boundary value problems, respectively

Here, the existence of the inverse matrix $\textbf{A} ^{-1}$ is guaranteed by the condition (1.2).

Theorem 2.3. Assume that $\boldsymbol{A}$ and $\boldsymbol{g}$ are twice continuously differentiable. Then the solution $\boldsymbol{u}$ of the system (1.1) can be decomposed as $\boldsymbol{u} = \boldsymbol{R}+ \boldsymbol{L}_L+ \boldsymbol{L}_R$, where $\boldsymbol{R}$ and $\boldsymbol{L} = \boldsymbol{L}_L+ \boldsymbol{L}_R$ satisfy

Proof. A detailed proof can be found in ^[42]. □

2.2.   Uniform Shishkin mesh

Let $N$ be an integer divisible by $2(\ell+1)$. We define the transition points

where $\sigma$ is a user-chosen constant with $\sigma = \mathcal{O}(1)$. In general, this parameter is chosen as $\sigma \geq k+1$ where $k$ is the order of polynomials in the approximation space. Then we divide each of the intervals $\Omega_s: = [\lambda_s, \lambda_{s+1}]$ and $_s \Omega: = [1-\lambda_{s+1}, 1-\lambda_s]$, $s = 0, \dots, \ell$ into $\cfrac{N}{2(\ell+1)}$ subintervals of equal mesh size

An example of a piecewise-uniform Shishkin mesh with $N = 32$ elements for a system of $\ell = 3$ reaction-diffusion equations is shown in Figure 1.

We next define the nodes recursively as

We denote the mesh and a partition of the domain $\varOmega$ by $I_n = [x_{n-1}, x_n], \quad n = 1, \dots, N$ and $\mathcal{T}_N = \{I_n:n = 1, \dots, N\}$, respectively. For $I_n\in \mathcal{T}_N$, the outward unit normal $\textbf{n}_{I_n}$ on $I_n$ is defined as $\textbf{n}_{I_n}(x_n) = 1$ and $\textbf{n}_{I_n}(x_{n-1}) = -1$; for simplicity, we use $\textbf{n}$ instead of $\textbf{n}_{I_n}$.

We use the following basic notations in the sequel. By $L^2(\Omega)$, we denote the space of square integrable functions on $\Omega$ with the norm $\Vert u\Vert^2_{L^2(\Omega)} = \int_{\Omega }u^2(x)\, d x$ and sometimes, we will use the abbreviation $\Vert \cdot \Vert = \Vert u\Vert^2_{L^2(\Omega)}$. The standard Sobolev space is denoted by $H^k(\Omega)$ with the norm $\Vert\cdot\Vert_{k, \Omega}$ and semi-norm $|\cdot|_{k, \Omega}$ given as

We define the norm for a vector-valued function $\textbf{u}$ as

For each interval $I_n$, the broken Sobolev space is defined by

and the corresponding norm and semi-norm

For the future reference we use the following notations

3.   The WG-FEM method

This section is devoted to introduce novel concepts such as weak functions and weak derivatives from which we define our method for the problem (1.1). For the rest of the paper, we denote by $\mathbb{P}_k(I_n)$ the set of polynomials defined on $I_n$ with degree at most $k$. The space of weak functions $\mathcal{W}(I_n)$ on $I_n$ is defined by

Here, a weak function $u = \{u_0, u_b\}$ has two components and the first component $u_0$ represents the value of $u$ in $(x_{n-1}, x_n)$ and $u_b$ is interpreted as the value of $u$ on $\partial I_n = \{x_{n-1}, x_n\}$. From now on, we assume that $k = 2$ unless otherwise mentioned.

Let $S_{N}(I_n)$ be a local weak Galerkin (WG) finite element space given by

where $\mathbb{P}_0(\partial I_n)$ stands for constant polynomials on $\partial I_n$. We remark that the results can be extended to $\mathbb{P}_k$ elements when $k > 2$. However, in this case {some } additional compatibility conditions of the data will be required in order to have (2.8).

Next, we define a global WG finite element space $S_N$ that comprises of weak functions $u = \{u_0, u_b\}$ such that $u_0|_{I_n}\in \mathbb{P}_k(I_n)$ and $u_b|_{x_n}$ is the constant for $n = 0, \dots, N$.

Let $S_N^0$ be the subspace of $S_N$ with zero boundary conditions, that is,

The weak derivative $d_{w, I_n} u \in \mathbb{P}_{{k-1}}(I_n)$ of a function $u \in S_N(I_n)$ is defined to be the solution of the following equation

where

and

The discrete weak derivative $d_{w} u$ of the weak function $u = \{u_0, u_b\}$ on the finite element space $S_{N}$ is defined by

Our WG-FEM scheme for the system of singularly perturbed reaction-diffusion problems (1.1) is given as follows.

 

Here, the bilinear and the linear forms are defined by, for any $u_i^N = \{u_{i0}, u_{ib}\}$,

where $\varrho_n \geq 0, n = 1, \dots N$ is the penalization parameter associated with the node $x_n$ defined as follows:

Choosing a penalty parameter in stabilized numerical methods is an important issue in uniform convergence estimates. Usually, the penalization parameter will depend on the perturbation parameters. For example, $\varrho_n = \varepsilon_\ell^2 h_n^{-1}$ is taken in the WG finite element schemes ^[22,26,29]. However, a uniform convergence rate can not be attained for a penalization constant depending on the perturbation parameters.

4.   Stability of the method

In the following analysis, we will recall the following multiplicative trace inequality and the inverse inequality.

We introduce the $\mathcal{E}$-weighted energy norm $\Vert|\cdot \Vert|$ in $[S_N^0]^\ell$ as follows: for $\textbf{v} = (v_1^N, \dots, v_\ell ^N)^T = (\{v_{10}, v_{1b}\}, \dots, \{v_{\ell 0}, v_{\ell b}\})^T \in [S_N^0]^\ell$,

where $\eta$ is the coercivity constant of $\textbf{A}$.

We also introduce the discrete $H^1$ energy-like norm $\Vert|\cdot\Vert|_{\varepsilon}$ in $S_N^\ell+H^1(\Omega)^\ell$ defined as

where ${ v_{i0}^\prime}$ is the ordinary derivative of a functions ${v_{i0}(x)}$.

We show that the norms $\Vert|\cdot\Vert|$ and $\Vert|\cdot\Vert|_\varepsilon$ defined by (4.3) and (4.4), respectively are equivalent in the {weak Galerkin} finite element space $[S_N^0]^\ell$ in the next lemma.

Lemma 4.1. Let $\textbf{v}_N\in [S_N^0]^\ell$. Then there are two positive constant $C_l$ and $C_s$ such that

Proof. For any $v^N_i = \{ v_{i0}, v_{ib}\}\in S_N^0$, by the definition of weak derivative (3.3) and integration by parts we arrive at

Choosing $w = d_w v_i^N$ in the above Eq (4.6) yields

Summing up the above equation over all $I_n\in \mathcal{T}_N$, and using the inverse inequality (4.2), we obtain

Therefore, we have

From the penalty parameter (3.6), we have

To see this, the minimal possible $h_n = H_0 = \cfrac{2(\ell+1)\lambda_1}{N}$ implies that $\cfrac{\varepsilon_1 h_n^{-1}}{\varrho_n} = \cfrac{\alpha}{2(\ell +1) \sigma} = :C$. Hence using (4.8), we obtain

which together with (4.7) implies that

On the other hand, taking $w = v_{i0}^\prime$ in the Eq (4.6) yields

Summing up the above equation over all $I_n\in \mathcal{T}_N$, using the inverse inequality (4.2), we have

Therefore, we have

With the help of (4.8), we result in

We obtain the desired result (4.5) in view of the inequalities (4.9) and (4.11) and the definition of the norms $\Vert|\cdot\Vert$ and $\Vert|\cdot \Vert|_\varepsilon$. Thus we complete the proof. □

We next show that the coercivity of the bilinear form $a(\cdot, \cdot)$ on $[S_N^0]^\ell$ in the energy norm $\Vert|\cdot\Vert|$ defined by (4.3).

Lemma 4.2. Let $\textbf{v}_N \in [S_N^0]^\ell$. Then there holds

Proof. Using the coercivity (2.4) of the reaction matrix $\bf{A}$, we have

The proof is completed. □

In light of Lemma 4.2, we deduce that

which in turn implies the problem (3.4) has a unique solution. The existence follows from the uniqueness.

As a result of Lemma 4.1 and Lemma 4.2, we conclude that the bilinear form $a(\cdot, \cdot)$ is also coercive in the energy like norm $\Vert|\cdot \Vert|_\varepsilon$ defined by (4.4).

Lemma 4.3. Let $\textbf{v}_N, \textbf{w}_N \in [S_N^0]^\ell$. Then there exist positive constants $C_c$ and $C_e$ such that

5.   Error analysis in the energy norm

In this section, we study the error analysis of the proposed numerical scheme applied to the problem (1.1) in the energy norm associated with the bilinear form. We will show that the WG-FEM solution converges uniformly in the energy norm with respect to the perturbation parameters. For the uniform convergence analysis on Shishkin mesh, we will use a special interpolation operator given in ^[11]. On each interval $I_n$, we introduce the set of $k+1$ nodal functional $N_\ell$ defined as follows: for any $v\in C(I_n)$

A local interpolation $\mathcal{I} :H^1(I_n)\to \mathbb{P}_{{k}}(I_n)$ is now defined by

The local interpolation operator $\mathcal{I}$ can used for constructing a continuous global interpolation.

Since $\mathcal{I} v|_{I_n}$ is continuous on $I_n$ and is in the $H^1(I_n)$ space, we denote $\mathcal{I} v|_{\partial I_n}$ by $\mathcal{I} v|_{I_n}$, for simplicity. Form this fact we observe that for any $v\in H^1(I_n)$ we have

Lemma 5.1. ^[11] For any $w\in H^{k+1}(I_n), I_n \in \mathcal{T}_N$, the interpolation $\mathcal{I} w$ defined by (5.1) has the following estimates:

where $h_n$ is the length of element $I_n$ and $C$ is independent of $h_n,$ and {$\varepsilon_i, \; i = 1, \dots, \ell.$}

Lemma 5.2. Let $\mathcal{I} \boldsymbol{R}$ and $\mathcal{I} \boldsymbol{L}$ be the interpolations of the regular part $\boldsymbol{R}$ and the layer part $\boldsymbol{L}$ of the solution {$\boldsymbol{u} \in H^{k+1}(\varOmega)$} on the piecewise-uniform Shishkin mesh, respectively. Assume also that $\varepsilon_\ell \ln N \leq \ell \alpha /(2(\ell+1)\sigma)$ and let $\varOmega_\lambda = [\lambda_\ell, 1-\lambda_\ell]$. Then, we have $\mathcal{I} \boldsymbol{u} = \mathcal{I} \boldsymbol{R}+ \mathcal{I} \boldsymbol{L}$ and the following interpolation estimates are satisfied for $i = 1, \dots, \ell$

where $\varepsilon^{1/2}: = \varepsilon_1^{1/2}+\dots, +\varepsilon_\ell^{1/2}$. Furthermore, the following estimates hold true

Proof. The linearity of the interpolation implies that $\mathcal{I} \textbf{u} = \mathcal{I}(\textbf{R}+ \textbf{L}) = \mathcal{I} \textbf{R}+ \mathcal{I} \textbf{L}.$ Applying the estimate (5.3), the bounds for the derivatives of regular components $R_i$ of the solution in Lemma 2.1 and using (2.8), we obtain

This completes the proof of estimates (5.5).

Using the fact that $\mathcal{B}_{\varepsilon_i}^\alpha(x)\leq \mathcal{B}_{\varepsilon_\ell}^\alpha (x)$ for $i = 1, \dots, \ell$ and $\lambda_\ell = \cfrac{\sigma \varepsilon_\ell}{\alpha} \ln N$, we have

The $L^2$- norm estimate of the layer part of the solution on the sub-interval $\varOmega_\lambda$ follows from

Hence, from the above inequalities we have

Thus, we complete the proof of the estimate (5.8).

We also have

This proves the estimate (5.9).

Due to (5.3) of Lemma 5.1 and the bounds for derivatives (2.9), we obtain at once

Thus, the estimate (5.6) is proved.

For the proof of (5.7) we follow ^[11]. An inverse estimate yields that

We will derive a bound for $\Vert \mathcal{I} L_i\Vert_{L^2(\varOmega_\lambda)}$. For the interval $I_n = (x_{n-1}, x_n)$, we have the estimate for the local nodal functional $N_m(L_i)$ as

The local representation

implies that

where we use the fact $\Vert\phi_m \Vert_{L^2(I_n)}\leq CN^{-1}$. Summing up over all $I_n \subset \varOmega_\lambda$ yields that

Since the mesh size on $\varOmega_\lambda$ is $H_\ell = H_{\ell+1}$, the term in the parenthesis on the right hand side of the above inequality can be written as

Integrating the above inequality on $I_n\subset \varOmega_\lambda$ and using the fact that $H_\ell = \mathcal{O}(N^{-1})$, we have

Summing up the above inequality for $n = \frac{\ell N}{2(\ell+1)}+1, \dots, \frac{(\ell+2) N}{2(\ell+1)}-1$ leads to

It remains to bound on the last interval $(x_{\frac{(\ell+2) N}{2(\ell+1)}-1}, x_{\frac{(\ell+2) N}{2(\ell+1)}})$. From the inequality (5.13), we have

These two last estimates give the desired estimate. Thus the estimate (5.7) is proved.

From (5.7) and (5.8), we get

which completes the proof of (5.10).

Using the triangle inequality and (5.7) and (5.9), we have

Similarly, using the inverse estimate, we get

Hence, we complete the proof of (5.11).

By (5.3) and (2.10), we have for $l = 1, 2$,

which shows (5.12). Thus we complete the proof.

□

The exact solution of problem (1.1) does not satisfy the WG-FEM scheme (3.4) and hence the WG-FEM lacks of consistency. Consequently, inconsistency leads to loss of the classical Galerkin orthogonality. As a result, we follow different techniques from the ones used in the standard finite element procedure to derive the error estimates.

Now Strang's second lemma provides a quasi-optimal bound for $\Vert| \textbf{u}- \textbf{u}_N\Vert|_\varepsilon$.

Theorem 5.3. Let $\boldsymbol{u}$ and $\boldsymbol{u}_N$ be the solutions of problems (1.1) and (3.4) respectively. Then there exists a positive constant $C$ independent of $N$ and $\varepsilon_i$ such that

where $a(\cdot, \cdot)$ is the bilinear form given by (3.5).

First, we will establish some error equations which will be needed in the error analysis below.

Lemma 5.4. Let $\boldsymbol{u} = (u_1, \dots, u_\ell)$ be the solution of the problem (1.1). Then for any $\boldsymbol{v}_N = (v_1^N, \dots, v_\ell^N) = (\{v_{10}, v_{1b}\}, \dots, \{v_{\ell 0}, v_{\ell b}\})\in [S_N^0]^\ell$, we have

where

Proof. For any $\textbf{v}_N\in [S_N^0]^\ell$, using the commutative property (5.2) of the interpolation operator we have

Using the definition of the weak derivative (3.3) and integration by parts, we have

From the definition of the interpolation and integration by parts, we obtain

which implies that

We infer from the Eqs (5.19)–(5.21) that

Summing up the Eq (5.22) over all interval $I_n\in \mathcal{T}_N$, we find

Using integration by parts, one can show that

Summing up the above equation over all interval $I_n\in \mathcal{T}_N$, we get

where we used the fact that $\big\langle u_i^{\prime }, \textbf{n} v_{ib}\big\rangle = 0.$ Finally, by plugging the Eq (5.24) into (5.23), we arrive at the desired result (5.15).

Lastly, the Eq (5.18) clearly holds. We complete the proof.

□

The following lemma will be useful in the error analysis.

Lemma 5.5. Assume that $\textbf{u} = (u_1, \dots, u_\ell), \; \mathit{\text{with}}\; u_i \in H^{k+1}(\varOmega)$ is the solution of the problem (1.1). Then we have the following estimate

where $\theta_i = u_i- \mathcal{I} u_i$ for $i = 1, \dots, \ell.$

Proof. From the trace inequality (4.1), we can write

It remains to estimate $\Vert \theta_i^\prime \Vert_{L^2(I_n)}$ and $\Vert \theta_i^{\prime \prime} \Vert _{L^2(I_n)}$, individually. From the estimate (5.5), one has

With the help of the estimate (5.11) and (5.12) one can show that

With the help of the above estimates, the fact that $\sigma\geq k+1$, and the triangle inequality, one can conclude that

and

The desired result follows from combining the above estimates and the mesh size $h_n$. Thus, we complete the proof. □

Lemma 5.6. Assume that $u_i \in H^{k+1}(\varOmega)$ and the penalization parameter $\varrho_n$ is given by (3.6). If $\sigma \geq k+1$, then we have

where $T(\boldsymbol{u}, \boldsymbol{v}_N) = \sum_{i = 1}^\ell T_1(u_i, v_i^N)+T_2(\boldsymbol{u}, \boldsymbol{v}_N)$ and $C$ is independent of $N$ and $\varepsilon_i, \; i = 1, \dots, \ell$.

Proof. It follows from Cauchy-Schwarz inequality, Lemma 5.5 and the penalization parameter (3.6) that

As a result

We next bound the term $T_2(\textbf{u}, \textbf{v}_N)$. We need to estimate $\Vert u_i- \mathcal{I} u_i \Vert, \quad i = 1, \dots, \ell$. Using the estimates (5.5)–(5.8) of Lemma 5.2 and Cauchy-Schwarz inequality, taking $\sigma \geq k+1$, we get

The above estimate (5.30) and Cauchy-Schwarz inequality yield the following bound

From the above estimate, we have

From the estimates (5.29) and (5.31), we have the desired result. Thus we complete the proof. □

Theorem 5.7. Let $\boldsymbol{u} = (u_1, \dots, u_\ell) = \boldsymbol{R}+ \boldsymbol{L}\; \mathit{\text{with}}\; u_i\in H^{k+1}(\varOmega)$ be the solution of the problem (1.1) and assume that the conditions of Lemma 5.2 hold with $\sigma\geq k+1$. Then, the following estimates hold true:

where $C$ is independent of $N$ and $\varepsilon_i, i = 1, \dots, \ell$.

Proof. Since $\theta_i^R: = R_i- \mathcal{I} R_i$ and $\theta_i^L: = L_i- \mathcal{I} L_i$ are continuous on $\varOmega$, we get $s(\theta_i^R, \theta_i^R) = s(\theta_i^L, \theta_i^L) = 0$ for $i = 1, \dots, \ell$. Then we have

In the light of the interpolation errors (5.5) and (2.8), we obtain for $i = 1, \dots, \ell$

which together with (5.33) yields

Using the inequalities (5.11), (5.12) and (5.30), we have

which together with (5.34) gives the desired result

where we have used $\sigma \geq k+1$. The proof is completed. □

We next estimate the consistency error $\sup_{\textbf{w}_N\in [S_N^0]^\ell } \cfrac{\vert a(\textbf{u}, \textbf{w}_N)-L(\textbf{w}_N)\vert}{\Vert|\textbf{w}_N \Vert|_\varepsilon}$.

Lemma 5.8. Assume that $\boldsymbol{u} = (u_1, \dots, u_\ell), u_i \in H^{k+1}(\Omega), i = 1, \dots, \ell$ is the solution of (1.1). If $\sigma \geq k+1$, then the following estimate holds true:

where $C$ is independent of $N$, and $\varepsilon_i, \; i = 1, \dots, \ell$.

Proof. Using the definition of bilinear form (3.4), the fact that $s(u_i, w_i^N) = s(u_i- \mathcal{I} u_i, w_i^N) = 0$ for $i = 1, \dots, \ell$ and Lemma 5.4, we have

where $T(\textbf{u}, \textbf{w}_N) = T_1(\textbf{u}, \textbf{w}_N)+T_2(\textbf{u}, \textbf{w}_N)$ and $T_1(\textbf{u}, \textbf{w}_N) = \sum_{i = 1}^\ell T_1(u_i, w_i^N)$ with $T_1(u_i, w_i^N)$ and $T_2(\textbf{u}, \textbf{w}_N)$ are given in (5.17) and (5.18), respectively. By Lemma 5.6, if $\sigma \geq k+1$, the first term on the right side of the above equation can be estimated as

For the second term, we use the continuity property (4.14) of bilinear form $a(\cdot, \cdot)$ and again the fact that $s(u_i - \mathcal{I} u_i, w_i^N) = 0, i = 1, \dots, \ell$ and we obtain

Appealing the estimates (5.5), (5.11), (5.12), (5.30) and using (2.8), if $\sigma \geq k+1$ we obtain

From (5.35) and (5.36), we arrive at

which is the desired result. We complete the proof. □

Theorem 5.9. Assume that $\boldsymbol{u} = (u_1, \dots, u_\ell), u_i \in H^{k+1}(\Omega), i = 1, \dots, \ell$ is the exact solution and $\boldsymbol{u}_N\in [S_N^0]^\ell$ is the WG-FEM solution given by (3.4) on the uniform Shishkin mesh for the problem (1.1), respectively. If $\sigma \geq k+1$, then we have the following estimate

where $C$ is independent of $N$ and $\varepsilon_i, i = 1, \dots, \ell$.

Proof. Using Theorem 5.7, if $\sigma\geq k+1$ we have

Hence, we obtain

Note that $\mathcal{I} \textbf{R}$ and $\mathcal{I} \textbf{L}$ are both in $[S_N^0]^\ell$, the set of piecewise discontinuous polynomials {of degree at most $k$}, and that $\mathcal{I} \textbf{R} + \mathcal{I} \textbf{L} \in [S_N^0]^\ell$. Take $\textbf{v}_N = (\mathcal{I} R_1+ \mathcal{I} L_1, \dots, \mathcal{I} R_\ell + \mathcal{I} L_\ell) \in S_N^\ell$. Invoking Theorem 5.3 and Lemma 5.8, the desired result follows.

□

6.   Balanced norm estimates

As stated in the introduction, error estimates in the corresponding energy norm of FEMs not adequate. The reason arises from the fact that the energy norm of the boundary layer functions $\exp(\frac{-x}{\varepsilon_\ell})$ and $\exp(\frac{-(1-x)}{\varepsilon_\ell})$ are of order $\mathcal{O}(\varepsilon_\ell^{1/2})$. Therefore, the error estimates in the energy norm is not much strong than the $L^2$-norm if $\varepsilon_\ell \ll 1$. A stronger norm obtained by scaling of the coefficient of the $H^1$-seminorm captures correctly the boundary layers. This norm is called the balanced norm defined as follows. For $\textbf{v} = (v_1^N, \dots, v_\ell ^N)^T = (\{v_{10}, v_{1b}\}, \dots, \{v_{\ell 0}, v_{\ell b}\})^T \in [S_N^0]^\ell$,

where $s^b(\textbf{u}_N, \textbf{v}_N)$ is given by

Here, the penalization parameter $\varrho_n^b$ is now defined as

where $\varepsilon = \sum_{i = 1}^\ell \varepsilon_i$.

We note that the error bound $N^{-(k+1)}$ independent of $\varepsilon^{1/2}$ in Theorem 5.9 comes from the estimate of the $L^2$-norm of $\textbf{u}-\mathcal{I} \textbf{u}$ in the energy norm error estimates. These terms can be handled by replacing the special interpolation operator $\mathcal{I}$ defined by (5.1) with a projection operator $Q_h:H^1(I_n)\to S_N$ defined as follows.

Let $P_h:L^2(I_n)\to \mathbb{P}_k(I_n)$ be the local weighted $L^2$-projection restricted to interval $I_n$ defined by

This weighted $L^2$-projection is well-defined because we assume that the diagonal elements are positive and the reaction coefficient matrix is strongly diagonally dominant matrix. With the aid of the Bramble-Hilbert lemma, one can show that for $i = 1, \dots, \ell$

We introduce the projection operator $Q_{h}:H^1(I_n)\to S_N$ such that

Clearly, $Q_h u_i\in S_N^0$ if $u_i\in H_0^1(I_n)$ for $i = 1, \dots, \ell$. By (6.5), we have

The following trace and inverse inequalities will be used in the forthcoming analysis ^[43]. For any function $\phi \in H(I_n),$ we have

We would like to derive similar estimates as Lemma 5.2 for the projection operator $Q_0$ which is essentially the generalized $L^2$-projection. The following lemma will serve this purpose.

Lemma 6.1. Assume that the conclusions of Lemma 5.2 hold. Then we have the following error estimates for the operator $Q_0$ on the uniform Shishkin mesh.

Proof. It is known that the $L^2$-projection $Q_0$ is $L^\infty$-stable ^[44]. Therefore, by the triangle inequality we have

which proves (6.10). Using this inequality, we get

where we used Lemma 5.2 and the fact that $h_n = \mathcal{O}(N^{-1})$ in $\varOmega_\lambda$. It follows from (6.7) that

The first term on the right side of the above inequality can be bounded as

Next, we estimate the layer parts on $\varOmega\setminus\varOmega_\lambda$.

From (6.15) and (6.16), we get

which completes the proof of (6.12).

With the help of an inverse inequality on $\varOmega_\lambda$, we obtain

because $\Vert \mathcal{I} u_i-u_i\Vert _{L^2(I_n)}$ and $\Vert Q_0 u_i-u_i\Vert _{L^2(I_n)}$ are both of order $\mathcal{O}(N^{-(k+1)})$ on $\varOmega_\lambda$. By Lemma 5.2 and the above estimate, we arrive at

On the other hand, from Lemma 5.2, we have

Combining (6.17) and (6.18) gives the desired result (6.13). Finally, using an inverse estimate we obtain at once

which proves (6.14). Thus, we complete the proof. □

We derive the following error equations involving the projection $Q_h$ which are similar to ones in Lemma 5.4. To this end, we still need another special projection operator defined as follows. We refer interested readers to ^[45] for details.

Lemma 6.2. ^[45] For $u_i\in H^1(\varOmega)$, there is a projection operator $\pi_h u_i\in H^1(0, 1)$, restricted on element $I_n$, $\pi_h u_i \in \mathbb{P}_{k+1}(I_n)$ satisfies

Lemma 6.3. Let $\boldsymbol{u} = (u_1, \dots, u_\ell)$ be the solution of the problem (1.1). Then for any $\boldsymbol{v}_N = (v_1^N, \dots, v_\ell^N) = (\{v_{10}, v_{1b}\}, \dots, \{v_{\ell 0}, v_{\ell b}\})\in [S_N^0]^\ell$, we have

where

Proof. Using the definition of the operator $Q_h$, the weak derivative (3.3), integration by parts and (6.19), we have

where $v_{n} = v(x_{n})$ and $v_{n-1} = v(x_{n-1})$ for a function $v$. This implies that

Following the same procedures as in the energy norm estimates, we prove (6.21). Clearly, we have the Eq (6.24). We complete the proof. □

Lemma 6.4. Assume that $\textbf{u} = (u_1, \dots, u_\ell), \;\mathit{\text{with}}\; u_i \in H^{k+1}(\varOmega)$ is the solution of the problem (1.1) and the penalization parameter $\varrho_n^b$ is given by (6.3). Then we have

where $C$ is independent of $N$ and $\varepsilon_i, \; i = 1, \dots, \ell$, $T^b(\boldsymbol{u}, \boldsymbol{v}_N) = T_1^b(\boldsymbol{u}, \boldsymbol{v}_N)+T_2^b(\boldsymbol{u}, \boldsymbol{v}_N)$, and $s^b(Q_h u_i, v_N)$ is given by (6.2).

Proof. Note that $T_2^b(\textbf{u}, \textbf{v}_N) = 0$ due to the definition of the projection $Q_h$. By the inverse estimate (6.9), Lemma 5.5 and Lemma 6.2, we obtain at once

where $\xi_i = u_i-\pi_h u_i$ and $\theta_i = u_i- \mathcal{I} u_i$ for $i = 1, \dots, \ell.$

İmitating the arguments in the energy norm estimates and using the above fact, one can prove that

It follows from Cauchy–Schwarz inequality, the trace inequality (6.8) and Lemma 6.1 that

Here, we used the fact that $\varepsilon^{-1}\varrho_n^b = \varrho_n$. Therefore, we complete the proof. □

The main result of this section is the following theorem.

Theorem 6.5. Assume that $\boldsymbol{u} = (u_1, \dots, u_\ell), u_i \in H^{k+1}(\Omega), i = 1, \dots, \ell$ is the exact solution and $\boldsymbol{u}_N = \{u_1^N, \dots, u_\ell^N\}\in [S_N^0]^\ell$ is the WG-FEM solution given by (3.4) on the uniform Shishkin mesh for the problem (1.1), respectively. If $\sigma \geq k+1$, then we have the following improved balanced error estimate

where $C$ is independent of $N$ and $\varepsilon_i, \; i = 1, \dots, \ell$.

Proof. From Lemma 6.1 and Lemma 6.4, we obtain at once

where we have used that $s^b(u_i, u_i) = 0$. Imitating the analyses in the energy norm estimates and using again Lemma 6.4, we have

Therefore, we obtain

Next, using the triangle inequality and combining the above estimate with (6.28) yield

The proof is now completed. □

7.   Numerical experiment

We present various numerical experiments to show the performance of the WG-FEM in this section. All the integration was calculated by using $5$-point Gauss-Legendre quadrature integral formula.

Example 7.1. Consider the following coupled system of reaction-diffusion problem with constant coefficients

where $\mathcal{E} = diag(\varepsilon_1^2, \varepsilon_2^2)$ with $0 < \varepsilon_1\leq \varepsilon_2 < < 1$, \quad $\boldsymbol{g} = (g_1, g_2)^T, \quad \boldsymbol{A} = \begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix}$ and $g_1, g_2$ are chosen such that

is the exact solution $\boldsymbol{u}(x) = (u^1(x), u^2(x))$ of the system of reaction-diffusion problem (7.1). Note that $R_i(x), i = 1, 2$ is constant and (2.8) holds. We know that the solution has exponential layers of width $\mathcal{O}(\varepsilon_2 |\ln \varepsilon_2|)$ at $x = 0$ and $x = 1$, while only $u^1(x)$ has an additional sublayer of width $\mathcal{O}(\varepsilon_1 |\ln \varepsilon_1|)$. We take $\rho > 1/2$, $\alpha = 0.99$ and $\sigma = 3$ for this problem.

We applied the WG-FEM (3.4) for solving the problem (7.1). The numerical errors $\textbf{e}: = \textbf{u}- \textbf{u}_N$ are computed in the energy norm by

for a fixed $\varepsilon_1, \varepsilon_2$ and $N$. We report the numerical experiments for the uniform error calculated by

in Table 1. The order of convergence $r_{\varepsilon}$ is computed using mesh levels $(N_1, ||| \textbf{e}^{N_1}|||_{\varepsilon})$ and $(N_2, ||| \textbf{e}^{N_2}|||_{\varepsilon})$:

Table 1 shows that the energy norm error estimates exhibit $k$-order convergence which agrees perfectly with the theoretical error estimates.

In order to pay attention to the dependency of the energy norm on the parameters, we compute the energy norm estimates for a fixed $\varepsilon_1$ and different values of $\varepsilon_2$. For instance, we first fixed $\varepsilon_1 = 10^{-10}$ and take different values of $\varepsilon_2 = 10^{-4}, \dots, 10^{-9}$. The results are presented in Table 2, Table 4, Figures 2a and 2b. These results verify that the method is robust on the uniform Shishkin mesh and the order of convergence is of $\mathcal{O}(\varepsilon^{1/2} (N^{-1}\ln N)^k)$, where $\varepsilon^{1/2} = \varepsilon_1^{1/2}+\varepsilon_2^{1/2}$ using the linear $k = 1$ and quadratic $k = 2$ element functions, which is in excellent agreement with the main result of Theorem 5.9. Moreover, we infer from Table 2 that $\cfrac{|||u-u_N|||_{\varepsilon_{{j}}}}{|||u-u_N|||_{\varepsilon_{{j+2}}}}\approx \sqrt{\frac{{{10^{-j}}}}{{{10^{-(j+2)}}}}}$ for $\varepsilon_{{j}} = \{{10^{-10}}, {{10^{-j}}}\}, \quad {j = 4, \dots, 9, }$ where $|||u-u_N|||^2_{\varepsilon_{{j}}} = (10^{-10})^2 \Vert d_w e_1^N \Vert^2+(10^{-j})^2 \Vert d_w e_2^N \Vert^2+ \eta \sum_{i = 1}^2 \Vert e_{i0}\Vert^2+\sum_{i = 1}^2 s(e_i^N, e_i^N)$.

This implies that the errors might be affected by a term involving $\sqrt{\varepsilon_2}$. We observe almost linear convergence up to a logarithmic factor using the linear elements and almost quadratic convergence using the quadratic elements if $N$ gets larger. Hence, for larger $N\geq 64$, the rate of convergence is of order $\mathcal{O}(\varepsilon_2^{1/2}(N^{-1}\ln N)^k)$ which agrees with the theory indicated by Theorem 5.9. Table 2 shows that the errors and the order of convergence are dominated by the term $N^{-(k+1)}$ when $N$ and $\varepsilon$ are smaller. We also observe that if $\varepsilon_2$ decreases for a fixed $\varepsilon_1$, the energy norm error estimates get smaller. These observations suggest that the main result of Theorem 5.9 is sharp.

On the other hand, we compute the numerical errors $\textbf{e}: = \textbf{u}- \textbf{u}_N$ with respect to the balanced norm by

for a fixed $\varepsilon_1, \varepsilon_2$ and $N$. We list the uniform balanced error bounds $\textbf{e}^{N, b}$ calculated as before in Table 3. We also report the numerical results in the balanced norm in Table 4 and we notice that the error estimates in the balanced norm remain almost unchanged as $\varepsilon_2$ decreases for a fixed $\varepsilon_1$ unlike the estimates in the energy norm. This confirms the theory stated in Theorem 6.5. We have plotted the balanced norm error estimates for a fixed $\varepsilon$ and varying $\varepsilon_2$ on log-log scale in Figures 2c and 2d for a viewable illustrations. Evidently, the errors stay almost constant while the parameters vary and behave like

This confirms the result of Theorem 6.5 up to a root of $\ln N$.

Example 7.2. We next consider the problem (1.1) with variable coefficients

We take $\rho = 3/4$, $\alpha = 0.80$ and $\sigma = 3$. The exact solution is unknown. Hence a finer mesh constructed as below is used for estimating the numerical errors.

We compute the errors $\textbf{e} = \textbf{u}_N- \textbf{u}_{2N}$ where $\textbf{u}_{2N}$ is our numerical solution computed on a mesh consisting of the initial uniform Shishkin mesh and the midpoints $x_{n+1/2} = \frac{x_n+x_{n+1}}{2}, n = 0, \dots, N-1$. Therefore we calculate

for a fixed $\varepsilon_1, \varepsilon_2, \varepsilon_3$ and $N$. The numerical results are listed in Tables 5 and 6 for the uniform errors in the energy and balanced norms, respectively

where $\textbf{e}_{\varepsilon_1, \varepsilon_2, \varepsilon_3 }^{N, b}$ is defined as before and the order for convergence is calculated by (7.2). The results clearly suggest that the $k$-order uniform convergence in the energy norm, which is in good agreement with the main result of Theorem 5.9. The errors in the balanced norm behave like $\mathcal{O}(N^{-1}\ln N)^k$ which agrees with the results of Theorem 6.5 up to a square root of $\ln N$. As before, we observe from Table 7 and Figures 3a and 3b that the energy norm estimates depend on $\varepsilon^{1/2} = \varepsilon_1^{1/2}+\varepsilon_2^{1/2}+\varepsilon_3^{1/2}$ and errors change decreasingly as $\varepsilon\to 0$ while the balanced norm estimates do not depend on the parameters and the errors remain almost unchanged as seen from Table 8 and Figures 3c and 3d.

8.   Conclusions

In this paper, we studied the WG-FEM for system of SPPs of reaction-diffusion type in which the equations have diffusion parameters of the different magnitudes on a piecewise uniform Shishkin mesh. With the help of a special interpolation operator, we derived optimal and uniform error bounds in the energy and the balanced norms up to a logarithmic factor. The proposed WG-FEM uses the procedure of elimination of the interior unknowns from the discrete linear system and thus the method is comparable with the classical FEM. We will investigate sharper error bounds in balanced norm and extend these results to to high dimensional problem on a tensor product meshes in the future work.

Acknowledgements

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

Conflict of interest

The authors declare that they have no conflict of interest.