∎

11institutetext: Corresponding author: Jin Zhang 22institutetext: School of
Mathematics and Statistics, Shandong Normal University, Jinan 250014, China
22email<EMAIL_ADDRESS>33institutetext: Chunxiao
Zhang44institutetext: School of Mathematics and Statistics, Shandong Normal
University, Jinan 250014, China
44email<EMAIL_ADDRESS>
# Uniform convergence of finite element method on Bakhvalov-type mesh for a
2-D singularly perturbed convection-diffusion problem with exponential layers
††thanks: The current research was partly supported by NSFC (11771257),
Shandong Provincial NSF (ZR2021MA004).

Jin Zhang  Chunxiao Zhang

(Received: date / Accepted: date)

###### Abstract

On Bakhvalov-type mesh, uniform convergence analysis of finite element method
for a 2-D singularly perturbed convection-diffusion problem with exponential
layers is still an open problem. Previous attempts have been unsuccessful. The
primary challenges are the width of the mesh subdomain in the layer adjacent
to the transition point, the restriction of the Dirichlet boundary condition,
and the structure of exponential layers. To address these challenges, a novel
analysis technique is introduced for the first time, which takes full
advantage of the characteristics of interpolation and the connection between
the smooth function and the layer function on the boundary. Utilizing this
technique in conjunction with a new interpolation featuring a simple
structure, uniform convergence of optimal order $k+1$ under an energy norm can
be proven for finite element method of any order $k$. Numerical experiments
confirm our theoretical results.

###### Keywords:

Singularly perturbed Convection–diffusion Bakhvalov-type mesh Finite element
method Uniform convergence

###### MSC:

65N12 65N30

##  1 Introduction

Take the following elliptic boundary value problem into consideration:

(1) $\displaystyle-\varepsilon\Delta u-\boldsymbol{b}\cdot\nabla u+cu=$
$\displaystyle f\quad$ $\displaystyle\text{in}\quad\Omega=(0,1)^{2},$
$\displaystyle u=$ $\displaystyle 0\quad$
$\displaystyle\text{on}\quad\partial\Omega,$

where $0<\varepsilon\ll 1$ is a small positive constant, and $\boldsymbol{b}$
denotes $(b_{1}(x,y),b_{2}(x,y))$. Let $b_{1}$, $b_{2}$, and $f$ be sufficient
smooth functions satisfying

(2) $b_{1}(x,y)\geq\beta_{1}>0,\ b_{2}(x,y)\geq\beta_{2}>0,\
c(x,y)+\frac{1}{2}\nabla\cdot\boldsymbol{b}(x,y)\geq\gamma>0\ \text{on}\
\Omega,$

with some positive constants $\beta_{1}$, $\beta_{2}$ and $\gamma$. Under the
given conditions, equation (1) possesses a unique solution in the space
$H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$ for all $f\in L^{2}(\Omega)$
Roo1Sty2:2008-Robust . This solution exhibits exponential boundary layers of
width $\mathcal{O}(\varepsilon\ln(1/\varepsilon))$ at $x=0$ and $y=0$ as well
as a corner layer at $(0,0)$. Since the diffusion parameter $\varepsilon$ can
be arbitrarily small, this problem is characterized as singularly perturbed
and is convection-dominated.

Studying such singularly perturbed convection-diffusion problems has
significant theoretical and practical value. In order to obtain the required
accuracy of the numerical solution, various uniformly convergent numerical
methods with respect to the small parameter $\varepsilon$ were proposed (see
Roo1Sty2:2008-Robust ). Among them, the combination of finite element method
and layer-adapted meshes is particularly effective. Shishkin-type mesh and
Bakhvalov-type mesh are two popular types of layer-adapted meshes
Lin:2009-Layer . Shishkin-type mesh has a simple structure and has been widely
used to study convergence (see Li1Nav2:1998-Uniformly ; Lin1:2000-Uniform ;
Zhang1:2003-Finite ; Roo1Sch2:2015-Convergence ; Lin1Sty2:2012-Balanced ;
Zhang1Liu2:2016-Optimal ; Liu1Zhang2:2018-Uniform  and their references
therein). In contrast to Shishkin-type mesh, Bakhvalov-type mesh has a complex
structure as it is graded in the mesh subdomain used to resolve layers.
Nevertheless, its transition point is independent of the mesh parameter $N$,
which is an essential property in certain cases. Furthermore, Bakhvalov-type
mesh performs numerically better than Shishkin-type mesh, particularly when
higher-order finite element schemes are employed.

Up to this point, uniform convergence analysis for singularly perturbed
problems on Bakhvalov-type mesh is still an open problem, which is attributed
to the intricate structure of the mesh (see (Roo1Sty2:2015-Some, , Question
4.1) for further information). In an effort to deal with the problems raised
in (Roo1Sty2:2015-Some, , Question 4.1), several convergence studies utilizing
finite element method have been carried out. In the case of 1D, Roos
Roo1:2006-Error  introduced a quasi-interpolation for linear finite element
method and proved uniform convergence of optimal order, while this approach
cannot be generalized to higher-order finite element schemes or higher-
dimensional cases. Subsequently, Zhang zhang1Liu2:2020-Optimal  developed a
novel interpolation with a simplified structure suitable for finite element
method of any order. In the case of 2D, Zhang zhang1Liu2:2023-Convergence
established uniform convergence for a singularly perturbed convection-
diffusion problem with parabolic layers by extending the novel interpolation
in zhang1Liu2:2020-Optimal  to the case of 2D. However, when considering a 2-D
singularly perturbed convection-diffusion problem with exponential layers,
previous research, including Roo1Sch2:2012-Analysis  and
zhang1Liu2:2023-Convergence , both failed to achieve uniform convergence. The
main difficulties come from the special width of the Bakhvalov-type mesh
subdomain in the layer adjacent to the transition point, the restriction of
the Dirichlet boundary condition, and the structure of exponential layers.
These difficulties result in suboptimal results when using the standard error
analysis methods to estimate some convection terms of the layer part, and this
situation does not arise in parabolic layers.

This article aims to present a novel analysis technique for solving the 2-D
singularly perturbed convection-diffusion problem with exponential layers. To
be specific, for the smooth part in particular regions, we convert the
corresponding estimates in the two-dimensional case into several types of one-
dimensional interpolation error estimates. Furthermore, we make full use of
the relationship between the smooth function and the layer function on the
boundary. For the layer part, we apply a new interpolation with a simplified
structure that were proposed in zhang1Liu2:2023-Convergence  and adopt the
standard error analysis approach. Our scheme is demonstrated to uniformly
converge to an optimal order of $k+1$. To the best of our knowledge, this is
the first proof of uniform convergence for a typical problem (1) on Bakhvalov-
type mesh in the two-dimensional setting. It is worth mentioning that the
novel analysis technique proposed in this paper provides a new path for the
future convergence analysis.

The structure of this article is as follows: In Section 2, we present the
regularity of the solution to (1), define a Bakhvalov-type mesh, and introduce
the $k$-th order finite element method. New interpolation and preliminary
results are displayed in Section 3. Detailed derivations of uniform
convergence under an energy norm can be found in Section 4. Section 5 carries
out the numerical experiments that illustrate our theoretical results.

Assume that $D$ is any measurable subset of $\Omega$. $(\cdot,\cdot)_{D}$,
$\|\cdot\|_{L^{1}(D)}$, $\|\cdot\|_{\infty,D}$, $\|\cdot\|_{D}$ and
$|\cdot|_{1,D}$ represent the standard inner product in $L^{2}(D)$, the
standard norms in $L^{1}(D)$, $L^{\infty}(D)$, $L^{2}(D)$ and the standard
seminorms in $H^{1}(D)$, respectively. When $D=\Omega$, drop $D$ from these
notations for the sake of clarity. In addition, all constants used in this
paper, including the generic constant $C$ and the fixed constant $C_{i}$, are
all positive and unaffected by the singular perturbation $\varepsilon$ as well
as the mesh parameter N.

##  2 Regularity, Bakhvalov-type mesh and finite element method

###  2.1 Regularity of the solution

In the process of analysis, we assume that $k\geq 1$ is a fixed integer.

###### Assumption 1

We can decompose the solution to (1) into the following parts:

(3) $u=S+E_{1}+E_{2}+E_{12}\quad\forall(x,y)\in\overline{\Omega},$

where $S$ is the smooth part, $E_{1}$ is the exponential layer part at $x=0$,
$E_{2}$ is another exponential layer part at $y=0$, and $E_{12}$ is the corner
layer part.

For any $(x,y)\in\bar{\Omega}$, the following bounds hold:

