# Enabling On-Demand Cyber-Physical Control Applications with UAV Access
Points

Igor Donevski, Jimmy Jessen Nielsen This work has been submitted to IEEE for
possible publication. Copyright may be transferred without notice, after which
this version may no longer be accessible. Department of Electronic Systems,
Aalborg University, Denmark  e-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Achieving cyber-physical control over a wireless channel requires satisfying
both the timeliness of a single packet and preserving the latency reliability
across several consecutive packets. To satisfy those requirements as an
ubiquitous service requires big infrastructural developments, or flexible on-
demand equipment such as UAVs. To avoid the upfront cost in terms of finance
and energy, this paper analyzes the capability of UAV access points (UAVAPs)
to satisfy the requirements for cyber-physical traffic. To investigate this,
we perform a Gilbert-Eliott burst-error analysis that is analytically derived
as a combination of two separate latency measurement campaigns and provide an
upper-bound analysis of the UAVAP system. The analysis is centered around a
UAVAP that uses its LTE connection to reach the backhaul, while providing
service to ground nodes (GNs) with a Wi-Fi access point (AP). Thus, we combine
both measurement campaigns to analyze the plausibility of the described setup
in casual, crowded or mixed network settings.

###### Index Terms:

Unmanned Aerial Vehicle, Drones, Access Point, Low Latency, Survival Time,
Burst Errors

##  I Introduction

Connected robotics is a prime driver in the development of the sixth
generation (6G) of mobile networks [1]. This comes as an extension to one of
the main service expected from the fifth generation of mobile networks (5G)
named Ultra-Reliable and Low-Latency Communication (URLLC) [2]. Offering URLLC
is difficult due to the combined requirement of achieving latency in the order
of few milliseconds with reliability of around $1-10^{-5}$. However, the URLLC
use case is concerned with latency-reliability of each separate packet which
does not capture the operation of multi-packet operations such as cyber-
physical control applications. These types of applications transmit regular
updates that are not gravely affected by a single untimely packet (a packet
that violates the latency requirements). Therefore, cyber-control applications
are dependent on the time on which they can operate nominally without having
received an eligible packet, referred to as _survival time_ [3]. Therefore, in
this paper we focus on analyzing experimentally derived insights on the
potential application of Unmanned Aerial Vehicles (UAVs), colloquially
referred to as drones. The analysis is targeted at improving communication
system reliability through reducing the likelihood of error bursts such as
consecutive untimely packets for ground nodes (GNs).

###  I-A Goals & Motivation

Figure 1: An illustration of the scenario investigating the probability of
lengthy error bursts in the use if a UAVAP.

UAV Access Points (UAVAPs) are a perfect adaptable solution for addressing
gaps in wireless communication coverage, mainly as a benefit to their mobility
as a modular mobile service equipment [4]. This aims to save large quantities
of money and energy through avoiding constant operation of Base Stations for
this particular service [5, 6]. Moreover, the adaptable horizontal position
and altitude, improve the chances of having a Line of Sight both towards base
station (BS) towers mounted on top of built up terrain, as well as GNs [7, 8].
Thus, UAVAPs are an intuitive solution for providing reliable service to GNs
with low latencies on both the GN access link (GNAL) and BS access link
(BSAL). Therefore, this work intends to address an interest in modeling the
error burst probability for a full stack of BSAL and GNAL, where both
connections are offered by widely available technologies such as LTE and Wi-Fi
respectively.

Despite the high likelihood for LoS on both BSAL and GNAL, both interfaces
might provide significant delays when the network conditions get crowded, and
channel competition is high. Thus, the goal is to frame the issue of an UAVAP
enabling ground based closed-loop end-to-end control applications in different
network conditions, by focusing on two key questions. Does UAVAP provided
wireless access offer improved latency-reliability in relaxed and crowded GNAL
and/or BSAL? What are the shortcomings of the UAVAP wireless access and how
can those be addressed?

The paper organization is done in the following manner. Section II discusses
significant works in the fields of cyber-physical control, burst error
analysis, and experimental UAVAP systems. Section III introduces the model for
burst errors in UAVAP systems. Section IV showcases statistical properties of
the considered measurement campaigns. Section V showcases the burst error
performance of a UAVAP system. Finally, Section VI summarizes the results of
the analysis, addresses concerns, and discusses possible future directions.

##  II Related Work

The experimental and theoretical communities are generally concerned with
several complementary problems of UAVAP systems, particularly for addressing
reliability concerns. The work of [9] does a throughput investigation for an
autonomously repositioning UAV in a highly controlled scenario. We previously
investigated how the UAV position can be exploited when having a fixed
coverage area [8], and upgraded the analysis for wireless network slicing
[10]. Nonetheless, analyses that focus on the positioning of a UAVAP are
generally concerned with optimizing the wireless conditions of the GNAL, and
rely on some strong assumptions for the BSAL [6].

There are experimental efforts to address some more common problems that stem
from the design of cellular systems for UAVAPs. As covered in the overview in
[11], the two central issues that appear in such systems are the higher inter-
cell interference from distant BSs and the issue of antenna propagation design
that is tilted toward the ground users. Although these issues add to lessen
the reliability of the system, they introduce frequent handover between the
close and distant BSs, and are covered in several different measurement
campaigns, and are well elaborated in the following works [12, 13, 14].
However, there was a lack of experimentally derived latency for UAVAPs, which
was thus covered in our previous work [15].

The survival time of cyber-physical control applications varies depending on
the application and can be generally inferred from the report in [3]. In
accord, the work of [16] observes that the safety of autonomously guided road
vehicle systems is directly reliant on the number of consecutive errors. Our
previous work in [17] analyzed the efficiency of multi-connectivity for
servicing cyber-physical control traffic. As such, this paper addresses a gap
in the literature for burst error analysis for UAVAPs by building a model for
analyzing temporal correlations in each interface, for two different sets of
network conditions.

##  III Burst Errors in UAVAP systems

In an end-to-end link between a GN and a central controller, data is
transferred periodically with an occurrence of $T_{s}$. Packet timeliness must
be satisfied (packet latency must not exceed some value $\theta\leq T_{s}$ )
with the goal of keeping all actions taken by the robot fresh and not
outdated. Therefore, just like in the URLLC use case we look at the latency-
reliability function (a CDF of the system latency) [18]. Moreover, since our
communication occurs on two different links we identify the total UAVAP
latency limit as the sum of both $\sum_{i}\theta_{i}<\theta$, where
$i\in\\{\text{Wi-Fi},\text{LTE}\\}$ for the GNAL and BSAL respectively.

The total latency budget $\theta$ needs to be split for each interface. For a
latency of $l_{i}$ on interface $i$, and its latency deadline $\theta_{i}$ we
define the latency-reliability function as:

$F_{i}(\theta_{i})=\Pr(l_{i}\leq\theta_{i}).$ (1)

Thus, we define a packet error in contrary to traditional cases where
$\theta\rightarrow\infty$ but as the case where a single of both interfaces is
untimely as in:

$P_{e}^{(i)}=1-F_{i}(\theta_{i}).$ (2)

Treating each link individually like this, means that we judge the UAVAP
system by its upper bound of performance where neither of the links has
violated its timeliness, regardless of the final collective $\theta$. This
allows for bad states of each interface to be identified quickly and treated
individually in the two hops of communications that need to be nominal for a
packet to be timely.

Therefore, each untimely transmission on each interface of the UAVAP can be
considered as a transition to a failed state for that link of the hop. In this
way, the complexity of the collective wireless conditions such as congestion,
small/fast fading, noise and interference are collected and simplified to two
states ( a _good_ state G and a _bad_ state B, as illustrated in Fig. 2) of a
discrete time Markov chain. This model is the Gilbert-Elliott (GE) burst error
model that has two parameters $p_{i}$ and $r_{i}$ that are the probabilities
to transition from the good to the bad state and from the bad state to the
good one, respectively [17]. Thus to model the GE system accurately, we
require statistically relevant data that define the state transitions in (3)
accurately.

Figure 2: The GE system for modeling burst errors.

$\mathbf{P}_{i}=\begin{bmatrix}1-p_{i}&p_{i}\\\ r_{i}&1-r_{i}\end{bmatrix},$
(3)

with steady states:

$\displaystyle\pi_{i,\text{G}}=\frac{r_{i}}{p_{i}+r_{i}},$ (4)
$\displaystyle\pi_{i,\text{B}}=\frac{p_{i}}{r_{i}+p_{i}},$ (5)

and total probability of the system must satisfy
$\pi_{i,\text{G}}+\pi_{i,\text{B}}=1$.

Therefore, given the GE model for each interface, we can convert the total
end-to-end link as a four state GE model as shown in Fig. 3. Here, each
interface moving over to a bad state results into a bad overall state and
compromises the total link.

Figure 3: The adapted GE system for modeling burst errors across a UAVAP with
independent GNAL and BSAL.

##  IV Measurement Campaigns

In this section we showcase and compare two different testing campaigns on Wi-
Fi and LTE latencies. The two sets of probabilities that are contained in the
subsections that follow are taken as, representative of direct crowded links
in an NLoS scenario (Subsection IV-A), relaxed radio conditions in a purely
LoS scenario (Subsection IV-B), and relaxed radio conditions for a direct GN
to BS link (Subsection IV-C).

###  IV-A Crowded Scenario

Figure 4: Empirical latency CDFs for Relaxed, Crowded, and Direct-relax.

The initial trace was taken in a challenging environment that was crowded with
competing devices, and is thus referred as the _Crowd_ scenario. This was a
trace of latency measurements of sending small (128 bytes) UDP packets every
100 ms between a server and a node that connects to an LTE BS, or Wi-Fi AP, at
the Aalborg University campus [17]. These measurements were taken in a fully
NLOS (indoor) environment with physically-static devices that had a varying
channel competition during the day. Nonetheless, the measurements were taken
in good-SNR radio conditions, which enhances the differences between LTE and
Wi-Fi protocol operation in the presence of competing traffic. Most
importantly, the medium access protocol (MAC) varies on whether it needs to
operate in some licensed spectrum (for LTE), or it needs to operate in a
highly contentious unlicensed spectrum (for Wi-Fi). In Fig. 4 we can see that
Wi-Fi can provide 5 ms delays for 90% of its packets. However, if an
application allows for more relaxed timeliness constraint, 99% of packets
arrive before the 80 ms deadline. On the other hand, it is very hard to
achieve very low latency for LTE, even for a very small portion of packets.
Therefore there is little difference in the LTE latency deadlines, 36 ms and
40 ms, for 90% and 99% of packets, respectively.

###  IV-B Relaxed Scenario

The second trace was taken in a more controlled scenario where the contention
for the wireless resources was relaxed and is thus referred to as the _Relax_
scenario. The measurements are taken outdoors, with full LOS links, in the 5G
Playground Carinthia testbed in the outskirts of Klagenfurt, Austria. This was
a full UAVAP implementation where the samples were taken while flying along a
pre-planned path with a constant speed of 1 m/s. Each 100 ms a small packet
(20 byte payload) was generated and transmitted from a ground user, sent
through the UAVAP to the LTE BS. The testbed was set in a way that is easy to
extract and separate the LTE and Wi-Fi behavior[15].

From the statistical performance of Wi-Fi and LTE In Fig. 4 we can notice that
both interfaces retain the same statistical behaviour as in the Crowd scenario
but with different timing performances. In general, the measured latencies
were slightly lower due to the smaller physical distance of the server from
the LTE BS. The Wi-Fi can provide 5 ms timeliness to 95% of the packets, while
the 99% packet pass-through rate is reached when the timeliness constraint is
set to be 200 ms. The LTE connection here has a significantly wider latency
gap for 90% and 99% of packets, 18 ms and 23 ms, respectively.

###  IV-C Direct-relaxed Scenario

Within the same _Relax_ testing scenario in Klagenfurt, a separate measurement
trace was taken only for the GN that directly connects to the BS. This was to
provide realistic and proportional measurement set to compare against the
UAVAP setup. In this case, the GN had mixed LOS and NLOS LTE link with the BS
along the same outdoor path and speed taken in the previous trace. This is
common outdoor environment where reliable implementation of cyber-physical
applications is needed. Finally, we note that this measurement is strongly
correlated with the LTE${}_{\text{relax}}$. This correlation can be exploited
to generate a synthetic GE model for a LTE${}_{\text{direct-crowded}}$

Table I: $p$ and $r$ transition probabilities used in evaluation. | Crowd |  | Relax
---|---|---|---
| $p$ | $r$ |  | $p$ | $r$
Wi-Fi | 0.0515 | 0.9468 |  | 0.0001 | 0.0845
LTE | 0.0178 | 0.2577 |  | 0.0127 | 0.8356

##  V Numerical Analysis & Results

