∎

11institutetext: Corresponding author: Jin Zhang 22institutetext: School of
Mathematics and Statistics, Shandong Normal University, Jinan 250014, China
22email<EMAIL_ADDRESS>33institutetext: Xiaoqi Ma 44institutetext:
School of Mathematics and Statistics, Shandong Normal University, Jinan
250014, China
44email<EMAIL_ADDRESS>
# Uniform convergence of optimal order of a local discontinuous Galerkin
method on a Shishkin mesh under a balanced norm

Xiaoqi Ma  Jin Zhang  Wenchao Zheng

(Received: date / Accepted: date)

###### Abstract

This article investigates a local discontinuous Galerkin (LDG) method for one-
dimensional and two-dimensional singularly perturbed reaction-diffusion
problems on a Shishkin mesh. During this process, due to the inability of the
energy norm to fully capture the behavior of the boundary layers appearing in
the solutions, a balanced norm is introduced. By designing novel numerical
fluxes and constructing special interpolations, optimal convergences under the
balanced norm are achieved in both 1D and 2D cases. Numerical experiments
support the main theoretical conclusions.

###### Keywords:

Singular perturbation Local discontinuous Galerkin Shishkin mesh Balanced norm

###### MSC:

65N30 65N50

##  1 Introduction

Singular perturbation problems are popularly used in various fields, such as
fluid mechanics, energy development and electronic science R2012-finite ;
m2021-local ; z2013-point . Their typical feature is that one or more layers
usually appear in their solutions. To fully resolve the layers and derive the
uniform convergence with respect to perturbation parameters, layer-adapted
meshes were introduced in the 1960s B1969-Towards  and gradually became an
active research field Ma1Zha2:2023-S ; miller1996-fitted . Due to its simple
structure, Shishkin meshes have attracted the attention of many researchers,
see Xie1Zhu2Zho3:2013-U ; Zha1Liu2:2017-S ; miller1996-fitted .

At present, many numerical methods have been well analyzed on Shishkin meshes,
such as the finite difference method, the finite element method and so on
Styn1Styn2:2018-Convection-diffusion ; R2008-Robust . However, these methods
cannot accurately calculate the gradient of the solution. To solve this
problem, the local discontinuous Galerkin (LDG) method C1998-local  is a good
choice. As a kind of discontinuous Galerkin (DG) method, the LDG method
inherits many advantages of the DG method C2012-Discontinuous , such as high-
order accuracy, good stability and flexibility on $hp$ adaptivity. These
advantages and its ease of calculating gradients make the LDG method more
suitable for handling those solutions with layers. Therefore, the LDG method
has been applied to different singularly perturbed problems in recent years
c2022-balanced ; z2013-point .

In this paper, we mainly discuss singularly perturbed reaction-diffusion
equations. For these reaction-diffusion problems, the authors in
m2021-convergence  have analyzed uniform convergence under an energy norm. But
the energy norm is too weak to fully capture the behavior of the layers
appearing in the solutions. For the sake of addressing this issue, researchers
have introduced a new norm–balanced norm, which is more suitable than the
energy norm for layers; see R2012-finite ; R2015-convergence . In
c2022-balanced , Cheng et al. derived a uniform convergence of optimal order
under the assumption that the smooth component of the solution belongs to the
finite element space. But the assumption is too strong and generally cannot be
established.

Our goal is to design a LDG method to avoid the aforementioned drawbacks. For
this purpose, we design a special numerical flux, and it is the key to our LDG
method. Furthermore, we also introduce a new projection, which consists of a
Gauß-Radau projection and a (weighted) local $L^{2}$ projection. Based on the
new numerical fluxes and projection, it is proven that the LDG solution is
convergent with optimal order under the balanced norm. Finally, the LDG method
and its convergence analysis are extended from one-dimensional to two-
dimensional.

The framework of the article is arranged as follows. First, we obtain the
uniform convergence of the LDG method for a one-dimensional singularly
perturbed problem in Section 2. In this progress, a LDG method is introduced
and a specific numerical flux is designed. Then we define a Shishkin mesh and
present some projections. Finally, the optimal convergence order is derived in
the balanced norm. In Section 3, we study the LDG method for a singular
perturbation problem in 2D, and the content distribution is the same as in
Section 2. Finally, a numerical experiment is presented to verify the
theoretical results. In the paper, $C$ is a general positive constant which is
independent of $\epsilon$ and $N$. Here $\epsilon$ is the perturbation
parameter and $N$ is the mesh parameter. Furthermore, let $k\geq 1$ be a fixed
positive integer.

##  2 LDG method for 1D case

We study the following singularly perturbed two-point boundary value problem:

(1) $-\epsilon
u^{\prime\prime}+b(x)u=f(x),\quad\text{$x\in\Omega=(0,1)$},\quad u(0)=u(1)=0,$

in which $0<\epsilon\ll 1$ is the perturbation parameter,
$b(x)\geqslant\beta^{2}>0$ and $f(x)$ are sufficiently smooth functions on
$\overline{\Omega}$. Here $\beta$ is some positive constant. The solution of
(1) is typically characterized by boundary layers at $x=0$ and $x=1$ of width
$O(\sqrt{\epsilon}\ln(1/\epsilon))$.

###  2.1 LDG method

First, we divide $\Omega=(0,1)$ into $0=x_{0}<x_{1}<\ldots<x_{N-1}<x_{N}=1.$
Define $[x_{i-1},x_{i}]$ as $I_{i}$ for $i=1,2,\ldots,N$, and the step size of
the interval $I_{i}$ as $h_{i}=x_{i}-x_{i-1}$. Then the discontinuous Sobolev
space $V_{N}$ is denoted as

$V_{N}=\\{v\in L^{2}(\Omega):\text{$v|_{I_{i}}\in\mathbb{P}_{k}(I_{i}),1\leq
i\leq N$ }\\},$

where $\mathbb{P}_{k}(I_{i})$ is the space of polynomials of degree at most
$k\geq 1$ on $I_{i}$. For $v\in V_{N}$, we introduce
$v^{\pm}_{i}=\lim\limits_{x\to x^{\pm}_{i}}v(x)$, and define a jump as:
$[[v]]_{i}=v^{-}_{i}-v^{+}_{i}$ for $i=1,2,\ldots,N-1$, $[[v]]_{0}=-v^{+}_{0}$
and $[[v]]_{N}=v^{-}_{N}$.

Next, the LDG method is introduced for the problem (1). Firstly, let us
rewrite it as the following first-order system,

$\displaystyle\epsilon^{-1}q=u^{\prime},\quad\text{in $\Omega$},$
$\displaystyle-q^{\prime}+b(x)u=f(x),\quad\text{in $\Omega$}.$

Let $\chi=(r,v)\in V_{N}\times V_{N}$ be any test function. Assume that
$\left<\cdot,\cdot\right>_{I_{i}}$ is the inner product in $L^{2}(I_{i})$.
Find $W=(Q,U)\in V_{N}\times V_{N}$ such that

(2)
$\displaystyle\epsilon^{-1}\left<Q,r\right>_{I_{i}}+\left<U,r^{\prime}\right>_{I_{i}}-\widehat{U}_{i}r^{-}_{i}+\widehat{U}_{i-1}r^{+}_{i-1}=0,$
(3)
$\displaystyle\left<Q,v^{\prime}\right>_{I_{i}}+\left<bU,v\right>_{I_{i}}-\widehat{Q}_{i}v^{-}_{i}+\widehat{Q}_{i-1}v^{+}_{i-1}=\left<f,v\right>_{I_{i}},$

where the numerical fluxes $\widehat{Q}$ and $\widehat{U}$ are defined by