(4a) $\displaystyle\left|\frac{\partial^{m+n}S}{\partial x^{m}\partial
y^{n}}(x,y)\right|$ $\displaystyle\leq C,$ (4b)
$\displaystyle\left|\frac{\partial^{m+n}E_{1}}{\partial x^{m}\partial
y^{n}}(x,y)\right|$ $\displaystyle\leq
C\varepsilon^{-m}e^{\frac{-\beta_{1}x}{\varepsilon}},$ (4c)
$\displaystyle\left|\frac{\partial^{m+n}E_{2}}{\partial x^{m}\partial
y^{n}}(x,y)\right|$ $\displaystyle\leq
C\varepsilon^{-n}e^{\frac{-\beta_{2}y}{\varepsilon}},$ (4d)
$\displaystyle\left|\frac{\partial^{m+n}E_{12}}{\partial x^{m}\partial
y^{n}}(x,y)\right|$ $\displaystyle\leq
C\varepsilon^{-(m+n)}e^{\frac{-\beta_{1}x-\beta_{2}y}{\varepsilon}},$

where $m$, $n$ are nonnegative integers with $0\leq m+n\leq k+1$.

We can refer to Lin1Sty2:2001-Asymptotic  and (Sty:2005-Steady, , Sec. 7) for
concrete information regarding Assumption 1.

###  2.2 Bakhvalov-type mesh

Bakhvalov mesh was first proposed in Bakh:1969-Optimization  to improve the
convergence order. Bakhvalov-type mesh, as an approximation of Bakhvalov mesh,
has been widely used to avoid solving nonlinear equations appearing on
Bakhvalov mesh Lin:2009-Layer .

In this study, we consider a Bakhvalov-type mesh introduced in Lin:2009-Layer
, which is defined by

