# Wavelet-based Edge Multiscale Finite Element Methods for Singularly
Perturbed Convection-Diffusion Equations

Shubin Fu Eastern Institute for Advanced Study, Eastern Institute of
Technology, Ningbo, Zhejiang 315200, P. R. China<EMAIL_ADDRESS>Eric
Chung Department of Mathematics, The Chinese University of Hong Kong, Hong
Kong Special Administrative Region<EMAIL_ADDRESS>Guanglian Li
Corresponding author. Department of Mathematics, The University of Hong Kong,
Pokfulam Road, Hong Kong Special Administrative Region<EMAIL_ADDRESS>
###### Abstract

We propose a novel efficient and robust Wavelet-based Edge Multiscale Finite
Element Method (WEMsFEM) motivated by [13, 22] to solve the singularly
perturbed convection-diffusion equations. The main idea is to first establish
a local splitting of the solution over a local region by a local bubble part
and local Harmonic extension part, and then derive a global splitting by means
of Partition of Unity. This facilitates a representation of the solution as a
summation of a global bubble part and a global Harmonic extension part, where
the first part can be computed locally in parallel. To approximate the second
part, we construct an edge multiscale ansatz space locally with hierarchical
bases as the local boundary data that has a guaranteed approximation rate
without higher regularity requirement on the solution. The key innovation of
this proposed WEMsFEM lies in a provable convergence rate with little
restriction on the mesh size or the regularity of the solution. Its
convergence rate with respect to the computational degree of freedom is
rigorously analyzed, which is verified by extensive 2-d and 3-d numerical
tests.

Keywords  Multiscale method, Convection-diffusion, Singularly perturbed,
Partition of Unity, Hierarchical bases

##  1 Introduction

This paper focuses on solving the singularly perturbed convection-diffusion
equation. Let $D\subset\mathbb{R}^{d}$ ($d=1,2,3$) be an open bounded domain
with a Lipschitz boundary. Then we seek a function $u\in V:=H^{1}_{0}(D)$ such
that

$\displaystyle-\epsilon\Delta u+b\cdot\nabla u$ $\displaystyle=f\quad\text{in
}D,$ (1.1) $\displaystyle u$ $\displaystyle=0\quad\text{on }\partial D,$

where the force term $f\in L^{2}(D)$, $0<\epsilon\ll 1$ is a singular
perturbation parameter and the velocity field $b\in
L^{\infty}(D;\mathbb{R}^{d})$. We assume that the velocity field $b$ is
incompressible, i.e., $\nabla\cdot b=0$. The focus of this paper is on the
convection-dominated regime with large Péclet number
$\mathrm{Pe}:=\left\|b\right\|_{L^{\infty}\left(D\right)}/\epsilon$.

Many important applications such as metamaterials, porous media and fluid
mechanics involves multiple scales. For example, metamaterials are periodic
artificial material with the size of each cell being on a subwavelength scale
that possesses desirable properties which are not available for the natural
material. Another example arises from singularly perturbed convection-
diffusion problem (1.1), whose solution has high oscillation in the form of
boundary layer or internal layers. Due to this disparity of scales, the
classical numerical treatment becomes prohibitively expensive and even
intractable for many multiscale applications. Nonetheless, motivated by the
broad spectrum of practical applications, a large number of multiscale model
reduction techniques, e.g., Multiscale Finite Element Methods (MsFEMs),
Heterogeneous Multiscale Methods (HMMs), Variational Multiscale Methods (VMM),
flux norm approach, Generalized Multiscale Finite Element Methods (GMsFEMs)
and Localized Orthogonal Decomposition (LOD), have been proposed in the
literature [2, 7, 10, 11, 18, 19, 26] over the last few decades. They have
achieved great success in the efficient and accurate simulation of problems
involving multiple scales.

The singularly perturbed convection-diffusion problem (1.1) is numerically
challenging due to the presence of boundary layers or internal layers with
width of $\mathcal{O}(\epsilon)$, which are usually exponential. To obtain a
reasonable macroscale solution, one has to resolve all these microscale
information using a very fine mesh, which is computational extremely
infeasible. To reduce the computational complexity, one category of numerical
methods is based upon stabilization, albeit at the loss of accuracy especially
near the layers [3, 8, 21]. Another category is to incorporate the information
of the layers in the form of local multiscale basis functions into the trial
space [8, 23, 28]. These methods usually achieve good accuracy near the
layers.

In this paper, we develop a novel Wavelet-based Edge Multiscale Finite Element
Method (WEMsFEM) [13, 14, 15, 22] which relies an edge multiscale ansatz space
that has a good approximation property with very low requirement on the
solution, and thereby bypass the requirement on a very fine mesh. This edge
multiscale ansatz space is constructed in each local domain by solving a
homogeneous convection-dominated diffusion equation with hierarchical bases as
Dirichlet boundary data and then use partition of unity functions to glue the
local basis functions together as global multiscale basis functions. WEMsFEM
shares a certain similarity with the Exponentially Convergent Multiscale
Finite Element Method (ExpMsFEM) [4, 5, 6] in that both methods utilize the
approximation over the coarse skeleton to define the local multiscale basis
functions. The most striking differences are twofold. Firstly, WEMsFEM can be
regarded as a variant of the Partition of Unity Finite Element Methods (PUFEM)
[27], which decomposes the global approximation space into summation of local
approximation spaces, while ExpMsFEM uses edge localization and coupling to
communicate local and global approximations. Secondly, the local multiscale
space in WEMsFEM is established by wavelets, and it inherits the intrinsic
hierarchical structure of the wavelets. Consequently, there is no need for a
further model reduction or deriving important modes by means of Singular Value
Decomposition (SVD).

###  1.1 Main contributions

  1. 1.

One key ingredient of WEMsFEM is to transfer the approximation properties over
the coarse skeleton to its interior by means of the transposition method [25].
This is rigorously justified for Problem (1.1) in Section A.

  2. 2.

To approximate the solution $u$ to Problem (1.1), we construct an edge
multiscale ansatz space with level parameter $\ell$ which encodes key
microscale information of the solution such that a guaranteed approximation
rate is achieved, which is established in Lemma 5.2.

  3. 3.

The converge rate of WEMsFEMs is proved rigorously without extra regularity
requirement on the solution $u$ or demanding grids, see Theorem 5.2.

  4. 4.

We conduct extensive numerical experiments both in 2-d and 3-d with large
Péclet number and even with high oscillation in the diffusion coefficients,
and all experiments demonstrate fast convergence rate of WEMsFEMs with respect
to the level parameter $\ell$.

###  1.2 Organization

We summarize in Section 2 the basic results on singularly perturbed
convection-diffusion problem (1.1), and demonstrate the main idea of our
approach. Then we present the basic notation and approximation properties of
hierarchical bases in Section 3.2 for the sake of completeness, which is
followed by the construction of multiscale ansatz space $V_{\text{ms},\ell}$.
Based upon this, we propose in Section 4 the novel Wavelet-based Edge
Multiscale Finite Element Methods (WEMsFEMs), which are the key findings of
this paper. Their theoretical and numerical performance are presented in
Sections 5 and 6. Finally, we conclude the paper with several remarks in
Section 7.

##  2 Preliminaries

We recap in this section important properties and numerical challenges of
Problem (1.1), and its finite element approximation.

###  2.1 Weak formulation and a priori  estimate

We define the bilinear form $a(\cdot,\cdot)$ on $V\times V$ associated to
(1.1) by

$\displaystyle a(u,v):=\epsilon\int_{D}\nabla u\cdot\nabla
v\,\mathrm{d}x+\int_{D}(b\cdot\nabla u)\,v\,\mathrm{d}x\qquad\text{for all
}u,v\in V.$ (2.1)

Since $\nabla\cdot b=0$, an integration by parts implies that the bilinear
form $a(\cdot,\cdot)$ is $V$-elliptic, i.e.,

$\displaystyle a(v,v)$ $\displaystyle=\epsilon\left\|\nabla
v\right\|_{L^{2}\left(D\right)}^{2}\qquad\text{ for all }v\in V.$ (2.2)

Nevertheless, this $V$-ellipticity vanishes when $\epsilon\to 0$, which makes
the standard finite element method unstable.

An application of the Poincarè-Friedrichs inequality leads to the boundedness
of the bilinear form,

$\displaystyle|a(u,v)|\leq{\rm C}(\rm{D},\rm{b})\left\|\nabla
u\right\|_{L^{2}\left(D\right)}\left\|\nabla v\right\|_{L^{2}\left(D\right)}.$

Here, the positive constant ${\rm C}(\rm{D},\rm{b})$ may depend on the
diameter of the domain $D$ and the $L^{\infty}$-norm of the velocity $b$.

The weak formulation for problem (1.1) is to find $u\in V$ such that

$\displaystyle a(u,v)=(f,v)_{D}\quad\text{for all }v\in V.$ (2.3)

Furthermore, this weak solution $u$ is Hölder continuous,

######  Proposition 2.1 (Hölder estimate).

For any $\alpha\in[0,1)$, there holds

$\displaystyle\|u\|_{C^{\alpha}(D)}\leq\mathrm{C}(\alpha)\left(\mathrm{Pe}\left\|\nabla
u\right\|_{L^{2}\left(D\right)}+\epsilon^{-1}\left\|f\right\|_{L^{2}\left(D\right)}\right).$
(2.4)

Here, the positive constant $\mathrm{C}(\alpha)$ is independent of $\epsilon$.

###### Proof.

Note that Problem (1.1) is equivalent to

$\displaystyle-\Delta u=\epsilon^{-1}\left(-{b}\cdot\nabla u+f\right)$
$\displaystyle\qquad\text{ in }D$ , $\displaystyle u=0$
$\displaystyle\qquad\text{ on }\partial D$ .

On the one hand, the right hand side can be bounded by

$\displaystyle\left\|-\epsilon^{-1}{b}\cdot\nabla
u+\epsilon^{-1}f\right\|_{L^{2}\left(D\right)}\leq\mathrm{Pe}\left\|\nabla
u\right\|_{L^{2}\left(D\right)}+\epsilon^{-1}\left\|f\right\|_{L^{2}\left(D\right)}.$
(2.5)

On the other hand, the Hölder estimate for Laplacian [16] leads to

$\displaystyle\|u\|_{C^{\alpha}(D)}\leq\mathrm{C}(\alpha)\left\|-\epsilon^{-1}{b}\cdot\nabla
u+\epsilon^{-1}f\right\|_{L^{2}\left(D\right)}$

for some constant $\mathrm{C}(\alpha)$ depends only on the domain $D$ and its
boundary $\partial D$. Together with (2.5), we obtain the desired assertion. ∎

###  2.2 Finite Element discretization

To discretize problem (1.1), let $\mathcal{T}^{H}$ be a regular partition of
the domain $D$ into finite elements with a mesh size $H$. For ease of
presentation, we denote the mesh Péclet number $\rm{Pe}_{H,b,\epsilon}$ of the
coarse mesh $\mathcal{T}^{H}$ by

$\rm{Pe}_{H,b,\epsilon}:=H\|b\|_{L^{\infty}(D)}/\epsilon.$

The vertices of $\mathcal{T}^{H}$ are denoted by $\\{O_{i}\\}_{i=1}^{N}$, with
$N$ being the total number of coarse nodes. The coarse neighborhood associated
with the node $O_{i}$ is denoted by

$\omega_{i}:=\bigcup\\{K_{j}\in\mathcal{T}^{H}:\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ O_{i}\in\overline{K}_{j}\\}.$ (2.6)