Combining latency measurements is a difficult feat due to the vast variability
of each testing setup, unpredictable environment, unforeseen biases, and
traffic proportions that result in incompatible latency distributions.
However, the GE model can abstract the complexity of the aforementioned
challenges by setting an adequate $\theta$ that well represents the
proportions of each respective testing setup. With this, we want to achieve a
thorough analysis of the burst error performance in any combination of
crowding on each links.

Since both measurement campaigns concern vastly different properties and BSAL
networks, the timing constraints had to be scaled, so that the results can be
comparable. We have done this by considering the equi-performance points for
both Wi-Fi and LTE. In detail, these are the delay values for which both
interfaces, LTE and Wi-Fi provide the same reliability within their respective
measurement setup, or simply, the point where the CDF lines cross. For the
crowded scenario the delay separating a bad versus a good state of the channel
is $\theta_{\text{crowd}}=38.25\,\text{ms}$. For the relaxed scenario we set
the cutoff delay per hop to $\theta_{\text{relax}}=22.5\,\text{ms}$. Since the
interfaces produce similar behaviour in both measurement campaigns, the
crossing point is the only relevant point that is correlated between both
campaigns. Setting the delay values like this gives a balanced weight to both
interfaces given the two different testing environments, and allows us to
analyze bursty losses across multiple packets. This isolates the GE models for
each interface and measurement campaign separately, given the different
distances to the server. In the case of the two different measurement
campaigns, the baseline latency difference is 0.7 times larger for the _Crowd_
setup. Therefore, the resulting $p$ and $r$ values are presented in Tab. I.

To compare the performance of the UAVAP we use the LTE-direct-relaxed
connection for a GN outdoors and on the ground. Since the Relax and the LTE-
direct measurements are taken in the exact same testing conditions, we used
the following reasoning for the latency deadline $\theta_{\text{direct-
relaxed}}=2\cdot\theta_{\text{relax}}=45\leavevmode\nobreak\ \text{ms}$. This
is a best case scenario for the $\theta_{\text{direct-relaxed}}$ as it is
equivalent to the upper performance bound where not one, but both interfaces
are in a bad state in the UAVAP setup. The resulting GE model values for this
measurement scenario are $p_{\text{direct-relaxed}}=0.1457$ and
$r_{\text{direct-relaxed}}=0.7857$. Following this, we used the
LTE${}_{\text{crowd}}$ measurements to synthesize the $p$ and $r$ values for a
potential LTE${}_{\text{Direct-crowd}}$ by following the proportional
correlation between LTE${}_{\text{relax}}$ and LTE${}_{\text{Direct-relax}}$.
This resulted in the $p_{\text{direct-crowd}}=0.2042$ and $r_{\text{direct-
crowd}}=0.2423$ values for the fictional LTE${}_{\text{Direct-crowd}}$
dataset.

Figure 5: Burst error length analysis for all combinations of possible
environments _Crowd_ and _Relax_ for each interface.

An analysis between both sets of measurements are that the LTE $p$
probabilities are almost equal, with $p_{\text{LTE relax}}$ having a slightly
smaller probability to move over to the bad state. Moreover, there is a larger
discrepancy between the LTE $r$ values in both environments, where
$r_{\text{LTE relax}}$ has a much larger probability to return to the good
state once it has performed poorly. This makes a big difference in the steady
state probability of being in the bad channel state of 0.0150 for relaxed,
versus 0.0646 for the crowded. There are clear differences between the
performance of the LTE and Wi-Fi interfaces, however, as per the previous
section the issue of using UAVAP is that the negative effects of both
interfaces are stacked when relaying.

There are more distinct differences between both sets of Wi-Fi measurements.
While the probability for the Wi-Fi to go to the bad state are generally small
$p_{\text{Wi-Fi relax}}$ is orders of magnitude lower than in the crowded
scenario. However, the probability of the Wi-Fi connection to recover in the
crowded scenario is much higher than in the relaxed one. This is the case
because outages in the relaxed scenario are due to rare but persistent channel
outages. Therefore, the outage probability of the relaxed Wi-Fi has a value of
0.0015, versus the outage probability of crowded Wi-Fi of 0.0516.

In Fig. 5 we show results of a Monte-Carlo simulation, of $10^{9}$ samples,
for all possible combinations of the measurement campaigns. It is common that
cyber-physical control applications for the GN to allow for error bursts of 3
consecutive packet errors [17, 3] (equivalent to 300 ms). With such an
implementation in mind we can notice that all UAVAP implementations outperform
a direct-relaxed LTE link. Moreover, the usefulness of the UAVAP is not very
strong when the LTE connection is in the crowded mode. However, this would
mean that the LTE link for the direct implementation.

In Fig. 6 we expand the investigation to longer burst length tolerances, up to
30 samples (equivalent to 3 s). We can see that the low probability of Wi-Fi
to exit a bad state has a bad influence on the system. This is particularly
visible in the LTE${}_{\text{relax}}$-Wi-Fi${}_{\text{Relax}}$ scenario that
rarely went into a bad state, but remains there for a very long time. This
”flattens” the burst length performance when poor SNR conditions, or other
systematic failures that affect Wi-Fi occur.

Figure 6: Burst error length analysis for all combinations of possible
environments _Crowd_ and _Relax_ for each interface.

##  VI Conclusions

In this work we analyzed the burst packet error probability of UAVAPs for
providing on-demand communications services for the purpose of cyber-physical
applications. We did this by using the GE 2-state model for simulating
temporal correlation in wireless channels from two separate sets of latency
measurements of LTE and Wi-Fi. The two different sets of latency measurements
were subjected to different latency deadline values as a way to demonstrate
the timeliness sensitivity of packets. The two sets, one taken in a busy
spectrum and another in a relaxed spectrum are used to provide a more
generalized picture from the initial measurements. The results show that using
a UAVAP in order to avoid burst packet errors was always superior to
transmitting to a distant LTE BS directly. The advantages can vary depending
on the crowding of the radio resources.

We need to note that even though this analysis contains varied measurement
data, further investigation is needed to verify the usefulness of UAVAPs for
ubiquitous deployment. In particular, the LTE connection to the main grid can
drastically vary on the deployed environment. Finally, this work inspires
future studies where the GN exploits both interfaces, such as in the case of
multi-connectivity, to exploit the advantages of both the direct and the UAVAP
links.

## Acknowledgment

The work was supported by the European Union’s research and innovation
programme under the Marie Sklodowska-Curie grant agreement No. 812991
”PAINLESS” within the Horizon 2020 Program.

## References

  * [1] W. Saad, M. Bennis, and M. Chen, “A vision of 6g wireless systems: Applications, trends, technologies, and open research problems,” _IEEE Network_ , vol. 34, no. 3, pp. 134–142, Oct. 2020.
  * [2] P. Popovski, Č. Stefanović, J. J. Nielsen, E. De Carvalho, M. Angjelichinoski, K. F. Trillingsgaard, and A.-S. Bana, “Wireless access in ultra-reliable low-latency communication (URLLC),” _IEEE Transactions on Communications_ , vol. 67, no. 8, pp. 5783–5801, 2019.
  * [3] 3GPP, “Technical specification TS 22.104 V16.5.0. 5G; Service requirements for cyber-physical control applications in vertical domains,” 2021.
  * [4] N. Babu, I. Donevski, A. Valcarce, P. Popovski, J. J. Nielsen, and C. B. Papadias, “Fairness based energy-efficient 3d path planning of a portable access point: A deep reinforcement learning approach,” _IEEE Open Journal of the Communications Society_ , pp. 1–1, 2022.
  * [5] M. Virgili, N. Babu, M. Javidsharifi, I. Valiulahi, C. Masouros, A. J. Forsyth, T. Kerekes, and C. B. Papadias, “Cost-efficient design of an energy-neutral uav-based mobile network,” _IEEE Transactions on Communications_ , pp. 1–1, 2022.
  * [6] A. Fotouhi, H. Qiang, M. Ding, M. Hassan, L. G. Giordano, A. Garcia-Rodriguez, and J. Yuan, “Survey on UAV Cellular Communications: Practical Aspects, Standardization Advancements, Regulation, and Security Challenges,” _IEEE Communications Surveys Tutorials_ , vol. 21, no. 4, pp. 3417–3442, Mar. , 2019.
  * [7] M. Mozaffari, W. Saad, M. Bennis, Y. Nam, and M. Debbah, “A Tutorial on UAVs for Wireless Networks: Applications, Challenges, and Open Problems,” _IEEE Communications Surveys Tutorials_ , vol. 21, no. 3, pp. 2334–2360, Mar. , 2019.
  * [8] I. Donevski and J. J. Nielsen, “Dynamic Standalone Drone-Mounted Small Cells,” in _Proc. of European Conference on Networks and Communications (EuCNC)_ , Dubrovnik, Croatia, Sep. , 2020, pp. 342–347.
  * [9] R. Gangula, O. Esrafilian, D. Gesbert, C. Roux, F. Kaltenberger, and R. Knopp, “Flying rebots: First results on an autonomous uav-based lte relay using open airinterface,” in _2018 IEEE 19th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_. IEEE, 2018, pp. 1–5.
  * [10] I. Donevski, J. J. Nielsen, and P. Popovski, “Standalone Deployment of a Dynamic Drone Cell for Wireless Connectivity of Two Services,” in _Proc. of IEEE Wireless Communications and Networking Conference (WCNC)_ , no. TBP, Nanjing, China , Apr. , 2021.
  * [11] X. Lin, R. Wiren, S. Euler, A. Sadam, H.-L. Määttänen, S. Muruganathan, S. Gao, Y.-P. E. Wang, J. Kauppi, Z. Zou _et al._ , “Mobile network-connected drones: Field trials, simulations, and design insights,” _IEEE Vehicular Technology Magazine_ , vol. 14, no. 3, pp. 115–125, 2019.
  * [12] A. Fakhreddine, C. Bettstetter, S. Hayat, R. Muzaffar, and D. Emini, “Handover challenges for cellular-connected drones,” in _Proceedings of the 5th Workshop on Micro Aerial Vehicle Networks, Systems, and Applications_ , 2019, pp. 9–14.
  * [13] H. C. Nguyen, R. Amorim, J. Wigard, I. Z. Kovács, T. B. Sørensen, and P. E. Mogensen, “How to ensure reliable connectivity for aerial vehicles over cellular networks,” _IEEE Access_ , vol. 6, pp. 12 304–12 317, 2018\.
  * [14] B. Van der Bergh, A. Chiumento, and S. Pollin, “Lte in the sky: Trading off propagation benefits with interference costs for aerial nodes,” _IEEE Communications Magazine_ , vol. 54, no. 5, pp. 44–50, 2016.
  * [15] I. Donevski, C. Raffelsberger, M. Sende, A. Fakhreddine, and J. J. Nielsen, “An experimental analysis on drone-mounted access points for improved latency-reliability,” in _Proceedings of the 7th ACM Workshop on Micro Aerial Vehicle Networks, Systems, and Applications_ , 2021.
  * [16] P. M. de Sant Ana, N. Marchenko, P. Popovski, and B. Soret, “Wireless control of autonomous guided vehicle using reinforcement learning,” in _Proc. IEEE Global Communications Conference (GLOBECOM)_ , dec 2020.
  * [17] I. Donevski, I. Leyva-Mayorga, J. J. Nielsen, and P. Popovski, “Performance trade-offs in cyber–physical control applications with multi-connectivity,” _Frontiers in Communications and Networks_ , vol. 2, 2021.
  * [18] J. J. Nielsen, R. Liu, and P. Popovski, “Ultra-reliable low latency communication using interface diversity,” _IEEE Transactions on Communications_ , vol. 66, no. 3, pp. 1322–1334, 2018.


# Convergence analysis of a weak Galerkin finite element method on a Shishkin
mesh for a singularly perturbed fourth-order problem in 2D

Shicheng Liu<EMAIL_ADDRESS>Xiangyun Meng<EMAIL_ADDRESS>Qilong
Zhai<EMAIL_ADDRESS>School of Mathematics, Jilin University, Changchun
130012, China  School of Mathematics and Statistics, Beijing Jiaotong
University, Beijing 100044, China

###### Abstract

In this paper, we apply the weak Galerkin (WG) finite element method to solve
the singularly perturbed fourth-order boundary value problem in 2D domain. A
Shishkin mesh is used to ensure that the method exhibits uniform convergence,
regardless of the singular perturbation parameter. Asymptotically optimal
order error estimate in a $H^{2}$ discrete norm is established for the
corresponding WG solutions. Numerical tests are presented to verify the theory
of convergence.

###### keywords:

Weak Galerkin finite element method, fourth-order differential equation,
singularly perturbed, Shishkin mesh.

###### MSC:

[2020] 65N15 , 65N30 , 35B25

††journal: Journal of Computational and Applied Mathematics

##  1 Introduction

This paper is concerned with a singularly perturbed fourth-order boundary
value problem in a square region $\Omega$. The problem is solved by the weak
Galerkin (WG) method with a Shishkin mesh, as follows:

$\displaystyle\varepsilon^{2}\Delta^{2}u-\Delta u$
$\displaystyle=f,\quad\text{in}~{}\Omega,$ (1.1) $\displaystyle u$
$\displaystyle=0,\quad\text{on}~{}\partial\Omega,$ (1.2)
$\displaystyle\frac{\partial u}{\partial\mathbf{n}}$
$\displaystyle=0,\quad\text{on}~{}\partial\Omega,$ (1.3)

with a positive parameter $\varepsilon$ satisfying $0<\varepsilon\ll 1$ and
$f\in L^{2}(\Omega)$. This problem is used for thin elastic plates clamped in
tension. $f$ represents the transverse load, $\varepsilon$ symbolizes the
ratio of bending stiffness to tensile stiffness of the plate, and the
function’s solution, denoted as $u$, represents the displacement of the plate.
This problem emerges in the investigation of the linearization of the fourth-
order perturbation associated with the fully nonlinear Monge-Ampère
equation[2, 3].

A variational formulation for the fourth-order equation (1.1) with the
boundary conditions (1.2) and (1.3) seeks $u\in H^{2}(\Omega)$ such that

$\varepsilon^{2}(\Delta u,\Delta v)+(\nabla u,\nabla v)=(f,v),\quad\forall
v\in H^{2}_{0}(\Omega),$ (1.4)

where $H^{2}_{0}(\Omega)$ is a subspace of the Sobolev space $H^{2}(\Omega)$
consisting of functions with vanishing value and normal derivative on
$\partial\Omega$.

Research into numerical methods for singularly perturbed differential
equations commenced in the early 1970s, with the frontier of research
continuously expanding ever since. Bakhvalov made an important early
contribution to the optimization of numerical methods by means of special
meshes [1] in 1969. In the early 1990s, G.I. Shishkin proposed piecewise-
equidistant meshes [20, 26]. These meshes, characterized by their very simple
structure, are usually easy to analyze. Shishkin meshes for various problems
and numerical methods have been studied since and they are still popular.
Numerous scholars have explored the numerical analysis of analogous problems
using various finite element methods, such as the mixed finite element method
in [10], the hp finite element method in [7], the continuous interior penalty
finite element method in [11], the upwind finite volume element method in
[36], the conforming finite element method in [24]. And also some papers
consider finite element methods on quasi-uniform meshes: a continuous interior
penalty finite element method in [3] and nonconforming finite element method
in [5, 13, 15, 19, 23, 31, 34, 38].

In this paper, we employ the WG method with a Shishkin mesh to investigate the
convergence behavior of a fourth-order boundary value problem that exhibits
singular perturbation. The WG method proves to be an effective numerical
technique for the partial differential equations(PDEs). The initial proposal
for its application in solving second-order elliptic problems was made by
Junping Wang and Xiu Ye in [28]. The core concept involves establishing
distinct basis functions for the interior and boundary of each partitioned
element, and substituting the traditional differential operator with a
discretized weak differential operator. The WG method has been applied to
Stokes equations [29, 32, 37], elasticity equations [4, 14, 27, 39], Maxwell’s
equations [22], biharmonic equations [8, 44], Navier-Stokes equations [16, 18,
43], Brinkman equations [21, 35, 40], the multigrid approach [6], the
incompressible flow [45], the maximum principle [17, 30], the post-processing
technique [33] and so on. For singular perturbed value problems, the WG method
has also yielded some results, such as the singularly perturbed convection-
diffusion problems for WG in 1D [42, 46] and 2D [41], the singularly perturbed
biharmonic equation for WG in uniform mesh [9].

This paper is organized as follows. In Section 2, we introduce the Shishkin
mesh and the assumptions associated. In Section 3, we give the definitions of
the weak Laplacian operator and weak gradient operator. We also present WG
finite element schemes for the singularly perturbed value problem. In Section
4, we introduce some local $L^{2}$ projection operators and give some
approximation properties. In Section 5, we establish error estimates for the
WG scheme in a $H^{2}$-equivalent discrete norm. And in Section 6, we report
the results of two numerical experiments.

##  2 Preliminaries and notations

To solve problem (1.1)-(1.3), we suppose the following assumption holds in
[19], which involves structuring the solution and decomposing $u$ into smooth
and layered components.

######  Assumption 2.1.

The solution $u$ to the singularly perturbed fourth-order boundary value
problem (1.1)-(1.3) can be expressed as the sum of its smooth and layered
components, as follows:

$\displaystyle u=S+\sum_{l=1}^{4}E_{l}+E_{12}+E_{23}+E_{34}+E_{41}.$

In this decomposition, $S$ represents a smooth function, each $E_{l}$
corresponds to a boundary layer component along the sides of $\bar{\Omega}$ in
anti-clockwise order, and the remaining components are corner layer parts
along the corners of $\bar{\Omega}$ in anti-clockwise order. Furthermore,
there exists a constant $C$ for all points $(x,y)\in\bar{\Omega}$, which is
independent of $x$ and $y$. This constant satisfies the following conditions
for $0\leq i+j\leq k+1$,

$\displaystyle\left|\frac{\partial^{i+j}S}{\partial x^{i}\partial
y^{j}}\right|\leq C,$ $\displaystyle\left|\frac{\partial^{i+j}E_{1}}{\partial
x^{i}\partial y^{j}}\right|\leq C\varepsilon^{1-j}e^{-\frac{y}{\varepsilon}},$
$\displaystyle\left|\frac{\partial^{i+j}E_{4}}{\partial x^{i}\partial
y^{j}}\right|\leq C\varepsilon^{1-i}e^{-\frac{x}{\varepsilon}},$
$\displaystyle\left|\frac{\partial^{i+j}E_{41}}{\partial x^{i}\partial
y^{j}}\right|\leq
C\varepsilon^{1-i-j}e^{-\frac{x}{\varepsilon}}e^{-\frac{y}{\varepsilon}}.$

Moreover, the other components of the decomposition are bounded in a similar
manner.

In order to solve the layer structure in the solution of problems (1.1)-(1.3),
a well-suited layer-adapted Shishkin mesh be considered. This mesh is refined
in the layers. For a comprehensive discussion on the construction of Shishkin
meshes, please refer to [25].

Consider a positive integer $N\geq 4$ that is divisible by $4$. We introduce a
mesh transition parameter $\lambda$ to determine the location at which the
mesh switches from coarse to fine. This parameter is defined by

$\lambda=min\left\\{\alpha\varepsilon lnN,\frac{1}{4}\right\\},$ (2.1)

where $\alpha$ is a positive constant, selected to be $k+1$ for the subsequent
analysis.

Create a piecewise equidistant mesh for the interval $[0,1]$ by dividing it as
follows: divide $[0,\lambda]$ into $N/4$ subintervals, $[\lambda,1-\lambda]$
into $N/2$ subintervals, and $[1-\lambda,1]$ into $N/4$ subintervals. The
Shishkin mesh for the problem (1.1)-(1.3) is formed by taking the tensor
product of two such one-dimensional meshes, as illustrated in Figure 1. The
fine meshwidth denoted as $h$ and the coarse meshwidth denoted as $H$ in the
Shishkin mesh represented by $\mathcal{T}_{N}$ exhibit the following
characteristics:

$h=\frac{4\lambda}{N}\leq C\varepsilon N^{-1}lnN\quad\text{and}\quad
H=2\frac{1-2\lambda}{N}\leq CN^{-1}$ (2.2)

for some constant $C$. The domain $\Omega$ is divided into some subdomains,
see Figure 1. Denote $\Omega_{l1}\cup\Omega_{s}\cup\Omega_{l3}$ by
$\Omega_{r}^{1}$ and $\Omega_{l2}\cup\Omega_{s}\cup\Omega_{l4}$ by
$\Omega_{r}^{2}$.

Figure 1: A rectangular Shishkin mesh with $N=8$ and dissection of $\Omega$.

##  3 Weak differential operator and WG scheme

To propose the weak Galer-kin method, we introduce some key concepts. Consider
a element $T$ belonging to the partition $\mathcal{T}_{N}$ with a boundary
denoted by $\partial T$. Let the set of all edges in $\mathcal{T}_{N}$ be
represented as $\mathcal{E}_{N}$. We define a weak function
$v=\\{v_{0},v_{b},\mathbf{v}_{g}\\}$ on the element $T$, where $v_{0}\in
L^{2}(T)$, $v_{b}\in H^{\frac{1}{2}}(\partial T)$, and
$\mathbf{v}_{g}\cdot\mathbf{n}\in H^{\frac{1}{2}}(\partial T)$, with
$\mathbf{n}$ representing the outward normal direction on $\partial T$.
Furthermore, the first and second components, $v_{0}$ and $v_{b}$, correspond
to values of $v$ in the interior and on the boundary of $T$. The third
component $\mathbf{v}_{g}$ is employed to approximate the gradient of $v$
along the boundary of $T$. It’s important to note that each edge
$e\in\mathcal{E}_{N}$ has a unique value for $v_{b}$ and $\mathbf{v}_{g}$.
Additionally, it’s worth mentioning that $v_{b}$ and $\mathbf{v}_{g}$ may not
necessarily be associated with the trace of $v_{0}$ and $\nabla v_{0}$ on
$\partial T$.

For any integer $k\geq 3$, we establish a local discrete weak function space
for any element $T\in\mathcal{T}_{N}$ denoted as $W_{k}(T)$:

$W_{k}(T)=\left\\{v=\\{v_{0},v_{b},\mathbf{v}_{g}\\}:v_{0}\in\mathbb{Q}_{k}(T),v_{b}\in\mathbb{P}_{k}(e),\mathbf{v}_{g}\in[\mathbb{P}_{k}(e)]^{2},e\in\partial
T\right\\},$

where $e$ is the edge of $\partial T$, and $\mathbb{Q}_{k}$ is the space of
polynomials which are of degree not exceeding $k$ with respect to each one of
the variables $x$ and $y$. By extending $W_{k}(T)$ to encompass all element
$T\in\mathcal{T}_{N}$, we introduce the definition of a weak Galerkin space:

$V_{N}=\left\\{v=\\{v_{0},v_{b},\mathbf{v}_{g}\\}:v|_{T}\in W_{k}(T),\forall
T\in\mathcal{T}_{N}\right\\}.$

Let $V_{N}^{0}$ represent the subspace of $V_{h}$ where the traces vanish:

$V_{N}^{0}=\left\\{v=\\{v_{0},v_{b},\mathbf{v}_{g}\\}\in
V_{h},v_{b}|_{e}=0,\mathbf{v}_{g}\cdot\mathbf{n}|_{e}=0,e\subset\partial
T\cap\partial\Omega\right\\}.$

For any $v=\\{v_{0},v_{b},\mathbf{v}_{g}\\}$ and a fixed integer $k\geq 3$, we
define a discrete weak Laplacian operator $\Delta_{w,k}$ as a unique
polynomial $\Delta_{w,k}v\in\mathbb{Q}_{k}(T)$ on $T$ satisfying the following
equation:

$\displaystyle\left(\Delta_{w,k}v,\varphi\right)_{T}=\left(v_{0},\Delta\varphi\right)_{T}-\left\langle
v_{b},\nabla\varphi\cdot\mathbf{n}\right\rangle+\left\langle\mathbf{v}_{g}\cdot\mathbf{n},\varphi\right\rangle,\quad\forall\varphi\in\mathbb{Q}_{k}(T),$
(3.1)

where $\mathbf{n}$ is the unit outward normal vector to $\partial T$.
Likewise, a discrete weak gradient $\nabla_{w,k}$ is defined on $T$ as a
unique polynomial $\nabla_{w,k}v\in[\mathbb{Q}_{k}(T)]^{2}$ satisfying:

$\displaystyle\left(\nabla_{w,k}v,\mathbf{q}\right)_{T}=-\left(v_{0},\nabla\cdot\mathbf{q}\right)_{T}+\left\langle
v_{b},\mathbf{q}\cdot\mathbf{n}\right\rangle_{\partial
T},\quad\forall\mathbf{q}\in\left[\mathbb{Q}_{k}(T)\right]^{2}.$ (3.2)

For simplicity, when there is no confusion, we drop the subscript $k$ in the
notations $\Delta_{w,k}$ and $\nabla_{w,k}$ for the discrete weak Laplacian
and the discrete weak gradient. Additionally, we introduce the following
notations:

$\displaystyle\left(\Delta_{w}v,\Delta_{w}w\right)_{\mathcal{T}_{N}}:$
$\displaystyle=\sum_{T\in\mathcal{T}_{N}}\left(\Delta_{w}v,\Delta_{w}w\right)_{T},$
$\displaystyle\left(\nabla_{w}v,\nabla_{w}w\right)_{\mathcal{T}_{N}}:$
$\displaystyle=\sum_{T\in\mathcal{T}_{N}}\left(\nabla_{w}v,\nabla_{w}w\right)_{T}.$