(5) $\displaystyle x_{i}=\left\\{\begin{aligned}
&-\frac{\sigma\varepsilon}{\beta_{1}}\ln(1-2(1-\varepsilon)i/N)\quad&&\text{for}\quad
i=0,1,\dots,N/2,\\\ &1-2(1-x_{N/2})(N-i)/N\quad&&\text{for}\quad
i=N/2+1,\dots,N,\end{aligned}\right.$ $\displaystyle
y_{j}=\left\\{\begin{aligned}
&-\frac{\sigma\varepsilon}{\beta_{2}}\ln(1-2(1-\varepsilon)j/N)\quad&&\text{for}\quad
j=0,1,\dots,N/2,\\\ &1-2(1-y_{N/2})(N-j)/N\quad&&\text{for}\quad
j=N/2+1,\dots,N,\end{aligned}\right.$

where $\sigma\geq k+1$ and $N$ is an even positive integer. Transition points
are $x_{N/2}=-\frac{\sigma\varepsilon}{\beta_{1}}\ln\varepsilon$ and
$y_{N/2}=-\frac{\sigma\varepsilon}{\beta_{2}}\ln\varepsilon$, indicating a
shift in mesh from coarse to fine. Also, one has
$e^{-\beta_{1}x_{N/2}/\varepsilon}=\varepsilon^{\sigma}$ and
$e^{-\beta_{1}x_{N/2-1}/\varepsilon}=N^{-\sigma}$.

By connecting the mesh points $\left\\{(x_{i},y_{j})|i,j=0,1,\dots,N\right\\}$
with lines parallel to the x-axis and the y-axis, we obtain a rectangulation
of $\Omega$ denoted as $\mathcal{T}_{N}$ ( see Figure 1). We define
$K_{i,j}=[x_{i},x_{i+1}]\times[y_{j},y_{j+1}]$ as a specific mesh rectangle,
where $i,j=0,1,\dots,N-1$, and $K\in\mathcal{T}_{N}$ as a general mesh
rectangle. The lengths of $K_{i,j}$ in the $x$\- and $y$-directions are
represented by $h_{x,i}:=x_{i+1}-x_{i}$ and $h_{y,j}:=y_{j+1}-y_{j}$,
respectively.

Figure 1: Rectangulation of $\Omega$

###### Assumption 2

In this paper, suppose that $\varepsilon\leq N^{-1}$, because it is not a
practical restriction.

In the following lemma, we will present some important properties of
Bakhvalov-type mesh (5) (refer to (zhang1Liu2:2020-Optimal, , Lemma 3)).

###### Lemma 1

Let Assumption 2 hold, then for Bakhvalov-type mesh (5), one can obtain

(6) $\displaystyle C_{1}\varepsilon N^{-1}\leq h_{x,0}\leq C_{2}\varepsilon
N^{-1},$ (7) $\displaystyle h_{x,0}\leq h_{x,1}\leq\dots\leq h_{x,N/2-2},$ (8)
$\displaystyle\frac{1}{4}\sigma\varepsilon\leq
h_{x,N/2-2}\leq\sigma\varepsilon,$ (9)
$\displaystyle\frac{1}{2}\sigma\varepsilon\leq h_{x,N/2-1}\leq 2\sigma
N^{-1},$ (10) $\displaystyle N^{-1}\leq h_{x,i}\leq 2N^{-1}\quad N/2\leq i\leq
N-1,$ (11) $\displaystyle C_{3}\sigma\varepsilon\ln N\leq x_{N/2-1}\leq
C_{4}\sigma\varepsilon\ln N,\quad x_{N/2}\geq
C\sigma\varepsilon|\ln\varepsilon|,$ (12) $\displaystyle
h_{x,i}^{\rho}e^{-\beta_{1}x/\varepsilon}\leq
C\varepsilon^{\rho}N^{-\rho}\quad 0\leq i\leq N/2-2\quad\text{and}\quad
0\leq\rho\leq\sigma.$

Bounds for $h_{y,j}$ $0\leq j\leq N-1$, are analogous.

###  2.3 Finite element method

The variational form of (1) can be described as

(13) $\left\\{\begin{aligned} &\text{find}\ u\in H_{0}^{1}(\Omega)\ \text{such
that for all}\ v\in H_{0}^{1}(\Omega)\\\ &a(u,v)=(f,v),\end{aligned}\right.$

where

(14) $a(u,v):=\varepsilon(\nabla u,\nabla v)+(-\boldsymbol{b}\cdot\nabla
u,v)+(cu,v).$

Define a finite element space $V^{N}\subset H_{0}^{1}(\Omega)$ on
$\mathcal{T}_{N}$:

$V^{N}:=\left\\{v^{N}\in C(\overline{\Omega}):\ v^{N}|_{\partial\Omega}=0\
\text{and}\ v^{N}|_{K}\in\mathcal{Q}_{k}(K)\quad\forall
K\in\mathcal{T}_{N}\right\\},$

with $\mathcal{Q}_{k}(K)=\text{span}\left\\{x^{i}y^{j}:\ 0\leq i,j\leq
k\right\\}$.

Then the $k$-th order finite element method is

(15) $\left\\{\begin{aligned} &\text{find}\ u^{N}\in V^{N}\ \text{such that
for all}\ v\in V^{N}\\\ &a(u^{N},v^{N})=(f,v^{N}),\end{aligned}\right.$

where

(16) $a(u^{N},v^{N}):=\varepsilon(\nabla u^{N},\nabla
v^{N})+(-\boldsymbol{b}\cdot\nabla u^{N},v^{N})+(cu^{N},v^{N}).$

Condition (2) gives the coercivity

(17) $a(v^{N},v^{N})\geq C\|v^{N}\|_{\varepsilon}^{2}\quad\forall v^{N}\in
V^{N},$

here the energy norm $\|\cdot\|_{\varepsilon}$ is defined as

(18)
$\|v\|_{\varepsilon}^{2}:=\left(\varepsilon|v|_{1}^{2}+\|v\|^{2}\right)^{1/2}.$

The following Galerkin orthogonality property can be obtained by (13) and
(15):

(19) $a(u-u^{N},v^{N})=0\quad\forall v^{N}\in V^{N}.$

Furthermore, by employing the Lax-Milgram lemma, existence and uniqueness of
the finite element solution $u^{N}$ can be well demonstrated
Bre1Sco2:2008-Mathematical .

##  3 New interpolation, interpolation errors and some preliminary results

In this section, we will introduce a novel interpolation that were proposed in
zhang1Liu2:2023-Convergence . Interpolation errors and preliminary results are
also presented. To begin with, some notations need to be defined: For $k\geq
1$, set $x_{i}^{s}:=x_{i}+(s/k)h_{x,i}$ and $y_{j}^{t}:=y_{j}+(t/k)h_{y,j}$,
where $i,j=0,1,\dots,N-1$ and $s,t=0,1,\dots,k-1$. Besides, set
$x_{N}^{0}=x_{N}$ and $y_{N}^{0}=y_{N}$ to ensure notational uniformity. For
any $v\in C^{0}(\overline{\Omega})$, its standard Lagrange interpolation
$v^{I}\in V^{N}$ can be expressed as

(20) $\displaystyle v^{I}(x,y)=$
$\displaystyle\sum_{i=0}^{N-1}\sum_{s=0}^{k-1}\left(\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}v(x_{i}^{s},y_{j}^{t})\theta_{i,j}^{s,t}(x,y)+v(x_{i}^{s},y_{N}^{0})\theta_{i,N}^{s,0}(x,y)\right)$
$\displaystyle+\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}v(x_{N}^{0},y_{j}^{t})\theta_{N,j}^{0,t}(x,y)+v(x_{N}^{0},y_{N}^{0})\theta_{N,N}^{0,0}(x,y),$

here $\theta_{i,j}^{s,t}(x,y)\in V^{N}$ is the piecewise tensor-product
interpolation basis function of degree $k$ connected to the point
$(x_{i}^{s},y_{j}^{t})$. Write $\theta_{i,j}^{s,t}(x,y)$ as
$\theta_{i,j}^{s,t}$ for simplicity.

Recall (3), we can decompose the new interpolation $\Pi u$ as

(21) $\Pi u=S^{I}+\pi_{1}E_{1}+\pi_{2}E_{2}+\pi_{12}E_{12},$

where $S^{I}$ represents the standard Lagrange interpolation of $S$, and

(22)
$\displaystyle\pi_{1}E_{1}=E_{1}^{I}-\sum_{s=0}^{k-1}\left(\sum_{t=1}^{k-1}E_{1}(x_{N/2-1}^{s},y_{0}^{t})\theta_{N/2-1,0}^{s,t}-\sum_{j=1}^{N-1}\sum_{t=0}^{k-1}E_{1}(x_{N/2-1}^{s},y_{j}^{t})\theta_{N/2-1,j}^{s,t}\right),$
(23)
$\displaystyle\pi_{2}E_{2}=E_{2}^{I}-\sum_{t=0}^{k-1}\left(\sum_{s=1}^{k-1}E_{2}(x_{0}^{s},y_{N/2-1}^{t})\theta_{0,N/2-1}^{s,t}+\sum_{i=1}^{N-1}\sum_{s=0}^{k-1}E_{2}(x_{i}^{s},y_{N/2-1}^{t})\theta_{i,N/2-1}^{s,t}\right),$
(24)
$\displaystyle\pi_{12}E_{12}=E_{12}^{I}-\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}E_{12}(x_{N/2-1}^{s},y_{N/2-1}^{t})\theta_{N/2-1,N/2-1}^{s,t}.$

###### Remark 1

The concept of $\pi_{1}E_{1}$ involves subtracting some intractable terms that
fail to achieve uniform convergence from the standard Lagrange interpolation
$E_{1}^{I}$, and at the same time, ensuring $\Pi u$ satisfies the homogeneous
Dirichlet boundary condition. $\pi_{2}E_{2}$ and $\pi_{12}E_{12}$ are
similarly constructed.

According to (Apel:1999-Anisotropic, , Theorem 2.7), we have the following
interpolation errors.

###### Lemma 2

Let $p\in(1,\infty]$ and $K_{i,j}\in\mathcal{T}_{N}$ with
$i,j=0,1,\cdots,N-1$. For any $w\in W^{k+1,p}(\Omega)$, its standard Lagrange
interpolation $w^{I}$ at the vertices of $K_{i,j}$ satisfies

$\displaystyle\|w-w^{I}\|_{L^{p}(K_{i,j})}\leq
C\sum_{m+n=k+1}h_{x,i}^{m}h_{y,j}^{n}\|\frac{\partial^{k+1}w}{\partial
x^{m}\partial y^{n}}\|_{L^{p}(K_{i,j})},$
$\displaystyle\|(w-w^{I})_{x}\|_{L^{p}(K_{i,j})}\leq
C\sum_{m+n=k}h_{x,i}^{m}h_{y,j}^{n}\|\frac{\partial^{k+1}w}{\partial
x^{m+1}\partial y^{n}}\|_{L^{p}(K_{i,j})},$
$\displaystyle\|(w-w^{I})_{y}\|_{L^{p}(K_{i,j})}\leq
C\sum_{m+n=k}h_{x,i}^{m}h_{y,j}^{n}\|\frac{\partial^{k+1}w}{\partial
x^{m}\partial y^{n+1}}\|_{L^{p}(K_{i,j})},$

here $m$ and $n$ denote non-negative integers.

Next, some preliminary results will be given.

###### Lemma 3

Let Assumptions 1 and 2 hold. Let $\sigma\geq k+1$. Denote $E_{i}^{I}$ as the
standard Lagrange interpolations of $E_{i}\,(i=1,2,12)$, respectively. Then,
the following interpolation error estimates can be obtained on Bakhvalov-type
mesh (5).

(25) $\displaystyle\|E_{1}-E_{1}^{I}\|+\|E_{2}-E_{2}^{I}\|\leq CN^{-(k+1)},$
(26) $\displaystyle\|E_{12}-E_{12}^{I}\|\leq C\varepsilon
N^{-k}+C\varepsilon^{1/2}N^{-(k+1)}+CN^{-(1+2\sigma)},$ (27)
$\displaystyle\|E_{1}-E_{1}^{I}\|_{\varepsilon}+\|E_{2}-E_{2}^{I}\|_{\varepsilon}+\|E_{12}-E_{12}^{I}\|_{\varepsilon}\leq
CN^{-k}.$

###### Proof

Detailed demonstrations can refer to (zhang1Liu2:2023-Convergence, , Lemmas
3.2 and 3.3).

We also require the following inequalities for the subsequent analysis:
$\forall K\in\mathcal{T}_{N}$, for $i=1,2,12$, we have

(28) $\displaystyle\|E_{i}^{I}\|_{\infty,K}\leq\|E_{i}\|_{\infty,K},$ (29)
$\displaystyle\|\pi_{i}E_{i}\|_{\infty,K}\leq\|E_{i}\|_{\infty,K},$

where $E_{i}^{I}$ denote the standard Lagrange interpolation and
$\pi_{i}E_{i}$ are defined in (22)-(24).

##  4 Uniform convergence

In this section, we will prove uniform convergence of optimal order for
problem (1) on Bakhvalov-type mesh (5).

For the sake of briefness, set

$\displaystyle\times[y_{N-1},y_{N}]$ $\displaystyle=:K_{a},$
$\displaystyle[x_{N-1},x_{N}]\times[y_{N/2-1},y_{N/2}]$
$\displaystyle=:K_{b},$ $\displaystyle\Omega\textbackslash(K_{a}\cup K_{b})$
$\displaystyle=:\Omega_{0}.$

Let $v^{N}:=\Pi u-u^{N}$. Recall (16), (17), (19), (3) and (21), one can
obtain

(30) $\displaystyle C\|v^{N}\|_{\varepsilon}^{2}\leq a(v^{N},v^{N})=a(\Pi
u-u,v^{N})$ $\displaystyle=$ $\displaystyle\varepsilon\int_{\Omega}\nabla(\Pi
u-u)\nabla v^{N}\mathrm{d}x\mathrm{d}y+\int_{\Omega}c(\Pi
u-u)v^{N}\mathrm{d}x\mathrm{d}y$
$\displaystyle-\int_{K_{a}}\boldsymbol{b}\cdot\nabla(\Pi
u-u)v^{N}\mathrm{d}x\mathrm{d}y-\int_{K_{b}}\boldsymbol{b}\cdot\nabla(\Pi
u-u)v^{N}\mathrm{d}x\mathrm{d}y$
$\displaystyle-\sum_{i=1,2,12}\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{i}E_{i}-E_{i})v^{N}\mathrm{d}x\mathrm{d}y-\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(S^{I}-S)v^{N}\mathrm{d}x\mathrm{d}y$
$\displaystyle=:$ $\displaystyle I+II+III+IV+V+VI.$

Estimates for $I-VI$ are presented in the following lemmas.

###### Lemma 4

Let Assumptions 1 and 2 hold. Let $\sigma\geq k+1$. Then one has

(31) $|I+II|\leq CN^{-k}\|v^{N}\|_{\varepsilon}.$

###### Proof

Consider (3) and (21), $\Pi u-u$ can be split into

$\Pi
u-u=(\pi_{1}E_{1}-E_{1})+(\pi_{2}E_{2}-E_{2})+(\pi_{12}E_{12}-E_{12})+(S^{I}-S).$

Triangle inequality yields

(32) $\|\Pi
u-u\|^{2}\leq\|\pi_{1}E_{1}-E_{1}\|^{2}+\|\pi_{2}E_{2}-E_{2}\|^{2}+\|\pi_{12}E_{12}-E_{12}\|^{2}+\|S^{I}-S\|^{2}.$

Triangle inequality, Lemma 3, Lemma 1 and (22) yield