The overlapping constant $C_{\mathrm{ov}}$ is defined by

$\displaystyle
C_{\mathrm{ov}}:=\max\limits_{K\in\mathcal{T}^{H}}\\#\\{O_{i}:K\subset\omega_{i}\text{
for }i=1,2,\cdots,N\\}.$ (2.7)

We assume that $\epsilon\ll 1$ is a small parameter, and that the coarse mesh
$\mathcal{T}^{H}$ fails to resolve the small parameter, i.e.,
$\rm{Pe}_{H,b,\epsilon}\gg 1$.

Over the coarse mesh $\mathcal{T}^{H}$, let $V_{H}$ be the conforming
piecewise linear finite element space,

$V_{H}:=\\{v\in\mathcal{C}(D):V|_{T}\in\mathbb{P}_{1}\text{ for all
}T\in\mathcal{T}^{H}\\},$

where $\mathbb{P}_{1}$ denotes the space of linear polynomials. Then the
classical Galerkin approximation of Problem (2.3) reads as finding $u_{H}\in
V_{H}$, satisfying

$\displaystyle a(u_{H},v_{H})=(f,v_{H})_{D}\quad\text{ for all }v_{H}\in
V_{H}.$ (2.8)

A straightforward calculation leads to the following estimate

$\displaystyle\left\|\nabla(u-u_{H})\right\|_{L^{2}\left(D\right)}$
$\displaystyle\lesssim(1+\rm{Pe}_{H,b,\epsilon})\min\limits_{w\in
V_{H}}\left\|\nabla(u-w)\right\|_{L^{2}\left(D\right)}.$

This implies that if the solution $u$ has sufficient regularity $u\in
H^{2}(D)$, then a standard interpolation property leads to

$\displaystyle\left\|\nabla(u-u_{H})\right\|_{L^{2}\left(D\right)}$
$\displaystyle\lesssim(1+\rm{Pe}_{H,b,\epsilon})H\|u\|_{H^{2}(D)}.$

Here, and throughout this paper, the notation $A\lesssim B$ means $A\leq cB$
for some benign constant that does not depend on the small parameter
$\epsilon$ and the mesh size $H$. On each coarse neighborhood $\omega_{i}$,
the local bilinear form $a_{\omega_{i}}(\cdot,\cdot)$ is defined analogous to
(2.1) by replacing the global domain $D$ with $\omega_{i}$.

Even though the hidden constant in the estimate above is independent of the
small parameter $\epsilon$, $\|u\|_{H^{2}(D)}$ relies on $\epsilon$.
Furthermore, there is no convergence unless $u\in H^{1+\delta}(D)$ for some
$\delta>0$.

###  2.3 Main objective

The main goal of this paper is to construct a multiscale ansatz space
$V_{\text{ms},\ell}$ with a non-negative level parameter $\ell\in\mathbb{N}$
that incorporates certain a priori information of the exact solution $u$, such
that the following quantity

$\displaystyle\min\limits_{w\in
V_{\text{ms},\ell}}\left\|\nabla(u-w)\right\|_{L^{2}\left(D\right)}$ (2.9)

has a certain decay rate with respect to the number of basis functions in
$V_{\text{ms},\ell}$, and in the meanwhile, this does not require further
regularity on the solution $u$.

##  3 Edge multiscale ansatz space $V_{\text{ms},\ell}$

We introduce in this section the methodology for the construction of
$V_{\text{ms},\ell}$ to obtain a guaranteed decay rate in (2.9) without extra
regularity requirement on the solution $u$.

###  3.1 Local-Global splitting

To start, note that the solution $u$ from (1.1) satisfies the following
equation

$\displaystyle\mathcal{L}_{i}u$ $\displaystyle=f\quad$ $\displaystyle\text{ in
}\omega_{i},$

which can be split into the summation of two parts, namely

$\displaystyle u|_{\omega_{i}}=u^{i,\mathrm{I}}+u^{i,\mathrm{II}}.$ (3.1)

Here, the two components $u^{i,\mathrm{I}}$ and $u^{i,\mathrm{II}}$ are
respectively given by