Let us introduce a stabilizer, which is a bilinear form for any
$u_{N}=\\{u_{0},u_{b},\mathbf{u}_{g}\\}$ and
$v=\\{v_{0},v_{b},\mathbf{v}_{g}\\}$ in the space $V_{h}$. It is defined as
follows:

$\displaystyle\begin{aligned}
s\left(u_{h},v\right)=\sum_{T\in\mathcal{T}_{N}}\Big{(}\varepsilon^{2}h^{-1}&\left\langle\nabla
u_{0}-\mathbf{u}_{g},\nabla v_{0}-\mathbf{v}_{g}\right\rangle_{\partial T}\\\
&+\left(\varepsilon^{2}h^{-2}H^{-1}+H^{-1}\right)\left\langle
u_{0}-u_{b},v_{0}-v_{b}\right\rangle_{\partial T}\Big{)},\end{aligned}$ (3.3)

where $h$ and $H$ is defined as (2.2).

Algorithm 1 WG Algorithm

To obtain a numerical approximation for (1.1)-(1.3), we seek
$u_{N}=\\{u_{0},u_{b},\mathbf{u}_{g}\\}\in V_{N}$ that satisfies the following
equation:

$\displaystyle\varepsilon^{2}(\Delta_{w}u_{N},\Delta_{w}v)_{N}+(\nabla_{w}u_{N},\nabla_{w}v)_{N}+s(u_{N},v)=(f,v),$
(3.4)

for any $v=\\{v_{0},v_{b},\mathbf{v}_{g}\\}\in V_{N}^{0}$.

The following lemma provides a valuable result regarding the finite element
space $V^{0}_{N}$.

######  Lemma 3.1.

For any $v\in V^{0}_{N}$, let ${|||}v{|||}$ be given as follows:

$\displaystyle{|||}v{|||}^{2}=\varepsilon^{2}(\Delta_{w}v,\Delta_{w}v)_{\mathcal{T}_{N}}+(\nabla_{w}v,\nabla_{w}v)_{\mathcal{T}_{N}}+s(v,v),$
(3.5)

Then, ${|||}\cdot{|||}$ defines a norm in $V^{0}_{N}$.

###### Proof.

We shall only confirm the positivity property for ${|||}\cdot{|||}$. Consider
$v=\\{v_{0},v_{b},\mathbf{v}_{g}\\}\in V^{0}_{N}$ with the assumption that
${|||}v{|||}=0$. This implies $\Delta_{w}v=0$ and $\nabla_{w}v=0$ in each
element $T$, while also satisfying $v_{0}=v_{b}$ and $\nabla
v_{0}=\mathbf{v}_{g}$ on $\partial T$. Next, we will proof that $\Delta
v_{0}=0$ in each element $T$. To this end, for any
$\varphi\in\mathbb{Q}_{k}(T)$, employing the definition (3.1) and the fact
that $\Delta_{w}v=0$, we obtain

$\begin{split}0&=\left(\Delta_{w}v,\varphi\right)_{T}\\\
&=\left(v_{0},\Delta\varphi\right)_{T}-\left\langle
v_{b},\nabla\varphi\cdot\mathbf{n}\right\rangle_{\partial
T}+\left\langle\mathbf{v}_{g}\cdot\mathbf{n},\varphi\right\rangle_{\partial
T}\\\ &=\left(\Delta v_{0},\varphi\right)_{T}+\left\langle
v_{0}-v_{b},\nabla\varphi\cdot\mathbf{n}\right\rangle_{\partial
T}+\left\langle\mathbf{v}_{g}\cdot\mathbf{n}-\nabla
v_{0}\cdot\mathbf{n},\varphi\right\rangle_{\partial T}\\\ &=\left(\Delta
v_{0},\varphi\right)_{T},\end{split}$ (3.6)

where we have applied the fact that $v_{0}-v_{b}=0$ and $\nabla
v_{0}-\mathbf{v}_{g}=0$ in the final equality. The identity (3.6) indicates
that $\Delta v_{0}=0$ in each element $T$. Together with the conditions
$v_{0}=v_{b}$ and $\nabla v_{0}=\mathbf{v}_{g}$ on $\partial T$, we conclude
that $v$ is a globally smooth harmonic function on $\Omega$. Considering the
boundary conditions $v_{b}=0$ and $\mathbf{v}_{g}\cdot\mathbf{n}=0$, we infer
that the unique solution is $v\equiv 0$ on $\Omega$. This concludes the proof.
∎

##  4 Local $L^{2}$ projection operators and approximation properties

In this section, we introduce some projection operators for each element
$T\in\mathcal{T}_{N}$. Consider the $L^{2}$ projection operator
$\mathcal{Q}_{0}$, which projects onto $\mathbb{Q}_{k}(T)$. Additionally, for
each edge $e\in\partial T$, we consider the $L^{2}$ projection operators
$\mathcal{Q}_{b}$ and $\mathbf{Q}_{g}$ onto local polynomial spaces
$\mathbb{P}_{k}(e)$ and $[\mathbb{P}_{k}(e)]^{2}$, respectively. We define a
projection $\mathcal{Q}_{N}$ of $u$ into the finite element space $V_{N}$ such
that on each element $T$

$\displaystyle\mathcal{Q}_{N}u=\left\\{\mathcal{Q}_{0}u,\mathcal{Q}_{b}u,\mathbf{Q}_{g}(\nabla
u)\right\\}.$

Furthermore, let $\mathbf{Q}_{N}$ represent the local $L^{2}$ projection onto
$[\mathbb{Q}_{k}(T)]^{2}$. The following lemma demonstrates that the weak
Laplacian $\Delta_{w}$ is the polynomial projection of the classical Laplacian
$\Delta$.

######  Lemma 4.1.

On each element $T\in\mathcal{T}_{N}$, for any $v\in H^{2}(T)$,

$\displaystyle\Delta_{w}(\mathcal{Q}_{N}v)=\mathcal{Q}_{0}(\Delta v).$ (4.1)

###### Proof.

For any $\tau\in\mathbb{Q}_{k}(T)$, we obtain that

$\displaystyle\left(\Delta_{w}\mathcal{Q}_{N}u,\tau\right)_{T}$
$\displaystyle=\left(\mathcal{Q}_{0}u,\Delta\tau\right)_{T}-\left\langle\mathcal{Q}_{b}u,\nabla\tau\cdot\mathbf{n}\right\rangle_{\partial
T}+\left\langle\mathbf{Q}_{g}(\nabla
u)\cdot\mathbf{n},\tau\right\rangle_{\partial T}$
$\displaystyle=\left(u,\Delta\tau\right)_{T}-\left\langle
u,\nabla\tau\cdot\mathbf{n}\right\rangle_{\partial T}+\left\langle\nabla
u\cdot\mathbf{n},\tau\right\rangle_{\partial T}$ $\displaystyle=\left(\Delta
u,\tau\right)_{T}$ $\displaystyle=\left(\mathcal{Q}_{0}\Delta
u,\tau\right)_{T}.$

This concludes the proof. ∎

A similar lemma holds for the weak gradient $\nabla_{w}$, as indicated in the
subsequent lemma.

######  Lemma 4.2.

On each element $T\in\mathcal{T}_{N}$, for any $v\in H^{2}(T)$,

$\displaystyle\nabla_{w}(\mathcal{Q}_{N}u)=\mathbf{Q}_{N}(\nabla u).$ (4.2)

###### Proof.

For any $\mathbf{q}\in[\mathbb{Q}_{k}(T)]^{2}$, from the weak gradient
definition and integration by parts, we get

$\displaystyle\left(\nabla_{w}\mathcal{Q}_{N}u,\mathbf{q}\right)_{T}$
$\displaystyle=-\left(\mathcal{Q}_{0}u,\nabla\cdot\mathbf{q}\right)_{T}+\left\langle\mathcal{Q}_{b}u,\mathbf{q}\cdot\mathbf{n}\right\rangle_{\partial
T}$ $\displaystyle=-\left(u,\nabla\cdot\mathbf{q}\right)_{T}+\left\langle
u,\mathbf{q}\cdot\mathbf{n}\right\rangle_{\partial T}$
$\displaystyle=\left(\nabla u,\mathbf{q}\right)_{T}$
$\displaystyle=\left(\mathbf{Q}_{N}\nabla u,\mathbf{q}\right)_{T}.$

This corresponds to identity (4.2). ∎

Now, we introduce notation that will be employed in the following lemmas.
Consider any element $T$ in the partition $\mathcal{T}_{N}$. We define
$\partial T_{1}$ and $\partial T_{2}$ as the sets of element edges that are
parallel to the $x$ and $y$ axes, respectively. The following lemmas are
employed in the convergence analysis, and readers are directed to [12, chapter
3.1] for a detailed proof process.

######  Lemma 4.3.

Consider $\phi\in H^{k+1}(T)$ with $k\geq 3$. Let $\mathcal{Q}_{0}\phi$ denote
the $L^{2}$-projection of $\phi$ onto $\mathbb{Q}_{k}(T)$. Then the following
inequality estimate holds,

$\displaystyle\|\phi-\mathcal{Q}_{0}\phi\|_{\partial T_{i}}\leq
C\left(h_{j}^{k+\frac{1}{2}}\|\partial_{j}^{k+1}\phi\|_{T}+h_{i}^{k+1}h_{j}^{-\frac{1}{2}}\|\partial_{i}^{k+1}\phi\|_{T}+h_{j}^{\frac{1}{2}}h_{i}^{k}\|\partial_{i}^{k}\partial_{j}\phi\|_{T}\right),$
(4.3) $\displaystyle\|\partial_{i}(\phi-\mathcal{Q}_{0}\phi)\|_{\partial
T_{i}}\leq
C\left(h_{i}^{k}h_{j}^{-\frac{1}{2}}\|\partial_{i}^{k+1}\phi\|_{T}+h_{i}^{k-1}h_{j}^{\frac{1}{2}}\|\partial_{i}^{k}\partial_{j}\phi\|_{T}+h_{j}^{k-\frac{1}{2}}\|\partial_{j}^{k}\partial_{i}\phi\|_{T}\right),$
(4.4) $\displaystyle\|\partial_{j}(\phi-\mathcal{Q}_{0}\phi)\|_{\partial
T_{i}}\leq
C\left(h_{j}^{k-\frac{1}{2}}\|\partial_{j}^{k+1}\phi\|_{T}+h_{i}^{k}h_{j}^{-\frac{1}{2}}\|\partial_{i}^{k}\partial_{j}\phi\|_{T}\right),$
(4.5)

where $i,j\in\\{1,2\\}$, and $i\neq j$.

######  Lemma 4.4.

Let $v\in\mathbb{Q}_{k}(T)$ with $k\geq 3$ such that the following
inequalities holds,

$\displaystyle\|v\|_{\partial T_{i}}\leq$ $\displaystyle
Ch_{j}^{-\frac{1}{2}}\|v\|_{T},$ (4.6)
$\displaystyle\|\partial_{i}v\|_{T}\leq$ $\displaystyle Ch_{i}^{-1}\|v\|_{T},$
(4.7)

where $C$ is a constant only depends on $k$ and $i\in\\{1,2\\}$, $i\neq j$.

By applying Lemma 4.3 and Lemma 4.4, we can deduce the following estimates
which are valuable for the convergence analysis of the WG finite element
schemes (3.4).

######  Lemma 4.5.

Let $k\geq 3$, $u\in H^{k+1}(\Omega)$. There exists a constant $C$ such that
the following estimates hold true,

$\displaystyle\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\Delta
u-\mathcal{Q}_{0}\Delta u\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}\leq
C\left(N^{-(k-\frac{3}{2})}+\varepsilon^{-1}N^{-(k-\frac{3}{2})}\ln^{k-\frac{3}{2}}N+\varepsilon^{-1}N^{\frac{1}{2}-\alpha}\right),$
(4.8)
$\displaystyle\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\nabla\left(\Delta
u-\mathcal{Q}_{0}\Delta u\right)\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\leq
C\left(N^{-(k-\frac{5}{2})}+\varepsilon^{-2}N^{-(k-\frac{5}{2})}\ln^{k-\frac{5}{2}}N+\varepsilon^{-2}N^{\frac{3}{2}-\alpha}\right),$
(4.9) $\displaystyle\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\nabla
u-\mathbf{Q}_{N}\nabla u\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}\leq
C\left(N^{-(k-\frac{1}{2})}+N^{-(k-\frac{1}{2})}\ln^{k-\frac{1}{2}}N+N^{\frac{1}{2}-\alpha}\right),$
(4.10)
$\displaystyle\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\nabla(\mathcal{Q}_{0}u-u)\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\leq
C\left(N^{-(k-\frac{1}{2})}+N^{-(k-\frac{1}{2})}\ln^{k-\frac{1}{2}}N+N^{\frac{3}{2}-\alpha}\right),$
(4.11)
$\displaystyle\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\mathcal{Q}_{0}u-u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\leq C\left(N^{-(k+\frac{1}{2})}+\varepsilon
N^{-(k+\frac{1}{2})}\ln^{k+\frac{1}{2}}N+\varepsilon
N^{\frac{1}{2}-\alpha}\right).$ (4.12)