(33) $\displaystyle\|\pi_{1}E_{1}-E_{1}\|^{2}$
$\displaystyle\leq\|\pi_{1}E_{1}-E_{1}^{I}\|^{2}+\|E_{1}^{I}-E_{1}\|^{2}\leq
CN^{-2(k+1)},$

where we have used

$\displaystyle\|\pi_{1}E_{1}-E_{1}^{I}\|^{2}$ $\displaystyle\leq
CN^{-2\sigma}\sum_{s=0}^{k-1}\left(\sum_{t=1}^{k-1}\|\theta_{N/2-1,0}^{s,t}\|^{2}+\sum_{j=1}^{N-1}\sum_{t=0}^{k-1}\|\theta_{N/2-1,j}^{s,t}\|^{2}\right)$
$\displaystyle\leq CN^{-2\sigma}\sum_{j=0}^{N-1}(h_{x,N/2-1}h_{y,j})$
$\displaystyle\leq CN^{-2\sigma-1}.$

In a similar manner as $\|\pi_{1}E_{1}-E_{1}\|^{2}$, we can obtain

(34) $\displaystyle\|\pi_{2}E_{2}-E_{2}\|^{2}\leq CN^{-2(k+1)},$ (35)
$\displaystyle\|\pi_{12}E_{12}-E_{12}\|^{2}\leq
C\varepsilon^{2}N^{-2k}+C\varepsilon N^{-2(k+1)}+CN^{-4\sigma-2}.$

(36) $\|S^{I}-S\|^{2}\leq CN^{-2(k+1)},$

can be deduced by Lemma 2, Lemma 1 and Assumption 1.

Collecting (33)-(36), we can prove

(37) $\|\Pi u-u\|\leq CN^{-(k+1)}.$

$\|\Pi u-u\|_{\varepsilon}$ can be analyzed similarly as above, then one has

(38) $\|\Pi u-u\|_{\varepsilon}\leq CN^{-k}.$

In conclusion,

(39) $|I+II|\leq C\|\Pi u-u\|_{\varepsilon}\|v^{N}\|_{\varepsilon}\leq
CN^{-k}\|v^{N}\|_{\varepsilon}.$

In the next lemma, a novel analysis technique will be established to get an
optimal result for III and IV.

###### Lemma 5

Let Assumptions 1 and 2 hold. Let $\sigma\geq k+1$. Then one can obtain

(40) $|III+IV|\leq
CN^{-(k+1)}\|v^{N}\|_{\varepsilon}+C\varepsilon^{\sigma-1/2}N^{-1/2}\|v^{N}\|_{\varepsilon}.$

###### Proof

First, we will analyze $\int_{K_{a}}b_{1}(\Pi
u-u)_{x}v^{N}\mathrm{d}x\mathrm{d}y$, which is the main focus of this article.
Standard error estimates, i.e., studying
$\int_{K_{a}}b_{1}(S^{I}-S)_{x}v^{N}\mathrm{d}x\mathrm{d}y$ and
$\int_{K_{a}}b_{1}(\pi_{i}E_{i}-E_{i})_{x}v^{N}\mathrm{d}x\mathrm{d}y\
(i=1,2,12)$ separately, would result in the inability of

(41)
$|\int_{K_{a}}b_{1}\sum_{s=0}^{k-1}E_{1}(x_{N/2-1}^{s},y_{N}^{0})(\theta_{N/2-1,N}^{s,0})_{x}v^{N}\mathrm{d}x\mathrm{d}y|\leq
Ch_{x,N/2-1}^{-1/2}h_{y,N-1}^{1/2}N^{-\sigma}\|v^{N}\|_{\varepsilon}$

to achieve uniform convergence, where $\frac{1}{2}\sigma\varepsilon\leq
h_{x,N/2-1}\leq 2\sigma N^{-1}$ and $N^{-1}\leq h_{y,N-1}\leq 2N^{-1}$.

To avoid analyzing (41) and to get an optimal result, we propose a new
analysis technique. Firstly, according to the boundary condition $\Pi u=0\
{\text{o}n}\ \partial\Omega$ and (21)-(24), we derive that $\Pi u|_{K_{a}}$
satisfies:

$\left\\{\begin{aligned} &\Pi u(x_{N/2-1}^{s},y_{N}^{0})=0\quad
s=0,1,\dots,k-1,\\\ &\Pi u(x_{N/2}^{0},y_{N}^{0})=0,\\\ &\Pi
u(x_{N/2}^{0},y_{N-1}^{t})=u(x_{N/2}^{0},y_{N-1}^{t})\quad t=0,1,\dots,k-1,\\\
&\Pi
u(x_{N/2-1}^{s},y_{N-1}^{t})=u(x_{N/2-1}^{s},y_{N-1}^{t})-E_{1}(x_{N/2-1}^{s},y_{N-1}^{t})\quad
s,t=0,1,\dots,k-1,\end{aligned}\right.$

and its specific form is

$\displaystyle\Pi u=$
$\displaystyle\sum_{t=0}^{k-1}u(x_{N/2}^{0},y_{N-1}^{t})\theta_{N/2,N-1}^{0,t}+\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}u(x_{N/2-1}^{s},y_{N-1}^{t})\theta_{N/2-1,N-1}^{s,t}$
$\displaystyle-\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}E_{1}(x_{N/2-1}^{s},y_{N-1}^{t})\theta_{N/2-1,N-1}^{s,t}.$

Then, decompose $\int_{K_{a}}b_{1}(\Pi u-u)_{x}v^{N}\mathrm{d}x\mathrm{d}y$ as
the smooth part $\mathcal{B}_{S}$ and the layer part $\mathcal{B}_{E}$:

$\int_{K_{a}}b_{1}(\Pi
u-u)_{x}v^{N}\mathrm{d}x\mathrm{d}y=\int_{K_{a}}b_{1}\left((\Pi
u)_{x}-u_{x}\right)v^{N}\mathrm{d}x\mathrm{d}y=\mathcal{B}_{S}+\mathcal{B}_{E},$

where

$\displaystyle\mathcal{B}_{S}:=$
$\displaystyle\int_{K_{a}}b_{1}\left[\sum_{t=0}^{k-1}S(x_{N/2}^{0},y_{N-1}^{t})(\theta_{N/2,N-1}^{0,t})_{x}\right.$
$\displaystyle\left.+\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}S(x_{N/2-1}^{s},y_{N-1}^{t})(\theta_{N/2-1,N-1}^{s,t})_{x}-S_{x}(x,y)\right]v^{N}\mathrm{d}x\mathrm{d}y,$

and

$\displaystyle\mathcal{B}_{E}:=$
$\displaystyle\int_{K_{a}}b_{1}\left[\sum_{t=0}^{k-1}E(x_{N/2}^{0},y_{N-1}^{t})(\theta_{N/2,N-1}^{0,t})_{x}+\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}E_{2}(x_{N/2-1}^{s},y_{N-1}^{t})(\theta_{N/2-1,N-1}^{s,t})_{x}\right.$
$\displaystyle\left.+\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}E_{12}(x_{N/2-1}^{s},y_{N-1}^{t})(\theta_{N/2-1,N-1}^{s,t})_{x}-E_{x}(x,y)\right]v^{N}\mathrm{d}x\mathrm{d}y.$

After conducting the aforementioned steps, it is evident that there is no
further requirement to address the intractable term (41). Then the primary
difficulty turns into analyzing $\mathcal{B}_{S}$. To obtain the optimal
convergence order of $\mathcal{B}_{S}$, we convert its corresponding estimates
in the two-dimensional case into several types of one-dimensional
interpolation error estimates. The specific conversion is as follows:

The piecewise tensor-product interpolation basis function
$\theta_{i,j}^{s,t}(x,y)$ can be expressed as the product of two univariate
basis functions: $\theta_{i}^{s}(x)\times\theta_{j}^{t}(y)$, where
$\theta_{i}^{s}(x)$ and $\theta_{j}^{t}(y)$ represent the basis function at
point $(x_{i}^{s},y_{j}^{t})$ in the $x$\- and $y$-directions, respectively.
In light of this, we can obtain