$\left\\{\begin{aligned} \mathcal{L}_{i}u^{i,\mathrm{I}}&=f\quad&&\text{ in
}\omega_{i}\\\ u^{i,\mathrm{I}}&=0\quad&&\text{ on
}\partial\omega_{i},\end{aligned}\right.$ (3.2)

and

$\left\\{\begin{aligned} \mathcal{L}_{i}u^{i,\mathrm{II}}&=0\quad&&\text{ in
}\omega_{i}\\\ u^{i,\mathrm{II}}&=u\quad&&\text{ on
}\partial\omega_{i}.\end{aligned}\right.$ (3.3)

This implies that $u^{i,\mathrm{I}}$ contains the local information of the
source term $f$, and $u^{i,\mathrm{II}}$ encodes the global information over
the coarse skeleton $\mathcal{E}^{H}:=\cup\\{\partial
K:K\in\mathcal{T}^{H}\\}$. Note that $u^{i,\mathrm{I}}$ can be solved in each
coarse neighborhood $\omega_{i}$ for $i=1,2,\cdots,N$ in parallel.

Similar to Proposition 2.1, we obtain, for any $\alpha\in[0,1)$, there holds

$\displaystyle\|u^{i,\mathrm{II}}\|_{C^{\alpha}(\omega_{i})}\leq\mathrm{C}(\alpha)\mathrm{Pe}\left\|\nabla
u^{i,\mathrm{II}}\right\|_{L^{2}\left(\omega_{i}\right)}.$ (3.4)

Here, the positive constant $\mathrm{C}(\alpha)$ is independent of $\epsilon$
and $H$, and it can change value from context to context.

Next, we introduce a partition of unity $\\{\chi_{i}\\}_{i=1}^{N}$ subordinate
to the cover $\\{\omega_{i}\\}_{i=1}^{N}$ to induce a global decomposition
from (3.1). Recall that $O_{i}$ denotes the ith coarse node in the coarse mesh
$\mathcal{T}^{H}$. We assume this Partition of Unity
$\\{\chi_{i}\\}_{i=1}^{N}$ satisfies the following properties:

$\left\\{\begin{aligned} {\text{supp}(\chi_{i})}&\subset\omega_{i}\\\
\sum\limits_{i=1}^{N}\chi_{i}&=1\text{ in }D\\\
\|\chi_{i}\|_{L^{\infty}(\omega_{i})}&\leq C_{\infty}\\\
\|\nabla\chi_{i}\|_{L^{\infty}(\omega_{i})}&\leq
C_{\text{G}}H^{-1}\end{aligned}\right.$

for some positive constants $C_{\infty}$ and $C_{\text{G}}$, which is named
the $(C_{\mathrm{ov}},C_{\infty},C_{\text{G}})$ partition of unity [27]. For
example, the standard bilinear basis functions on the coarse mesh
$\mathcal{T}^{H}$ are the partition of unity function with
$C_{\infty}=C_{\text{G}}=1$. Note that there are many other alternatives for
the partition of unity functions besides using this bilinear basis functions.
For instance, one can utilize the flat top type of partition of unity
functions, cf. [17].

Consequently, we can represent $u$ as a summation of local parts,

$\displaystyle
u=\left(\sum_{i=1}^{N}\chi_{i}\right)u=\sum_{i=1}^{N}\left(\chi_{i}u|_{\omega_{i}}\right).$

Inserting the local splitting (3.1), we derive

$\displaystyle
u=\sum_{i=1}^{N}\chi_{i}\left(u^{i,\mathrm{I}}+u^{i,\mathrm{II}}\right)$
$\displaystyle=\sum_{i=1}^{N}\chi_{i}u^{i,\mathrm{I}}+\sum_{i=1}^{N}\chi_{i}u^{i,\mathrm{II}}$
$\displaystyle:=u^{\mathrm{I}}+u^{\mathrm{II}}.$ (3.5)

Here, the first term

$\displaystyle u^{\mathrm{I}}:=\sum_{i=1}^{N}\chi_{i}u^{i,\mathrm{I}}$ (3.6)

denotes the global bubble part which contains information of the source term
$f$, and the second term

$\displaystyle u^{\mathrm{II}}:=\sum_{i=1}^{N}\chi_{i}u^{i,\mathrm{II}}.$
(3.7)

is the global Harmonic extension part that encodes global information over the
coarse skeleton $\mathcal{E}^{H}$.

Since the global bubble part $u^{\mathrm{I}}$ can be retrieved locally, our
main objective (2.9) is reduced to seeking a good approximative space
$V_{\text{ms},\ell}$ for $u^{\mathrm{II}}$ in the sense that the following
quantity has a certain decay rate,

$\displaystyle\min\limits_{w\in
V_{\text{ms},\ell}}\left\|\nabla(u^{\mathrm{II}}-w)\right\|_{L^{2}\left(D\right)}.$
(3.8)

To obtain this approximative space $V_{\text{ms},\ell}$, we aim to design
global multiscale basis functions taking the same form as (3.7). Since the
partition of unity $\\{\chi_{i}\\}_{i=1}^{N}$ is known, then we only need to
design local multiscale basis functions to approximate $u^{i,\mathrm{II}}$
(3.3). Note that $u^{i,\mathrm{II}}$ is completely determined by the boundary
data $u|_{\partial\omega_{i}}$, hence the question is reduced to seeking a
good approximation to $u|_{\partial\omega_{i}}$. Note also that if we can
obtain a proper approximation to $u|_{\partial\omega_{i}}$, this approximation
property can transfer into the local neighborhood $\omega_{i}$ thanks to the a
priori estimate established in Section A.

In short, we need to approximate $u|_{\partial\omega_{i}}$, which, however,
has no guaranteed regularity higher than $C^{1}(\partial\omega_{i})$ as stated
in Proposition 2.1. This motivates us to introduce hierarchical bases, since
their approximation properties do not rely on the regularity of the target
function. We present in Section 3.2 the basic definition of hierarchical bases
with level parameter $\ell\in\mathbb{N}^{+}$ on a unit interval, and then
define the approximative space $V_{\text{ms},\ell}$ in Section 3.3.

###  3.2 Hierarchical subspace splitting over $I=:[0,1]$

In this section, we introduce the hierarchical bases on the unit interval
$I:=[0,1]$, which facilitate hierarchically splitting the space $L^{2}(I)$
[9].

$x_{0,0}$$x_{0,1}$$x_{1,0}$$x_{1,1}$$x_{1,2}$$x_{2,0}$$x_{2,2}$$x_{2,4}$$x_{2,1}$$x_{2,3}$$W_{0}$$W_{1}$$W_{2}$$\oplus$$\oplus$$x_{1,0}$$x_{1,1}$$x_{1,2}$$x_{2,0}$$x_{2,1}$$x_{2,2}$$x_{2,3}$$x_{2,4}$$V_{1}$$V_{2}$
Figure 1: An illustration of hierarchical bases for $\ell=0,1,2$.

Let the level parameter and the mesh size be $\ell$ and $h_{\ell}:=2^{-\ell}$
with $\ell\in\mathbb{N}$, respectively. Then the grid points on level $\ell$
are

$x_{\ell,j}=j\times h_{\ell},\qquad 0\leq j\leq 2^{\ell}.$

We can define the basis functions on level $\ell$ by

$\psi_{\ell,j}(x)=\left\\{\begin{aligned} &1-|x/h_{\ell}-j|,&&\text{ if
}x\in[(j-1)h_{\ell},(j+1)h_{\ell}]\cap[0,1],\\\ &0,&&\text{
otherwise.}\end{aligned}\right.$

Define the set on each level $\ell$ by

$B_{\ell}:=\Bigg{\\{}j\in\mathbb{N}\Bigg{|}\begin{aligned}
&j=1,\cdots,2^{\ell}-1,j\,\rm{is\,odd},&&\rm{if}\,\ell>0\\\
&j=0,1,&&\rm{if}\,\ell=0\end{aligned}\Bigg{\\}}.$

The subspace of level $\ell$ is

$W_{\ell}:=\text{span}\\{\psi_{\ell,j}:\quad j\in B_{\ell}\\}.$

We denote $V_{\ell}$ as the subspace in $L^{2}(I)$ up to level $\ell$, which
is defined by the direct sum of subspaces

$V_{\ell}:=\oplus_{m\leq\ell}W_{m}.$

Consequently, this yields the hierarchical structure of the subspace
$V_{\ell}$, namely,

$V_{0}\subset V_{1}\subset\cdots\subset V_{\ell}\subset V_{\ell+1}\cdots$

Furthermore, the following hierarchical decomposition of the space $L^{2}(I)$
holds

$L^{2}(I)=\lim_{\ell\to\infty}\oplus_{m\leq\ell}W_{m}.$

Note that one can derive the hierarchical decomposition of the space
$L^{2}(I^{d})$ for $d>1$ by means of tensor product. Note further that we will
use the subspace $V_{\ell}$ to approximate the restriction of the exact
solution $u$ on each coarse edge.

Let $I_{\ell}:C(I)\to V_{\ell}$ be nodal interpolation for each level
$\ell\geq 0$, then there holds [29],

$\displaystyle\|v-I_{\ell}v\|_{L^{2}(I)}$ $\displaystyle\lesssim
2^{-\ell}\|v\|_{C(I)}.$ (3.9)

###  3.3 Edge multiscale ansatz space

Now we are ready to define the approximative space $V_{\text{ms},\ell}$.

To begin, we first define the linear space over the boundary of each coarse
neighborhood $\partial\omega_{i}$. Let the level parameter $\ell\in\mathbb{N}$
be fixed, and let $\Gamma_{i}^{k}$ with $k=1,2,3,4$ be a partition of
$\partial\omega_{i}$ with no mutual intersection, i.e.,
$\cup_{k=1}^{4}\overline{\Gamma_{i}^{k}}=\partial\omega_{i}$ and
$\Gamma_{i}^{k}\cap\Gamma_{i}^{k^{\prime}}=\emptyset$ if $k\neq k^{\prime}$.
Furthermore, we denote $V_{i,\ell}^{k}\subset C(\Gamma_{i}^{k})$ as the linear
space spanned by hierarchical bases up to level $\ell$ on each coarse edge
$\Gamma_{i}^{k}$, then the local edge space $V_{i,\ell}$ defined over
$\partial\omega_{i}$ is the smallest linear space having $V_{i,\ell}^{k}$ as a
subspace. Let $\\{\psi_{i,\ell}^{j}\\}_{j=1}^{2^{\ell+2}}$ and
$\\{x_{i,\ell}^{j}\\}_{j=1}^{2^{\ell+2}}$ be the nodal basis functions and the
associated nodal points for $V_{i,\ell}$, then we can represent the local edge
space $V_{i,\ell}$ by

$\displaystyle V_{i,\ell}:=\text{span}\left\\{\psi_{i,\ell}^{j}:1\leq j\leq
2^{\ell+2}\right\\}.$ (3.10)

We depict the grid points $\\{x_{i,\ell}^{j}\\}_{j=1}^{2^{\ell+2}}$ of the
hierarchical bases over $\partial\omega_{i}$ in Figure 2 for level parameter
$\ell$ ranging from 0 to 2.

$\omega_{i}$$\ell=2$$\ell=1$$\ell=0$ Figure 2: Grid points for $\ell=0,1,2$
over $\partial\omega_{i}$.

Next, we introduce the local multiscale space over each coarse neighborhood
$\omega_{i}$, which is defined by

$\displaystyle\mathcal{L}^{-1}_{i}(V_{i,\ell}):=\text{span}\left\\{\mathcal{L}^{-1}_{i}(\psi_{i,\ell}^{j}):1\leq
j\leq 2^{\ell+2}\right\\}.$ (3.11)

Here, $\mathcal{L}^{-1}_{i}(\psi_{i,\ell}^{j}):=v\in H^{1}(\omega_{i})$,
satisfying

$\left\\{\begin{aligned} \mathcal{L}_{i}v&:=-\epsilon\Delta v+b\cdot\nabla
v=0\quad\mbox{in }\omega_{i},\\\ v&=\psi_{i,\ell}^{j}\quad\mbox{on
}\partial\omega_{i}.\end{aligned}\right.$ (3.12)

By this construction, the dimension of $\mathcal{L}^{-1}_{i}(V_{i,\ell})$
equals $2^{\ell+2}$.

Finally, the edge multiscale ansatz space is defined by the Partition of Unity
$\\{\chi_{i}\\}_{i=1}^{N}$,

$\displaystyle
V_{\text{ms},\ell}:=\text{span}\left\\{\chi_{i}\mathcal{L}^{-1}_{i}(\psi_{i,\ell}^{j}):\,\,1\leq
i\leq N\text{ and }1\leq j\leq 2^{\ell+2}\right\\}.$ (3.13)

##  4 Wavelet-based Edge Multiscale Finte Element Method (WEMsFEM)

We present in this section the Wavelet-based Edge Multiscale Finte Element
Method (WEMsFEM) to solve (1.1) based upon the global splitting (3.5).

To get a good approximation to $u^{\mathrm{II}}$, we solve for
$u_{\text{ms},\ell}^{\mathrm{II}}\in V_{\text{ms},\ell}$, s.t.,

$\displaystyle
a\left(u_{\text{ms},\ell}^{\mathrm{II}},v_{\text{ms},\ell}\right)=(f,v_{\text{ms},\ell})_{D}-a(u^{\mathrm{I}},v_{\text{ms},\ell})\quad\text{
for all }v_{\text{ms},\ell}\in V_{\text{ms},\ell}.$ (4.1)

Then the multiscale solution $u_{\text{ms},\ell}$ is defined as the summation
of these two parts,

$\displaystyle
u_{\text{ms},\ell}:=u^{\mathrm{I}}+u_{\text{ms},\ell}^{\mathrm{II}}.$ (4.2)

Our main algorithm is summarized in Algorithm 1.

######  Remark 4.1 (Local decomposition (3.1)).

The local decomposition (3.1) is different from the one in [13] in that we
have to solve $u^{i,\mathrm{I}}$. This arises from estimate (5.2), and hence
$u^{\mathrm{I}}$ is not negligible.

Data: The level parameter $\ell\in\mathbb{N}$; coarse neighborhood
$\omega_{i}$ and its four coarse edges $\Gamma_{i}^{k}$ with $k=1,2,3,4$; the
subspace $V_{\ell,i}^{k}\subset L^{2}(\Gamma_{i}^{k})$ up to level $\ell$ on
each coarse edge $\Gamma_{i}^{k}$.

Result: $u_{\text{ms},\ell}$

1 Construct the local edge space $V_{i,\ell}$ (3.10);

2 Calculate the local multiscale space $\mathcal{L}^{-1}_{i}(V_{i,\ell})$
(3.11);

3 Construct the global multiscale space $V_{\text{ms},\ell}$ (3.13);

4 Solve for the local bubble function $u^{i,\mathrm{I}}$ from (3.2) for
$i=1,\cdots,N$;

5 Calculate the global bubble function $u^{\mathrm{I}}$ from (3.6);

6 Solve for $u_{\text{ms},\ell}^{\mathrm{II}}\in V_{\text{ms},\ell}$ from
(2.8);

Obtain $u_{\text{ms},\ell}$ from (4.2).

Algorithm 1 Wavelet-based Edge Multiscale Finte Element Method (WEMsFEM)

To analyze the convergence of Algorithm 1, we next introduce the global
interpolation operator $\mathcal{P}_{\ell}$ of level $\ell$: $V\cap C(D)\to
V_{\text{ms},\ell}$. Since the edge multiscale ansatz space
$V_{\text{ms},\ell}$ is generated by the local multiscale space
$\mathcal{L}^{-1}(V_{i,\ell})$ by means of the partition of unity (3.13), we
only need to define local interpolation operator $\mathcal{P}_{i,\ell}$ of
level $\ell$: $C(\partial\omega_{i})\to\mathcal{L}^{-1}(V_{i,\ell})$.

Recall that $\\{\psi_{i,k}^{j}\\}_{j=1}^{2^{k+2}}$ and
$\\{x_{i,k}^{j}\\}_{j=1}^{2^{k+2}}$ are the hierarchical bases and their nodal
points on $\partial\omega_{i}$ of level $k$, then
$\mathcal{P}_{i,\ell}:C(\partial\omega_{i})\to\mathcal{L}^{-1}(V_{i,\ell})$ is
defined by

$\displaystyle\mathcal{P}_{i,\ell}v:=\sum\limits_{k=0}^{\ell}\sum\limits_{j=1}^{2^{k+2}}v(x_{i,k}^{j})\mathcal{L}_{i}^{-1}(\psi_{i,k}^{j}).$
(4.3)

Note that $\mathcal{P}_{i,\ell}v|_{\partial\omega_{i}}$ is the interpolation
by hierarchical bases up to level $\ell$, this is because

$\displaystyle\mathcal{P}_{i,\ell}v|_{\partial\omega_{i}}$
$\displaystyle=\sum\limits_{k=0}^{\ell}\sum\limits_{j=1}^{2^{k+2}}v(x_{i,k}^{j})\mathcal{L}_{i}^{-1}(\psi_{i,k}^{j})|_{\partial\omega_{i}}$
$\displaystyle=\sum\limits_{k=0}^{\ell}\sum\limits_{j=1}^{2^{k+2}}v(x_{i,k}^{j})(\psi_{i,k}^{j}).$

Note that any $v\in V\cap C(D)$ can be expressed by

$\displaystyle v=\sum_{i=1}^{N}\chi_{i}v|_{\omega_{i}},$

then the global interpolation $\mathcal{P}_{\ell}$ of level $\ell$: $V\cap
C(D)\to V_{\text{ms},\ell}$ can be defined by means of the local projection
(4.3),

$\displaystyle\mathcal{P}_{\ell}(v):=\sum_{i=1}^{N}\chi_{i}(\mathcal{P}_{i,\ell}v|_{\partial\omega_{i}}).$
(4.4)

##  5 Error estimate

We present in this section the convergence analysis for Algorithm 1, which is
divided into two steps. First, we derive the approximation properties of the
edge multiscale ansatz space $V_{\text{ms},\ell}$ defined in (3.13), which is
based upon establishing the local approximation properties of the local
multiscale space $\mathcal{L}^{-1}_{i}(V_{i,\ell})$ in each coarse
neighborhood $\omega_{i}$. Second, we provide a quasi-optimal result by Céa
type of estimate in Theorem 5.2.

Recall that $\\{\chi_{i}\\}_{i=1}^{N}$ is a
$(C_{\mathrm{ov}},C_{\infty},C_{\text{G}})$ partition of unity subordinate to
the cover $\\{\omega_{i}\\}_{i=1}^{N}$, then a slight modification of [27,
Theorem 2.1] leads to the following result,

######  Theorem 5.1 (Local-global approximation).

Let $v_{i}\in H^{1}(\omega_{i})$ for $i=1,\cdots,N$, and let

$\displaystyle v=\sum_{i=1}^{N}\chi_{i}v_{i}\in H^{1}(D).$

Then there holds

$\displaystyle\|v\|_{L^{2}(D)}$
$\displaystyle\leq\sqrt{C_{\mathrm{ov}}}C_{\infty}\left(\sum_{i=1}^{N}\|v_{i}\|_{L^{2}(\omega_{i})}^{2}\right)^{1/2},$
$\displaystyle\|\nabla v\|_{L^{2}(D)}$
$\displaystyle\leq\sqrt{2C_{\mathrm{ov}}}\left(\sum_{i=1}^{N}C_{\text{G}}^{2}H^{-2}\|v_{i}\|_{L^{2}(\omega_{i})}^{2}+\|\chi_{i}\nabla
v_{i}\|_{L^{2}(\omega_{i})}^{2}\right)^{1/2}.$ (5.1)

######  Remark 5.1.

The global estimate in (5.1) is slightly different from [27, Theorem 2.1]
since the local estimate $\|\chi_{i}\nabla v_{i}\|_{L^{2}(\omega_{i})}$ is
utilized instead of $\|\nabla v_{i}\|_{L^{2}(\omega_{i})}$. The main reason is
that our proposed local multiscale space $\mathcal{L}_{i}^{-1}(V_{i,\ell})$ do
not have approximation rate in the sense that the following quantity

$\displaystyle\min\limits_{w\in\mathcal{L}_{i}^{-1}(V_{i,\ell})}\left\|\nabla(u^{i,\mathrm{II}}-w)\right\|_{L^{2}\left(\omega_{i}\right)}$

may not admit a certain decay rate with respect to the level parameter $\ell$.
However, we can prove the decay rate of the following modified quantity,

$\displaystyle\min\limits_{w\in\mathcal{L}_{i}^{-1}(V_{i,\ell})}\left\|\chi_{i}\nabla(u^{i,\mathrm{II}}-w)\right\|_{L^{2}\left(\omega_{i}\right)}.$

######  Lemma 5.1 (A priori estimate for $u^{\mathrm{II}}$).

Let $u^{i,\mathrm{II}}$ and $u^{\mathrm{II}}$ be defined in (3.3) and (3.5)
and let $f\in L^{2}(D)$. Then it holds

$\displaystyle\|\nabla u^{i,\mathrm{II}}\|_{L^{2}(D)}$
$\displaystyle\lesssim\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(\omega_{i}\right)}+\|\nabla
u\|_{L^{2}(\omega_{i})}$ $\displaystyle\|\nabla u^{\mathrm{II}}\|_{L^{2}(D)}$
$\displaystyle\lesssim\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(D\right)}+\|\nabla
u\|_{L^{2}(D)}.$

###### Proof.

On the one hand, multiplying (3.2) by $u^{i,\mathrm{I}}$ and then taking its
integral over $\omega_{i}$ lead to

$\displaystyle\epsilon\left\|\nabla
u^{i,\mathrm{I}}\right\|_{L^{2}\left(\omega_{i}\right)}^{2}=\int_{\omega_{i}}fu^{i,\mathrm{I}}\,\mathrm{d}x.$

Notice that $u^{i,\mathrm{I}}=0$ on $\partial\omega_{i}$, therefore, an
application of the Poincaré inequality reveals

$\displaystyle\left\|u^{i,\mathrm{I}}\right\|_{L^{2}\left(\omega_{i}\right)}+H\left\|\nabla
u^{i,\mathrm{I}}\right\|_{L^{2}\left(\omega_{i}\right)}\lesssim\rm{Pe}_{H,b,\epsilon}H\left\|f\right\|_{L^{2}\left(\omega_{i}\right)}.$

Together with the local splitting (3.1), this implies

$\displaystyle\left\|\nabla
u^{i,\mathrm{II}}\right\|_{L^{2}\left(\omega_{i}\right)}$
$\displaystyle=\left\|\nabla(u-u^{i,\mathrm{I}})\right\|_{L^{2}\left(\omega_{i}\right)}\leq\left\|\nabla
u^{i,\mathrm{I}}\right\|_{L^{2}\left(\omega_{i}\right)}+\left\|\nabla
u\right\|_{L^{2}\left(\omega_{i}\right)}$
$\displaystyle\lesssim\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(\omega_{i}\right)}+\|\nabla
u\|_{L^{2}(\omega_{i})}.$

On the other hand, using Theorem 5.1, we obtain

$\displaystyle\left\|u^{\mathrm{I}}\right\|_{L^{2}\left(D\right)}+H\left\|\nabla
u^{\mathrm{I}}\right\|_{L^{2}\left(D\right)}\lesssim\rm{Pe}_{H,b,\epsilon}H\left\|f\right\|_{L^{2}\left(D\right)}.$
(5.2)

Then the global splitting (3.5) indicates

$\displaystyle\left\|\nabla
u^{\mathrm{II}}\right\|_{L^{2}\left(D\right)}=\left\|\nabla(u-u^{\mathrm{I}})\right\|_{L^{2}\left(D\right)}\leq\left\|\nabla
u\right\|_{L^{2}\left(D\right)}+\left\|\nabla
u^{\mathrm{I}}\right\|_{L^{2}\left(D\right)},$

which, together with (5.2), leads to the desired assertion. ∎

Next, we present the approximation properties of the global interpolation
operator $\mathcal{P}_{\ell}$ (4.4),

######  Lemma 5.2 (Approximation properties of $\mathcal{P}_{\ell}$).

Assume that $f\in L^{2}(D)$ and the level parameter $\ell\in\mathbb{N}$ is
non-negative. Let $u\in V$ be the solution to Problem (1.1) and the global
Harmonic extension $u^{\mathrm{II}}$ be defined in (3.7). Then there holds

$\displaystyle\left\|u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}}\right\|_{L^{2}\left(D\right)}$
$\displaystyle\lesssim(\rm{Pe}_{H,b,\epsilon}+1)^{2}\epsilon^{1/2}2^{-\ell}\left(\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(D\right)}+\left\|\nabla
u\right\|_{L^{2}\left(D\right)}\right),$ (5.3)
$\displaystyle\left\|\nabla(u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}})\right\|_{L^{2}\left(D\right)}$
$\displaystyle\lesssim(\rm{Pe}_{H,b,\epsilon}+1)^{5/2}H^{-1}\epsilon^{1/2}2^{-\ell}\left(\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(D\right)}+\left\|\nabla
u\right\|_{L^{2}\left(D\right)}\right).$ (5.4)