(4) $\displaystyle\widehat{Q}_{i}=\left\\{\begin{aligned}
&Q^{+}_{0}+\lambda_{0}U^{+}_{0},\quad&&\text{$i=0$},\\\
&Q^{+}_{i},\quad&&\text{$i=1,2,\ldots,N-1$},\\\
&Q^{-}_{N}-\lambda_{N}U^{-}_{N},\quad&&\text{$i=N$},\end{aligned}\right.$
$\displaystyle\widehat{U}_{i}=\left\\{\begin{aligned}
&0,\quad&&\text{$i=0,N$},\\\
&U^{-}_{\frac{3}{4}N}-\lambda_{q}[[Q]]_{\frac{3}{4}N},\quad&&i=\frac{3}{4}N,\\\
&U^{-}_{i},\quad&&\text{$i=1,2,\ldots,\frac{3}{4}N-1,\frac{3}{4}N+1,\ldots,N-1$},\\\
\end{aligned}\right.$

where we choose $\lambda_{0}=\lambda_{N}=\epsilon^{\frac{1}{2}}$ and
$\lambda_{q}=\epsilon^{-\frac{1}{2}}$.

###### Remark 1

We design a new numerical flux $\widehat{U}_{i}$, which is the key to
convergence analysis of the LDG method in the balanced norm. Unlike
c2022-balanced , an additional term involved with $[[Q]]_{\frac{3}{4}N}$
appears in $\widehat{U}_{\frac{3}{4}N}$. By utilizing this new numerical flux,
the optimal convergence order can be derived without the impossible condition
in c2022-balanced .

###### Remark 2

The selection of numerical flux (4) in this paper is a variant of the LDG
method, see (Di1Ern2:2012-M, , Section 4.4.2) and Ngu1Per2:2009-motified  for
more details. Due to the numerical flux $\widehat{U}_{\frac{3}{4}N}$ depends
on $[[Q]]_{\frac{3}{4}N}$, we cannot use (2) to solve for $Q$ in terms of $U$
element-by-element on $I_{\frac{3}{4}N}\cup I_{\frac{3}{4}N+1}$. This
precludes the local elimination of $Q$ on $I_{\frac{3}{4}N}\cup
I_{\frac{3}{4}N+1}$, thereby increasing the computational cost of the
approximation method. However, in this case, we can derive an optimal
convergence order under the balanced norm.

Define $\left<w,v\right>=\sum\limits_{i=1}^{N}\left<w,v\right>_{I_{i}}$ and
express the scheme (2) and (3) by a compact form: Find $W=(Q,U)\in V_{N}\times
V_{N}$ such that

$B(W;\chi)=\left<f,v\right>,\quad\forall\chi=(r,v)\in V_{N}\times V_{N},$

where

(5) $\displaystyle B(W;\chi)=$
$\displaystyle\epsilon^{-1}\left<Q,r\right>+\left<bU,v\right>+\left<U,r^{\prime}\right>-\sum_{i=1}^{N-1}U^{-}_{i}[[r]]_{i}+\left<Q,v^{\prime}\right>-\sum_{i=0}^{N-1}Q^{+}_{i}[[v]]_{i}$
$\displaystyle-(Qv)_{N}^{-}+\sum_{i\in\\{0,N\\}}\lambda_{i}[[u]]_{i}[[v]]_{i}+\lambda_{q}[[Q]]_{\frac{3}{4}N}[[r]]_{\frac{3}{4}N}.$

###  2.2 Regularity of the solution

###### Lemma 1

(Referene (miller1996-fitted, , Lemma 6.2)) For $i=0,1,2,3,\cdots,k+2$,
suppose the solution $u$ of (1) allows the decomposition
$u(x)=\bar{u}(x)+u_{\epsilon,1}(x)+u_{\epsilon,2}(x)$ with
$x\in\overline{\Omega}$, then

$|\frac{d^{i}\bar{u}(x)}{dx^{i}}|\leq
C,\quad|\frac{d^{i}u_{\epsilon,1}(x)}{dx^{i}}|\leq
C\epsilon^{-\frac{i}{2}}e^{-\frac{\beta
x}{\sqrt{\epsilon}}},\quad|\frac{d^{i}u_{\epsilon,2}(x)}{dx^{i}}|\leq
C\epsilon^{-\frac{i}{2}}e^{-\frac{\beta(1-x)}{\sqrt{\epsilon}}}.$

Furthermore, for $i=0,1,2,\cdots,k+1$, we have
$q=\bar{q}+q_{\epsilon,1}+q_{\epsilon,2}=\epsilon\frac{d\bar{u}}{dx}+\epsilon\frac{du_{\epsilon,1}}{dx}+\epsilon\frac{du_{\epsilon,2}}{dx}$
and

$|\frac{d^{i}\bar{q}(x)}{dx^{i}}|\leq
C\epsilon,\quad|\frac{d^{i}q_{\epsilon,1}(x)}{dx^{i}}|\leq
C\epsilon^{-\frac{i-1}{2}}e^{-\frac{\beta
x}{\sqrt{\epsilon}}},\quad|\frac{d^{i}q_{\epsilon,2}(x)}{dx^{i}}|\leq
C\epsilon^{-\frac{i-1}{2}}e^{-\frac{\beta(1-x)}{\sqrt{\epsilon}}}.$

###  2.3 Shishkin mesh

Now a Shishkin mesh is introduced. Let $N\geq 4$ be a multiple of $4$ and the
transition point be

$\displaystyle\tau_{t}=\frac{\sigma\sqrt{\epsilon}\ln
N}{\beta}\leq\frac{1}{4},$

where we take $\sigma\geq k+1$. Divide $[\tau_{t},1-\tau_{t}]$ into $N/2$
equidistant subintervals and in addition, the $N/4$ intervals are evenly
distributed on both the intervals $[0,\tau_{t}]$ and $[1-\tau_{t},1]$.
Therefore, the Shishkin mesh are presented by

(6) $x_{i}=\left\\{\begin{split}&4\frac{\sigma\sqrt{\epsilon}}{\beta}t_{i}\ln
N,\quad&&\text{for $i=0,1,\ldots,\frac{1}{4}N$},\\\
&\tau_{t}+2(1-2\tau_{t})(t_{i}-\frac{1}{4}),&&\text{for
$i=\frac{1}{4}N+1,\ldots,\frac{3}{4}N$},\\\ &1-\frac{4\sigma\sqrt{\epsilon}\ln
N}{\beta}(1-t_{i}),\quad&&\text{for
$i=\frac{3}{4}N+1,\ldots,N$},\end{split}\right.$

in which $t_{j}\equiv\frac{j}{N}$ for $j=0,1,\ldots,N$. It is obvious that on
the layer-adapted mesh (6),

$C\sqrt{\epsilon}N^{-1}\ln N\leq h_{j}\leq CN^{-1}$

for $1\leq j\leq N$ hold.

###  2.4 Projections

For the later analytical proof, we design a special projection that is
obtained by combining the local (or weighted local) $L^{2}$ projection and the
Gauß-Radau projection.

  * •

Local $L^{2}$ projection $\pi$. For any $w\in L^{2}(I_{i})$, $\pi
w\in\mathbb{P}_{k}(I_{i})$ is denoted as:

$\displaystyle\left<\pi
w,v\right>_{I_{i}}=\left<w,v\right>_{I_{i}},\quad\forall
v\in\mathbb{P}_{k}(I_{i}),$

for $I_{i}=(x_{i-1},x_{i})$, $i=1,2,\ldots,N$.

  * •

Weighted local $L^{2}$ projection $\pi_{b}$. For any $w\in L^{2}(I_{i})$,
$\pi_{b}w\in\mathbb{P}_{k}(I_{i})$ such that

$\displaystyle\left<b\pi_{b}w,v\right>_{I_{i}}=\left<bw,v\right>_{I_{i}},\quad\forall
v\in\mathbb{P}_{k}(I_{i}),$

hold for each cell $I_{i}=(x_{i-1},x_{i})$, $i=1,2,\dots,N$, where the weight
function $b(x)$ satisfies $b(x)\geq\beta^{2}>0$ for some constant $\beta$. In
the special case of $b=1$, $\pi_{b}\equiv\pi$.

  * •

Gauß-Radau projection $\pi^{\pm}$. For any $w\in H^{1}(I_{i})$,
$\pi^{\pm}w\in\mathbb{P}_{k}(I_{i})$ is defined as:

$\displaystyle\left<\pi^{\pm}w,v\right>_{I_{i}}=\left<w,v\right>_{I_{i}},\quad\forall
v\in\mathbb{P}_{k-1}(I_{i}),$
$\displaystyle(\pi^{+}w)^{+}_{i-1}=w^{+}_{i-1},\quad(\pi^{-}w)^{-}_{i}=w^{-}_{i},$

for $I_{i}=(x_{i-1},x_{i})$, $i=1,2,\ldots,N$.

The existence and uniqueness of Gauß-Radau projection is referred to A1999-An
. Following (c2021-local, , Lemma 5.1), we obtain the local stabilities and
local approximation properties in Lemma 2.

###### Lemma 2

Let $||v||_{I_{i}}=||v||_{L^{2}(I_{i})}$ and $||v||_{L^{\infty}(I_{i})}$ be
the usual $L^{2}$ and $L^{\infty}$ norms on $I_{i}$,

$\displaystyle||\pi w||_{I_{i}}\leq C||w||_{I_{i}},$
$\displaystyle||\pi_{b}w||_{I_{i}}\leq C||w||_{I_{i}},$
$\displaystyle||\pi^{-}w||_{I_{i}}\leq
C\left(||w||_{I_{i}}+h_{i}^{1/2}|w_{i}|\right),$
$\displaystyle||\pi^{+}w||_{I_{i}}\leq
C\left(||w||_{I_{i}}+h_{i}^{1/2}|w_{i-1}|\right),$
$\displaystyle||\mathbb{R}w||_{L^{\infty}(I_{i})}\leq
C||w||_{L^{\infty}(I_{i})},\quad\mathbb{R}=\pi,\pi_{b},\pi^{\pm},$
$\displaystyle||w-\mathbb{R}w||_{L^{m}(I_{i})}\leq
Ch_{i}^{k+1}||w^{(k+1)}||_{L^{m}(I_{i})},\quad
k=2,\infty,\quad\mathbb{R}=\pi,\pi_{b},\pi^{\pm},$

where $C>0$ is a constant independent of $h_{i}$ and $w$.

Now we provide the following projections:

$\displaystyle(P^{-}_{N}u)|_{I_{i}}=\left\\{\begin{aligned}
&\pi^{-}u,\quad&&\text{for $i=1,2,\ldots,N/4,3N/4+1,\ldots,N$},\\\
&\pi_{b}u,\quad&&\text{for $i=N/4+1,\ldots,3N/4$},\\\ \end{aligned}\right.$
$\displaystyle(P^{+}_{N}q)|_{I_{i}}=\left\\{\begin{aligned}
&\pi^{+}q,\quad&&\text{for $i=2,3,\cdots,N$},\\\ &\pi q,\quad&&\text{for
$i=1$}.\\\ \end{aligned}\right.$

Assume that $||v||^{2}=\sum\limits_{i=1}^{N}||v||_{I_{i}}^{2}$ and
$||v||_{I_{i}}^{2}=\left<v,v\right>_{I_{i}}$.

###### Lemma 3

Let $\epsilon\leq CN^{-1}$ and $\sigma\geq k+1$. Then on the layer-adapted
mesh (6), there is

(7) $\displaystyle\qquad||u-P^{-}_{N}u||_{[0,x_{N/4}]\cup[x_{3N/4},1]}\leq
C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1},$ (8)
$\displaystyle\qquad||q-P^{+}_{N}q||\leq C\epsilon^{\frac{3}{4}}(N^{-1}\ln
N)^{k+1},$ (9) $\displaystyle\qquad||u-P^{-}_{N}u||_{L^{\infty}(I_{i})}\leq
CN^{-(k+1)},\quad i=N/4+1,\cdots,3N/4,$ (10)
$\displaystyle\qquad||q-P^{+}_{N}q||_{L^{\infty}(I_{i})}\leq
C\epsilon^{\frac{1}{2}}(N^{-1}\ln N)^{k+1},\quad
i=1,\cdots,N/4,3N/4+1,\cdots,N,$ (11)
$\displaystyle\qquad||q-P^{+}_{N}q||_{L^{\infty}(I_{i})}\leq C(\epsilon
N^{-(k+1)}+\epsilon^{\frac{1}{2}}N^{-\sigma}),\quad i=N/4+1,\cdots,3N/4.$

###### Proof

According to the arguments in c2022-balanced  and Lemma 2, this lemma can be
derived without any difficulties.

###  2.5 Convergence analysis

Recall that ${w}=(q,u)$ is the solution of eqrefeq1.1. According to (5), we
define a natural norm by

(12)
$|||w|||^{2}_{E}=B(w;w)=\epsilon^{-1}||q||^{2}+||b(x)^{\frac{1}{2}}u||^{2}+\sum_{j\in\\{0,N\\}}\lambda_{j}[[u]]_{j}^{2}+\lambda_{q}[[q]]^{2}_{\frac{3}{4}N}.$

Thena more powerful balanced norm is introduced, which is denoted as

(13)
$|||w|||^{2}_{B}=\epsilon^{-3/2}||q||^{2}+||b(x)^{\frac{1}{2}}u||^{2}+\sum_{j\in\\{0,N\\}}[[u]]_{j}^{2}+\epsilon^{-1}[[q]]^{2}_{\frac{3}{4}N}.$

###### Theorem 2.1

Let $\epsilon\leq CN^{-1}$ and $\sigma\geq k+1$ on the mesh (6). Suppose that
${W}=(Q,U)\in V_{N}\times V_{N}$ is the solution of (2) and (3). And suppose
that ${w}=(q,u)$ is the solution of (1) meeting Lemma 1. Then define ${e=w-W}$
and

$\displaystyle|||e|||_{E}\leq C\left(\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+\epsilon^{\frac{1}{2}}N^{-k}+N^{-(k+1)}\right),$
$\displaystyle|||e|||_{B}\leq C\left((N^{-1}\ln
N)^{k+1}+\epsilon^{\frac{1}{4}}N^{-k}\right),$

where $C>0$ is a constant.

###### Proof

First, we divide the error $e$ into ${e}=(q-Q,u-U)=\mathbb{\eta-\chi}$ with

$\displaystyle\mathbb{\eta}=(\eta_{q},\eta_{u})=(q-P^{+}_{N}q,u-P^{-}_{N}u),$
$\displaystyle\mathbb{\chi}=(\chi_{q},\chi_{u})=(Q-P^{+}_{N}q,U-P^{-}_{N}u)\in
V_{N}\times V_{N}.$

From Lemma 1, the consistency of numerical fluxes and (5), we have

(14) $|||\chi|||_{E}^{2}=B(\chi;\chi)=B(\eta;\chi)\equiv\sum_{j=1}^{8}S_{j}.$

Here,

$\displaystyle S_{1}=\epsilon^{-1}\left<\eta_{q},\chi_{q}\right>,\quad$
$\displaystyle S_{2}=\left<\eta_{u},\chi^{\prime}_{q}\right>,$ $\displaystyle
S_{3}=\left<\eta_{q},\chi^{\prime}_{u}\right>,\quad$ $\displaystyle
S_{4}=\left<b\eta_{u},\chi_{u}\right>,$ $\displaystyle
S_{5}=-\sum_{i=0}^{N-1}(\eta_{q})^{+}_{i}[[\chi_{u}]]_{i}-(\eta_{q})^{-}_{N}(\chi_{u})^{-}_{N},\quad$
$\displaystyle S_{6}=-\sum_{i=1}^{N-1}(\eta_{u})^{-}_{i}[[\chi_{q}]]_{i},$
$\displaystyle
S_{7}=\lambda_{q}[[\eta_{q}]]_{\frac{3}{4}N}[[\chi_{q}]]_{\frac{3}{4}N},\quad$
$\displaystyle
S_{8}=\sum_{i\in\\{0,N\\}}\lambda_{i}[[\eta_{u}]]_{i}[[\chi_{u}]]_{i}.$

Below, we will estimate them in sequence.

From Hölder inequalities and (8),

$\displaystyle S_{1}\leq
C(\epsilon^{-\frac{1}{2}}||\eta_{q}||)(\epsilon^{-\frac{1}{2}}||\chi_{q}||)\leq
C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1}|||\chi|||_{E}.$

Then from the orthogonality of the projection $P_{N}^{-}u$ and $P_{N}^{+}q$,
we have $S_{2}=S_{3}=0$.

For $S_{4}$, we divide it into the following two parts,

$S_{4}=\left<b\eta_{u},\chi_{u}\right>_{[0,x_{N/4}]\cup[x_{3N/4},1]}+\left<b\eta_{u},\chi_{u}\right>_{[x_{N/4},x_{3N/4}]},$

where from Hölder inequalities and (7),

$\displaystyle\left<b\eta_{u},\chi_{u}\right>_{[0,x_{N/4}]\cup[x_{3N/4},1]}$
$\displaystyle\leq C||\eta_{u}||_{[0,x_{N/4}]\cup[x_{3N/4},1]}|||\chi|||_{E}$
$\displaystyle\leq C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1}|||\chi|||_{E},$

and the definition of $\pi_{b}$,
$\left<b\eta_{u},\chi_{u}\right>_{[x_{N/4},x_{3N/4}]}=0$.

From Hölder inequalities, (10) and
$\lambda_{0}=\lambda_{N}=\epsilon^{\frac{1}{2}}$,

$\displaystyle S_{5}$
$\displaystyle=(\eta_{q})^{+}_{0}(\chi_{u})^{+}_{0}-(\eta_{q})^{-}_{N}(\chi_{u})^{-}_{N}$
$\displaystyle\leq
C\left(\sum_{j\in\\{0,N\\}}\lambda_{i}^{-1}[[\eta_{q}]]^{2}_{i}\right)^{1/2}\left(\sum_{j\in\\{0,N\\}}\lambda_{i}[[\chi_{u}]]^{2}_{i}\right)^{1/2}$
$\displaystyle\leq
C\epsilon^{-\frac{1}{4}}||\eta_{q}||_{L^{\infty}(\Omega)}|||\chi|||_{E}$
$\displaystyle\leq C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1}|||\chi|||_{E}.$

From Hölder inequalities, the trace inequality, (9),
$\lambda_{q}=\epsilon^{-\frac{1}{2}}$ and $\epsilon\leq CN^{-1}$, one has

$\displaystyle S_{6}$
$\displaystyle=|-\sum_{i=N/4+1}^{3N/4-1}(\eta_{u})^{-}_{i}[[\chi_{q}]]_{i}-(\eta_{u})^{-}_{3N/4}[[\chi_{q}]]_{\frac{3}{4}N}|$
$\displaystyle\leq
C\sum_{i=N/4+1}^{3N/4-1}||\eta_{u}||_{L^{\infty}(I_{i})}h_{i}^{\frac{1}{2}}||\chi_{q}||_{I_{i}\cup
I_{i+1}}+C\lambda_{q}^{-\frac{1}{2}}||\eta_{u}||_{L^{\infty}(I_{3N/4})}(\sqrt{\lambda_{q}}[[\chi_{q}]]_{\frac{3}{4}N})$
$\displaystyle\leq
C\left(\epsilon\sum_{i=N/4+1}^{3N/4-1}h_{i}^{-1}||\eta_{u}||^{2}_{L^{\infty}(I_{i})}\right)^{\frac{1}{2}}|||\chi|||_{E}+C\epsilon^{\frac{1}{4}}N^{-(k+1)}|||\chi|||_{E}$
$\displaystyle\leq
C\left(\epsilon^{\frac{1}{2}}N^{-k}+C\epsilon^{\frac{1}{4}}N^{-(k+1)}\right)|||\chi|||_{E}.$

According to (10), (11) and $\lambda_{q}=\epsilon^{-\frac{1}{2}}$, we obatin

$\displaystyle S_{7}\leq
C\left(\lambda_{q}[[\eta_{q}]]^{2}_{\frac{3}{4}N}\right)^{\frac{1}{2}}\left(\lambda_{q}[[\chi_{q}]]^{2}_{\frac{3}{4}N}\right)^{\frac{1}{2}}\leq
C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1}|||\chi|||_{E}.$