###### Proof.

To derive (4.8), we estimate $\sum_{T\in\mathcal{T}{N}}\|\Delta
u-\mathcal{Q}_{0}\Delta u\|_{\partial T_{i}}$ by breaking down the function
$u$ in Assumption 2.1. Each term in the decomposition will be considered
individually. To begin, we can apply inequality (4.3) to obtain

$\displaystyle\left(\sum_{T\in\Omega_{c}}\|\Delta S-\mathcal{Q}_{0}\Delta
S\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle
Ch^{k-\frac{3}{2}}\cdot N\cdot h\leq CN^{-(k-\frac{3}{2})},$
$\displaystyle\left(\sum_{T\in\Omega_{s}}\|\Delta S-\mathcal{Q}_{0}\Delta
S\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle
CH^{k-\frac{3}{2}}\cdot N\cdot H\leq CN^{-(k-\frac{3}{2})},$
$\displaystyle\left(\sum_{T\in\Omega_{l}}\|\Delta S-\mathcal{Q}_{0}\Delta
S\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle CN\cdot
h_{1}^{\frac{1}{2}}h_{2}^{\frac{1}{2}}\left(h_{j}^{k-\frac{3}{2}}+h_{j}^{-\frac{1}{2}}h_{i}^{k-1}+h_{j}^{\frac{1}{2}}h_{i}^{k-2}\right)$
$\displaystyle\leq$ $\displaystyle CN^{-(k-\frac{3}{2})},$

where we have used the fact that $h_{i}=h$ in the domain $\Omega_{c}$ and
$h_{i}=H$ in the domain $\Omega_{s}$, for $i=1,2$. Next, we provide estimates
only for the sets of element edges parallel to the $x$ axis, as the orther
part follows a similar way. Considering the boundary layer $E_{1}$, we obtain

$\displaystyle\left(\sum_{T\in\Omega_{c}}\|\Delta E_{1}-\mathcal{Q}_{0}\Delta
E_{1}\|_{\partial T_{1}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle
Ch^{k-\frac{3}{2}}\cdot\varepsilon^{-k}\left(\int_{\Omega_{c}}e^{-\frac{2y}{\varepsilon}}\,dxdy\right)^{\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle C\varepsilon^{-k}\cdot
h^{k-\frac{3}{2}}\cdot\varepsilon\ln^{\frac{1}{2}}N$ $\displaystyle\leq$
$\displaystyle C\varepsilon^{-\frac{1}{2}}N^{-(k-\frac{3}{2})}\ln^{k-1}N,$
$\displaystyle\left(\sum_{T\in\Omega_{l}}\|\Delta E_{1}-\mathcal{Q}_{0}\Delta
E_{1}\|_{\partial T_{1}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle
C\left(\right.\varepsilon^{2-k}h_{2}^{k-\frac{3}{2}}+\varepsilon
h_{1}^{k-1}h_{2}^{-\frac{1}{2}}+h_{1}^{k-2}h_{2}^{\frac{1}{2}}+\varepsilon^{-k}h_{2}^{k-\frac{3}{2}}$
$\displaystyle+\varepsilon^{-1}h_{1}^{k-1}h_{2}^{-\frac{1}{2}}+\varepsilon^{-2}h_{1}^{k-2}h_{2}^{\frac{1}{2}}\left.\right)\left(\int_{\Omega_{l}}e^{-\frac{2y}{\varepsilon}}\,dxdy\right)^{\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-\frac{3}{2}}N^{-(k-\frac{3}{2})}\ln^{k-\frac{3}{2}}N\cdot\varepsilon^{\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-1}N^{-(k-\frac{3}{2})}\ln^{k-\frac{3}{2}}N,$

where we have used the inequality (4.3). As for the region $\Omega_{r}^{2}$
that remains in the partition, by applying inequality (4.6), we derive

$\displaystyle\left(\sum_{T\in\Omega_{r}^{2}}\|\Delta
E_{1}-\mathcal{Q}_{0}\Delta E_{1}\|_{\partial
T_{1}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle
C\sum_{T\in\Omega_{r}^{2}}\|\Delta E_{1}\|_{\partial
T_{1}}+\sum_{T\in\Omega_{r}^{2}}h_{2}^{-\frac{1}{2}}\|\mathcal{Q}_{0}\Delta
E_{1}\|_{T}$ $\displaystyle\leq$ $\displaystyle
C(\varepsilon+\varepsilon^{-1})\left(N\int_{0}^{1}e^{-\frac{2y}{\varepsilon}}\,dx\right)^{\frac{1}{2}}+H^{-\frac{1}{2}}\left\|\Delta
E_{1}\right\|_{\Omega_{r}^{2}}$ $\displaystyle\leq$ $\displaystyle
C\left(\varepsilon^{-1}N^{\frac{1}{2}-\alpha}+H^{-\frac{1}{2}}(\varepsilon^{\frac{3}{2}}N^{-\alpha}+\varepsilon^{-\frac{1}{2}}N^{-\alpha})\right)$
$\displaystyle\leq$ $\displaystyle C\varepsilon^{-1}N^{\frac{1}{2}-\alpha}.$

A similar bound can be readily obtained for $E_{2},E_{3}$, and $E_{4}$. Let’s
focus on estimating $E_{41}$ for the concer layers, as the other concer layers
in the decomposition from Assumption 2.1 follow a similar way. By applying
inequalities (4.3) and (4.6), we arrive at

$\displaystyle\left(\sum_{T\in\Omega_{c}}\|\Delta E_{41}-\mathcal{Q}_{0}\Delta
E_{41}\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle
C\varepsilon^{-k}h^{k-\frac{3}{2}}\cdot\left(\int_{\Omega_{c}}e^{-\frac{2x}{\varepsilon}}e^{-\frac{2y}{\varepsilon}}\,dxdy\right)^{\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-\frac{1}{2}}N^{-(k-\frac{3}{2})}\ln^{k-\frac{3}{2}}N.$
$\displaystyle\left(\sum_{T\in\Omega_{r}^{1}}\|\Delta
E_{41}-\mathcal{Q}_{0}\Delta E_{41}\|_{\partial
T_{1}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle
C\sum_{T\in\Omega_{r}^{1}}\left(\left\|\Delta E_{41}\right\|_{\partial
T_{1}}+h_{2}^{-\frac{1}{2}}\left\|\mathcal{Q}_{0}\Delta
E_{41}\right\|_{T}\right)$ $\displaystyle\leq$ $\displaystyle
C\varepsilon^{-1}\left(N\int_{\lambda}^{1-\lambda}e^{-\frac{2x}{\varepsilon}}e^{-\frac{2y}{\varepsilon}}\,dx\right)^{\frac{1}{2}}+h^{-\frac{1}{2}}\left\|\Delta
E_{41}\right\|_{\Omega_{r}^{1}}$ $\displaystyle\leq$ $\displaystyle
C\left(\varepsilon^{-1}\cdot\varepsilon^{\frac{1}{2}}N^{\frac{1}{2}-\alpha}+h^{-\frac{1}{2}}\cdot
N^{-\alpha}\right)$ $\displaystyle\leq$ $\displaystyle
C\varepsilon^{-\frac{1}{2}}N^{\frac{1}{2}-\alpha},$
$\displaystyle\left(\sum_{T\in\Omega_{r}^{2}}\|\Delta
E_{41}-\mathcal{Q}_{0}\Delta E_{41}\|_{\partial
T_{1}}^{2}\right)^{\frac{1}{2}}\leq$ $\displaystyle
C\sum_{T\in\Omega_{r}^{2}}\left(\left\|\Delta E_{41}\right\|_{\partial
T_{1}}+h_{2}^{-\frac{1}{2}}\left\|\mathcal{Q}_{0}\Delta
E_{41}\right\|_{T}\right)$ $\displaystyle\leq$ $\displaystyle
C\varepsilon^{-1}\left(N\int_{0}^{1}e^{-\frac{2x}{\varepsilon}}e^{-\frac{2y}{\varepsilon}}\,dx\right)^{\frac{1}{2}}+H^{-\frac{1}{2}}\left\|\Delta
E_{41}\right\|_{\Omega_{r}^{2}}$ $\displaystyle\leq$ $\displaystyle
C\left(\varepsilon^{-1}\cdot\varepsilon^{\frac{1}{2}}N^{\frac{1}{2}-\alpha}+H^{-\frac{1}{2}}\cdot
N^{-\alpha}\right)$ $\displaystyle\leq$ $\displaystyle
C\varepsilon^{-\frac{1}{2}}N^{\frac{1}{2}-\alpha}.$

By combining the aforementioned proofs, we establish inequality (4.8).
Likewise, for inequalities (4.9), (4.10), (4.11), and (4.12), the proof
follows a similar way as above. Therefore, we omit the detailed explanation. ∎

##  5 Error estimate

In this section, the objective is to provide error estimates for the WG
solution $u_{N}$ obtained from (3.4).

###  5.1 Error equation

We introduce notation used in error analysis to represent the error between
the finite element solution and the $L^{2}$ projection of the exact solution,
as follows

$\displaystyle e_{N}=\mathcal{Q}_{N}u-u_{N}=\\{e_{0},e_{b},\mathbf{e}_{g}\\}.$

The convergence analysis relies on the error equation, and in the following
lemma, we will establish an equation that the error $e_{N}$ satisfies.

######  Lemma 5.1.

Suppose $u$ and $u_{N}=\\{u_{0},u_{b},\mathbf{u}_{g}\\}\in V_{N}$ represent
the solutions of (1.1) and (3.4), respectively. Then for any $v\in V_{N}^{0}$,
we have

$\displaystyle
a(e_{N},v)=l_{1}(u,v)-l_{2}(u,v)+l_{3}(u,v)+s(\mathcal{Q}_{N}u,v),$ (5.1)

where

$\displaystyle
l_{1}(u,v)=\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}\left\langle\Delta
u-\mathcal{Q}_{0}\Delta u,\left(\nabla
v_{0}-\mathbf{v}_{g}\right)\cdot\mathbf{n}\right\rangle_{\partial T},$
$\displaystyle
l_{2}(u,v)=\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}\left\langle\nabla\left(\Delta
u-\mathcal{Q}_{0}\Delta
u\right)\cdot\mathbf{n},v_{0}-v_{b}\right\rangle_{\partial T},$ $\displaystyle
l_{3}(u,v)=\sum_{T\in\mathcal{T}_{N}}\left\langle\left(\nabla
u-\mathbf{Q}_{N}\nabla
u\right)\cdot\mathbf{n},v_{0}-v_{b}\right\rangle_{\partial T}.$

###### Proof.

Using the definition of weak Laplacian (3.2), integration by parts and Lemma
4.1, for any $v\in V_{N}^{0}$, we yield

$\displaystyle\left(\Delta_{w}\mathcal{Q}_{N}u,\Delta_{w}v\right)_{T}$
$\displaystyle=\left(v_{0},\Delta\left(\Delta_{w}\mathcal{Q}_{N}u\right)\right)_{T}+\left\langle\mathbf{v}_{g}\cdot\mathbf{n},\Delta_{w}\mathcal{Q}_{N}u\right\rangle_{\partial
T}-\left\langle
v_{b},\nabla\left(\Delta_{w}\mathcal{Q}_{N}u\right)\cdot\mathbf{n}\right\rangle_{\partial
T}$ $\displaystyle=\left(\Delta
v_{0},\Delta_{w}\mathcal{Q}_{N}u\right)_{T}+\left\langle
v_{0},\nabla\left(\Delta_{w}\mathcal{Q}_{N}u\right)\cdot\mathbf{n}\right\rangle_{\partial
T}-\left\langle\nabla
v_{0}\cdot\mathbf{n},\Delta_{w}\mathcal{Q}_{N}u\right\rangle_{\partial T}$
$\displaystyle\quad+\left\langle\mathbf{v}_{g}\cdot\mathbf{n},\Delta_{w}\mathcal{Q}_{N}u\right\rangle_{\partial
T}-\left\langle
v_{b},\nabla\left(\Delta_{w}\mathcal{Q}_{N}u\right)\cdot\mathbf{n}\right\rangle_{\partial
T}$ $\displaystyle=\left(\Delta
v_{0},\Delta_{w}\mathcal{Q}_{N}u\right)_{T}+\left\langle
v_{0}-v_{b},\nabla\left(\Delta_{w}\mathcal{Q}_{N}u\right)\cdot\mathbf{n}\right\rangle_{\partial
T}-\left\langle\left(\nabla
v_{0}-\mathbf{v}_{g}\right)\cdot\mathbf{n},\Delta_{w}\mathcal{Q}_{N}u\right\rangle_{\partial
T}$ $\displaystyle=\left(\Delta v_{0},\mathcal{Q}_{0}\Delta
u\right)_{T}+\left\langle v_{0}-v_{b},\nabla\left(\mathcal{Q}_{0}\Delta
u\right)\cdot\mathbf{n}\right\rangle_{\partial T}-\left\langle\left(\nabla
v_{0}-\mathbf{v}_{g}\right)\cdot\mathbf{n},\mathcal{Q}_{0}\Delta
u\right\rangle_{\partial T}$ $\displaystyle=\left(\Delta u,\Delta
v_{0}\right)_{T}+\left\langle v_{0}-v_{b},\nabla\left(\mathcal{Q}_{0}\Delta
u\right)\cdot\mathbf{n}\right\rangle_{\partial T}-\left\langle\left(\nabla
v_{0}-\mathbf{v}_{g}\right)\cdot\mathbf{n},\mathcal{Q}_{0}\Delta
u\right\rangle_{\partial T},$

which implies that

$\begin{split}\sum_{T\in\mathcal{T}_{h}}\left(\Delta u,\Delta
v_{0}\right)_{T}=&\left(\Delta_{w}\mathcal{Q}_{N}u,\Delta_{w}v\right)_{N}-\sum_{T\in\mathcal{T}_{h}}\left\langle
v_{0}-v_{b},\nabla\left(\mathcal{Q}_{0}\Delta
u\right)\cdot\mathbf{n}\right\rangle_{\partial T}\\\
&+\sum_{T\in\mathcal{T}_{h}}\left\langle\left(\nabla
v_{0}-\mathbf{v}_{g}\right)\cdot\mathbf{n},\mathcal{Q}_{0}\Delta
u\right\rangle_{\partial T}.\end{split}$ (5.2)

Likewise, we deduce from integration by parts and Lemma 4.2 the following

$\displaystyle\left(\nabla_{w}\mathcal{Q}_{N}u,\nabla_{w}v\right)_{T}$
$\displaystyle=-\left(v_{0},\nabla\cdot\nabla_{w}\mathcal{Q}_{N}u\right)_{T}+\left\langle
v_{b},\nabla_{w}\mathcal{Q}_{N}u\cdot\mathbf{n}\right\rangle_{\partial T}$
$\displaystyle=\left(\nabla
v_{0},\nabla_{w}\mathcal{Q}_{N}u\right)_{T}+\left\langle
v_{b}-v_{0},\nabla_{w}\mathcal{Q}_{N}u\cdot\mathbf{n}\right\rangle_{\partial
T}$ $\displaystyle=\left(\nabla v_{0},\mathbf{Q}_{N}\nabla
u\right)_{T}+\left\langle v_{b}-v_{0},\mathbf{Q}_{N}\nabla
u\cdot\mathbf{n}\right\rangle_{\partial T}$ $\displaystyle=\left(\nabla
v_{0},\nabla u\right)_{T}-\left\langle v_{0}-v_{b},\left(\mathbf{Q}_{N}\nabla
u\right)\cdot\mathbf{n}\right\rangle_{\partial T},$

which implies that

$\begin{split}\sum_{T\in\mathcal{T}_{h}}\left(\nabla u,\nabla
v_{0}\right)_{T}=\left(\nabla_{w}\mathcal{Q}_{N}u,\nabla_{w}v\right)_{N}+\sum_{T\in\mathcal{T}_{h}}\left\langle
v_{0}-v_{b},\mathbf{Q}_{N}\nabla u\cdot\mathbf{n}\right\rangle_{\partial
T}.\end{split}$ (5.3)

Test equation (1.1) with the vector $v_{0}$ of
$v=\\{v_{0},v_{b},\mathbf{v}g\\}\in V_{N}^{0}$, we find

$\displaystyle\varepsilon^{2}\left(\Delta^{2}u,v_{0}\right)-\left(\Delta
u,v_{0}\right)=\left(f,v_{0}\right).$

Using the boundary conditions that $\mathbf{v}_{g}\cdot\mathbf{n}$ and $v_{b}$
vanish on $\partial\Omega$, along with integration by parts, we derive

$\displaystyle\varepsilon^{2}\left(\Delta^{2}u,v_{0}\right)-\left(\Delta
u,v_{0}\right)=$
$\displaystyle\sum_{T\in\mathcal{T}_{h}}\Big{[}\varepsilon^{2}\left(\Delta
u,\Delta v_{0}\right)_{T}-\varepsilon^{2}\left\langle\Delta u,\nabla
v_{0}\cdot\mathbf{n}\right\rangle_{\partial T}+\left(\nabla u,\nabla
v_{0}\right)_{T}$ $\displaystyle+\varepsilon^{2}\left\langle\nabla(\Delta
u)\cdot\mathbf{n},v_{0}\right\rangle_{\partial T}-\left\langle\nabla
u\cdot\mathbf{n},v_{0}\right\rangle_{\partial T}\Big{]}$ $\displaystyle=$
$\displaystyle\sum_{T\in\mathcal{T}_{h}}\Big{[}\varepsilon^{2}\left(\Delta
u,\Delta v_{0}\right)_{T}-\varepsilon^{2}\left\langle\Delta u,\left(\nabla
v_{0}-\mathbf{v}_{g}\right)\cdot\mathbf{n}\right\rangle_{\partial
T}+\left(\nabla u,\nabla v_{0}\right)_{T}$
$\displaystyle+\varepsilon^{2}\left\langle\nabla(\Delta
u)\cdot\mathbf{n},v_{0}-v_{b}\right\rangle_{\partial T}-\left\langle\nabla
u\cdot\mathbf{n},v_{0}-v_{b}\right\rangle_{\partial T}\Big{]}.$

Upon applying the previously mentioned equation together with (5.2) and (5.3),
we get

$\displaystyle\varepsilon^{2}\left(\Delta_{w}\mathcal{Q}_{N}u,\Delta_{w}v\right)_{N}+\left(\nabla_{w}\mathcal{Q}_{N}u,\nabla_{w}v\right)_{N}=\left(f,v_{0}\right)+l_{1}(u,v)-l_{2}(u,v)+l_{3}(u,v).$

By adding $s\left(\mathcal{Q}_{N}u,v\right)$ to both sides of the above
equation, we arrive at

$a\left(\mathcal{Q}_{N}u,v\right)=\left(f,v_{0}\right)+l_{1}(u,v)-l_{2}(u,v)+l_{3}(u,v)+s\left(\mathcal{Q}_{N}u,v\right).$
(5.4)

Subtracting (3.4) from (5.4) yields the error equation as follows

$\displaystyle
a(e_{N},v)=l_{1}(u,v)-l_{2}(u,v)+l_{3}(u,v)+s(\mathcal{Q}_{N}u,v),$

for all $v\in V_{N}^{0}$. This completes the derivation of (5.1). ∎

###  5.2 Error estimate

The following theorem is the estimate for the error function $e_{N}$ in the
triple-bar norm (3.5), which is an $H^{2}$-equivalent norm in $V_{N}^{0}$.

######  Theorem 5.1.

Consider the weak Galerkin finite element solution from (3.4), denoted as
$u_{N}\in V_{N}$. Assuming that $u\in H^{k+1}(\Omega)$, it follows that there
exists a constant $C$ such that

${|||}e_{N}{|||}\leq C\left(N^{-(k-1)}\ln^{k-\frac{3}{2}}N\right).$ (5.5)

###### Proof.

Upon substituting $v=e_{N}$ into the error equation (5.1), we derive the
following equation,

${|||}e_{N}{|||}^{2}=l_{1}(u,e_{N})-l_{2}(u,e_{N})+l_{3}(u,e_{N})+s(\mathcal{Q}_{N}u,e_{N}).$
(5.6)

Using the Cauchy-Schwarz inequality, the meshwidth characteristics (2.2) and
the inequality (4.8), we derive

$\begin{split}l_{1}\left(u,e_{N}\right)&=\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}\left\langle\Delta
u-\mathcal{Q}_{0}\Delta u,\left(\nabla
e_{0}-\mathbf{e}_{g}\right)\cdot\mathbf{n}\right\rangle_{\partial T}\\\
&\leq\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h\|\Delta
u-\mathcal{Q}_{0}\Delta u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{-1}\|\nabla
e_{0}-\mathbf{e}_{g}\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}\\\ &\leq
C\varepsilon
h^{\frac{1}{2}}\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\Delta
u-\mathcal{Q}_{0}\Delta u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}{|||}e_{N}{|||}\\\ &\leq C\varepsilon
h^{\frac{1}{2}}\left(N^{-(k-\frac{3}{2})}+\varepsilon^{-1}N^{-(k-\frac{3}{2})}\ln^{k-\frac{3}{2}}N+\varepsilon^{-1}N^{\frac{1}{2}-\alpha}\right){|||}e_{N}{|||}\\\
&\leq C\left(\varepsilon
N^{-(k-1)}+\varepsilon^{\frac{1}{2}}N^{-(k-1)}\ln^{k-1}N+\varepsilon^{\frac{1}{2}}N^{\frac{1}{2}-\alpha}\right){|||}e_{N}{|||}.\end{split}$
(5.7)

From both the Cauchy-Schwarz inequality and (2.2) and the inequality (4.9), it
can be deduced that

$\begin{split}l_{2}\left(u,e_{N}\right)&=\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}\left\langle\nabla\left(\Delta
u-\mathcal{Q}_{0}\Delta
u\right)\cdot\mathbf{n},e_{0}-e_{b}\right\rangle_{\partial T}\\\
&\leq\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{2}H\|\nabla\left(\Delta
u-\mathcal{Q}_{0}\Delta u\right)\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{-2}H^{-1}\|e_{0}-e_{b}\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\\\ &\leq C\varepsilon
hH^{\frac{1}{2}}\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\nabla\left(\Delta
u-\mathcal{Q}_{0}\Delta u\right)\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}{|||}e_{N}{|||}\\\ &\leq C\varepsilon
hH^{\frac{1}{2}}\left(N^{-(k-\frac{5}{2})}+\varepsilon^{-2}N^{-(k-\frac{5}{2})}\ln^{k-\frac{5}{2}}N+\varepsilon^{-2}N^{\frac{3}{2}-\alpha}\right){|||}e_{N}{|||}\\\
&\leq C\left(\varepsilon
N^{-(k-1)}+N^{-(k-1)}\ln^{k-\frac{3}{2}}N+N^{1-\alpha}\right){|||}e_{N}{|||}.\end{split}$
(5.8)

Similarly, it follows from the Cauchy-Schwarz inequality and (2.2) and (4.10)
that

$\begin{split}l_{3}\left(u,e_{N}\right)&=\sum_{T\in\mathcal{T}_{N}}\left\langle\left(\nabla
u-\mathbf{Q}_{N}\nabla
u\right)\cdot\mathbf{n},e_{0}-e_{b}\right\rangle_{\partial T}\\\
&\leq\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}H\|\nabla
u-\mathbf{Q}_{N}\nabla u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{N}}H^{-1}\|e_{0}-e_{b}\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\\\ &\leq
CH^{\frac{1}{2}}\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\nabla
u-\mathbf{Q}_{N}\nabla u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}{|||}e_{N}{|||}\\\ &\leq
CH^{\frac{1}{2}}\left(N^{-(k-\frac{1}{2})}+N^{-(k-\frac{1}{2})}\ln^{k-\frac{1}{2}}N+N^{\frac{1}{2}-\alpha}\right){|||}e_{N}{|||}\\\
&\leq
C\left(N^{-k}+N^{-k}\ln^{k-\frac{1}{2}}N+N^{-\alpha}\right){|||}e_{N}{|||}.\end{split}$
(5.9)

In the same way, considering $s(\mathcal{Q}_{N}u,e_{N})$, it follows from the
Cauchy-Schwarz inequality and (2.2) and (4.11)-(4.12) that

$\displaystyle\bigg{|}\sum_{i=1,2}\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{-1}\left\langle\nabla\mathcal{Q}_{0}u-\mathbf{Q}_{g}\nabla
u,\nabla e_{0}-\mathbf{e}_{g}\right\rangle_{\partial T_{i}}\bigg{|}$
$\displaystyle\leq\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{-1}\|\nabla\mathcal{Q}_{0}u-\nabla
u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{-1}\|\nabla
e_{0}-\mathbf{e}_{g}\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}$
$\displaystyle\leq C\varepsilon
h^{-\frac{1}{2}}\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\nabla\mathcal{Q}_{0}u-\nabla
u\|_{\partial T_{i}}^{2}\right)^{\frac{1}{2}}{|||}e_{N}{|||}$
$\displaystyle\leq C\varepsilon
h^{-\frac{1}{2}}\left(N^{-(k-\frac{1}{2})}+N^{-(k-\frac{1}{2})}\ln^{k-\frac{1}{2}}N+N^{\frac{3}{2}-\alpha}\right){|||}e_{N}{|||}$
$\displaystyle\leq
C\left(\varepsilon^{\frac{1}{2}}N^{-(k-1)}+\varepsilon^{\frac{1}{2}}N^{-(k-1)}\ln^{k-1}N+\varepsilon^{\frac{1}{2}}N^{2-\alpha}\right){|||}e_{N}{|||},$
(5.10)
$\displaystyle\bigg{|}\sum_{i=1,2}\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{-2}H^{-1}\left\langle\mathcal{Q}_{0}u-\mathcal{Q}_{b}u,e_{0}-e_{b}\right\rangle_{\partial
T_{i}}\bigg{|}$
$\displaystyle\leq\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{-2}H^{-1}\|\mathcal{Q}_{0}u-u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{N}}\varepsilon^{2}h^{-2}H^{-1}\|e_{0}-e_{b}\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}$ $\displaystyle\leq C\varepsilon
h^{-1}H^{-\frac{1}{2}}\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\mathcal{Q}_{0}u-u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}{|||}e_{N}{|||}$ $\displaystyle\leq
C\varepsilon h^{-1}H^{-\frac{1}{2}}\left(N^{-(k+\frac{1}{2})}+\varepsilon
N^{-(k+\frac{1}{2})}\ln^{k+\frac{1}{2}}N+\varepsilon
N^{\frac{1}{2}-\alpha}\right)$ $\displaystyle\leq
C\left(N^{-(k-1)}+\varepsilon N^{-(k-1)}\ln^{k-\frac{1}{2}}N+\varepsilon
N^{2-\alpha}\right),$ (5.11)