###### Proof.

Let ${e}_{\mathrm{}}:=u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}}$ be
the global error, then the property of the partition of unity of
$\\{\chi_{i}\\}_{i=1}^{N}$, together with (3.7), leads to

${e}_{\mathrm{}}=\sum\limits_{i=1}^{N}\chi_{i}{e}^{i}_{\mathrm{}}\qquad\text{
with
}\qquad{e}^{i}_{\mathrm{}}:=u^{i,\mathrm{II}}-\mathcal{P}_{i,\ell}u^{i,\mathrm{II}}.$

Our proof is composed of three steps.

Step 1. Estimate the local error ${e}^{i}_{\mathrm{}}$ over each local
boundary $\partial\omega_{i}$. Using the nodal interpolation error for
hierarchical bases (3.9) and a scaling argument, we obtain

$\displaystyle\|{e}^{i}_{\mathrm{}}\|_{L^{2}(\partial\omega_{i})}$
$\displaystyle=\|u^{i,\mathrm{II}}-\mathcal{P}_{i,\ell}u^{i,\mathrm{II}}\|_{L^{2}(\partial\omega_{i})}$
$\displaystyle\leq 2^{-\ell}H\|u^{i,\mathrm{II}}\|_{C(\partial\omega_{i})}.$

Then an application of the Hölder estimate (3.4) of the local solution
$u^{i,\mathrm{II}}$ yields,

$\displaystyle\|{e}^{i}_{\mathrm{}}\|_{L^{2}(\partial\omega_{i})}$
$\displaystyle\leq
2^{-\ell}\mathrm{C}(\alpha)\rm{Pe}_{H,b,\epsilon}\left\|\nabla
u^{i,\mathrm{II}}\right\|_{L^{2}\left(\omega_{i}\right)}.$ (5.5)

Step 2. Estimate the local error ${e}^{i}_{\mathrm{}}$ over each coarse
neighborhood $\omega_{i}$ mainly by the transposition method established in
Section A.

Note that each local error ${e}^{i}_{\mathrm{}}$ satisfies the following
equation,

$\left\\{\begin{aligned}
\mathcal{L}_{i}{e}^{i}_{\mathrm{}}:=-\epsilon\Delta{e}^{i}_{\mathrm{}}+b\cdot\nabla{e}^{i}_{\mathrm{}}&=0&&\text{
in }\omega_{i},\\\ v&={e}^{i}_{\mathrm{}}&&\text{ on
}\partial\omega_{i}.\end{aligned}\right.$

Theorem A.1 indicates that the local error ${e}^{i}_{\mathrm{}}$ inside of
$\omega_{i}$ can be bounded by its boundary data,

$\displaystyle\|{e}^{i}_{\mathrm{}}\|_{L^{2}(\omega_{i})}$
$\displaystyle\leq{\rm C}_{{\rm
weak}}\epsilon^{1/2}(1+{\rm{Pe}_{H,b,\epsilon}})\|{e}^{i}_{\mathrm{}}\|_{L^{2}(\partial\omega_{i})}$
$\displaystyle\left\|\chi_{i}\nabla{e}^{i}_{\mathrm{}}\right\|_{L^{2}\left(\omega_{i}\right)}$
$\displaystyle\leq{\rm C}_{{\rm
weak}}H^{-1}\epsilon^{1/2}(1+{\rm{Pe}_{H,b,\epsilon}})^{3/2}\|{e}^{i}_{\mathrm{}}\|_{L^{2}(\partial\omega_{i})}.$

Then together with (5.5), this leads to

$\displaystyle\|{e}^{i}_{\mathrm{}}\|_{L^{2}(\omega_{i})}$
$\displaystyle\leq{\rm C}_{{\rm
weak}}\epsilon^{1/2}(1+{\rm{Pe}_{H,b,\epsilon}})^{2}2^{-\ell}\mathrm{C}(\alpha)\left\|\nabla
u^{i,\mathrm{II}}\right\|_{L^{2}\left(\omega_{i}\right)}$ (5.6)
$\displaystyle\left\|\chi_{i}\nabla{e}^{i}_{\mathrm{}}\right\|_{L^{2}\left(\omega_{i}\right)}$
$\displaystyle\leq{\rm C}_{{\rm
weak}}H^{-1}\epsilon^{1/2}(1+{\rm{Pe}_{H,b,\epsilon}})^{5/2}2^{-\ell}\mathrm{C}(\alpha)\left\|\nabla
u^{i,\mathrm{II}}\right\|_{L^{2}\left(\omega_{i}\right)}.$ (5.7)

Step 3. Estimate the global error by summation of local error. Using Theorem
5.1, we obtain

$\displaystyle\left\|{e}_{\mathrm{}}\right\|_{L^{2}\left(D\right)}\leq\sqrt{C_{\mathrm{ov}}}C_{\infty}\left(\sum_{i=1}^{N}\|{e}^{i}_{\mathrm{}}\|_{L^{2}(\omega_{i})}^{2}\right)^{1/2},$

which, together with (5.6) and the overlapping condition (2.7) leads to

$\displaystyle\left\|{e}_{\mathrm{}}\right\|_{L^{2}\left(D\right)}\leq
C_{\text{weak}}\mathrm{C}(\alpha)\sqrt{C_{\mathrm{ov}}}C_{\infty}(\rm{Pe}_{H,b,\epsilon}+1)^{2}\epsilon^{1/2}2^{-\ell}\left(\sum_{i=1}^{N}\|\nabla
u^{i,\mathrm{II}}\|_{L^{2}(\omega_{i})}^{2}\right)^{1/2}.$

Lemma 5.1 implies

$\displaystyle\left\|{e}_{\mathrm{}}\right\|_{L^{2}\left(D\right)}\leq
C_{\text{weak}}\mathrm{C}(\alpha)\sqrt{C_{\mathrm{ov}}}C_{\infty}(\rm{Pe}_{H,b,\epsilon}+1)^{2}\epsilon^{1/2}2^{-\ell}\left(\sum_{i=1}^{N}\left(\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(\omega_{i}\right)}+\|\nabla
u\|_{L^{2}(\omega_{i})}\right)^{2}\right)^{1/2}.$