$\displaystyle\sum_{t=0}^{k-1}S(x_{N/2}^{0},y_{N-1}^{t})(\theta_{N/2,N-1}^{0,t})_{x}+\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}S(x_{N/2-1}^{s},y_{N-1}^{t})(\theta_{N/2-1,N-1}^{s,t})_{x}-S_{x}(x,y)$
$\displaystyle=$
$\displaystyle\sum_{t=0}^{k-1}\theta_{N-1}^{t}(y)\left(S(x_{N/2}^{0},y_{N-1}^{t})(\theta_{N/2}^{0}(x))_{x}+\sum_{s=0}^{k-1}S(x_{N/2-1}^{s},y_{N-1}^{t})(\theta_{N/2-1}^{s}(x))_{x}-S_{x}(x,y_{N-1}^{t})\right)$
$\displaystyle+\left(\sum_{t=0}^{k-1}\theta_{N-1}^{t}(y)S_{x}(x,y_{N-1}^{t})+\theta_{N}^{0}(y)S_{x}(x,y_{N}^{0})-S_{x}(x,y)\right)-\theta_{N}^{0}(y)S_{x}(x,y_{N}^{0})$
$\displaystyle=$
$\displaystyle\sum_{t=0}^{k-1}\theta_{N-1}^{t}(y)\left(S_{x}^{I}(x,y_{N-1}^{t})-S_{x}(x,y_{N-1}^{t})\right)+\left(S_{x}^{I}(x,y)-S_{x}(x,y)\right)-\theta_{N}^{0}(y)S_{x}(x,y_{N}^{0}),$

here $S_{x}^{I}(x,y_{N-1}^{t})-S_{x}(x,y_{N-1}^{t})$ denotes the one-
dimensional interpolation of variable $x$, and $S_{x}^{I}(x,y)-S_{x}(x,y)$
denotes the one-dimensional interpolation of variable $y$.

The follwing one-dimensional interpolation errors ((Cia:2002-Finite, , Theorem
3.1.4)) will be employed for the next estimates:

(42) $\displaystyle\|v-v^{I}\|_{W^{l,q}}(I_{i})\leq
h_{x,i}^{k+1-l+1/q-1/p}|v|_{W^{k+1,p}(I_{i})}\quad\forall v\in
W^{k+1,p}(I_{i}),$ (43) $\displaystyle\|v-v^{I}\|_{W^{l,q}}(J_{j})\leq
h_{y,j}^{k+1-l+1/q-1/p}|v|_{W^{k+1,p}(J_{j})}\quad\forall v\in
W^{k+1,p}(J_{j}),$

where $I_{i}$ denotes $[x_{i},x_{i+1}]$ with $i=0,1,\dots,N-1$ and $J_{j}$
denotes $[y_{j},y_{j+1}]$ with $j=0,1,\dots,N-1$, $l=0,1$ and $1\leq
p,q\leq\infty$.

Now, we present a detailed demonstration for $\int_{K_{a}}b_{1}(\Pi
u-u)_{x}v^{N}\mathrm{d}x\mathrm{d}y$, which can be decomposed as

(44) $\displaystyle\int_{K_{a}}b_{1}(\Pi u-u)_{x}v^{N}\mathrm{d}x\mathrm{d}y$
$\displaystyle=:$
$\displaystyle\mathcal{I}_{1}+\mathcal{I}_{2}+\mathcal{I}_{3}+\mathcal{I}_{4}+\mathcal{I}_{5}+\mathcal{I}_{6}+\mathcal{I}_{7},$

where

(45)
$\displaystyle\mathcal{I}_{1}=\int_{K_{a}}b_{1}\sum_{t=0}^{k-1}\theta_{N-1}^{t}(y)\left(S_{x}^{I}(x,y_{N-1}^{t})-S_{x}(x,y_{N-1}^{t})\right)v^{N}\mathrm{d}x\mathrm{d}y,$
(46)
$\displaystyle\mathcal{I}_{2}=\int_{K_{a}}b_{1}\left(S_{x}^{I}(x,y)-S_{x}(x,y)\right)v^{N}\mathrm{d}x\mathrm{d}y,$
(47)
$\displaystyle\mathcal{I}_{3}=-\int_{K_{a}}b_{1}\theta_{N}^{0}(y)S_{x}(x,y_{N}^{0})v^{N}\mathrm{d}x\mathrm{d}y,$
(48)
$\displaystyle\mathcal{I}_{4}=\int_{K_{a}}b_{1}\sum_{t=0}^{k-1}E(x_{N/2}^{0},y_{N-1}^{t})(\theta_{N/2,N-1}^{0,t})_{x}v^{N}\mathrm{d}x\mathrm{d}y,$
(49)
$\displaystyle\mathcal{I}_{5}=\int_{K_{a}}b_{1}\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}E_{2}(x_{N/2-1}^{s},y_{N-1}^{t})(\theta_{N/2-1,N-1}^{s,t})_{x}v^{N}\mathrm{d}x\mathrm{d}y,$
(50)
$\displaystyle\mathcal{I}_{6}=\int_{K_{a}}b_{1}\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}E_{12}(x_{N/2-1}^{s},y_{N-1}^{t})(\theta_{N/2-1,N-1}^{s,t})_{x}v^{N}\mathrm{d}x\mathrm{d}y,$
(51)
$\displaystyle\mathcal{I}_{7}=-\int_{K_{a}}b_{1}E_{x}(x,y)v^{N}\mathrm{d}x\mathrm{d}y.$

Estimates for $\mathcal{I}_{1}-\mathcal{I}_{7}$ are discussed in detail below.

From Hölder inequalities, (42), Assumption 1 and Lemma 1, one has

(52) $|\mathcal{I}_{1}|\leq
CN^{-k}\sum_{t=0}^{k-1}\|\frac{\mathrm{d}^{k+1}S(x,y_{N-1}^{t})}{\mathrm{d}x^{k+1}}\|_{K_{a}}\|v^{N}\|_{K_{a}}\leq
CN^{-(k+1)}\|v^{N}\|_{\varepsilon}.$

Hölder inequalities, (43), Assumption 1 and Lemma 1 give

(53) $|\mathcal{I}_{2}|\leq
CN^{-(k+1)}\|\frac{\mathrm{d}^{k+1}S_{x}(x,y)}{\mathrm{d}y^{k+1}}\|_{K_{a}}\|v^{N}\|_{K_{a}}\leq
CN^{-(k+2)}\|v^{N}\|_{\varepsilon}.$

Green’s formula yields

$\displaystyle\mathcal{I}_{3}=$
$\displaystyle\int_{K_{a}}((b_{1})_{x}v^{N}+b_{1}v_{x}^{N})\theta_{N}^{0}(y)S(x,y_{N}^{0})\mathrm{d}x\mathrm{d}y-\int_{y_{N-1}}^{y_{N}}b_{1}\theta_{N}^{0}(y)S(x_{N/2},y_{N}^{0})v^{N}(x_{N/2},y)\mathrm{d}y$
$\displaystyle+\int_{y_{N-1}}^{y_{N}}b_{1}\theta_{N}^{0}(y)S(x_{N/2-1},y_{N}^{0})v^{N}(x_{N/2-1},y)\mathrm{d}y.$

Since $u=0$ on the boundary $\partial\Omega$, the smooth function
$S(x,y_{N}^{0})$ can be seen as the layer function $-E(x,y_{N}^{0})$, which
holds

$\|S(x,y_{N}^{0})\|_{K_{a}}=\|E(x,y_{N}^{0})\|_{K_{a}}\leq
C\varepsilon^{1/2}N^{-\sigma}.$

Thus,

(54)
$|\int_{K_{a}}((b_{1})_{x}v^{N}+b_{1}v_{x}^{N})\theta_{N}^{0}(y)S(x,y_{N}^{0})\mathrm{d}x\mathrm{d}y|\leq
C\|S(x,y_{N}^{0})\|_{K_{a}}\|v_{x}\|_{K_{a}}\leq
CN^{-\sigma}\|v^{N}\|_{\varepsilon}.$

Besides,

(55)
$\displaystyle|\int_{y_{N-1}}^{y_{N}}b_{1}\theta_{N}^{0}(y)S(x_{N/2},y_{N}^{0})v^{N}(x_{N/2},y)\mathrm{d}y|+|\int_{y_{N-1}}^{y_{N}}b_{1}\theta_{N}^{0}(y)S(x_{N/2-1},y_{N}^{0})v^{N}(x_{N/2-1},y)\mathrm{d}y|$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{\sigma}\|v^{N}\|_{\varepsilon}+CN^{-\sigma-1/2}\ln^{1/2}N\|v^{N}\|_{\varepsilon},$

where we have used $|S(x_{N/2},y_{N}^{0})|=|-E(x_{N/2},y_{N}^{0})|\leq
C\varepsilon^{\sigma}$,
$|S(x_{N/2-1},y_{N}^{0})|=|-E(x_{N/2-1},y_{N}^{0})|\leq CN^{-\sigma}$ and