From Hölder inequalities, (9), (10) and
$\lambda_{0}=\lambda_{N}=\epsilon^{\frac{1}{2}}$, one obtains

$\displaystyle S_{8}$
$\displaystyle=\sum_{i\in\\{0,N\\}}\lambda_{i}[[\eta_{u}]]_{i}[[\chi_{u}]]_{i}$
$\displaystyle\leq
C\left(\sum_{i\in\\{0,N\\}}\lambda_{i}[[\eta_{u}]]^{2}_{i}\right)^{\frac{1}{2}}\left(\sum_{i\in\\{0,N\\}}\lambda_{i}[[\chi_{u}]]^{2}_{i}\right)^{\frac{1}{2}}$
$\displaystyle\leq
C\epsilon^{\frac{1}{4}}||\eta_{u}||_{L^{\infty}(\Omega)}|||\chi|||_{E}$
$\displaystyle\leq C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1}|||\chi|||_{E}.$

By gathering all the above estimates, we get

$\displaystyle B(\eta;\chi)\leq C\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}|||\chi|||_{E}.$

Using (14), there is

$|||\chi|||^{2}_{E}=B(\chi;\chi)=B(\eta;\chi)\leq
C\left(\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+\epsilon^{\frac{1}{2}}N^{-k}\right)|||\chi|||_{E}.$