Together with the overlapping condition (2.7), we derive

$\displaystyle\left\|{e}_{\mathrm{}}\right\|_{L^{2}\left(D\right)}\leq\sqrt{2}C_{\text{weak}}\mathrm{C}(\alpha)C_{\mathrm{ov}}C_{\infty}(\rm{Pe}_{H,b,\epsilon}+1)^{2}\epsilon^{1/2}2^{-\ell}\left(\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(D\right)}+\|\nabla
u\|_{L^{2}(D)}\right).$

This proves (5.3).

Next, we prove (5.4). Using again Theorem 5.1, we obtain

$\displaystyle\|\nabla{e}_{\mathrm{}}\|_{L^{2}(D)}$
$\displaystyle\leq\sqrt{2C_{\mathrm{ov}}}\left(\sum_{i=1}^{N}C_{\text{G}}^{2}H^{-2}\|{e}^{i}_{\mathrm{}}\|_{L^{2}(\omega_{i})}^{2}+\|\chi_{i}\nabla{e}^{i}_{\mathrm{}}\|_{L^{2}(\omega_{i})}^{2}\right)^{1/2}$
(5.8)

Analogously, we can derive the following upper bound for the second term by
(5.7) and Lemma 5.1,

$\displaystyle\|\chi_{i}\nabla{e}^{i}_{\mathrm{}}\|_{L^{2}(\omega_{i})}\lesssim(1+\rm{Pe}_{H,b,\epsilon})^{5/2}H^{-1}\epsilon^{1/2}2^{-\ell}\left(\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(\omega_{i}\right)}+\|\nabla
u\|_{L^{2}(\omega_{i})}\right).$

Combining with (5.8) and (5.6), we have proved (5.4), and this completes our
proof. ∎

Finally, we investigate the error between the exact solution $u$ to Problems
(1.1) and the multiscale solution $u_{{\rm ms},\ell}$, which is derived from
first establishment of the quasi-optimality of $u_{{\rm ms},\ell}$ that is
analogous to the classical Céa’s lemma, and then apply the approximation
properties of the edge multiscale ansatz space stated in Lemma 5.2.

######  Theorem 5.2 (Error estimate for Algorithm 1).

Let $\ell\in\mathbb{N}$ be non-negative. Let $u\in V$ be the solution to
Problem (1.1) and let $u_{\text{ms},\ell}$ be the multiscale solution defined
in (4.2). Then there holds

$\displaystyle\left\|\nabla(u-u_{{\rm
ms},\ell})\right\|_{L^{2}\left(D\right)}\lesssim\rm{Pe}_{H,b,\epsilon}^{3}H^{-1}\epsilon^{1/2}2^{-\ell}\left(\rm{Pe}_{H,b,\epsilon}\left\|f\right\|_{L^{2}\left(D\right)}+\|\nabla
u\|_{L^{2}(D)}\right).$ (5.9)

###### Proof.

We obtain from the global decomposition (3.5), the definition of $u_{{\rm
ms},\ell}^{\mathrm{II}}$ (2.8), and (4.2),

$\displaystyle\epsilon\left\|\nabla(u-u_{{\rm
ms},\ell})\right\|_{L^{2}\left(D\right)}^{2}$
$\displaystyle=a(u^{\mathrm{II}}-u_{{\rm
ms},\ell}^{\mathrm{II}},u^{\mathrm{II}}-u_{{\rm ms},\ell}^{\mathrm{II}}).$

Since $V_{\text{ms},\ell}\subset V$, this implies

$\displaystyle\epsilon\left\|\nabla(u-u_{{\rm
ms},\ell})\right\|_{L^{2}\left(D\right)}^{2}$
$\displaystyle=a(u^{\mathrm{II}}-u_{{\rm
ms},\ell}^{\mathrm{II}},u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}}).$

The definition of the bilinear from $a(\cdot,\cdot)$ (2.1), together with the
Cauchy-Schwarz inequality, further leads to

$\displaystyle\qquad\epsilon\left\|\nabla(u^{\mathrm{II}}-u^{\mathrm{II}}_{{\rm
ms},\ell})\right\|_{L^{2}\left(D\right)}^{2}$
$\displaystyle=\epsilon\int_{D}\nabla(u^{\mathrm{II}}-u_{{\rm
ms},\ell}^{\mathrm{II}})\cdot\nabla(u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}})\mathrm{d}x+\int_{D}b\dot{\nabla}(u^{\mathrm{II}}-u_{{\rm
ms},\ell}^{\mathrm{II}})(u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}})\mathrm{d}x$
$\displaystyle\leq{\epsilon}\left\|\nabla(u^{\mathrm{II}}-u_{{\rm
ms},\ell}^{\mathrm{II}})\right\|_{L^{2}\left(D\right)}\left\|\nabla(u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}})\right\|_{L^{2}\left(D\right)}+\|b\|_{L^{\infty}(D)}\left\|\nabla(u^{\mathrm{II}}-u_{{\rm
ms},\ell}^{\mathrm{II}})\right\|_{L^{2}\left(D\right)}\left\|u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}}\right\|_{L^{2}\left(D\right)}.$

Consequently, we derive

$\displaystyle\left\|\nabla(u-u_{{\rm
ms},\ell})\right\|_{L^{2}\left(D\right)}\leq\left\|\nabla(u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}})\right\|_{L^{2}\left(D\right)}+\mathrm{Pe}\left\|u^{\mathrm{II}}-\mathcal{P}_{\ell}u^{\mathrm{II}}\right\|_{L^{2}\left(D\right)}.$

Finally, an application of (5.3) and (5.4) yields the desired assertion. ∎

######  Remark 5.2 (Convergence with respect to $\epsilon$ in (5.9)).

The error estimate in (5.9) indicates plausibly half order convergence with
respect to the parameter $\epsilon$ when the element Péclet number
$\rm{Pe}_{H,b,\epsilon}$ is fixed. However, this is not observed from the
numerical test. One possibly explanation is that the upper bound of the term
$\|\nabla u\|_{L^{2}(D)}$ depends on $\epsilon$. When $d=1$, it is proved that
$\|\nabla u\|_{L^{2}(D)}\leq\mathcal{O}(\epsilon^{-1/2})$ [21, Proposition 8].

##  6 Numerical tests

In this section, several numerical experiments are presented to illustrate the
computational performance of Algorithm 1. In our experiments, we set the
computational domain to be $D:=[0,1]^{d}$ for $d=2$ and $3$ and the constant
force is employed, namely $f:=1$. Let $\mathcal{T}^{H}$ be a regular quasi-
uniform rectangular mesh over $D$ with maximal mesh size $H$ and let
$\mathcal{T}^{h}$ be a regular quasi-uniform rectangular mesh over each coarse
element $T\in\mathcal{T}^{H}$ with maximal mesh size $h$. To estimate the
accuracy of the numerical solution $u_{\text{ms},\ell}$ (4.2), we calculate
the relative $L^{2}(D)$-error and the relative $H^{1}(D)$-error defined by

$e_{L^{2}}:=\frac{\|u_{\text{ms},\ell}-u_{\text{ref}}\|_{L^{2}(D)}}{\|u_{\text{ref}}\|_{L^{2}(D)}},\quad
e_{H^{1}}:=\sqrt{\frac{a(u_{\text{ms},\ell}-u_{\text{ref}},u_{\text{ms},\ell}-u_{\text{ms},\ell})}{a(u_{\text{ref}},u_{\text{ref}})}},$

where $u_{\text{ref}}$ is the reference solution computed using classical
$\mathbb{Q}^{1}$ conforming Galerkin over $\mathcal{T}^{h}$ with
$h:=\sqrt{2}/2^{10}$ for $d=2$ and $h:=\sqrt{2}/2^{6}$ for $d=3$.

### Example 1

We take the perturbation parameter $\epsilon:=10^{-2}$ and the velocity field
is assumed to be

$\displaystyle b:=\alpha\left[\sin(k\pi x)\cos(k\pi y),-\cos(k\pi x)\sin(k\pi
y)\right]^{T},$

which has a cellular structure with several eddies and separatrices. We
consider two pairs of parameter,

$\displaystyle(\alpha,k):=(2,24)\text{ and }(\alpha,k):=(8,48).$

$H$ | $\ell=0$ | $\ell=1$ | $\ell=2$
---|---|---|---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
$\sqrt{2}/8$ | 0.83% | 4.44% | 0.32% | 2.03% | 0.19% | 1.06%
$\sqrt{2}/16$ | 0.26% | 4.48% | 0.07% | 1.84% | 0.12% | 0.81%
$\sqrt{2}/32$ | 0.28% | 7.19% | 0.03% | 2.09% | 0.0061% | 0.42%
$\sqrt{2}/64$ | 0.11% | 6.28% | 0.01% | 1.08% | 0.0012% | 0.20%
Table 1: Convergence history of Algorithm 1 for Example 1 with
$(\alpha,k)=(2,24)$.

The convergence history of Algorithm 1 with different coarse mesh size $H$ and
different level parameter $\ell$ is displayed in Tables 1 and 2, respectively.
In specific, we depict the reference solution and the multiscale solution
$u_{\text{ms},\ell}$ from Algorithm 1 with $H=\sqrt{2}/16$ and $\ell=1$ in
Figures 3 and 4.

We observe that both relative errors $e_{L^{2}}$ and $e_{H^{1}}$ decrease as
level parameter $\ell$ increases as expected, errors corresponding to the case
$(\alpha,k)=(8,48)$ are generally larger than $(\alpha,k)=(2,24)$ since the
former case has a larger Péclet number and thus results in a more oscillating
solution as exhibited in Figures 3 and 4. We also observe from Figures 3 and 4
that the multiscale solution $u_{\text{ms},\ell}$ with $H=\sqrt{2}/16$ and
$\ell=1$ approximates the reference solution well and recover multiscale
details.

(a) $u_{\text{ref}}$ (b) $u_{\text{ms},\ell}$

Figure 3: Comparison of reference solution $u_{\text{ref}}$ with $u_{\text{ms},\ell}$ from Algorithm 1 with $H=\sqrt{2}/16$ and $\ell=1$ for Example 1 with $(\alpha,k)=(2,24)$. $H$ | $\ell=0$ | $\ell=1$ | $\ell=2$
---|---|---|---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
$\sqrt{2}/8$ | 0.91% | 3.79% | 0.32% | 1.09% | 0.30% | 0.91%
$\sqrt{2}/16$ | 1.20% | 3.46% | 1.22% | 2.48% | 0.83% | 1.73%
$\sqrt{2}/32$ | 2.18% | 5.47% | 0.82% | 2.58% | 0.37% | 1.43%
$\sqrt{2}/64$ | 1.09% | 8.62% | 0.05% | 3.57% | 0.013% | 0.74%
Table 2: Convergence history of Algorithm 1 for Example 1 with
$(\alpha,k)=(8,48)$.