(56) $\displaystyle\begin{aligned}
|\int_{y_{N-1}}^{y_{N}}v^{N}(x_{N/2},y)\mathrm{d}y|\leq
CN\|v^{N}\|_{L^{1}(K_{N/2,N-1})}\leq C\|v^{N}\|_{\varepsilon},\end{aligned}$
$\displaystyle\begin{aligned}
|\int_{y_{N-1}}^{y_{N}}v^{N}(x_{N/2-1},y)\mathrm{d}y|&\leq
C|\int_{y_{N-1}}^{y_{N}}\int_{x_{0}}^{x_{N/2-1}}v_{x}^{N}(x,y)\mathrm{d}x\mathrm{d}y|\\\
&\leq C\|v_{x}^{N}\|_{L^{1}([x_{0},x_{N/2-1}]\times[y_{N-1},y_{N}])}\\\ &\leq
CN^{-1/2}\ln^{1/2}N\|v^{N}\|_{\varepsilon}.\end{aligned}$

(54)-(55) give

(57) $|\mathcal{I}_{3}|\leq CN^{-\sigma}\|v^{N}\|_{\varepsilon}.$

By $|E(x_{N/2},y_{N-1}^{t})|\leq C\varepsilon^{\sigma}$,
$|E_{2}(x_{N/2-1}^{s},y_{N-1}^{t})|\leq C\varepsilon^{\sigma}$,
$|E_{12}(x_{N/2-1}^{s},y_{N-1}^{t})|\leq C\varepsilon^{\sigma}N^{-\sigma}$,
inverse inequalities and Lemma 1, we have

(58) $\displaystyle|\mathcal{I}_{4}|\leq
C\varepsilon^{\sigma}(h_{x,N/2-1}^{-1}h_{x,N/2-1}^{1/2}h_{y,N-1}^{1/2})\|v^{N}\|_{K_{a}}\leq
C\varepsilon^{\sigma-1/2}N^{-1/2}\|v^{N}\|_{\varepsilon},$ (59)
$\displaystyle|\mathcal{I}_{5}|\leq
C\varepsilon^{\sigma}(h_{x,N/2-1}^{-1}h_{x,N/2-1}^{1/2}h_{y,N-1}^{1/2})\|v^{N}\|_{K_{a}}\leq
C\varepsilon^{\sigma-1/2}N^{-1/2}\|v^{N}\|_{\varepsilon},$ (60)
$\displaystyle|\mathcal{I}_{6}|\leq
C\varepsilon^{\sigma}N^{-\sigma}(h_{x,N/2-1}^{-1}h_{x,N/2-1}^{1/2}h_{y,N-1}^{1/2})\|v^{N}\|_{K_{a}}\leq
C\varepsilon^{\sigma-1/2}N^{-\sigma-1/2}\|v^{N}\|_{\varepsilon}.$

Green’s formula yields

$\displaystyle\mathcal{I}_{7}=$
$\displaystyle\int_{K_{a}}((b_{1})_{x}v^{N}+b_{1}v_{x}^{N})E(x,y)\mathrm{d}x\mathrm{d}y-\int_{y_{N-1}}^{y_{N}}b_{1}E(x_{N/2},y)v^{N}(x_{N/2},y)\mathrm{d}y$
$\displaystyle+\int_{y_{N-1}}^{y_{N}}b_{1}E(x_{N/2-1},y)v^{N}(x_{N/2-1},y)\mathrm{d}y,$

where

(61)
$|\int_{K_{a}}((b_{1})_{x}v^{N}+b_{1}v_{x}^{N})E(x,y)\mathrm{d}x\mathrm{d}y|\leq
CN^{-\sigma}\|v^{N}\|_{\varepsilon},$

and

(62)
$\displaystyle|\int_{y_{N-1}}^{y_{N}}b_{1}E(x_{N/2},y)v^{N}(x_{N/2},y)\mathrm{d}y|+|\int_{y_{N-1}}^{y_{N}}b_{1}E(x_{N/2-1},y)v^{N}(x_{N/2-1},y)\mathrm{d}y|$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{\sigma}|\int_{y_{N-1}}^{y_{N}}v^{N}(x_{N/2},y)\mathrm{d}y|+CN^{-\sigma}|\int_{y_{N-1}}^{y_{N}}v^{N}(x_{N/2-1},y)\mathrm{d}y|$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{\sigma}\|v^{N}\|_{\varepsilon}+CN^{-\sigma-1/2}\ln^{1/2}N\|v^{N}\|_{\varepsilon}.$

(61)-(62) give

(63) $|\mathcal{I}_{7}|\leq CN^{-\sigma}\|v^{N}\|_{\varepsilon}.$

Substituting (52)-(53), (57), (58)-(60) and (63) into (44), we prove

(64) $|\int_{K_{a}}b_{1}(\Pi u-u)_{x}v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1)}\|v^{N}\|_{\varepsilon}+C\varepsilon^{\sigma-1/2}N^{-1/2}\|v^{N}\|_{\varepsilon}.$

Next, we are going to deal with $\int_{K_{a}}b_{2}(\Pi
u-u)_{y}v^{N}\mathrm{d}x\mathrm{d}y$, which can be divided into

(65) $\displaystyle\int_{K_{a}}b_{2}(\Pi u-u)_{y}v^{N}\mathrm{d}x\mathrm{d}y$
$\displaystyle=$
$\displaystyle\int_{K_{a}}b_{2}(S^{I}-S)_{y}v^{N}\mathrm{d}x\mathrm{d}y+\int_{K_{a}}b_{2}(\pi_{1}E_{1}-E_{1})_{y}v^{N}\mathrm{d}x\mathrm{d}y+\int_{K_{a}}b_{2}(E_{2}^{I}-E_{2})_{y}v^{N}\mathrm{d}x\mathrm{d}y$
$\displaystyle+\int_{K_{a}}b_{2}(E_{12}^{I}-E_{12})_{y}v^{N}\mathrm{d}x\mathrm{d}y.$

(66) $|\int_{K_{a}}b_{2}(S^{I}-S)_{y}v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1)}\|v^{N}\|_{\varepsilon},$

is straightforward by Hölder inequalities, Lemma 2, Assumption 1 and Lemma 1.

From Hölder inequalities, inverse inequalities, Lemma 1 and (29), we have

(67)
$|\int_{K_{a}}b_{2}(\pi_{1}E_{1}-E_{1})_{y}v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN\|E_{1}\|_{\infty,K_{a}}\|v^{N}\|_{L^{1}(K_{a})}\leq
CN^{-\sigma}\|v^{N}\|_{\varepsilon}.$

Under a similar analysis as
$\int_{K_{a}}b_{2}(\pi_{1}E_{1}-E_{1})_{y}v^{N}\mathrm{d}x\mathrm{d}y$, we can
get

(68)
$\displaystyle|\int_{K_{a}}b_{2}(E_{2}^{I}-E_{2})_{y}v^{N}\mathrm{d}x\mathrm{d}y|\leq
C\varepsilon^{\sigma}\|v^{N}\|_{\varepsilon},$ (69)
$\displaystyle|\int_{K_{a}}b_{2}(E_{12}^{I}-E_{12})_{y}v^{N}\mathrm{d}x\mathrm{d}y|\leq
C\varepsilon^{\sigma}N^{-\sigma}\|v^{N}\|_{\varepsilon}.$

Substituting (66)-(69) into (65), we prove

(70) $|\int_{K_{a}}b_{2}(\Pi u-u)_{y}v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1)}\|v^{N}\|_{\varepsilon}.$

In conclusion, combined (64) with (70), we obtain

(71) $|III|=|-\int_{K_{a}}\boldsymbol{b}\cdot\nabla(\Pi
u-u)v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1)}\|v^{N}\|_{\varepsilon}+C\varepsilon^{\sigma-1/2}N^{-1/2}\|v^{N}\|_{\varepsilon}.$

Derivations of IV is similar to that of III, thus we have

(72) $|IV|=|-\int_{K_{b}}\boldsymbol{b}\cdot\nabla(\Pi
u-u)v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1)}\|v^{N}\|_{\varepsilon}+C\varepsilon^{\sigma-1/2}N^{-1/2}\|v^{N}\|_{\varepsilon}.$

###### Remark 2

We briefly summarize the novel analysis technique. Firstly, the specific form
of the interpolation $\Pi u$ on $K_{a}$ is described. Then, decompose
$\int_{K_{a}}b_{1}(\Pi u-u)_{x}v^{N}\mathrm{d}x\mathrm{d}y$ as the smooth part
$\mathcal{B}_{S}$ and the layer part $\mathcal{B}_{E}$. For $\mathcal{B}_{S}$,
we convert the corresponding estimates in the two-dimensional case into
several types of one-dimensional interpolation error estimates, see (45)-(46)
and (52)-(53). Furthermore, we take full advantage of the boundary condition
where the sum of the smooth function and the layer function is zero, so that
analysis of the smooth function on the boundary can be transformed into
studying the layer function, see (47) and (54)-(55). For $\mathcal{B}_{E}$, we
employ the standard error estimats, see (48)-(51) and (58)-(63).