Therefore, it is straightforward to derive

(15) $\displaystyle|||\chi|||_{E}\leq C\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+C\epsilon^{\frac{1}{2}}N^{-k},$ (16) $\displaystyle|||\chi|||_{B}\leq
C\epsilon^{-\frac{1}{4}}|||\chi|||_{E}\leq C(N^{-1}\ln
N)^{k+1}+C\epsilon^{\frac{1}{4}}N^{-k}.$

Furthermore, from Lemma 3, 12 and 13, there are

(17) $\displaystyle|||\eta|||_{E}\leq C\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+CN^{-(k+1)},$ (18) $\displaystyle|||\eta|||_{B}\leq C(N^{-1}\ln
N)^{k+1}.$

Combining (15), (16), (17) and (18), Theorem 2.1 can be obtained.

###### Remark 3

For Theorem 2.1, when we choose $\epsilon\leq CN^{-4}$,

$|||e|||_{E}\leq C\left(N^{-(k+2)}(\ln
N)^{k+1}+N^{-(k+1)}\right),\quad|||e|||_{B}\leq C(N^{-1}\ln N)^{k+1}.$

Unlike c2022-balanced , the optimal convergence order can be proved without
the condition that the smooth component of the solution belongs to the finite
element space.

##  3 LDG method for two-dimensional case

We consider a two-dimensional singularly perturbed problem:

(19) $\displaystyle-\epsilon\Delta u+b(x,y)u=f(x,y),\quad\text{in
$\Omega=(0,1)^{2}$},$ $\displaystyle u=0,\quad\text{on $\partial\Omega$}.$

where $0<\epsilon\ll 1$ is a perturbation parameter, $b(x,y)\geqslant
2\beta^{2}>0$ and $f(x,y)$ are sufficient smooth functions on
$\overline{\Omega}$. It is widely known that there is a unique solution to
this problem, characterized by the presence of a boundary layer with a width
of $\mathcal{O}(\sqrt{\epsilon}\ln(1/\epsilon))$ along $\partial\Omega$.

###  3.1 LDG method

The two-dimensional Shishkin mesh $\\{(x_{i},y_{j}),i,j=1,2,\cdots,N\\}$ is
constructed by the tensor product of the Shishkin mesh in 1D (6) in both
vertical and horizontal directions. Here $x_{i}$ is presented in (6) and
$y_{j}$ is defined in a same way.

Then for $i,j=1,\cdots,N$, we define $I_{i,x}=(x_{i-1},x_{i})$ and
$J_{j,y}=(y_{j-1},y_{j})$. Let $h_{i,x}=x_{i}-x_{i-1}$ and
$h_{j,y}=y_{j}-y_{j-1}$. In addition, suppose that

$\Omega_{N}=\\{\text{$\kappa_{ij}$: $\kappa_{ij}=I_{i,x}\times J_{j,y}$ for
$i,j=1,2,3,\cdots,N$}\\}$

is a rectangle partition of the domain $\Omega$. We define the $k$-degree
discontinuous finite element space by:

$R_{N}=\\{v\in
L^{2}(\Omega):\text{$v|_{\kappa_{ij}}\in\mathbb{Q}_{k}(\kappa_{ij})$
},\kappa_{ij}\in\Omega_{N}\\},$

where $\mathbb{Q}_{k}(\kappa_{ij})$ is the space of polynomials on
$\kappa_{ij}$, with a maximum degree of $k\geq 1$ in each variable. For $v\in
R_{N}$ and $y\in J_{j,y}$, $j=1,2,\ldots,N$, $v^{\pm}_{i,y}=\lim\limits_{x\to
x^{\pm}_{i}}v(x,y)$ and $v^{\pm}_{x,j}=\lim\limits_{y\to y^{\pm}_{j}}v(x,y)$
are used to express the traces evaluated from the four directions. Then we
define the jumps as:

$\displaystyle[[v]]_{i,y}=v^{-}_{i,y}-v^{+}_{i,y}\text{ for
$i=1,2,\cdots,N-1$},\;$ $\displaystyle[[v]]_{0,y}=-v^{+}_{0,y},\;$
$\displaystyle[[v]]_{N,y}=v^{-}_{N,y}.$
$\displaystyle[[v]]_{x,j}=v^{-}_{x,j}-v^{+}_{x,j}\text{ for
$j=1,2,\cdots,N-1$},\;$ $\displaystyle[[v]]_{x,0}=-v^{+}_{x,0},\;$
$\displaystyle[[v]]_{x,N}=v^{-}_{x,N}.$

Below, we shall introduce the LDG method for (19). Firstly, we rewrite the
problem (19) as follows:

$\displaystyle\epsilon^{-1}p=u_{x},\quad$ $\displaystyle\text{in $\Omega$},$
$\displaystyle\epsilon^{-1}q=u_{y},\quad$ $\displaystyle\text{in $\Omega$},$
$\displaystyle-p_{x}-q_{y}+bu=f,\quad$ $\displaystyle\text{in $\Omega$}.$

Suppose that $Z=(v,s,r)\in R_{N}\times R_{N}\times R_{N}$ is any test
function, and $(\cdot,\cdot)_{D}$ is the $L^{2}$ inner product in
$D\subset\mathbb{R}^{2}$. Suppose that $\left<\cdot,\cdot\right>_{I}$ is the
inner product in $L^{2}(I)$ with $I\subset\mathbb{R}$. Find $T=(U,P,Q)\in
R_{N}\times R_{N}\times R_{N}$ such that the variational formulas hold in
$\kappa_{ij}$,

(20)
$\displaystyle\epsilon^{-1}(P,s)_{\kappa_{ij}}+(U,s_{x})_{\kappa_{ij}}-\left<\widehat{U}_{i,y},s^{-}_{i,y}\right>_{J_{j,y}}+\left<\widehat{U}_{i-1,y},s^{+}_{i-1,y}\right>_{J_{j,y}}=0,$
$\displaystyle\epsilon^{-1}(Q,r)_{\kappa_{ij}}+(U,r_{y})_{\kappa_{ij}}-\left<\widehat{U}_{x,j},r^{-}_{x,j}\right>_{I_{i,x}}+\left<\widehat{U}_{x,j-1},r^{+}_{x,j-1}\right>_{I_{i,x}}=0,$
$\displaystyle(P,v_{x})_{\kappa_{ij}}-\left<\widehat{P}_{i,y},v^{-}_{i,y}\right>_{J_{j,y}}+\left<\widehat{P}_{i-1,y},v^{+}_{i-1,y}\right>_{J_{j,y}}+(Q,v_{y})_{\kappa_{ij}}-\left<\widehat{Q}_{x,j},v^{-}_{x,j}\right>_{I_{i,x}}$
$\displaystyle+\left<\widehat{Q}_{x,j-1},v^{+}_{x,j-1}\right>_{I_{i,x}}+(bU,v)_{\kappa_{ij}}=(f,v)_{\kappa_{ij}}.$

For $y\in J_{j,y},j=1,2,\ldots,N$, the numerical fluxes are defined by