(a) Reference FEM solution (b) WEMsFEM solution

Figure 4: Comparison of reference solution $u_{\text{ref}}$ with
$u_{\text{ms},\ell}$ from Algorithm 1 with $H:=\sqrt{2}/16$ and $\ell:=1$ for
Example 1 with $(\alpha,k):=(8,48)$.

### Example 2

This example is adopted from [3]. We take $\epsilon=10^{-3}$ and define the
velocity to be

$\displaystyle b:=\left[-\frac{\partial{g}}{\partial
y},\frac{\partial{g}}{\partial x}\right]^{T}\text{ with
}g(x,y):=\frac{1}{60\pi}\sin(5\pi x)\sin(6\pi y)+\frac{1}{20}(x+y).$

This velocity field also has a cellular structure with eddies and channels as
in Example 1. The convergence of $u_{\text{ms},\ell}$ with different coarse-
grid size $H$ and level parameter $\ell$ is reported in Table 3. Convergence
with respect to the level parameter $\ell$ can be observed in this example, we
also observe convergence with respect to $H$. Surprisingly, Algorithm 1 with
$\ell=0$ and $H=\sqrt{2}/8$ can deliver an accurate solution. The reference
solution $u_{\text{ref}}$ and multiscale solution $u_{\text{ms},\ell}$ with
$H:=\sqrt{2}/16$ and $\ell=0$ are displayed in Figure 5, which confirm the
data in Table 3. It can also be seen that $u_{\text{ms},\ell}$ with
$H:=\sqrt{2}/16$ and $\ell=0$ can capture the microscale feature of the
reference solution.

$H$ | $\ell=0$ | $\ell=1$ | $\ell=2$
---|---|---|---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
$\sqrt{2}/8$ | 3.20% | 7.93% | 0.19% | 2.07% | 0.03% | 0.42%
$\sqrt{2}/16$ | 0.61% | 6.01% | 0.05% | 1.37% | 0.0054% | 0.22%
$\sqrt{2}/32$ | 0.10% | 2.91% | 0.0090% | 0.55% | 0.0011% | 0.10%
$\sqrt{2}/64$ | 0.02% | 1.14% | 0.0017% | 0.20% | 0.0002% | 0.04%
Table 3: Convergence history of Algorithm 1 for Example 2.

(a) $u_{\text{ref}}$ (b) $u_{\text{ms},\ell}$

Figure 5: Comparison of reference solution $u_{\text{ref}}$ with
$u_{\text{ms},\ell}$ from Algorithm 1 with $H:=\sqrt{2}/16$ and $\ell=1$ for
Example 2.

### Example 3

A channelized flow field without eddies is considered with $\epsilon:=1$ and
$b:=200[\sin(48\pi y),0]^{T}.$ The numerical results are presented in Table 4,
and similar convergence behavior as Example 1 is observed. We depict the
reference solution $u_{\text{ref}}$ and $u_{\text{ms},\ell}$ from Algorithm 1
with $H:=\sqrt{2}/16$, $\ell=0$ in Figure 6. The errors displayed in Table 4
verify that Algorithm 1 can provide an accurate solution. In particular, even
if the level parameter $\ell=0$, $L^{2}(D)$-relative errors and
$H^{1}(D)$-relative errors are all below $0.56\%$ and $7.10\%$, respectively.

$H$ | $\ell=0$ | $\ell=1$ | $\ell=2$
---|---|---|---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
$\sqrt{2}/8$ | 0.47% | 2.67% | 0.04% | 0.68% | 0.02% | 0.45%
$\sqrt{2}/16$ | 0.16% | 1.95% | 0.03% | 1.31% | 0.02% | 1.10%
$\sqrt{2}/32$ | 0.25% | 4.48% | 0.08% | 2.49% | 0.01% | 0.54%
$\sqrt{2}/64$ | 0.56% | 7.10% | 0.02% | 1.42% | 0.0014% | 0.27%
Table 4: Convergence history of Algorithm 1 for Example 3.

(a) $u_{\text{ref}}$ (b) $u_{\text{ms},\ell}$

Figure 6: Comparison of the reference solution $u_{\text{ref}}$ with $u_{\text{ms},\ell}$ from Algorithm 1 with $H:=\sqrt{2}/16$ and $\ell=0$ for Example 3. $H$ | $\ell=0$ | $\ell=1$ | $\ell=2$
---|---|---|---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
$\sqrt{2}/8$ | 0.74% | 4.66% | 0.14% | 1.43% | 0.04% | 0.52%
$\sqrt{2}/16$ | 0.41% | 4.73% | 0.11% | 1.55% | 0.02% | 0.38%
$\sqrt{2}/32$ | 0.18% | 3.49% | 0.04% | 0.97% | 0.0047% | 0.20%
$\sqrt{2}/64$ | 0.06% | 2.06% | 0.0085% | 0.41% | 0.0010% | 0.08%
Table 5: Convergence history of Algorithm 1 for Example 4 with
$\epsilon_{1}=1$.

### Example 4

This example is adopted from [21]. In particular, we consider a diffusion
coefficient with multiple scales, such that we aim to solve

$\displaystyle-\nabla\cdot(A(x,y)\nabla u)+b\cdot\nabla u=f$

with the diffusion coefficient $A(x,y)$ being

$\displaystyle
A(x,y)=\epsilon\bigg{(}1+0.5\cos\bigg{(}\frac{2\pi}{\epsilon_{1}}x\bigg{)}\bigg{)}.$

Here, the oscillation parameter $\epsilon:=\frac{1}{128}$ and we take a
constant velocity $b:=(1,1)^{T}$. To guarantee the accuracy of reference
solution $u_{\text{ref}}$, the fine scale mesh size $h$ is set to be
$\sqrt{2}/2^{12}$.

(a) $u_{\text{ref}}$ (b) $u_{\text{ms},\ell}$

Figure 7: Comparison of reference solution $u_{\text{ref}}$ with $u_{\text{ms},\ell}$ from Algorithm 1 with $H:=\sqrt{2}/16$ and $\ell=1$ for Example 4 with $\epsilon_{1}=1$. $H$ | $\ell=0$ | $\ell=1$ | $\ell=2$
---|---|---|---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
| $e_{L^{2}}$
---
| $e_{H^{1}}$
---
$\sqrt{2}/8$ | 0.85% | 5.33% | 0.18% | 1.82% | 0.05% | 0.77%
$\sqrt{2}/16$ | 0.54% | 6.83% | 0.13% | 2.08% | 0.02% | 0.52%
$\sqrt{2}/32$ | 0.40% | 8.27% | 0.08% | 1.97% | 0.0075% | 0.36%
$\sqrt{2}/64$ | 0.21% | 6.51% | 0.02% | 1.37% | 0.0021% | 0.16%
Table 6: Convergence history of Algorithm 1 for Example 4 with
$\epsilon_{1}=1/64$.

(a) $u_{\text{ref}}$ (b) $u_{\text{ms},\ell}$

Figure 8: Comparison of reference solution $u_{\text{ref}}$ with
$u_{\text{ms},\ell}$ from Algorithm 1 with $H:=\sqrt{2}/16$ and $\ell=1$ for
Example 4 with $\epsilon_{1}=1/64$.

If $\epsilon_{1}=1/64$ and $H=1/16$, then $L^{2}(D)$-relative errors for MsFEM
and Adv-MsFEM proposed in [21] are 27% and 11%, and $H^{1}(D)$-relative errors
are 113% and 74%, which are not reliable, especially near the boundary layer.
To make our comparison relevant, we implement these two proposed methods and
obtain that $L^{2}(D)$-relative errors are 26% and 11%, $H^{1}(D)$-relative
errors are 108% and 73%, which agree with the data reported in [21].
Remarkably, Tables 5 and 6 clearly demonstrate that Algorithm 1 yields more
accuracy even with the level parameter $\ell=0$, both near boundary layer and
away from the boundary layer. Furthermore, we depict in Figures 7 and 8 the
reference solution $u_{\text{ref}}$ and multiscale solution
$u_{\text{ms},\ell}$ with $H:=\sqrt{2}/16$ and $\ell=0$. We observe that the
multiscale solution $u_{\text{ms},\ell}$ with $H:=\sqrt{2}/16$ and $\ell=0$
can capture the boundary layer accurately.

### Example 5

Example 5 is a 3-d example with $\epsilon:=1/16$, and the velocity field being

$b:=[g_{2}-g_{3},g_{3}-g_{1},g_{1}-g_{2}]^{T}\text{ with }$ $\displaystyle
g_{1}:=$ $\displaystyle\alpha k\cos(kx)\sin(ky)\sin(kz),$ $\displaystyle
g_{2}:=$ $\displaystyle\alpha k\sin(kx)\cos(ky)\sin(kz),$ $\displaystyle
g_{3}:=$ $\displaystyle\alpha k\sin(kx)\sin(ky)\cos(kz).$

Different combinations of $(\alpha,k)$ are tested and simulation results are
reported in Table 7. We only test Algorithm 1 with $\ell=0$ due to the huge
computational complexity involved for $d=3$.

$(\alpha,k)$ $H=\sqrt{2}/{8}$ $H=\sqrt{2}/{16}$ $e_{L^{2}}$ $e_{H^{1}}$
$T_{\text{solve}}$ $e_{L^{2}}$ $e_{H^{1}}$ $T_{\text{solve}}$ (1,4$\pi$) 0.92%
7.2% 0.8 0.29% 3.45% 21.5 (1,8$\pi$) 1.38% 6.05% 0.8 0.92% 11.06% 22.3
(4,4$\pi$) 3.72% 21.36% 0.8 1.55% 13.69% 22.1 (4,8$\pi$) 2.11% 7.33% 0.8 3.71%
31.66% 22.5

Table 7: Convergence history of Algorithm 1 and CPU time in seconds for
solving linear system in Example 5. Figure 9: Comparison of reference solution
$u_{\text{ref}}$ with $u_{\text{ms},\ell}$ from Algorithm 1 with $\ell=0$ for
Example 5 with $(\alpha,k)=(1,4\pi)$.

We can see the performance of Algorithm 1 can provide an accurate solution in
most cases. For example, the $L^{2}(D)$-relative error is merely $3.17\%$ when
$H=\sqrt{2}/{8}$ and $(\alpha,k)=(4,8\pi)$. Moreover, Figures 9 and 10
demonstrate that Algorithm 1 with $\ell=0$ can capture the microscale feature
of the reference solution. We note that the CPU time for applying a direct
solver for solving the linear system resulting from FEM is about 147 seconds
regardless of $\alpha$ and $k$. In comparison, the CPU time for Algorithm 1
are about 0.8 seconds and 22 seconds for $H=\sqrt{2}/{8}$ and
$H=\sqrt{2}/{16}$ respectively, and hence huge computational cost is saved.

Figure 10: Comparison of reference solution $u_{\text{ref}}$ with
$u_{\text{ms},\ell}$ from Algorithm 1 with $\ell=0$ for Example 5 with
$(\alpha,k)=(4,8\pi)$.

### Example 6

As the last example, we test Algorithm 1 for the challenging case with
$\epsilon\to 0$. In specific, we set $\epsilon:=10^{-9}$, $b:=[2,3]^{T}$, and
set $f$ such as the true solution is,