###### Lemma 6

Let Assumptions 1 and 2 hold. Let $\sigma\geq k+1$. Then we have

$|V+VI|\leq CN^{-k}\|v^{N}\|_{\varepsilon}.$

###### Proof

To begin with, we discuss term
$\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{1}E_{1}-E_{1})v^{N}\mathrm{d}x\mathrm{d}y$.
Green’s formula generates

(73)
$\displaystyle\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{1}E_{1}-E_{1})v^{N}\mathrm{d}x\mathrm{d}y$
$\displaystyle=$
$\displaystyle-\int_{\Omega_{0}}(\nabla\cdot\boldsymbol{b})(\pi_{1}E_{1}-E_{1})v^{N}\mathrm{d}x\mathrm{d}y-\int_{\Omega_{0}}(\pi_{1}E_{1}-E_{1})(\boldsymbol{b}\cdot\nabla
v^{N})\mathrm{d}x\mathrm{d}y$
$\displaystyle+\int_{\partial\Omega_{0}}(\boldsymbol{b}\cdot\boldsymbol{n})(\pi_{1}E_{1}-E_{1})v^{N}\mathrm{d}s.$

Employing Hölder inequalities and (33), we have

(74)
$|-\int_{\Omega_{0}}(\nabla\cdot\boldsymbol{b})(\pi_{1}E_{1}-E_{1})v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1)}\|v^{N}\|_{\varepsilon}.$

The reader is referred to (zhang1Liu2:2023-Convergence, , Lemma 4.1) for
detailed derivations of
$\int_{\Omega_{0}}(\pi_{1}E_{1}-E_{1})(\boldsymbol{b}\cdot\nabla
v^{N})\mathrm{d}x\mathrm{d}y$, from which one can get

(75) $|\int_{\Omega_{0}}(\pi_{1}E_{1}-E_{1})(\boldsymbol{b}\cdot\nabla
v^{N})\mathrm{d}x\mathrm{d}y|\leq CN^{-(k+1/2)}\|v^{N}\|_{\varepsilon}.$

Then take
$\int_{\partial\Omega_{0}}(\boldsymbol{b}\cdot\boldsymbol{n})(\pi_{1}E_{1}-E_{1})v^{N}\mathrm{d}s$
into consideration:

(76)
$\displaystyle|\int_{\partial\Omega_{0}}(\boldsymbol{b}\cdot\boldsymbol{n})(\pi_{1}E_{1}-E_{1})v^{N}\mathrm{d}s|$
$\displaystyle\leq$ $\displaystyle
C\|\pi_{1}E_{1}-E_{1}\|_{L^{\infty}(\partial\Omega_{0}\textbackslash\partial\Omega)}\|v^{N}\|_{L^{1}(\partial\Omega_{0}\textbackslash\partial\Omega)}$
$\displaystyle\leq$ $\displaystyle
C\|E_{1}\|_{L^{\infty}(\partial\Omega_{0}\textbackslash\partial\Omega)}\|v^{N}\|_{L^{1}(\partial\Omega_{0}\textbackslash\partial\Omega)}$
$\displaystyle\leq$ $\displaystyle
CN^{-\sigma}\left(|\int_{y_{N/2-1}}^{y_{N/2}}v^{N}(x_{N-1},y)\mathrm{d}y|+|\int_{y_{N-1}}^{1}v^{N}(x_{N/2},y)\mathrm{d}y|+|\int_{y_{N-1}}^{1}v^{N}(x_{N/2-1},y)\mathrm{d}y|\right)$
$\displaystyle+CN^{-\sigma}\left(|\int_{x_{N/2-1}}^{x_{N/2}}v^{N}(x,y_{N-1})\mathrm{d}x|+|\int_{x_{N-1}}^{1}v^{N}(x,y_{N/2})\mathrm{d}x|+|\int_{x_{N-1}}^{1}v^{N}(x,y_{N/2-1})\mathrm{d}x|\right)$
$\displaystyle\leq$ $\displaystyle CN^{-\sigma}\|v^{N}\|_{\varepsilon},$

here we have used Höler inequalities, (29), Lemma 1 and (56).

From (74)-(76), one can get

(77)
$|\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{1}E_{1}-E_{1})v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1/2)}\|v^{N}\|_{\varepsilon}.$

Similarly,

(78)
$|\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{2}E_{2}-E_{2})v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1/2)}\|v^{N}\|_{\varepsilon}.$

Next, we study
$\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{12}E_{12}-E_{12})v^{N}\mathrm{d}x\mathrm{d}y$,
which can be divided into the follwing terms using Green’s formula:

(79)
$\displaystyle\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{12}E_{12}-E_{12})v^{N}\mathrm{d}x\mathrm{d}y$
$\displaystyle=$
$\displaystyle-\int_{\Omega_{0}}(\nabla\cdot\boldsymbol{b})(\pi_{12}E_{12}-E_{12})v^{N}\mathrm{d}x\mathrm{d}y-\int_{\Omega_{0}}(\pi_{12}E_{12}-E_{12})(\boldsymbol{b}\cdot\nabla
v^{N})\mathrm{d}x\mathrm{d}y$
$\displaystyle+\int_{\partial\Omega_{0}}(\boldsymbol{b}\cdot\boldsymbol{n})(\pi_{12}E_{12}-E_{12})v^{N}\mathrm{d}s.$

According to Hölder inequalities and (35), one obtains

(80)
$|-\int_{\Omega_{0}}(\nabla\cdot\boldsymbol{b})(\pi_{12}E_{12}-E_{12})v^{N}\mathrm{d}x\mathrm{d}y|\leq(C\varepsilon
N^{-k}+C\varepsilon^{1/2}N^{-(k+1)}+CN^{-(2\sigma+1)})\|v^{N}\|_{\varepsilon}.$

Refer to (zhang1Liu2:2023-Convergence, , Lemma 4.2) for detailed
demonstrations of the second term in the right-hand side of (79):

(81) $|\int_{\Omega_{0}}(\pi_{12}E_{12}-E_{12})(\boldsymbol{b}\cdot\nabla
v^{N})\mathrm{d}x\mathrm{d}y|\leq CN^{-(k+1/2)}\|v^{N}\|_{\varepsilon}.$

With an approach similar to (76), we can derive that

(82)
$|\int_{\partial\Omega_{0}}(\boldsymbol{b}\cdot\boldsymbol{n})(\pi_{12}E_{12}-E_{12})v^{N}\mathrm{d}s|\leq
CN^{-\sigma}\|v^{N}\|_{\varepsilon}.$

Substituting (80)-(82) into (79), one can obtain

(83)
$|\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{12}E_{12}-E_{12})v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1/2)}\|v^{N}\|_{\varepsilon}.$

(83), combined with (77) and (78), proves

(84)
$|V|=|-\sum_{i=1,2,12}\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(\pi_{i}E_{i}-E_{i})v^{N}\mathrm{d}x\mathrm{d}y|\leq
CN^{-(k+1/2)}\|v^{N}\|_{\varepsilon}.$

At last, Hölder inequalities, Lemma 2, Lemma 1 and Assumption 1 yield

(85)
$|VI|=|-\int_{\Omega_{0}}\boldsymbol{b}\cdot\nabla(S^{I}-S)v^{N}\mathrm{d}x\mathrm{d}y|\leq
C\|\nabla(S^{I}-S)\|_{\Omega_{0}}\|v^{N}\|_{\Omega_{0}}\leq
CN^{-k}\|v^{N}\|_{\varepsilon}.$

Thus we are done.

Now we are in a position to convey the main conclusion of this article.

###### Theorem 4.1

Let Assumptions 1 and 2 hold. Let $u$ denote the exact solution to (1) and
$u^{N}$ denote the corresponding finite element solution to (15). Then we have

$\|u-u^{N}\|_{\varepsilon}\leq CN^{-k}.$

###### Proof

Triangle inequality generates

$\|u-u^{N}\|_{\varepsilon}\leq\|u-\Pi u\|_{\varepsilon}+\|\Pi
u-u^{N}\|_{\varepsilon},$

where from (38), one has

(86) $\|u-\Pi u\|_{\varepsilon}\leq CN^{-k},$

and from Lemmas 4, 5 and 6, one has

(87) $\|\Pi u-u^{N}\|_{\varepsilon}\leq CN^{-k}.$

Thus we are done.

##  5 Numerical experiments

Numerical experiments are carried out in this section to confirm our
theoretical results. All calculations were performed using Intel Visual
Fortran 11, and discrete problems were solved with the aid of the nonsymmetric
iterative solver GMRES; see, e.g., Ben1Gol2:2005-Numerical .

Consider the following singularly perturbed convection-diffusion equation:

(88) $\displaystyle-\varepsilon\Delta u-(2+x-y)u_{x}-(2-x+y)u_{y}+2u=$
$\displaystyle f\quad$ $\displaystyle\text{in}\quad\Omega=(0,1)^{2},$
$\displaystyle u=$ $\displaystyle 0\quad$
$\displaystyle\text{on}\quad\partial\Omega,$

where we select an appropriate function $f$ such that

(89) $u(x,y)=2\sin(\pi
x)(1-e^{-\frac{2x}{\varepsilon}})(1-y)^{2}(1-e^{-\frac{y}{\varepsilon}})$

is the exact solution to (1). This solution typically exhibits exponential
layers, as stated in Assumption 1.

Errors for $k=1$ are listed in Table 1 under the energy norm
$\|u-u^{N}\|_{\varepsilon}$, for $\varepsilon=10^{-4},10^{-5},\dots,10^{-8}$
and $N=8,16,32,64,128,256$; Errors for $k=2$ are listed in Table 2 under the
energy norm $\|u-u^{N}\|_{\varepsilon}$, for
$\varepsilon=10^{-4},10^{-5},\dots,10^{-8}$ and $N=8,16,32,64,128$. These data
indicate that the finite element solution converges uniformly to the exact
solution at an optimal rate of order $k+1$ under the energy norm, confirming
our main conclusion, i.e., Theorem 4.1.

Table 1: Errors of $\|u-u^{N}\|_{\varepsilon}$ and convergence order for $k=1$ $\varepsilon$ | $N$
---|---
| 8 | 16 | 32 | 64 | 128 | 256
$10^{-4}$ | 0.339E+00 | 0.167E+00 | 0.834E-01 | 0.418E-01 | 0.208E-01 | 0.104E-01
| 1.02 | 1.00 | 1.00 | 1.00 | 1.00 | —
$10^{-5}$ | 0.339E+00 | 0.167E+00 | 0.834E-01 | 0.418E-01 | 0.208E-01 | 0.104E-01
| 1.02 | 1.00 | 1.00 | 1.00 | 1.00 | —
$10^{-6}$ | 0.339E+00 | 0.167E+00 | 0.834E-01 | 0.418E-01 | 0.208E-01 | 0.104E-01
| 1.02 | 1.00 | 1.00 | 1.00 | 1.00 | —
$10^{-7}$ | 0.339E+00 | 0.167E+00 | 0.834E-01 | 0.418E-01 | 0.208E-01 | 0.104E-01
| 1.02 | 1.00 | 1.00 | 1.00 | 1.00 | —
$10^{-8}$ | 0.339E+00 | 0.167E+00 | 0.834E-01 | 0.418E-01 | 0.208E-01 | 0.104E-01
| 1.02 | 1.00 | 1.00 | 1.00 | 1.00 | —
Table 2: Errors of $\|u-u^{N}\|_{\varepsilon}$ and convergence order for $k=2$ $\varepsilon$ | $N$
---|---
| 8 | 16 | 32 | 64 | 128
$10^{-4}$ | 0.103E+00 | 0.257E-01 | 0.643E-02 | 0.162E-02 | 0.470E-03
| 2.00 | 2.00 | 1.99 | 1.79 | —
$10^{-5}$ | 0.103E+00 | 0.257E-01 | 0.643E-02 | 0.161E-02 | 0.402E-03
| 2.00 | 2.00 | 2.00 | 2.00 | —
$10^{-6}$ | 0.103E+00 | 0.257E-01 | 0.643E-02 | 0.161E-02 | 0.402E-03
| 2.00 | 2.00 | 2.00 | 2.00 | —
$10^{-7}$ | 0.103E+00 | 0.257E-01 | 0.643E-02 | 0.161E-02 | 0.402E-03
| 2.00 | 2.00 | 2.00 | 2.00 | —
$10^{-8}$ | 0.103E+00 | 0.257E-01 | 0.643E-02 | 0.161E-02 | 0.402E-03
| 2.00 | 2.00 | 2.00 | 2.00 | —

##  6 Conflict of interest statement

We declare that we have no conflict of interest.

## References

  * (1) T. Apel.  Anisotropic Finite Elements: Local Estimates and Applications.  Teubner Stuttgart, 1999.
  * (2) N. S. Bakhvalov.  On the optimization of the methods for solving boundary value problems in the presence of a boundary layer.  Zh. Vychisl. Mat. Mat. Fiz., 9(4):841–859, 1969.
  * (3) M. Benzi, G. H Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numer., 14:1–137, 2005.
  * (4) M. Brdar and H. Zarin.  A singularly perturbed problem with two parameters on a Bakhvalov-type mesh.  J. Comput. Appl. Math., 292:307–319, 2016\.
  * (5) S. C Brenner and L R. Scott.  The Mathematical Theory of Finite Element Methods, volume 3.  Springer, 2008.
  * (6) P. G Ciarlet.  Finite Element Method for Elliptic Problems.  SIAM, 2002.
  * (7) J. Li and I. M. Navon.  Uniformly convergent finite element methods for singularly perturbed elliptic boundary value problems: reaction-diffusion type.  Comput. Math. Appl., 35(3):57–70, 1998.
  * (8) 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.
  * (9) T. Linß.  Uniform superconvergence of a Galerkin finite element method on Shishkin-type meshes.  Numer. Methods for Partial Differential Equations, 16(5):426–440, 2000.
  * (10) T. Linß.  Layer-Adapted Meshes for Reaction-Convection-Diffusion Problems.  Springer, 2009.
  * (11) T. Linß and M. Stynes.  Asymptotic analysis and Shishkin-type decomposition for an elliptic convection-diffusion problem.  J. Math. Anal. Appl., 261(2):604–632, 2001.
  * (12) X. Liu and J. Zhang.  Uniform supercloseness of Galerkin finite element method for convection-diffusion problems with characteristic layers.  Comput. Math. Appl., 75(2):444–458, 2018\.
  * (13) X. Liu and J. Zhang.  Uniform convergence of optimal order for a finite element method on a Bakhvalov-type mesh for a singularly perturbed convection-diffusion equation with parabolic layers.  Numer. Algorithms, pages 1–20, 2023.
  * (14) H.-G. Roos and M. Schopf  Convergence and stability in balanced norms of finite element methods on Shishkin meshes for reaction-diffusion problems.  ZAMM J. Appl. Math., 95(6):551–565, 2015.
  * (15) H.-G. Roos.  Error estimates for linear finite elements on Bakhvalov-type meshes.  Appl. Math., 51:63–72, 2006.
  * (16) H.-G. Roos and M. Stynes.  Some open questions in the numerical analysis of singularly perturbed differential equations.  Comput. Methods Appl. Math., 15(4):531–550, 2015\.
  * (17) H.-G. Roos, M. Stynes, and L. Tobiska.  Robust Numerical Methods for Singularly Perturbed Differential Equations, volume 24.  Springer Science & Business Media, 2008.
  * (18) M. Stynes.  Steady-state convection-diffusion problems.  Acta Numer., 14:445–508, 2005.
  * (19) J. Zhang and X. Liu.  Convergence of a finite element method on a Bakhvalov-type mesh for singularly perturbed reaction-diffusion equation.  Appl. Math. Comput., 385:125403, 2020.
  * (20) J. Zhang and X. Liu.  Optimal order of uniform convergence for finite element method on Bakhvalov-type meshes.  J. Sci. Comput., 85(1):2, 2020.
  * (21) J. Zhang and X. Liu.  Convergence and supercloseness in a balanced norm of finite element methods on Bakhvalov-type meshes for reaction-diffusion problems.  J. Sci. Comput., 88(1):27, 2021.
  * (22) J. Zhang and X. Liu.  Convergence of a finite element method on a Bakhvalov-type mesh for a singularly perturbed convection-diffusion equation in 2D.  Numer. Methods for Partial Differential Equations, 39(2):1201–1219, 2023.
  * (23) J. Zhang, X. Liu, and M. Yang.  Optimal order $L^{2}$ error estimate of SDFEM on Shishkin triangular meshes for singularly perturbed convection-diffusion equations.  SIAM J. Numer. Anal., 54(4):2060–2080, 2016.
  * (24) J. Zhang and Y. Lv.  High-order finite element method on a Bakhvalov-type mesh for a singularly perturbed convection-diffusion problem with two parameters.  Appl. Math. Comput., 397:125953, 2021.
  * (25) Z. Zhang.  Finite element superconvergence on shishkin mesh for 2-D convection-diffusion problems.  Math. Comput., 72(243):1147–1177, 2003.
  * (26) H.-G. Roos and M. Schopf.  Analysis of finite element methods on Bakhvalov-type meshes for linear convection-diffusion problems in 2D.  Appl. Math., 57:97-108, 2012.