$\displaystyle\widehat{P}_{i,y}=\left\\{\begin{aligned}
&P^{+}_{0,y}+\lambda_{0,y}U^{+}_{0,y},\quad&&\text{$i=0$},\\\
&P^{+}_{i,y},\quad&&\text{$i=1,2,\cdots,N-1$},\\\
&P^{-}_{N,y}-\lambda_{N,y}U^{-}_{N,y},\quad&&\text{$i=N$},\end{aligned}\right.$
$\displaystyle\widehat{U}_{i,y}=\left\\{\begin{aligned}
&0,\quad&&\text{$i=0,N$},\\\
&U_{{\frac{3}{4}N},y}^{-}-\lambda_{P}[[P]]_{{\frac{3}{4}N},y},\quad&&\text{$i=\frac{3}{4}N$},\\\
&U^{-}_{i,y},\quad&&\text{$i=1,2,\cdots,\frac{3}{4}N-1,\frac{3}{4}N+1,\ldots,N-1$},\\\
\end{aligned}\right.$

where we take $\lambda_{0,y}=\lambda_{N,y}=\epsilon^{\frac{1}{2}}$ and
$\lambda_{P}=\epsilon^{-\frac{1}{2}}$. Similarly, when $x\in
I_{i,x}(i=1,2,\cdots,N)$, we define $\widehat{U}_{x,j}$ and
$\widehat{Q}_{x,j}$ for $j=1,2,\ldots,N$, where
$\widehat{U}_{x,\frac{3}{4}N}=U_{x,{\frac{3}{4}N}}^{-}-\lambda_{Q}[[Q]]_{x,{\frac{3}{4}N}}$
and $\lambda_{Q}=\epsilon^{-\frac{1}{2}}$.

Then the scheme (20) by a compact form can be expressed as: Find $T=(U,P,Q)\in
R_{N}\times R_{N}\times R_{N}$ to satisfy

$B(T;Z)=(f,v),\quad\forall Z=(v,s,r)\in R_{N}\times R_{N}\times R_{N},$

where

(21) $\displaystyle B(T;Z)=$
$\displaystyle(bU,v)+\epsilon^{-1}(P,s)+\epsilon^{-1}(Q,r)$
$\displaystyle+(U,s_{x})-\sum_{j=1}^{N}\sum_{i=1}^{N-1}\left<U^{-}_{i,y},[[s]]_{i,y}\right>_{J_{j,y}}+(U,r_{y})-\sum_{i=1}^{N}\sum_{j=1}^{N-1}\left<U^{-}_{x,j},[[r]]_{x,j}\right>_{I_{i,x}}$
$\displaystyle+(P,v_{x})-\sum_{j=1}^{N}\left\\{\sum_{i=0}^{N-1}\left<P^{+}_{i,y},[[v]]_{i,y}\right>_{J_{j,y}}+\left<P^{-}_{N,y},[[v]]_{N,y}\right>_{J_{j,y}}\right\\}$
$\displaystyle+(Q,v_{y})-\sum_{i=1}^{N}\left\\{\sum_{j=0}^{N-1}\left<Q^{+}_{x,j},[[v]]_{x,j}\right>_{I_{i,x}}+\left<Q^{-}_{x,N},[[v]]_{x,N}\right>_{I_{i,x}}\right\\}$
$\displaystyle+\sum_{j=1}^{N}\sum_{i\in\\{0,N\\}}\left<\lambda_{i,y}[[U]]_{i,y},[[v]]_{i,y}\right>_{J_{j,y}}+\sum_{i=1}^{N}\sum_{j\in\\{0,N\\}}\left<\lambda_{x,j}[[U]]_{x,j},[[v]]_{x,j}\right>_{I_{i,x}}$
$\displaystyle+\sum_{j=1}^{N}\left<\lambda_{P}[[P]]_{\frac{3}{4}N,y},[[s]]_{\frac{3}{4}N,y}\right>_{J_{j,y}}+\sum_{i=1}^{N}\left<\lambda_{Q}[[Q]]_{x,\frac{3}{4}N},[[r]]_{x,\frac{3}{4}N}\right>_{I_{i,x}}.$

###  3.2 Regularity of the solution

###### Lemma 4

(c2022-balanced ; Cl1Gr2:2005-motified ) Suppose that $f$ satisfies sufficient
compatibility conditions and $b$ is sufficiently smooth. The solution $u$ of
(19) admits the decomposition

$\displaystyle
u=S+\sum_{t=1}^{4}W_{t}+\sum_{t=1}^{4}Z_{t},\quad\text{$(x,y)\in\overline{\Omega}$},$

where $S$ is a smooth component, $W_{t}$ is a boundary layer component and
$Z_{t}$ is a corner layer component. Furthermore, for $0\leq i+j\leq k+2$, one
has

$\displaystyle|\frac{\partial^{i+j}S}{\partial x^{i}\partial y^{j}}|\leq C,$
$\displaystyle|\frac{\partial^{i+j}W_{1}}{\partial x^{i}\partial y^{j}}|\leq
C\epsilon^{-\frac{i}{2}}e^{-\frac{\beta x}{\sqrt{\epsilon}}},$
$\displaystyle|$ $\displaystyle\frac{\partial^{i+j}Z_{1}}{\partial
x^{i}\partial y^{j}}|\leq
C\epsilon^{-\frac{i+j}{2}}e^{-\frac{\beta(x+y)}{\sqrt{\epsilon}}},$

and so on for the remaining terms.

###  3.3 Projections

For the later analytical proof, we construct a new projection by combining the
local (or the weighted local) $L^{2}$ projection and the Gauß-Radau projection
m2021-local .

  * •

Local $L^{2}$ projection $\Pi$. For each cell $\kappa_{ij}\in\Omega_{N}$ and
$z\in L^{2}({\Omega})$, $\Pi z\in R_{N}$ is defined as:

$\displaystyle(\Pi z,v)_{\kappa_{ij}}=(z,v)_{\kappa_{ij}},\quad\forall
v\in\mathbb{Q}_{k}(\kappa_{ij}).$

  * •

Weighted local $L^{2}$ projection $\Pi_{b}$. Let $w\in C^{1}(\Omega_{N})$ and
$w\geq w_{0}>0$ be a general weight function. The weighted local $L^{2}$
projection $\Pi_{w}$ is defined as: For each cell $\kappa_{ij}\in\Omega_{N}$
and $z\in L^{2}(\Omega)$, $\Pi_{w}z\in R_{N}$ meets

$(w\Pi_{w}z,v)_{\kappa_{ij}}=(wz,v)_{\kappa_{ij}},\quad\forall
v\in\mathbb{Q}_{k}(\kappa_{ij})$

When $w=1$, the operator shall reduce to the local $L^{2}$ projection $\Pi$.
In this paper, we take $b(x)$ as a weight function and $b(x)\geq
2\beta^{2}>0$.

  * •

Gauß-Radau projection $\Pi^{\pm}$. First, we introduce the following notation:
For $l$, $m\geq 1$

$\mathbb{Q}_{l,m}=\left\\{\text{$\sum_{i=0}^{l}\sum_{j=0}^{m}a_{ij}x^{i}y^{j}$:
$a_{ij}\in\mathbb{R}$}\right\\}.$

For any $\kappa_{ij}\in\Omega_{N}$ and $z\in H^{1}(\kappa_{ij})$, we define
$\Pi^{-}_{x}z,\Pi^{-}_{y}z,\Pi^{+}_{x}z,\Pi^{+}_{y}z\in\mathbb{Q}_{k}(\kappa_{ij})$
as:

$\displaystyle\left\\{\begin{aligned}
&(\Pi^{-}_{x}z,w)_{\kappa_{ij}}=(z,w)_{\kappa_{ij}},\quad&&\text{$\forall
w\in\mathbb{Q}_{k-1,k}$},\\\
&\left<(\Pi^{-}_{x}z)^{-}_{i,y},w\right>_{J_{j,y}}=\left<z^{-}_{i,y},w\right>_{J_{j,y}},\quad&&\forall
w\in\mathbb{P}_{k}(J_{j,y}),\end{aligned}\right.$
$\displaystyle\left\\{\begin{aligned}
&(\Pi^{-}_{y}z,w)_{\kappa_{ij}}=(z,w)_{\kappa_{ij}},\quad&&\text{$\forall
w\in\mathbb{Q}_{k,k-1}$},\\\
&\left<(\Pi^{-}_{y}z)^{-}_{x,j},w\right>_{I_{i,x}}=\left<z^{-}_{x,j},w\right>_{I_{i,x}},\quad&&\forall
w\in\mathbb{P}_{k}(I_{i,x}),\end{aligned}\right.$
$\displaystyle\left\\{\begin{aligned}
&(\Pi^{+}_{x}z,w)_{\kappa_{ij}}=(z,w)_{\kappa_{ij}},\quad&&\text{$\forall
w\in\mathbb{Q}_{k-1,k}$},\\\
&\left<(\Pi^{+}_{x}z)^{+}_{i,y},w\right>_{J_{j,y}}=\left<z^{+}_{i,y},w\right>_{J_{j,y}},\quad&&\forall
w\in\mathbb{P}_{k}(J_{j,y}),\end{aligned}\right.$
$\displaystyle\left\\{\begin{aligned}
&(\Pi^{+}_{y}z,w)_{\kappa_{ij}}=(z,w)_{\kappa_{ij}},\quad&&\text{$\forall
w\in\mathbb{Q}_{k,k-1}$},\\\
&\left<(\Pi^{+}_{y}z)^{+}_{x,j},w\right>_{I_{i,x}}=\left<z^{+}_{x,j},w\right>_{I_{i,x}},\quad&&\forall
w\in\mathbb{P}_{k}(I_{i,x}).\end{aligned}\right.$

The existence and uniqueness of Gauß-Radau projection is referred to A1999-An
, and the corresponding approximate estimation refers to (m2021-local, , Lemma
4.3).

Now we design the following projection:

$\displaystyle(P^{-}u)|_{\kappa_{ij}}=\left\\{\begin{aligned}
&\Pi^{-}_{x}u,\quad&&\text{$i=1,\cdots,N/4,3N/4+1,\cdots,N-1,\quad
j=N/4+1,\cdots,3N/4$},\\\
&\Pi^{-}_{y}u,\quad&&\text{$i=N/4+1,\cdots,3N/4,\quad
j=1,\cdots,N/4,3N/4+1,\cdots,N-1$},\\\ &\Pi_{b}u,\quad&&\text{other else},\\\
\end{aligned}\right.$
$\displaystyle(P^{+}_{x}p)|_{\kappa_{ij}}=\left\\{\begin{aligned} &\Pi
p,\quad&&\text{$i=1,\quad j=1,\cdots,N$},\\\ &\Pi^{+}_{x}p,\quad&&\text{other
else},\\\ \end{aligned}\right.$
$\displaystyle(P^{+}_{y}q)|_{\kappa_{ij}}=\left\\{\begin{aligned} &\Pi
q,\quad&&\text{$i=1,\cdots,N,\quad j=1$},\\\ &\Pi^{+}_{y}q,\quad&&\text{other
else}.\\\ \end{aligned}\right.$

Let
$||z||^{2}=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}||z||_{\kappa_{ij}}^{2}$
and $||z||^{2}_{\kappa_{ij}}=(z,z)_{\kappa_{ij}}$. And we will give the
following estimates for $u$ and $p$, and the error estimates for $q$ can be
derived in a similar way.

###### Lemma 5

(m2021-local ) Let $\epsilon\leq CN^{-1}$ and $\sigma\geq k+1$, then on the
two-dimensional Shishkin mesh, one has

(22) $\displaystyle||u-P^{-}u||_{\Omega/\Omega_{22}}\leq
C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1},$ (23)
$\displaystyle||p-P^{+}_{x}p||\leq C\epsilon^{\frac{3}{4}}(N^{-1}\ln
N)^{k+1}+C\epsilon N^{-(k+1)}(\ln N)^{k+\frac{3}{2}},$ (24)
$\displaystyle||u-P^{-}u||_{L^{\infty}(\Omega)}\leq C(N^{-1}\ln N)^{k+1},$
(25) $\displaystyle||p-P^{+}_{x}p||_{L^{\infty}(\Omega)}\leq
C\epsilon^{\frac{1}{2}}(N^{-1}\ln N)^{k+1},$

where $\Omega_{22}:=[x_{N/4},x_{3N/4}]\times[y_{N/4},y_{3N/4}]$.

###### Proof

According to the arguments in m2021-local , the proofs of (22), (23), (24) and
(25) can be obtained easily.

###  3.4 Convergence analysis

Recall that $t=(u,p,q)$ is thesolution of (19). From (21), an natural norm is
denoed as:

(26) $\displaystyle|||t|||^{2}_{E}=B(t;t)$
$\displaystyle=\epsilon^{-1}||p||^{2}+\epsilon^{-1}||q||^{2}+||b(x)^{\frac{1}{2}}u||^{2}$
$\displaystyle+\sum^{N}_{j=1}\sum_{i\in\\{0,N\\}}(\lambda_{i,y},[[u]]_{i,y}^{2})_{J_{j,y}}+\sum^{N}_{i=1}\sum_{j\in\\{0,N\\}}(\lambda_{x,j},[[u]]_{x,j}^{2})_{I_{i,x}}$
$\displaystyle+\sum_{j=1}^{N}(\lambda_{P},[[p]]^{2}_{\frac{3}{4}N,y})_{J_{j,y}}+\sum_{i=1}^{N}(\lambda_{Q},[[q]]^{2}_{x,\frac{3}{4}N})_{I_{i,x}}.$

Then we denote the balanced norm by:

(27) $\displaystyle|||t|||^{2}_{B}$
$\displaystyle=\epsilon^{-3/2}||p||^{2}+\epsilon^{-3/2}||q||^{2}+||b(x)^{\frac{1}{2}}u||^{2}+\sum^{N}_{j=1}\sum_{i\in\\{0,N\\}}(1,[[u]]_{i,y}^{2})_{J_{j,y}}$
$\displaystyle+\sum^{N}_{i=1}\sum_{j\in\\{0,N\\}}(1,[[u]]_{x,j}^{2})_{I_{i,x}}+\sum_{j=1}^{N}(\epsilon^{-1},[[p]]^{2}_{\frac{3}{4}N,y})_{J_{j,y}}+\sum_{i=1}^{N}(\epsilon^{-1},[[q]]^{2}_{x,\frac{3}{4}N})_{I_{i,x}}.$

###### Theorem 3.1

Suppose that $\epsilon\leq CN^{-1}$ and $\sigma\geq k+1$ on the layer-adapted
mesh $\\{(x_{i},y_{j}),i,j=1,2,\ldots,N\\}$. Let $T=(U,P,Q)\in R_{N}\times
R_{N}\times R_{N}$ be the solution of (20), and $t=(u,p,q)$ is the solution of
(19). Moreover, we have

$\displaystyle|||e|||_{E}\leq C\left(\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+\epsilon^{\frac{1}{2}}N^{-k}(\ln
N)^{k+1}+\epsilon^{\frac{3}{4}}N^{-k}(\ln
N)^{k+\frac{3}{2}}+N^{-(k+1)}\right),$ $\displaystyle|||e|||_{B}\leq
C\left((N^{-1}\ln N)^{k+1}+\epsilon^{\frac{1}{4}}N^{-k}(\ln
N)^{k+1}+\epsilon^{\frac{1}{2}}N^{-k}(\ln N)^{k+\frac{3}{2}}\right),$

where we define ${e=t-T}$ and $C>0$ is a constant.

###### Proof

First, by using projection $P^{-}u$, $P^{+}_{x}p$ and $P^{+}_{y}q$, the error
is devided as ${e=t-T}={\eta-\chi}$ with

$\displaystyle\mathbb{\eta}=(\eta_{u},\eta_{p},\eta_{q})=(u-P^{-}u,p-P^{+}_{x}p,q-P^{+}_{y}q),$
$\displaystyle\mathbb{\chi}=(\chi_{u},\chi_{p},\chi_{q})=(U-P^{-}u,P-P^{+}_{x}p,Q-P^{+}_{y}q)\in
R_{N}\times R_{N}\times R_{N}.$

According to Lemma 4 and the consistency of numerical fluxes, there is

(28) $B(\chi;Z)=B(\eta;Z),\quad\forall Z=(u,p,q)\in R_{N}\times R_{N}\times
R_{N}.$

Set $Z=\chi$ in (28) and use (21), we derive

(29) $|||\chi|||_{E}^{2}=B(\chi;\chi)=B(\eta;\chi)\equiv\sum_{i=1}^{15}Y_{i},$

where

$\displaystyle Y_{1}=(b\eta_{u},\chi_{u}),\quad
Y_{2}=\epsilon^{-1}(\eta_{p},\chi_{p}),$ $\displaystyle
Y_{3}=(\eta_{u},\chi_{p,x}),\quad Y_{4}=(\eta_{p},\chi_{u,x}),$ $\displaystyle
Y_{5}=-\sum_{j=1}^{N}\sum_{i=1}^{N-1}\left<(\eta_{u})^{-}_{i,y},[[\chi_{p}]]_{i,y}\right>_{J_{j,y}},$
$\displaystyle
Y_{6}=-\sum_{j=1}^{N}\left\\{\sum_{i=0}^{N-1}\left<(\eta_{p})^{+}_{i,y},[[\chi_{u}]]_{i,y}\right>_{J_{j,y}}+\left<(\eta_{p})^{-}_{N,y},[[\chi_{u}]]_{N,y}\right>_{J_{j,y}}\right\\},$
$\displaystyle
Y_{7}=\sum_{j=1}^{N}\sum_{i\in\\{0,N\\}}\left<\lambda_{i,y}[[\eta_{u}]]_{i,y},[[\chi_{u}]]_{i,y}\right>_{J_{j,y}},$
$\displaystyle
Y_{8}=\sum_{j=1}^{N}\left<\lambda_{P}[[\eta_{p}]]_{\frac{3}{4}N,y},[[\chi_{p}]]_{\frac{3}{4}N,y}\right>_{J_{j,y}}.$

Because the proof of $q$ is similar to $p$, we just estimate the part of $p$.
Next, we will analyze the above parts in sequence.

By applying Hölder inequalities, the definition of the weighted local $L^{2}$
projection, (22) and (23),

$\displaystyle Y_{1}\leq
C||b^{\frac{1}{2}}\eta_{u}||_{\Omega/\Omega_{22}}||b^{\frac{1}{2}}\chi_{u}||\leq
C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1}|||\chi|||_{E}.$ $\displaystyle
Y_{2}\leq(\epsilon^{-1/2}||\eta_{p}||)(\epsilon^{-1/2}||\chi_{p}||)\leq
C\left(\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+\epsilon^{\frac{1}{2}}N^{-(k+1)}(\ln
N)^{k+\frac{3}{2}}\right)|||\chi|||_{E}.$

Here $\Omega_{22}:=[x_{N/4},x_{3N/4}]\times[y_{N/4},y_{3N/4}]$.

Furthermore, from the orthogonality of the projection $P^{-}u$ and
$P^{+}_{x}p$, we have $Y_{3}=Y_{4}=0$.

Let
$\Omega_{t}=\Omega-\left\\{\left([0,x_{N/4}]\times[y_{N/4},y_{3N/4}]\right)\cup\left([x_{3N/4},x_{N-1}]\times[y_{N/4},y_{3N/4}]\right)\right\\}$.
According to the trace inequality, (24), $\epsilon\leq CN^{-1}$ and
$\lambda_{P}=\epsilon^{-\frac{1}{2}}$, we get

$\displaystyle|Y_{5}|$ $\displaystyle\leq
C\sum_{\kappa_{ij}\in\Omega_{t},i\neq
3N/4}h_{j,y}^{\frac{1}{2}}||\eta_{u}||_{L^{\infty}(\kappa_{ij})}||\chi_{p}||_{L^{\infty}((I_{i}\cup
I_{i+1})\times
J_{j,y})}+\sum_{j=1}^{N}\left<(\eta_{u})^{-}_{3N/4,y},[[\chi_{p}]]_{3N/4,y}\right>_{J_{j,y}}$
$\displaystyle\leq C\sum_{\kappa_{ij}\in\Omega_{t},i\neq
3N/4}h_{j,y}^{\frac{1}{2}}h_{i,x}^{-\frac{1}{2}}||\eta_{u}||_{L^{\infty}(\kappa_{ij})}||\chi_{p}||+C\left(\sum_{j=1}^{N}\lambda_{P}^{-1}h_{j,y}||\eta_{u}||_{L^{\infty}(\kappa_{ij})}^{2}\right)^{\frac{1}{2}}|||\chi|||_{E}$
$\displaystyle\leq C\left(\epsilon^{\frac{1}{2}}N^{-k}(\ln
N)^{k+1}+\epsilon^{\frac{3}{4}}N^{-k}(\ln
N)^{k+\frac{3}{2}}+\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}\right)|||\chi|||_{E}.$

According to Hölder inequalities, (25), (24) and the condition
$\lambda_{0,y}=\lambda_{N,y}=\epsilon^{\frac{1}{2}}$,

$\displaystyle|Y_{6}|$
$\displaystyle=|-\sum_{j=1}^{N}\left(\left<(\eta_{p})_{0,y}^{+},[[\chi_{u}]]_{0,y}\right>_{J_{j,y}}+\left<(\eta_{p})_{N,y}^{-},[[\chi_{u}]]_{N,y}\right>_{J_{j,y}}\right)|$
$\displaystyle\leq\left(\sum_{j=1}^{N}\sum_{i\in\\{0,N\\}}\lambda_{i,y}^{-1}h_{j}||\eta_{p}||^{2}_{L^{\infty}(\Omega)}\right)^{\frac{1}{2}}|||\chi|||_{E}$
$\displaystyle\leq C\left(N\epsilon^{-\frac{1}{2}}N^{-1}\epsilon(N^{-1}\ln
N)^{2(k+1)}\right)^{\frac{1}{2}}|||\chi|||_{E}$ $\displaystyle\leq
C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1}|||\chi|||_{E}.$

And in a similar way,

$\displaystyle Y_{7}$ $\displaystyle\leq
C\left(\sum_{j=1}^{N}\sum_{i\in\\{0,N\\}}\lambda_{i,y}h_{j}||\eta_{u}||^{2}_{L^{\infty}(\Omega)}\right)^{\frac{1}{2}}|||\chi|||_{E}\leq
C\epsilon^{\frac{1}{4}}(N^{-1}\ln N)^{k+1}|||\chi|||_{E}.$

From Hölder inequalities, $\lambda_{P}=\epsilon^{-\frac{1}{2}}$ and (25), one
has

$\displaystyle Y_{8}\leq C\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}|||\chi|||_{E}.$

Combining (28) and (29), there is

$\displaystyle|||\chi|||^{2}_{E}=B(\eta;\chi)$ $\displaystyle\leq
C\left(\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+\epsilon^{\frac{1}{2}}N^{-k}(\ln N)^{k+1}\right.$
$\displaystyle\left.+C\epsilon^{\frac{3}{4}}N^{-k}(\ln
N)^{k+\frac{3}{2}}\right)|||\chi|||_{E}.$

Therefore, it is straightforward to derive

(30) $\displaystyle|||\chi|||_{E}\leq C\left(\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+\epsilon^{\frac{1}{2}}N^{-k}(\ln
N)^{k+1}+\epsilon^{\frac{3}{4}}N^{-k}(\ln N)^{k+\frac{3}{2}}\right),$ (31)
$\displaystyle|||\chi|||_{B}\leq C\epsilon^{-\frac{1}{4}}|||\chi|||_{E}\leq
C\left((N^{-1}\ln N)^{k+1}+\epsilon^{\frac{1}{4}}N^{-k}(\ln
N)^{k+1}+\epsilon^{\frac{1}{2}}N^{-k}(\ln N)^{k+\frac{3}{2}}\right).$

From Lemma 5, 26 and 27,

(32) $\displaystyle|||\eta|||_{E}\leq C\epsilon^{\frac{1}{4}}(N^{-1}\ln
N)^{k+1}+C\epsilon^{\frac{1}{2}}N^{-(k+1)}(\ln
N)^{k+\frac{3}{2}}+CN^{-(k+1)},$ (33) $\displaystyle|||\eta|||_{B}\leq
C\epsilon^{\frac{1}{4}}N^{-(k+1)}(\ln N)^{k+\frac{3}{2}}+C(N^{-1}\ln
N)^{k+1}.$

From (30), (31), (32) and (33), The proof of Theorem 3.1 is complete.

###### Remark 4

For Theorem 3.1, if we choose $\epsilon\leq CN^{-4}$, then

$|||e|||_{E}\leq CN^{-(k+2)}(\ln N)^{k+1}+CN^{-(k+1)},\quad|||e|||_{B}\leq
C(N^{-1}\ln N)^{k+1}.$

Unlike c2022-balanced , the optimal convergence order can be derived when the
smooth component of the solution does not belong to the finite element space.

##  4 Numerical experiments

In this section, we provide a numerical example to verify Theorem 2.1. The
following singularly perturbed reaction-diffusion problem is considered:

$-\epsilon u^{\prime\prime}(x)+u(x)=f(x)\quad\text{in $\Omega:=(0,1)$},\quad
u(0)=u(1)=0.$

Here $f(x)$ is selected to make the solution of the above problem satisfy

$u(x)=\frac{1-e^{-1/\sqrt{\epsilon}}}{e^{-x/\sqrt{\epsilon}}-e^{-(1-x)/\sqrt{\epsilon}}}-\cos(\pi
x).$

In this experiment, we choose
$\epsilon=10^{-4},10^{-6},10^{-8},10^{-10},10^{-12}$, $N=32,64,\ldots,1024$
and set $k=1,2,3$. On the layer-adapted mesh (6), we choose $\sigma=k+1$. Set
$e_{N}$ as the error of $|||e|||_{E}$ or $|||e|||_{B}$, where $N$ is the mesh
parameter. The convergence rates are calculated with the formula

$r_{p}=\frac{\log e_{N}-\log e_{2N}}{\log(2\ln N/\ln 2N)}.$

Table 1 shows the errors and convergence rates in the energy norm. The
convergence rates and errors under the balanced norm are displayed in Table 2.
Then the data provided by Table 1 and Table 2 support our main conclusion.
Furthermore, by comparing the data in Table 1 and Table 2 for each $\epsilon$,
it is obvious that the results are approximately
$|||e|||_{B}=\epsilon^{-1/4}|||e|||_{E}$. The example in 2D is not provided
because it is similar to the situation in 1D.

Table 1: The errors and convergence rates under the energy norm |  | $\epsilon=10^{-4}$ | $\epsilon=10^{-6}$ | $\epsilon=10^{-8}$ | $\epsilon=10^{-10}$ | $\epsilon=10^{-12}$
---|---|---|---|---|---|---
$k$ | $N$ | $|||e|||_{E}$ | $r_{p}$ | $|||e|||_{E}$ | $r_{p}$ | $|||e|||_{E}$ | $r_{p}$ | $|||e|||_{E}$ | $r_{p}$ | $|||e|||_{E}$ | $r_{p}$
$1$ | $32$ | 0.26E-1 | 1.59 | 0.82E-1 | 1.59 | 0.26E-2 | 1.59 | 0.82E-3 | 1.59 | 0.26E-3 | 1.59
| $64$ | 0.12E-1 | 1.75 | 0.37E-2 | 1.75 | 0.12E-2 | 1.75 | 0.37E-3 | 1.75 | 0.12E-3 | 1.75
| $128$ | 0.45E-2 | 1.85 | 0.14E-2 | 1.85 | 0.45E-3 | 1.85 | 0.14E-3 | 1.85 | 0.45E-4 | 1.85
| $256$ | 0.16E-2 | 1.92 | 0.51E-3 | 1.92 | 0.16E-3 | 1.92 | 0.51E-4 | 1.92 | 0.16E-4 | 1.92
| $512$ | 0.53E-3 | 1.95 | 0.17E-3 | 1.95 | 0.53E-4 | 1.95 | 0.17E-4 | 1.95 | 0.53E-5 | 1.95
| $1024$ | 0.17E-3 | — | 0.54E-4 | — | 0.17E-4 | — | 0.54E-5 | — | 0.17E-5 | —
$2$ | $32$ | 0.63E-2 | 2.44 | 0.20E-2 | 2.44 | 0.63E-3 | 2.44 | 0.20E-3 | 2.44 | 0.63E-4 | 2.44
| $64$ | 0.18E-2 | 2.65 | 0.57E-3 | 2.65 | 0.18E-3 | 2.65 | 0.57E-4 | 2.65 | 0.18E-4 | 2.65
| $128$ | 0.43E-3 | 2.80 | 0.14E-3 | 2.80 | 0.43E-4 | 2.80 | 0.14E-4 | 2.80 | 0.43E-5 | 2.80
| $256$ | 0.91E-4 | 2.88 | 0.29E-4 | 2.88 | 0.91E-5 | 2.88 | 0.29E-5 | 2.88 | 0.91E-6 | 2.88
| $512$ | 0.17E-4 | 2.94 | 0.55E-5 | 2.94 | 0.17E-5 | 2.94 | 0.55E-6 | 2.94 | 0.17E-6 | 2.94
| $1024$ | 0.31E-5 | — | 0.97E-6 | — | 0.31E-6 | — | 0.97E-7 | — | 0.31E-7 | —
$3$ | $32$ | 0.15E-2 | 3.30 | 0.48E-3 | 3.30 | 0.15E-3 | 3.30 | 0.48E-4 | 3.30 | 0.15E-4 | 3.30
| $64$ | 0.28E-3 | 3.56 | 0.89E-4 | 3.56 | 0.28E-4 | 3.56 | 0.89E-5 | 3.56 | 0.28E-5 | 3.56
| $128$ | 0.41E-4 | 3.74 | 0.13E-4 | 3.74 | 0.41E-5 | 3.74 | 0.13E-5 | 3.74 | 0.41E-6 | 3.74
| $256$ | 0.51E-5 | 3.86 | 0.16E-5 | 3.86 | 0.51E-6 | 3.86 | 0.16E-6 | 3.86 | 0.51E-7 | 3.86
| $512$ | 0.56E-6 | — | 0.18E-6 | — | 0.56E-7 | — | 0.18E-7 | — | 0.56E-8 | —
Table 2: The errors and convergence rates under the balanced norm |  | $\epsilon=10^{-4}$ | $\epsilon=10^{-6}$ | $\epsilon=10^{-8}$ | $\epsilon=10^{-10}$ | $\epsilon=10^{-12}$
---|---|---|---|---|---|---
$k$ | $N$ | $|||e|||_{B}$ | $r_{p}$ | $|||e|||_{B}$ | $r_{p}$ | $|||e|||_{B}$ | $r_{p}$ | $|||e|||_{B}$ | $r_{p}$ | $|||e|||_{B}$ | $r_{p}$
$1$ | $32$ | 0.25E+0 | 1.56 | 0.25E+0 | 1.56 | 0.25E+0 | 1.56 | 0.25E+0 | 1.56 | 0.25E+0 | 1.56
| $64$ | 0.11E+0 | 1.72 | 0.11E+0 | 1.72 | 0.11E+0 | 1.72 | 0.11E+0 | 1.72 | 0.11E+0 | 1.72
| $128$ | 0.44E-1 | 1.83 | 0.44E-1 | 1.83 | 0.44E-1 | 1.83 | 0.44E-1 | 1.83 | 0.44E-1 | 1.83
| $256$ | 0.16E-1 | 1.90 | 0.16E-1 | 1.90 | 0.16E-1 | 1.90 | 0.16E-1 | 1.90 | 0.16E-1 | 1.90
| $512$ | 0.53E-2 | 1.95 | 0.53E-2 | 1.95 | 0.53E-2 | 1.95 | 0.53E-2 | 1.95 | 0.53E-2 | 1.95
| $1024$ | 0.17E-2 | — | 0.17E-2 | — | 0.17E-2 | — | 0.17E-2 | — | 0.17E-2 | —
$2$ | $32$ | 0.61E-1 | 2.42 | 0.61E-1 | 2.42 | 0.61E-1 | 2.42 | 0.61E-1 | 2.42 | 0.61E-1 | 2.42
| $64$ | 0.18E-1 | 2.64 | 0.18E-1 | 2.64 | 0.18E-1 | 2.64 | 0.18E-1 | 2.64 | 0.18E-1 | 2.64
| $128$ | 0.43E-2 | 2.79 | 0.43E-2 | 2.79 | 0.43E-2 | 2.79 | 0.43E-2 | 2.79 | 0.43E-2 | 2.79
| $256$ | 0.90E-3 | 2.88 | 0.90E-3 | 2.88 | 0.90E-3 | 2.88 | 0.90E-3 | 2.88 | 0.90E-3 | 2.88
| $512$ | 0.17E-3 | 2.92 | 0.17E-3 | 2.92 | 0.17E-3 | 2.92 | 0.17E-3 | 2.92 | 0.17E-3 | 2.92
| $1024$ | 0.31E-4 | — | 0.31E-4 | — | 0.31E-4 | — | 0.31E-4 | — | 0.31E-4 | —
$3$ | $32$ | 0.15E-1 | 3.27 | 0.15E-1 | 3.27 | 0.15E-1 | 3.27 | 0.15E-1 | 3.27 | 0.15E-1 | 3.27
| $64$ | 0.28E-2 | 3.55 | 0.28E-2 | 3.55 | 0.28E-2 | 3.55 | 0.28E-2 | 3.55 | 0.28E-2 | 3.55
| $128$ | 0.41E-3 | 3.73 | 0.41E-3 | 3.73 | 0.41E-3 | 3.73 | 0.41E-3 | 3.73 | 0.41E-3 | 3.73
| $256$ | 0.51E-4 | 3.84 | 0.51E-4 | 3.84 | 0.51E-4 | 3.84 | 0.51E-4 | 3.84 | 0.51E-4 | 3.84
| $512$ | 0.55E-5 | — | 0.55E-5 | — | 0.55E-5 | — | 0.55E-5 | — | 0.55E-5 | —

##  5 Statements and Declarations

###  5.1 Funding

This research is supported by National Natural Science Foundation of China
(11771257) and Shandong Provincial Natural Science Foundation, China
(ZR2023YQ002).

###  5.2 Data availability statement

The authors confirm that the data supporting the findings of this study are
available within the article and its supplementary materials.

###  5.3 Conflict of interests

The authors declare that they have no conflict of interest.

## References

  * (1) R. Lin and M. Stynes.  A balanced finite element method for singularly perturbed reaction-diffusion problems.  SIAM J. Numer. Anal. 50(5), 2729–2743 (2012).
  * (2) H.G. Roos and M. Schopf.  Convergence and stability in balanced norms of finite element methods on Shishkin meshes for reaction-diffusion problems.  ZAMMZ. Angew. Math.Mech. 95(6), 551–565 (2015).
  * (3) Z. Xie, P. Zhu, and S. Zhou.  Uniform convergence of a coupled method for convection-diffusion problems in 2D Shishkin mesh.  Int. J. Numer. Anal. Model., 10(4):845–859, 2013.
  * (4) J. Zhang and X. Liu.  Supercloseness of the SDFEM on Shishkin triangular meshes for problems with exponential layers.  Adv. Comput. Math., 43(4):759–775, 2017.
  * (5) X. Ma and J. Zhang.  Supercloseness analysis of the nonsymmetric interior penalty Galerkin method for a singularly perturbed problem on Bakhvalov-type mesh.  Appl. Math. Lett., 144:108701, 7, 2023.
  * (6) D. Di Pietro and A. Ern.  Mathematical Aspects of Discontinuous Galerkin Methods, volume 69 of Mathématiques & Applications (Berlin) [Mathematics &Applications].  Springer, Heidelberg, 2012.
  * (7) Y. Cheng, L. Yan, and Y.J. Mei.  Balanced norm error estimate of the local discontinuous Galerkin method on layer-adapted meshes for reaction-diffusion problems.  Numer. Algorithms. pages 1–30, 2022.
  * (8) B. Cockburn, G.E. Karniadakis, and C.W. Shu.  Discontinuous Galerkin Methods. Theory, Computation and Applications.  Lect. Notes Comput. Sci. Eng. 11, Springer, Berlin, 2000.
  * (9) B. Cockburn and C.W. Shu.  The local discontinuous Galerkin method for time-dependent convection-diffusion systems.  SIAM J. Numer. Anal. 35 (1998), pp. 2440–2463.
  * (10) N. C. Nguyen, J. Peraire, and B. Cockburn.  An implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusion equations.  J. Comput. Phys., 228(9):3232–3254, 2009.
  * (11) Y.J. Mei, Y. Cheng, S. Wang, and Z. Xu.  Local discontinuous Galerkin method on layer-adapted meshes for singularly perturbed reaction-diffusion problems in two dimensions.  arXiv preprint arXiv:2103.01083, 2021.
  * (12) Y.J. Mei, S. Wang, Z. Xu, C. Song, and Y. Cheng.  Convergence analysis of the LDG method for singularly perturbed reaction-diffusion problems.  Symmetry. 13(12):2291, 2021.
  * (13) J.J.H. Miller, E. O’Riordan, and G.I. Shishkin.  Fitted Numerical Methods for Singular Perturbation Problems: error estimates in the maximum norm for linear problems in one and two dimensions.  Sci. World. J. 1996.
  * (14) Y. Cheng.  On the local discontinuous Galerkin method for singularly perturbed problem with two parameters.  J. Comput. Appl. Math., 392:113485, 2021\.
  * (15) N.S. Bakhvalov.  Towards optimization of methods for solving boundary value problems in the presence of boundary layers. Zh. Vychisl. Mat. Mat. Fiz. 9:841–859, 1969.
  * (16) M. Stynes and D. Stynes.  Convection-diffusion problems, volume 196 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI; Atlantic Association for Research in the Mathematical Sciences (AARMS), Halifax, NS, 2018.
  * (17) C. Clavero, J. L. Gracia, and E. O’Riordan.  A parameter robust numerical method for a two dimensional reaction-diffusion problem.  Math. Comp., 74(252):1743–1758, 2005.
  * (18) H. Roos, M. Stynes and L. Tobiska.  Robust Numerical Methods for Singularly Perturbed Differential Equations. Springer Series in Computational Mathematics, 2nd edn. Springer-Verlag, Berlin (2008).
  * (19) H. Zhu and Z. Zhang.  Pointwise error estimates for the LDG method applied to 1-d singularly perturbed reaction-diffusion problems.  Comput. Meth. Appl. Mat. 13(1):79–94, 2013\.
  * (20) X. Meng, C.W. Shu and B. Wu.  Optimal error estimates for discontinuous Galerkin methods based on upwind-biased fluxes for linear hyperbolic equations.  Math. Comp. 85, no. 299, 1225–1261, 2016.