and

$\displaystyle\bigg{|}\sum_{i=1,2}\sum_{T\in\mathcal{T}_{N}}H^{-1}\left\langle\mathcal{Q}_{0}u-\mathcal{Q}_{b}u,e_{0}-e_{b}\right\rangle_{\partial
T_{i}}\bigg{|}$
$\displaystyle\leq\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}H^{-1}\|\mathcal{Q}_{0}u-u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{N}}H^{-1}\|e_{0}-e_{b}\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}$ $\displaystyle\leq
CH^{-\frac{1}{2}}\sum_{i=1,2}\left(\sum_{T\in\mathcal{T}_{N}}\|\mathcal{Q}_{0}u-u\|_{\partial
T_{i}}^{2}\right)^{\frac{1}{2}}{|||}e_{N}{|||}$ $\displaystyle\leq
CH^{-\frac{1}{2}}\left(N^{-(k+\frac{1}{2})}+\varepsilon
N^{-(k+\frac{1}{2})}\ln^{k+\frac{1}{2}}N+\varepsilon
N^{\frac{1}{2}-\alpha}\right)$ $\displaystyle\leq C\left(N^{-k}+\varepsilon
N^{-k}\ln^{k+\frac{1}{2}}N+\varepsilon N^{1-\alpha}\right).$ (5.12)