$u(x,y)=2\sin
x\bigg{(}1-\exp\bigg{(}-\frac{2(1-x)}{\epsilon}\bigg{)}\bigg{)}y^{2}\bigg{(}1-\exp\bigg{(}-\frac{2(1-y)}{\epsilon}\bigg{)}\bigg{)}.$

$D_{11}$$D_{12}$$D_{21}$$D_{22}$ Figure 11: Dissection of $D$ and Shishkin
mesh.

Note that the true solution $u$ has exponential boundary layer near the upper
boundary and right boundary. Since the Péclet number $\mathrm{Pe}$ in this
case is too large, we have to adopt non-uniform mesh for both the coarse mesh
$\mathcal{T}^{H}$ and the fine mesh $\mathcal{T}^{h}$. For simplicity, we
utilize the Shishkin meshes [12, 24], where the domain $D$ is dissected into
four parts as depicted in Figure 11 such that $\bar{D}=D_{11}\cup D_{21}\cup
D_{12}\cup D_{22}$ with

$\displaystyle D_{11}=$
$\displaystyle[0,1-\lambda_{x}]\times[0,1-\lambda_{y}],\quad
D_{21}=[1-\lambda_{x},1]\times[0,1-\lambda_{y}],$ $\displaystyle D_{12}=$
$\displaystyle[0,1-\lambda_{x}]\times[1-\lambda_{y},1],\quad
D_{22}=[1-\lambda_{x},1]\times[1-\lambda_{y},1].$

Here, $\lambda_{x}=\lambda_{y}=2.5\epsilon\log(1024)$.

Each subdomain $D_{ij}$ is further divided into a mesh with $512\times 512$
rectangles so the number of elements in $\mathcal{T}^{h}$ is $1024\times
1024$. The coarse mesh $\mathcal{T}^{H}$ is also a Shishkin type of mesh
composed of $M\times M$ large rectangles for some $M\in\mathbb{N}_{+}$.

Figure 12: Convergence history (against exact solution) of Algorithm 1 with
$\ell=0$ with respect to the number of coarse elements $M$ in each axis for
Example 6.

(a) $u_{\text{ref}}$ (b) $u_{\text{ms},\ell}$

Figure 13: Comparison of the true solution $u_{\text{ref}}$ with
$u_{\text{ms},\ell}$ from Algorithm 1 with $\ell=0$ and $M=64$ for Example 6.

In our implementation, for neighborhoods $\omega_{i}\subset D_{11}$, only
partition of unity function is utilized, for other neighborhoods, level zero
wavelet-based multiscale basis functions are computed and incorporated into
the coarse space. Figure 12 shows relative errors of $u_{\text{ms}}$ with
respect to the exact solution for different $M$, convergence with respect to
number of coarse element $M$ is observed. Figure 13 compares the exact
solution and multiscale solution which further verifies high accuracy of
coarse-grid simulation. We note for this experiment, it just takes about 2.5
seconds to solve the linear system resulting from multiscale simulations
$(M=128)$, if one utilize FEM with $1024\times 1024$ fine elements, CPU time
for solving corresponding system with a direct solver in Matlab is 119
seconds, therefore Algorithm 1 is a promising candidate for solving singularly
perturbed convection-diffusion equation.

###  6.1 Sensitivity with respect to the mesh Péclet number
$\rm{Pe}_{H,b,\epsilon}$

We set the level parameter $\ell:=0$ and coarse mesh size $H:=\sqrt{2}/16$ in
Algorithm 1 and test it with different values of oscillation parameter
$\epsilon$ for Examples 1-4 to study its robustness with respect to the mesh
Péclet number $\rm{Pe}_{H,b,\epsilon}$. Experiments results are provided in
Figure 14. It is observed relative errors generally increase as mesh Péclet
number $\rm{Pe}_{H,b,\epsilon}$ becomes larger. However, for examples 2-4,
relative errors are quite stable for $\rm{Pe}_{H,b,\epsilon}=\mathcal{O}(10)$.
For example 1, although accuracy of Algorithm 1 have a relative stronger
dependence on the mesh Péclet number $\rm{Pe}_{H,b,\epsilon}$, the relative
errors are still very small even if $\rm{Pe}_{H,b,\epsilon}$ is as large as
250 despite $\ell=0$.

Figure 14: Performance of Algorithm 1 with $\ell=0$ and $H=\sqrt{2}/16$
against the mesh Péclet number $\rm{Pe}_{H,b,\epsilon}$. Here,
$(\alpha,k)=(8,48)$ in Example 1 and $\epsilon_{1}:=1/64$ in Example 4.

##  7 Conclusions

We develop a novel Wavelet-based Edge Multiscale Finite Element Method
(WEMsFEM) to solve singularly perturbed convection - diffusion problems. The
key feature of this approach is to introduce an edge multiscale ansatz space
which has guaranteed approximation properties without further regularity
requirement on the solution. The construction of this edge multiscale ansatz
space is motivated by the local and global splitting of the solution, and
their basis functions are solved locally using hierarchical bases as the
Dirichlet data. The convergence of WEMsFEM is rigorously justified, which is
confirmed numerically for a wide range of 2-d or 3-d convection-dominated
diffusion problems with high oscillation and boundary layers.

## Acknowledgements

GL acknowledges the support from Newton International Fellowships Alumni
following-on funding awarded by The Royal Society, Young Scientists fund
(Project number: 12101520) by NSFC and Early Career Scheme (Project number:
27301921), RGC, Hong Kong. SF’s research is supported by startup funding of
Eastern institute of technology, Ningbo and NSFC (Project number: 12301514).
The research of EC is partially supported by the Hong Kong RGC General
Research Fund (Project numbers: 14305222 and 14304021).

##  Appendix A Very-weak solutions to singularly perturbed convection-
diffusion problems

In this appendix, we derive an $L^{2}(\omega_{i})$-estimate for boundary
valued singularly perturbed convection-diffusion problem, which plays a
crucial role in the error analysis. Let $\omega_{i}$ be a coarse neighborhood
for any $i=1,\cdots,N$. For any $g\in L^{2}(\partial\omega_{i})$, we define
the following convection diffusion problem