Substituting (5.7)-(5.2) into (5.6) yields

$\displaystyle{|||}e_{N}{|||}^{2}\leq
C\left(N^{-(k-1)}\ln^{k-\frac{3}{2}}N\right){|||}e_{N}{|||}.$

which implies (5.5). This completes the proof of the theorem. ∎

##  6 Numerical Experiments

In this section we compute two numerical examples where we select $\alpha=k+1$
for creating the Shishkin mesh. Consider the singularly perturbed fourth-order
problem that seeks solution $u=u(x,y)$ satisfying

$\displaystyle\varepsilon^{2}\Delta^{2}u-\Delta u$
$\displaystyle=f,\quad\text{in}~{}\Omega,$ $\displaystyle
u=\partial_{\mathbf{n}}u$ $\displaystyle=0,\quad\text{on}~{}\partial\Omega,$

where $\Omega=(0,1)^{2}$, $f$ and $\varepsilon$ will choose later. Tables 1-6
display the errors ${|||}\mathcal{Q}_{N}u-u_{N}{|||}$ for several different
$\varepsilon$ and $N$.

######  Example 6.1.

Choose $f(x,y)$ such that the exact solution is $u(x,y)=g(x)g(y)$, where

$g(x)=\frac{1}{2}\left[\sin(\pi
x)+\frac{\pi\varepsilon}{1-e^{-1/\varepsilon}}\left(e^{-x/\varepsilon}+e^{(x-1)/\varepsilon}-1-e^{-1/\varepsilon}\right)\right].$

The numerical results are shown in Table 1-3.

In the case of $k=3$, the error and the order of convergence on Shishkin mesh
and uniform mesh are listed in Table 1 and Table 2, respectively. Compared
with the numerical results of the WG method for the problem on uniform mesh,
our results with Shishkin mesh are better and an $\varepsilon$-independent
asymptotically optimal order of convergence is achieved in all cases.

Table 1: Numerical results for Example 6.1 on Shishkin mesh. $\varepsilon$ | $N=8$ | $N=16$ | $N=32$ | $N=64$ | $N=128$
---|---|---|---|---|---
1e-00 | 1.01e-03 | 2.61e-04 | 6.58e-05 | 1.65e-05 | 4.12e-06
| 1.96 | 1.99 | 2.00 | 2.00 | –
1e-01 | 3.77e-03 | 1.06e-03 | 2.75e-04 | 6.94e-05 | 1.74e-05
| 1.83 | 1.95 | 1.99 | 2.00 | –
1e-02 | 1.17e-02 | 6.43e-03 | 3.03e-03 | 1.25e-03 | 4.59e-04
| 0.86 | 1.09 | 1.28 | 1.44 | –
1e-03 | 3.81e-03 | 2.08e-03 | 9.73e-04 | 4.00e-04 | 1.46e-04
| 0.87 | 1.10 | 1.28 | 1.45 | –
1e-04 | 1.22e-03 | 6.59e-04 | 3.08e-04 | 1.27e-04 | 4.64e-05
| 0.89 | 1.10 | 1.28 | 1.45 | –
1e-05 | 4.18e-04 | 2.09e-04 | 9.75e-05 | 4.01e-05 | 1.47e-05
| 1.00 | 1.10 | 1.28 | 1.45 | –
1e-06 | 2.09e-04 | 6.71e-05 | 3.09e-05 | 1.27e-05 | 4.64e-06
| 1.64 | 1.12 | 1.28 | 1.45 | –
1e-07 | 1.74e-04 | 2.44e-05 | 9.84e-06 | 4.01e-06 | 1.47e-06
| 2.84 | 1.31 | 1.29 | 1.45 | –
Table 2: Numerical results for Example 6.1 on uniform mesh. $\varepsilon$ | $N=8$ | $N=16$ | $N=32$ | $N=64$ | $N=128$
---|---|---|---|---|---
1e-00 | 1.01e-03 | 2.61e-04 | 6.58e-05 | 1.65e-05 | 4.12e-06
| 1.96 | 1.99 | 2.00 | 2.00 | –
1e-01 | 3.77e-03 | 1.06e-03 | 2.75e-04 | 6.94e-05 | 1.74e-05
| 1.83 | 1.95 | 1.99 | 2.00 | –
1e-02 | 3.87e-02 | 2.06e-02 | 8.03e-03 | 2.62e-03 | 7.53e-04
| 0.91 | 1.36 | 1.62 | 1.80 | –
1e-03 | 1.24e-02 | 1.60e-02 | 1.73e-02 | 1.44e-02 | 8.63e-03
| -0.37 | -0.12 | 0.26 | 0.74 | –
1e-04 | 1.30e-03 | 1.83e-03 | 2.57e-03 | 3.56e-03 | 4.74e-03
| -0.49 | -0.49 | -0.47 | -0.41 | –
1e-05 | 1.33e-04 | 1.85e-04 | 2.62e-04 | 3.70e-04 | 5.21e-04
| -0.47 | -0.50 | -0.50 | -0.49 | –
1e-06 | 1.96e-05 | 1.87e-05 | 2.62e-05 | 3.71e-05 | 5.25e-05
| 0.07 | -0.49 | -0.50 | -0.50 | –
1e-07 | 1.30e-05 | 2.24e-06 | 2.63e-06 | 3.71e-06 | 5.25e-06
| 2.54 | -0.23 | -0.50 | -0.50 | –

In the case of $k=4$, the numerical results are shown in Table 3, which has
yielded the asymptotically optimal order.

Table 3: Numerical results for Example 6.1 on Shishkin mesh. $\varepsilon$ | $N=8$ | $N=16$ | $N=32$ | $N=64$
---|---|---|---|---
1e-00 | 3.07e-05 | 3.90e-06 | 4.89e-07 | 6.12e-08
| 2.98 | 3.00 | 3.00 | –
1e-01 | 3.92e-04 | 5.35e-05 | 6.84e-06 | 8.61e-07
| 2.87 | 2.97 | 2.99 | –
1e-02 | 6.08e-03 | 2.56e-03 | 8.25e-04 | 2.11e-04
| 1.25 | 1.63 | 1.97 | –
1e-03 | 1.98e-03 | 8.29e-04 | 2.66e-04 | 6.77e-05
| 1.26 | 1.64 | 1.97 | –
1e-04 | 6.28e-04 | 2.63e-04 | 8.43e-05 | 2.15e-05
| 1.26 | 1.64 | 1.97 | –
1e-05 | 1.99e-04 | 8.32e-05 | 2.67e-05 | 6.79e-06
| 1.26 | 1.64 | 1.97 | –
1e-06 | 6.33e-05 | 2.63e-05 | 8.44e-06 | 2.24e-06
| 1.27 | 1.64 | 1.91 | –

######  Example 6.2.

Let $g(x)$ be as in Example 6.1 and set

$p(y)=2y(1-y^{2})+\varepsilon\left[ld(1-2y)-3\frac{q}{l}+\left(\frac{3}{l}-d\right)e^{-y/\varepsilon}+\left(\frac{3}{l}+d\right)e^{(y-1)/\varepsilon}\right],$

with $l=1-e^{-1/\varepsilon}$, $q=2-l$ and $d=1/(q-2\varepsilon l)$. Then
choose $f(x,y)$ such that the exact solution of (1.1) is $u(x,y)=g(x)p(y)$.
The numerical results are shown in Table 4-6.

The error and the order of convergence on Shishkin mesh and uniform mesh with
$k=3$ are listed in Table 4 and Table 5, respectively. And also we get better
results on Shishkin mesh than uniform mesh.

Table 4: Numerical results for Example 6.2 on Shishkin mesh. $\varepsilon$ | $N=8$ | $N=16$ | $N=32$ | $N=64$ | $N=128$
---|---|---|---|---|---
1e-00 | 1.66e-04 | 4.24e-05 | 1.07e-05 | 2.67e-06 | 6.67e-07
| 1.97 | 1.99 | 2.00 | 2.00 | –
1e-01 | 6.48e-03 | 1.82e-03 | 4.72e-04 | 1.19e-04 | 2.99e-05
| 1.83 | 1.94 | 1.99 | 2.00 | –
1e-02 | 2.10e-02 | 1.16e-02 | 5.44e-03 | 2.24e-03 | 8.25e-04
| 0.86 | 1.09 | 1.28 | 1.44 | –
1e-03 | 6.86e-03 | 3.74e-03 | 1.75e-03 | 7.20e-04 | 2.64e-04
| 0.87 | 1.10 | 1.28 | 1.45 | –
1e-04 | 2.18e-03 | 1.19e-03 | 5.55e-04 | 2.28e-04 | 8.35e-05
| 0.88 | 1.10 | 1.28 | 1.45 | –
1e-05 | 7.13e-04 | 3.76e-04 | 1.76e-04 | 7.22e-05 | 2.64e-05
| 0.92 | 1.10 | 1.28 | 1.45 | –
1e-06 | 2.87e-04 | 1.20e-04 | 5.56e-05 | 2.28e-05 | 8.36e-06
| 1.27 | 1.11 | 1.28 | 1.45 | –
1e-07 | 2.00e-04 | 4.01e-05 | 1.76e-05 | 7.29e-06 | 2.64e-06
| 2.32 | 1.19 | 1.27 | 1.46 | –
Table 5: Numerical results for Example 6.2 on uniform mesh. $\varepsilon$ | $N=8$ | $N=16$ | $N=32$ | $N=64$ | $N=128$
---|---|---|---|---|---
1e-00 | 1.66E-04 | 4.24E-05 | 1.07E-05 | 2.67E-06 | 6.67E-07
| 1.97 | 1.99 | 2.00 | 2.00 | –
1e-01 | 6.48E-03 | 1.82E-03 | 4.72E-04 | 1.19E-04 | 2.99E-05
| 1.83 | 1.94 | 1.99 | 2.00 | –
1e-02 | 6.96E-02 | 3.70E-02 | 1.44E-02 | 4.71E-03 | 1.35E-03
| 0.91 | 1.36 | 1.62 | 1.80 | –
1e-03 | 2.22E-02 | 2.88E-02 | 3.12E-02 | 2.60E-02 | 1.56E-02
| -0.37 | -0.12 | 0.26 | 0.74 | –
1e-04 | 2.35E-03 | 3.30E-03 | 4.63E-03 | 6.41E-03 | 8.54E-03
| -0.49 | -0.49 | -0.47 | -0.41 | –
1e-05 | 2.37E-04 | 3.34E-04 | 4.72E-04 | 6.66E-04 | 9.38E-04
| -0.49 | -0.50 | -0.50 | -0.49 | –
1e-06 | 2.85E-05 | 3.35E-05 | 4.73E-05 | 6.69E-05 | 9.45E-05
| -0.23 | -0.50 | -0.50 | -0.50 | –
1e-07 | 1.45E-05 | 3.62E-06 | 4.73E-06 | 6.69E-06 | 9.46E-06
| 2.00 | -0.39 | -0.50 | -0.50 | –

In the case of $k=4$, the numerical results are shown in Table 6, which has
also yielded the asymptotically optimal order, independent of $\varepsilon$.

Table 6: Numerical results for Example 6.2 on Shishkin mesh. $\varepsilon$ | $N=8$ | $N=16$ | $N=32$ | $N=64$
---|---|---|---|---
1e-00 | 3.84e-06 | 4.86e-07 | 6.09e-08 | 7.62e-09
| 2.98 | 3.00 | 3.00 | –
1e-01 | 6.93e-04 | 9.45e-05 | 1.21e-05 | 1.52e-06
| 2.87 | 2.97 | 2.99 | –
1e-02 | 1.09e-02 | 4.60e-03 | 1.48e-03 | 3.79e-04
| 1.25 | 1.63 | 1.97 | –
1e-03 | 3.57e-03 | 1.49e-03 | 4.79e-04 | 1.22e-04
| 1.26 | 1.64 | 1.97 | –
1e-04 | 1.13e-03 | 4.74e-04 | 1.52e-04 | 3.87e-05
| 1.26 | 1.64 | 1.97 | –
1e-05 | 3.58e-04 | 1.50e-04 | 4.81e-05 | 1.22e-05
| 1.26 | 1.64 | 1.97 | –
1e-06 | 1.14e-04 | 4.74e-05 | 2.53e-05 | 3.90e-06
| 1.26 | 0.90 | 2.70 | –

######  Remark 6.1.

The numerical examples show that we obtain the asymptotically optimal error
estimate. In fact, the $\ln^{k-\frac{3}{2}}N$ in the error results will have
an impact on the convergence order, and with the refinement of the Shishkin
mesh, the influence of $\ln^{k-\frac{3}{2}}N$ on the convergence order will
gradually become smaller. The convergence order is reduced by about half order
when $k=3$ and by one order in the case of $k=4$. Our numerical results also
confirm this well.

##  7 Conclusion

In this paper, we use the weak Galerkin finite element method to solve the
singularly perturbed fourth-order boundary value problem in 2D domain. By
constructing the Shishkin mesh which is suitable for the problem, we give the
corresponding numerical algorithm and the error estimate in the $H^{2}$
discrete norm. Compared with solving the problem on the uniform mesh, the WG
method obtains better results on the Shishkin mesh, which is verified by the
numerical results. Moreover, the results of our numerical experiments are
consistent with the error estimation theory, and the asymptotically optimal
convergence order is obtained.

## Statements and Declarations

Funding. This work was supported by the National Natural Science Foundation of
China (Grant No. 12101039, 12271208).

Data Availability. The code used in this work will be made available upon
request to the authors.

Competing Interests. The authors have no relevant financial or non-financial
interests to disclose.

## References

  * [1] N. S. Bahvalov, On the optimization of the methods for solving boundary value problems in the presence of a boundary layer, Ž. Vyčisl. Mat i Mat. Fiz., In Russian, 9 (1969), pp. 841–859.
  * [2] S. C. Brenner, T. Gudi, M. Neilan, and L.-y. Sung, $C^{0}$ penalty methods for the fully nonlinear monge-ampère equation, Math. Comp., 80 (2011), pp. 1979–1995.
  * [3] S. C. Brenner and M. Neilan, A $C^{0}$ interior penalty method for a fourth order elliptic singular perturbation problem, SIAM J. Numer. Anal., 49 (2011), pp. 869–892.
  * [4] G. Chen and X. Xie, A robust weak Galerkin finite element method for linear elasticity with strong symmetric stresses, Comput. Methods Appl. Math., 16 (2016), pp. 389–408.
  * [5] H. Chen, S. Chen, and L. Xiao, Uniformly convergent $C^{0}$-nonconforming triangular prism element for fourth-order elliptic singular perturbation problem, Numer. Methods Partial Differential Equations, 30 (2014), pp. 1785–1796.
  * [6] L. Chen, J. Wang, Y. Wang, and X. Ye, An auxiliary space multigrid preconditioner for the weak Galerkin method, Comput. Math. Appl., 70 (2015), pp. 330–344.
  * [7] P. Constantinou and C. Xenophontos, An $hp$ finite element method for a 4th order singularly perturbed boundary value problem in two dimensions, Comput. Math. Appl., 74 (2017), pp. 1565–1575.
  * [8] M. Cui, X. Ye, and S. Zhang, A modified weak Galerkin finite element method for the biharmonic equation on polytopal meshes, Commun. Appl. Math. Comput., 3 (2021), pp. 91–105.
  * [9] M. Cui and S. Zhang, On the uniform convergence of the weak Galerkin finite element method for a singularly-perturbed biharmonic equation, J. Sci. Comput., 82 (2020), pp. Paper No. 5, 15.
  * [10] S. Franz and H.-G. Roos, Robust error estimation in energy and balanced norms for singularly perturbed fourth order problems, Comput. Math. Appl., 72 (2016), pp. 233–247.
  * [11] S. Franz, H.-G. Roos, and A. Wachtel, A $C^{0}$ interior penalty method for a singularly-perturbed fourth-order elliptic problem on a layer-adapted mesh, Numer. Methods Partial Differential Equations, 30 (2014), pp. 838–861.
  * [12] E. H. Georgoulis, Discontinuous Galerkin methods on shape-regular and anisotropic meshes, PhD thesis, University of Oxford D. Phil. Thesis, 2003.
  * [13] J. Guzmán, D. Leykekhman, and M. Neilan, A family of non-conforming elements and the analysis of Nitsche’s method for a singularly perturbed fourth order problem, Calcolo, 49 (2012), pp. 95–125.
  * [14] G. Harper, J. Liu, S. Tavener, and B. Zheng, Lowest-order weak Galerkin finite element methods for linear elasticity on rectangular and brick meshes, J. Sci. Comput., 78 (2019), pp. 1917–1941.
  * [15] S. C. Hongru Chen, Uniformly convergent nonconforming element for 3D fourth order elliptic singular perturbation problem, J. Comput. Math., 32 (2014), pp. 687–695.
  * [16] X. Hu, L. Mu, and X. Ye, A weak Galerkin finite element method for the Navier-Stokes equations, J. Comput. Appl. Math., 362 (2019), pp. 614–625.
  * [17] W. Huang and Y. Wang, Discrete maximum principle for the weak Galerkin method for anisotropic diffusion problems, Commun. Comput. Phys., 18 (2015), pp. 65–90.
  * [18] X. Liu, J. Li, and Z. Chen, A weak Galerkin finite element method for the Navier-Stokes equations, J. Comput. Appl. Math., 333 (2018), pp. 442–457.
  * [19] X. Meng and M. Stynes, Convergence analysis of the Adini element on a Shishkin mesh for a singularly perturbed fourth-order problem in two dimensions, Adv. Comput. Math., 45 (2019), pp. 1105–1128.
  * [20] J. J. H. Miller, E. O’Riordan, and G. I. Shishkin, Fitted numerical methods for singular perturbation problems, World Scientific Publishing Co., Inc., River Edge, NJ, 1996.
  * [21] L. Mu, J. Wang, and X. Ye, A stable numerical algorithm for the Brinkman equations by weak Galerkin finite element methods, J. Comput. Phys., 273 (2014), pp. 327–342.
  * [22] L. Mu, J. Wang, X. Ye, and S. Zhang, A weak Galerkin finite element method for the Maxwell equations, J. Sci. Comput., 65 (2015), pp. 363–386.
  * [23] T. K. Nilssen, X.-C. Tai, and R. Winther, A robust nonconforming $H^{2}$-element, Math. Comp., 70 (2001), pp. 489–505.
  * [24] P. Panaseti, A. Zouvani, N. Madden, and C. Xenophontos, A $C^{1}$-conforming hp finite element method for fourth order singularly perturbed boundary value problems, Appl. Numer. Math., 104 (2016), pp. 81–97.
  * [25] H.-G. Roos, M. Stynes, and L. Tobiska, Robust numerical methods for singularly perturbed differential equations, vol. 24 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, second ed., 2008.
  * [26] G. I. Shishkin, Grid approximation of singularly perturbed elliptic and parabolic equations, Keldysh Institute of Applied Mathematics, Moscow, In Russian, 1990.
  * [27] C. Wang, J. Wang, R. Wang, and R. Zhang, A locking-free weak Galerkin finite element method for elasticity problems in the primal formulation, J. Comput. Appl. Math., 307 (2016), pp. 346–366.
  * [28] J. Wang and X. Ye, A weak Galerkin finite element method for second-order elliptic problems, J. Comput. Appl. Math., 241 (2013), pp. 103–115.
  * [29] J. Wang and X. Ye, A weak Galerkin finite element method for the stokes equations, Adv. Comput. Math., 42 (2016), pp. 155–174.
  * [30] J. Wang, X. Ye, Q. Zhai, and R. Zhang, Discrete maximum principle for the $P_{1}$-$P_{0}$ weak Galerkin finite element approximations, J. Comput. Phys., 362 (2018), pp. 114–130.
  * [31] L. Wang, Y. WU, and X. Xie, Uniformly stable rectangular elements for fourth order elliptic singular perturbation problems, Numer. Methods Partial Differential Equations, 29 (2013), pp. 721–737.
  * [32] R. Wang, X. Wang, Q. Zhai, and R. Zhang, A weak Galerkin finite element scheme for solving the stationary Stokes equations, J. Comput. Appl. Math., 302 (2016), pp. 171–185.
  * [33] R. Wang, R. Zhang, X. Zhang, and Z. Zhang, Supercloseness analysis and polynomial preserving recovery for a class of weak Galerkin methods, Numer. Methods Partial Differential Equations, 34 (2018), pp. 317–335.
  * [34] W. Wang, X. Huang, K. Tang, and R. Zhou, Morley-Wang-Xu element methods with penalty for a fourth order elliptic singular perturbation problem, Adv. Comput. Math., 44 (2018), pp. 1041–1061.
  * [35] X. Wang, Q. Zhai, and R. Zhang, The weak Galerkin method for solving the incompressible Brinkman flow, J. Comput. Appl. Math., 307 (2016), pp. 13–24.
  * [36] Y. Wang, Y. Li, and X. Meng, An upwind finite volume element method on a Shishkin mesh for singularly perturbed convection-diffusion problems, J. Comput. Appl. Math., 438 (2024), pp. Paper No. 115493, 20.
  * [37] Z. Wang, R. Wang, and J. Liu, Robust weak Galerkin finite element solvers for Stokes flow based on a lifting operator, Comput. Math. Appl., 125 (2022), pp. 90–100.
  * [38] P. Xie, D. Shi, and H. Li, A new robust $C^{0}$-type nonconforming triangular element for singular perturbation problems, Appl. Math. Comput., 217 (2010), pp. 3832–3843.
  * [39] S.-Y. Yi, A lowest-order weak Galerkin method for linear elasticity, J. Comput. Appl. Math., 350 (2019), pp. 286–298.
  * [40] Q. Zhai, R. Zhang, and L. Mu, A new weak Galerkin finite element scheme for the Brinkman model, Commun. Comput. Phys., 19 (2016), pp. 1409–1434.
  * [41] J. Zhang and X. Liu, Uniform convergence of a weak Galerkin finite element method on Shishkin mesh for singularly perturbed convection-diffusion problems in 2D, Appl. Math. Comput., 432 (2022), pp. Paper No. 127346, 12.
  * [42] J. Zhang and X. Liu, Uniform convergence of a weak Galerkin method for singularly perturbed convection-diffusion problems, Math. Comput. Simulation, 200 (2022), pp. 393–403.
  * [43] J. Zhang, K. Zhang, J. Li, and X. Wang, A weak Galerkin finite element method for the Navier-Stokes equations, Commun. Comput. Phys., 23 (2018), pp. 706–746.
  * [44] R. Zhang and Q. Zhai, A weak Galerkin finite element scheme for the biharmonic equations by using polynomials of reduced order, J. Sci. Comput., 64 (2015), pp. 559–585.
  * [45] T. Zhang and T. Lin, The weak Galerkin finite element method for incompressible flow, J. Math. Anal. Appl., 464 (2018), pp. 247–265.
  * [46] P. Zhu and S. Xie, A uniformly convergent weak Galerkin finite element method on Shishkin mesh for 1D convection-diffusion problem, J. Sci. Comput., 85 (2020), pp. Paper No. 34, 22.