$\left\\{\begin{aligned} \mathcal{L}_{i}v:=-\epsilon\Delta v+b\cdot\nabla
v&=0&&\text{ in }\omega_{i},\\\ v&=g&&\text{ on
}\partial\omega_{i}.\end{aligned}\right.$ (A.1)

Our goal is to derive the $L^{2}(\omega_{i})$-estimate of the solution $v$. To
this end, we employ a nonstandard variational form in the spirit of the
transposition method [25]. Denote the dual operator of $\mathcal{L}_{i}$ as
$\mathcal{L}_{i}^{*}$. Since $\nabla\cdot b=0$, we obtain
$\mathcal{L}_{i}^{*}=-\epsilon\Delta-b\cdot\nabla$. Let the test space
$X(\omega_{i})\subset H^{1}_{0}(\omega_{i})$ be defined by

$\displaystyle X(\omega_{i}):=\\{z:\mathcal{L}_{i}^{*}z\in
L^{2}(\omega_{i})\text{ and }z\in H^{1}_{0}(\omega_{i})\\}.$ (A.2)

This test space $X(\omega_{i})$ is endowed with the norm
$\|\cdot\|_{X(\omega_{i})}$:

$\forall z\in X(\omega_{i}):\|z\|_{X(\omega_{i})}^{2}=\int_{\omega_{i}}|\nabla
z|^{2}\,\mathrm{d}x+\int_{\omega_{i}}|\mathcal{L}_{i}^{*}z|^{2}\,\mathrm{d}x.$

Then we propose the following weak formulation corresponding to Problem (A.1):
seeking $v\in L^{2}(\omega_{i})$ such that

$\displaystyle\int_{\omega_{i}}v\mathcal{L}_{i}^{*}(z)\,\mathrm{d}x=-\epsilon\int_{\partial\omega_{i}}g\frac{\partial
z}{\partial n}\mathrm{d}s\quad\text{ for all }z\in X(\omega_{i}).$ (A.3)

The remaining of this section is devoted to proving the well posedness of the
nonstandard variational formulation (A.3) and deriving the a priori error
estimate. Our final result is the following,

######  Theorem A.1.

Given $g\in L^{2}(\partial\omega_{i})$. Let $v$ be the solution to (A.1), then
there exists a constant ${\rm C}_{{\rm weak}}$ independent of the parameter
$\epsilon$ such that

$\left\|v\right\|_{L^{2}\left(\omega_{i}\right)}\leq{\rm C}_{{\rm
weak}}\epsilon^{1/2}(1+{\rm{Pe}_{H,b,\epsilon}})\|g\|_{L^{2}(\partial\omega_{i})}.$

To prove Theorem A.1, one has to first derive the
$L^{2}(\partial\omega_{i})$-estimate of the normal trace $\frac{\partial
z}{\partial n}$ for any $z\in X(\omega_{i})$. This is established in the
following theorem:

######  Theorem A.2.

Let $w\in L^{2}(\omega_{i})$ and let $z\in X(\omega_{i})$ satisfy

$\left\\{\begin{aligned} \mathcal{L}_{i}^{*}(z):=-\epsilon\Delta
z-b\cdot\nabla z&=w&&\text{ in }\omega_{i},\\\ z&=0&&\text{ on
}\partial\omega_{i}.\end{aligned}\right.$ (A.4)

Then for some constant ${\rm C}_{{\rm weak}}$ independent of the parameter
$\epsilon$ and $H$, there holds

$\displaystyle\|\frac{\partial z}{\partial n}\|_{L^{2}(\partial\omega_{i})}$
$\displaystyle\leq{\rm C}_{{\rm
weak}}\epsilon^{-1/2}(1+{\rm{Pe}_{H,b,\epsilon}})\left\|w\right\|_{L^{2}\left(\omega_{i}\right)}.$

This constant ${\rm C}_{{\rm weak}}$ can change value from context to context.

###### Proof.

Note that $z\in X(\omega_{i})$ is the unique solution to the following weak
formulation

$\displaystyle\forall q\in
H^{1}_{0}(\omega_{i}):\int_{\omega_{i}}\mathcal{L}_{i}^{*}(z)q\;\,\mathrm{d}x=\int_{\omega_{i}}wq\,\mathrm{d}x.$

Taking $q:=z$ and applying integration by parts, we arrive at

$\displaystyle\epsilon\left\|\nabla
z\right\|_{L^{2}\left(\omega_{i}\right)}^{2}=\int_{\omega_{i}}wz\,\mathrm{d}x.$

Since $z|_{\partial\omega_{i}}=0$, then an application of the Poincaré
inequality leads to

$\displaystyle\left\|\nabla
z\right\|_{L^{2}\left(\omega_{i}\right)}\lesssim\rm{Pe}_{H,b,\epsilon}\left\|w\right\|_{L^{2}\left(\omega_{i}\right)}.$
(A.5)

Recall that the mesh Péclet number $\rm{Pe}_{H,b,\epsilon}$ of the coarse mesh
$\mathcal{T}^{H}$ is
$\rm{Pe}_{H,b,\epsilon}:=H\|b\|_{L^{\infty}(D)}/\epsilon$. Consequently,
plugging this estimate into (A.4) results in

$\displaystyle\left\|\Delta
z\right\|_{L^{2}\left(\omega_{i}\right)}\lesssim\epsilon^{-1}(1+\rm{Pe}_{H,b,\epsilon})\left\|w\right\|_{L^{2}\left(\omega_{i}\right)}.$
(A.6)

By [20, Page 20], we can obtain the following a priori estimate

$\displaystyle|{z}|_{H^{2}(\omega_{i})}$ $\displaystyle\lesssim\left\|\Delta
z\right\|_{L^{2}\left(\omega_{i}\right)}\lesssim\epsilon^{-1}(1+\rm{Pe}_{H,b,\epsilon})\left\|w\right\|_{L^{2}\left(\omega_{i}\right)}.$

Combining the preceding two estimates and applying interpolation between
$H^{1}(\omega_{i})$ and $H^{2}(\omega_{i})$ yield the $H^{3/2}(\omega_{i})$
regularity estimate

$\displaystyle\left\|z\right\|_{H^{3/2}\left(\omega_{i}\right)}\lesssim\epsilon^{-1/2}(1+\rm{Pe}_{H,b,\epsilon})\left\|w\right\|_{L^{2}\left(\omega_{i}\right)}.$
(A.7)

Since differentiation is continuous from $H^{3/2}(\omega_{i})$ to
$H^{1/2}(\omega_{i})$, by the trace theorem, we have

$\displaystyle\left\|\frac{\partial z}{\partial
n}\right\|_{L^{2}(\partial\omega_{i})}$ $\displaystyle\lesssim\|\nabla
z\|_{H^{1/2}(\omega_{i})}\lesssim\left\|z\right\|_{H^{3/2}\left(\omega_{i}\right)},$

which, together with (A.7), proves the desired assertion. ∎

Thanks to this a priori estimate presented in Theorem A.2, we are ready to
state the well posedness of the nonstandard variational formulation, cf.
(A.3),

######  Lemma A.1.

Let $v$ be the solution to problem (A.1) and let the test space
$X(\omega_{i})$ be defined in (A.2). Then the nonstandard variational form
(A.3) is well posed.

###### Proof.

To prove the well-posedness of the nonstandard variational form (A.3), we
introduce a bilinear form $c(\cdot,\cdot)$ on $L^{2}(\omega_{i})\times
L^{2}(\omega_{i})$ and a linear form $b(\cdot)$ on $L^{2}(\omega_{i})$,
defined by

$\displaystyle c(w_{1},w_{2})$
$\displaystyle:=\int_{\omega_{i}}w_{1}w_{2}\;\,\mathrm{d}x\quad\text{ for all
}w_{1},w_{2}\in L^{2}(\omega_{i})$ $\displaystyle r(w)$
$\displaystyle:=-\epsilon\int_{\partial\omega_{i}}\frac{\partial z}{\partial
n}w\;\mathrm{d}s\quad\text{ for all }w\in L^{2}(\omega_{i}),$

with $z$ being the unique solution to (A.4). It follows from Theorem A.2 that

$\displaystyle\|r\|:=\sup_{w\in
L^{2}(\omega_{i})}\frac{r(w)}{\|w\|_{L^{2}(\omega_{i})}}\leq{\rm C}_{{\rm
weak}}\epsilon^{1/2}(1+{\rm{Pe}_{H,b,\epsilon}})\left\|g\right\|_{L^{2}\left(\partial\omega_{i}\right)}.$
(A.8)

This yields well-posedness of the following variational problem: find $v\in
L^{2}(\omega_{i})$ such that

$\displaystyle c(v,w)=r(w)\quad\text{ for all }w\in L^{2}(\omega_{i}).$ (A.9)

The equivalence of problems (A.9) and (A.3) implies the desired well-posedness
of (A.3). ∎

Finally, we are ready to prove the main result stated in Theorem A.1:

###### Proof of Theorem A.1.

For all $w\in L^{2}(\omega_{i})$, we obtain from (A.8) and (A.9):

$\displaystyle\int_{\omega_{i}}vw\;\,\mathrm{d}x:=c(v,w)=r(w)\leq\text{C}_{\text{weak}}\epsilon^{1/2}(1+{\rm{Pe}_{H,b,\epsilon}})\left\|w\right\|_{L^{2}\left(\omega_{i}\right)}\|g\|_{\partial\omega_{i}}.$

Then this completes the proof by taking $w:=v$. ∎

Furthermore, we present the a priori estimate in energy norm to facilitate the
proofs in Section 5. This estimate is standard and is commonly known as the
Cacciappoli inquality, which can be found in many literature, e.g., [1, 22].

######  Corollary A.1.

Let $g\in L^{2}(\partial\omega_{i})$ and let $v$ be the solution to (A.1).
Then there exists a constant ${\rm C}_{{\rm weak}}$ independent of the
parameter $\epsilon$ such that

$\left\|\chi_{i}\nabla v\right\|_{L^{2}\left(\omega_{i}\right)}\leq{\rm
C}_{{\rm
weak}}H^{-1}\epsilon^{1/2}(1+{\rm{Pe}_{H,b,\epsilon}})^{3/2}\|g\|_{L^{2}(\partial\omega_{i})}.$

###### Proof.

Recall that $\\{\chi_{i}\\}_{i=1}^{N}$ is a
$(C_{\mathrm{ov}},C_{\infty},C_{\text{G}})$ partition of unity subordinate to
the cover $\\{\omega_{i}\\}_{i=1}^{N}$, which implies that $\chi_{i}=0$ on
$\partial\omega_{i}$. Multiplying (A.1) by $\chi_{i}^{2}v$ and applying
integration by parts, we arrive at

$\displaystyle\epsilon\int_{\omega_{i}}\chi_{i}^{2}|\nabla
v|^{2}\,\mathrm{d}x=-2\epsilon\int_{\omega_{i}}\nabla
v\cdot\nabla\chi_{i}\chi_{i}v\,\mathrm{d}x+\int_{\omega_{i}}b\cdot\nabla\chi_{i}\chi_{i}v^{2}\,\mathrm{d}x.$

Then an application of the Young’s inequality reveals the assertion. ∎

## References

  * [1] I. Babuska and R. Lipton.  Optimal local approximation spaces for generalized finite element methods with application to multiscale problems.  Multiscale Model. Simul., 9(1):373–406, 2011.
  * [2] L. Berlyand and H. Owhadi.  Flux norm approach to finite dimensional homogenization approximations with non-separated scales and high contrast.  Arch. Ration. Mech. Anal., 198(2):677–721, 2010.
  * [3] V. Calo, E. Chung, Y. Efendiev, and W. T. Leung.  Multiscale stabilization for convection-dominated diffusion in heterogeneous media.  Comput. Methods Appl. Mech. Engrg., 304:359–377, 2016.
  * [4] Y. Chen, T. Hou, and Y. Wang.  Exponential convergence for multiscale linear elliptic PDEs via adaptive edge basis functions.  Multiscale Model. Simul., 19(2):980–1010, 2021.
  * [5] Y. Chen, T. Hou, and Y. Wang.  Exponentially convergent multiscale method.  Commun. Appl. Math. Comput., 2023.
  * [6] Y. Chen, T. Hou, and Y. Wang.  Exponentially convergent multiscale methods for 2D high frequency heterogeneous Helmholtz equations.  Multiscale Model. Simul., 21(3):849–883, 2023.
  * [7] E. Chung, Y. Efendiev, and T. Hou.  Multiscale Model Reduction: Multiscale Finite Element Methods and Their Generalizations.  Springer.
  * [8] E. Chung, Y. Efendiev, and W. T. Leung.  Multiscale stabilization for convection–diffusion equations with heterogeneous velocity and diffusion coefficients.  Comput. Math. Appl., 79(8):2336–2349, 2020.
  * [9] I. Daubechies.  Ten Lectures on Wavelets, volume 61 of CBMS-NSF Regional Conference Series in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1992.
  * [10] W. E and B. Engquist.  The heterogeneous multiscale methods.  Commun. Math. Sci., 1(1):87–132, 2003.
  * [11] Y. Efendiev, J. Galvis, and T. Hou.  Generalized multiscale finite element methods.  J. Comput. Phys., 251:116–135, 2013.
  * [12] P. Farrell, A. Hegarty, J. M. Miller, E. O’Riordan, and G. I. Shishkin.  Robust computational techniques for boundary layers.  CRC Press, 2000.
  * [13] S. Fu, E. Chung, and G. Li.  Edge multiscale methods for elliptic problems with heterogeneous coefficients.  J. Comput. Phys., 396:228–242, 2019.
  * [14] S. Fu, E. Chung, and G. Li.  An edge multiscale interior penalty discontinuous Galerkin method for heterogeneous Helmholtz problems with large varying wavenumber.  J. Comput. Phys., 441:Paper No. 110387, 27, 2021.
  * [15] S. Fu, G. Li, R. Craster, and S. Guenneau.  Wavelet-based edge multiscale finite element method for Helmholtz problems in perforated domains.  Multiscale Model. Simul., 19(4):1684–1709, 2021.
  * [16] D. Gilbarg and N. Trudinger.  Elliptic Partial Differential Equations of Second Order.  Classics in Mathematics. Springer Berlin Heidelberg, 2001.
  * [17] M. Griebel and M. Schweitzer.  A particle-partition of unity method for the solution of elliptic, parabolic and hyperbolic PDE.  SIAM J. Sci. Comput., 22(3):853–890, 2000.
  * [18] T. Hou and X.-H. Wu.  A multiscale finite element method for elliptic problems in composite materials and porous media.  J. Comput. Phys., 134(1):169–189, 1997.
  * [19] T. Hughes, G. Feijóo, L. Mazzei, and J.-B. Quincy.  The variational multiscale method—a paradigm for computational mechanics.  Comput. Methods Appl. Mech. Engrg., 166(1-2):3–24, 1998.
  * [20] O. Ladyženskaja, V. Solonnikov, and N. Ural~ceva.  Linear and quasilinear equations of parabolic type.  Translated from the Russian by S. Smith. Translations of Mathematical Monographs, Vol. 23. American Mathematical Society, Providence, R.I., 1968.
  * [21] C. Le Bris, F. Legoll, and F. Madiot.  A numerical comparison of some multiscale finite element approaches for advection-dominated problems in heterogeneous media.  ESAIM: Mathematical Modelling and Numerical Analysis, 51(3):851–888, 2017.
  * [22] G. Li.  On the convergence rates of GMsFEMs for heterogeneous elliptic problems without oversampling techniques.  Multiscale Model. Simul., 17(2):593–619, 2019.
  * [23] G. Li, D. Peterseim, and M. Schedensack.  Error analysis of a variational multiscale stabilization for convection-dominated diffusion equations in two dimensions.  IMA J. Numer. Anal., 38(3):1229–1253, 2018.
  * [24] T. Linß and M. Stynes.  Numerical methods on shishkin meshes for linear convection–diffusion problems.  Comput. Methods Appl. Mech. Engrg., 190(28):3527–3542, 2001.
  * [25] J.-L. Lions and E. Magenes.  Non-homogeneous Boundary Value Problems and Applications. Vol. I.  Springer-Verlag, New York-Heidelberg, 1972.
  * [26] A. Målqvist and D. Peterseim.  Localization of elliptic multiscale problems.  Math. Comp., 83(290):2583–2603, 2014.
  * [27] J. Melenk and I. Babuška.  The partition of unity finite element method: basic theory and applications.  Comput. Methods Appl. Mech. Engrg., 139(1-4):289–314, 1996.
  * [28] P. Park and T. Hou.  Multiscale numerical methods for singularly perturbed convection-diffusion equations.  International Journal of Computational Methods, 1(01):17–65, 2004\.
  * [29] H. Yserentant.  On the multilevel splitting of finite element spaces.  Numer. Math., 49(4):379–412, 1986.