EXPONENTIALLY CONVERGENT MULTISCALE FINITE ELEMENT METHOD
YIFAN CHEN, THOMAS Y. HOU, YIXUAN WANG
Dedicated to Professor Stanley Osher’s 80th birthday with admiration and friendship
Abstract.
We provide a concise review of the exponentially convergent multiscale finite
element method (ExpMsFEM) for efficient model reduction of PDEs in heterogeneous
media without scale separation and in high-frequency wave propagation. ExpMsFEM
is built on the non-overlapped domain decomposition in the classical MsFEM while
enriching the approximation space systematically to achieve a nearly exponential con-
vergence rate regarding the number of basis functions. Unlike most generalizations of
MsFEM in the literature, ExpMsFEM does not rely on any partition of unity functions.
In general, it is necessary to use function representations dependent on the right-
hand side to break the algebraic Kolmogorov
푛
-width barrier to achieve exponential
convergence. Indeed, there are online and offline parts in the function representation
provided by ExpMsFEM. The online part depends on the right-hand side locally and
can be computed in parallel efficiently. The offline part contains basis functions that are
used in the Galerkin method to assemble the stiffness matrix; they are all independent
of the right-hand side, so the stiffness matrix can be used repeatedly in multi-query
scenarios.
1.
Introduction
Multiscale methods provide an efficient way to solve challenging PDEs. A few local
basis functions adapted to the problem are constructed offline to provide an effective
model reduction of the equation. One can then use the reduced model to compute the
solution online, possibly with different right-hand sides and in a way much faster than
solving the original equation. This property is beneficial in multi-query scenarios such
as optimal design and inverse problems. Moreover, multiscale methods are inevitable
for challenging problems in rough media and high-frequency wave propagation since
standard numerical methods suffer from a vast number of degrees of freedom. See
examples of the failure of finite element methods (FEMs) in elliptic equations with
rough coefficients [4] and the pollution effect in the Helmholtz equation [6].
In this paper, we present the framework of ExpMsFEM, the exponentially convergent
multiscale finite element method. It is a generalization of the classical MsFEM [22].
The main contribution of ExpMsFEM is the systematic improvement over MsFEM to
achieve exponentially convergent accuracy regarding the number of basis functions.
Also, unlike most generalizations of MsFEM in the literature, ExpMsFEM does not rely
on the partition of unity functions to connect local and global approximation spaces.
Instead, ExpMsFEM uses edge localization and coupling intrinsic to the non-overlapped
domain decomposition to communicate the local and global approximations.
In the literature, exponentially convergent multiscale methods have been pioneered
in the work of optimal basis [2] based on the partition of unity functions; see also the
developments in [45, 8, 9, 3, 44, 30, 31]. The work demonstrates the importance of
Caccioppoli’s inequality in establishing exponential convergence; more precisely, the
inequality implies the
low approximation complexity
of the restriction operator acting on
harmonic-type functions. The theory of ExpMsFEM is also based on some arguments
Date
: December 2022.
2010
Mathematics Subject Classification.
65N12, 65N15, 65N30, 31A35.
1
arXiv:2212.00823v2 [math.NA] 6 Feb 2023
2
YIFAN CHEN, THOMAS Y. HOU, YIXUAN WANG
using Caccioppoli’s inequality. Additionally, since no partition of unity functions is
used, technical tools such as
퐶
훼
estimates and trace theorems are needed to analyze
ExpMsFEM. We will comment on the similarity and differences between the optimal
basis work and ExpMsFEM at the end of the article.
This review is based on our previous work on exponentially convergent multiscale
methods for elliptic equations [11] and Helmholtz equations [12]. We focus on artic-
ulating the main ideas and the computational framework in the case of 2D stationary
problems with homogeneous boundary data. We provide references for the detailed
analysis in corresponding papers.
Organization.
In Section 2, we present the model problem that is the focus of this
article. In Section 3, we present the motivation and framework of the ExpMsFEM.
We provide numerical experiments to demonstrate the effectiveness of the ExpMsFEM
framework in Section 4. In Section 5, we discuss related literature, future possibilities,
and open questions.
2.
Model Problem
Consider the model problem in a bounded domain
Ω
⊂
R
푑
with a Lipschitz boundary
Γ
. Here,
푑
=
2
. For generality, the boundary can contain disjoint parts
Γ = Γ
1
∪
Γ
2
where
Γ
1
corresponds to the Dirichlet boundary conditions and
Γ
2
corresponds to the
Neumann and Robin type boundary conditions. The model equation is:
(2.1)
−∇·(
퐴
∇
푢
)+
푉푢
=
푓 ,
in
Ω
푢
=
0
,
on
Γ
1
퐴
∇
푢
·
휈
=
훽
푢,
on
Γ
2
.
Here,
퐴,푉,
훽
are functions in
퐿
∞
(
Ω
)
and can be rough, which makes the solution
oscillating and difficult to solve. The vector
휈
is the outer normal to the boundary.
In particular, when
푉
=
0
, the equation is the standard elliptic equation [11]. If
푉푢
=
−
푘
2
푢
and
푢
is a complex-valued function, one obtains the Helmholtz equation
[12] with wavenumber
푘
.
The weak formulation of (2.1) is given by
(2.2)
푎
(
푢,푣
)
:
=
(
퐴
∇
푢,
∇
푣
)
Ω
+(
푉푢,푣
)
Ω
−(
훽
푢,푣
)
Γ
2
=
(
푓 ,푣
)
Ω
,
∀
푣
∈ℋ(
Ω
)
,
where
(·
,
·)
푋
is the standard
퐿
2
inner product on the set
푋
. The space for
푣
is
ℋ(
Ω
)
:
=
{
푤
∈
퐻
1
(
Ω
)
:
푤
|
Γ
1
=
0
}
and the solution
푢
∈ ℋ(
Ω
)
. The energy norm
‖ · ‖
ℋ(
Ω
)
is
defined as
‖
푤
‖
2
ℋ(
Ω
)
:
=
(
퐴
∇
푤,
∇
푤
)
Ω
+|(
푉푤,푤
)
Ω
|
.
Here, we adopt an abuse of notation that the space can be real-valued or complex-
valued, depending on the context.
A generic assumption for
퐴
is
0
<
퐴
min
≤
퐴
(
푥
) ≤
퐴
max
<
∞
. We will present
more detailed assumptions on
푉,
훽
later in specific problems that our theory in [11, 12]
covers. Indeed, the theory can encompass the case for very general
푉
, provided
that
|(
푉푢,푢
)|
Ω
≤
푉
0
(
푢,푢
)
Ω
for some constant
푉
0
and the PDE satisfies good stability
estimates; see for example the rough Helmholtz example in [12]. In this review, we
mainly focus on the
conceptual algorithmic framework
of solving the equation (2.1) via
ExpMsFEM rather than a detailed analysis of the equation and the method.
3.
The ExpMsFEM Framework
In subsection 3.1, we discuss the general recipe for solving PDEs as a function
approximation problem. This motivates us to find accurate function representations
EXPONENTIALLY CONVERGENT MULTISCALE FINITE ELEMENT METHOD
3
to be used in the Galerkin method. We explain how ExpMsFEM manages to get
exponentially convergent representations in subsections 3.2, 3.3, 3.4 and 3.5.
3.1.
Solving PDEs as function approximation.
By the standard finite element theory
(e.g., [7]), when using the Galerkin method to solve (2.2), a key step is to find a func-
tion representation, or a space of basis functions that can approximate the solution
accurately. More precisely, suppose the space is
푆
, then, one usually wants
(3.1)
휂
(
푆
)
:
=
sup
푓
∈
퐿
2
(
Ω
)\{
0
}
inf
푣
∈
푆
푁
(
푓
)−
푣
ℋ(
Ω
)
‖
푓
‖
퐿
2
(
Ω
)
to be small. Here,
푁
:
푓
→
푢
is the solution operator
1
of (2.1).
For example, consider the elliptic equation with
푉
=
0
and
Γ
2
=
∅
. In such case, the
Galerkin method provides an optimal approximation of the solution in the space of basis
functions with respect to the energy norm [7, 11], due to the Galerkin orthogonality.
Therefore, a small
휂
(
푆
)
directly implies a small error in the solution. For the Helmholtz
equation, similar arguments hold based on the Gårding-type inequality, which leads
to the quasi-optimality of the solution; see, for example, [35, 12]. The failure of many
finite element methods in elliptic equations with rough coefficients [4] and Helmholtz’s
equations [6] is due to the poor approximation property.
휂
(
푆
)
is typically not small if
푆
is the standard finite element space, such as the space of tent functions.
Conceptually, ExpMsFEM finds an exponentially convergent function representa-
tion of the solution through the following three steps: (1) harmonic-bubble splitting,
(2) edge localization, (3) oversampling and exponentially convergent singular value
decomposition (SVD). We will detail the three steps and discuss relevant rigorous re-
sults at the end of subsections 3.2, 3.3, and 3.4. Then, we summarize the algorithm in
subsection 3.5.
3.2.
Harmonic-bubble splitting.
Consider a shape regular and uniform partition of
the domain
Ω
into finite elements with a mesh size
퐻
. The collection of elements is
denoted by
풯
퐻
=
{
푇
1
,푇
2
, ...,푇
푟
}
. Let
ℰ
퐻
=
{
푒
1
, 푒
2
, ..., 푒
푞
}
be the collection of edges in the
interior of
Ω
. We use
풩
퐻
=
{
푥
1
, 푥
2
, ..., 푥
푝
}
to denote the collection of interior nodes. We
also use
퐸
퐻
to denote the collection of interior edges as a set, i.e.,
퐸
퐻
=
–
푒
∈ℰ
퐻
푒
⊂
Ω
. A
more detailed explanation of the mesh structure can be found in [11, 12].
In each element
푇
∈풯
퐻
, we decompose the solution
푢
into
푢
=
푢
h
푇
+
푢
b
푇
such that
(3.2)
−∇·(
퐴
∇
푢
h
푇
)+
푉푢
h
푇
=
0
,
in
푇
푢
h
푇
=
푢,
on
휕
푇
\(
Γ
1
∪
Γ
2
)
푢
h
푇
=
0
,
on
휕
푇
∩
Γ
1
퐴
∇
푢
h
푇
·
휈
=
훽
푢
h
푇
,
on
휕
푇
∩
Γ
2
,
−∇·(
퐴
∇
푢
b
푇
)+
푉푢
b
푇
=
푓 ,
in
푇
푢
b
푇
=
0
,
on
휕
푇
\(
Γ
1
∪
Γ
2
)
푢
b
푇
=
0
,
on
휕
푇
∩
Γ
1
퐴
∇
푢
b
푇
·
휈
=
훽
푢
b
푇
,
on
휕
푇
∩
Γ
2
.
In short,
푢
h
푇
incorporates the interior boundary value of
푢
on the element, while
푢
b
푇
contains information of the right-hand side. All equations in (3.2) should be understood
in the standard weak sense as in (2.2).
We can further define a global decomposition
푢
=
푢
h
+
푢
b
, such that for each
푇
, it
holds that
푢
h
(
푥
)
=
푢
h
푇
(
푥
)
,
푢
b
(
푥
)
=
푢
b
푇
(
푥
)
when
푥
∈
푇
. Here, the component
푢
h
푇
(resp.
푢
h
)
1
Sometimes,
푁
is chosen to be the solution operator of the adjoint equation; for example see [35].
4
YIFAN CHEN, THOMAS Y. HOU, YIXUAN WANG
is called the local (resp. global)
harmonic part
,
푢
b
푇
(resp.
푢
b
) is the local (resp. global)
bubble part
, of the solution
푢
. Here, the harmonic part
푢
h
is not necessarily a harmonic
function due to the existence of
퐴
and
푉
, but it has a similar low complexity property
that a harmonic function has, due to the iterative argument of Caccioppoli’s inequality
first proposed in [2]. We will discuss this low complexity property in subsection 3.4.
Now, in the representation
푢
=
푢
h
+
푢
b
, the part
푢
b
can be directly computed by
solving local problems in parallel since the local boundary conditions are all known.
We are left to deal with the part
푢
h
.
Remark
3.1
.
We discuss several theoretical concerns and possible generalizations below:
•
A sufficient condition for the local components in (3.2) to be well-defined is that the
operator
푢
→−∇·(
퐴
∇
푢
)+
푉푢
(as well as the corresponding boundary conditions) is
elliptic in each local element, implied by the Poincaré inequality. In [11], we consider
elliptic equations with
푉
=
0
and
Γ
2
=
∅
, so this condition is satisfied. In [12], we
consider the Helmholtz equation where
푉
<
0
,
|
푉
|
=
푂
(
푘
2
)
and
Re
훽
=
0
,
Im
훽
=
푂
(
푘
)
. For such a case, the elliptic property is guaranteed when
퐻
=
푂
(
1
/
푘
)
.
•
For the global components
푢
h
,
푢
b
to be well-defined, we need the condition that the
solution
푢
is continuous. This can be guaranteed by the
퐶
훼
estimates of the equation
(2.1) under the assumptions mentioned earlier; see discussions in [11, 12].
•
We can generalize the above decomposition to PDEs with inhomogeneous boundary
conditions. To achieve so, we incorporate these boundary data into the equation for
푢
b
; see also Section 5.3 in [12] for a concrete example of problem with inhomogeneous
boundary data.
3.3.
Edge localization.
The next step is to find some local basis functions that accu-
rately approximate
푢
h
. ExpMsFEM uses the idea of
edge localization
to localize this
approximation task.
First, we define the “harmonic extension” operator
푄
퐸
퐻
that maps the edge values
̃
푢
h
=
푢
h
|
퐸
퐻
∈
퐻
1
/
2
(
퐸
퐻
)
to
푢
h
∈
퐻
1
(
Ω
)
, through the relation in the first set of equation
in (3.2). Here, we adopt the convention that if we write a tilde on the top of a function,
it is the restriction of this function on the edge set. We have that
푢
h
=
푄
퐸
퐻
̃
푢
h
=
푄
퐸
퐻
̃
푢
,
since
푢
h
and
푢
have the same edge values.
Then, let
퐶
(
퐸
퐻
)
be the space of continuous functions on
퐸
퐻
. We consider the edge
interpolation operator
퐼
퐻
:
퐻
1
/
2
(
퐸
퐻
)∩
퐶
(
퐸
퐻
)→
퐻
1
/
2
(
퐸
퐻
)∩
퐶
(
퐸
퐻
)
such that
퐼
퐻
̃
푢
=
’
푥
푖
∈풩
퐻
̃
푢
(
푥
푖
)
̃
휓
푖
where the edge function
̃
휓
푖
is linear on
퐸
퐻
and satisfies
̃
휓
푖
(
푥
푗
)
=
훿
푖푗
. Note that by the
convention of our notation we have
휓
푖
=
푄
퐸
퐻
̃
휓
푖
∈
퐻
1
(
Ω
)
. It is worth noting that
휓
′
푖
푠
are the basis functions used in the vanilla MsFEM.
With the interpolation operator, we can write
푄
퐸
퐻
̃
푢
=
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)+
’
푥
푖
∈풩
퐻
푢
(
푥
푖
)
휓
푖
.
Now, the residue
̃
푢
−
퐼
퐻
̃
푢
is zero at each interior node. This property allows us to
localize the residue to each edge. Indeed, by an abuse of notation, we can write
(3.3)
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)
=
’
푒
∈ℰ
퐻
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
,
where we equate the function
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
that is defined on
푒
to its zero extension to
퐸
퐻
, so that
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
∈
퐻
1
/
2
(
퐸
퐻
)
and thus
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
makes sense.
EXPONENTIALLY CONVERGENT MULTISCALE FINITE ELEMENT METHOD
5
Therefore, we localize the approximation task of
푢
h
to
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
, which is
defined for each edge
푒
.
Remark
3.2
.
Again, we discuss several theoretical concerns below:
•
Once the condition in Remark 3.1 is satisfied, the extension operator
푄
퐸
퐻
is well-
defined because the local equation is elliptic.
•
According to the comment in Remark 3.1, the solution
푢
is continuous, so the nodal
interpolation
퐼
퐻
̃
푢
is well-defined.
•
One can rigorously show that if we can approximate each local term with
‖
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
−
푤
푒
‖
ℋ(
Ω
)
≤
휖
푒
,
then the global approximation error satisfies
‖
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)−
’
푒
∈ℰ
퐻
푤
푒
‖
2
ℋ(
Ω
)
≤
퐶
mesh
’
푒
∈ℰ
퐻
휖
2
푒
,
where
퐶
mesh
is a constant dependent on the mesh structure only. In our previous
work [11, 12], we formalize the approximation in the edge space via the
퐻
1
/
2
00
(
푒
)
norm, which is equivalent to the
ℋ(
Ω
)
norm here after the extension by
푄
퐸
퐻
; see
Proposition 2.5 and Theorem 2.6 in [11]. In this review paper, we explain the ideas
using
푄
퐸
퐻
rather than
퐻
1
/
2
00
(
푒
)
, since the former is more concise in an algorithm-
focused exposition.
We call the step from local approximation to global approximation
edge coupling
.
3.4.
Exponentially convergent SVD.
Recall that by using the harmonic-bubble split-
ting and edge localization, we get the representation
(3.4)
푢
=
푢
h
+
푢
b
=
’
푒
∈ℰ
퐻
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
+
’
푥
푖
∈풩
퐻
푢
(
푥
푖
)
휓
푖
+
푢
b
.
ExpMsFEM then relies on oversampling and local SVD to get an exponentially con-
vergent approximation of each
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
. For each
푒
, consider an oversampling
domain
푤
푒
⊃
푒
. Any domain containing
푒
in the interior may be used, and as an
illustrative example, we set
휔
푒
=
ÿ
{
푇
∈풯
퐻
:
푇
∩
푒
≠
∅}
.
An illustration of this choice for a quadrilateral mesh is given in Figure 1.
푒
휔
푒
interior edge
edge connected to boundary
푒
휔
푒
Figure 1.
Illustration of oversampling domains. On the right, we use an
edge connected to the upper boundary as an illustrating example.
We can view
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
as the image of an operator acting on
푢
|
휔
푒
∈
퐻
1
(
휔
푒
)
.
We denote this operator by
푅
푒
such that
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
=
푄
퐸
퐻
푅
푒
(
푢
|
휔
푒
)
. Now, we
apply the harmonic-bubble splitting in subsection 3.2 to the domain
휔
푒
, which leads to
푢
|
휔
푒
=
푢
h
휔
푒
+
푢
b
휔
푒
. It follows that
(3.5)
푄
퐸
퐻
(
̃
푢
−
퐼
퐻
̃
푢
)|
푒
=
푄
퐸
퐻
푅
푒
푢
h
휔
푒
+
푄
퐸
퐻
푅
푒
푢
b
휔
푒
.
6
YIFAN CHEN, THOMAS Y. HOU, YIXUAN WANG
The term
푅
푒
푢
h
휔
푒
is a restriction of a harmonic part. As we mentioned at the beginning
of this article, one can prove that the restriction operator acting on harmonic-type
functions is of
low approximation complexity
. More precisely, consider the space of
harmonic parts in
휔
푒
, defined via
(3.6)
푈
(
휔
푒
)
:
=
{
푣
∈ℋ(
휔
푒
)
:
−∇·(
퐴
∇
푣
)+
푉푣
=
0
,
in
휔
푒
퐴
∇
푣
·
휈
=
훽
푣,
on
Γ
1
∩
휕휔
푒
}
.
The space is equipped with the norm
‖·‖
ℋ(
휔
푒
)
. Then, one can show that the left singular
values (in descending order) of the local operator
푄
퐸
퐻
푅
푒
:
(
푈
(
휔
푒
)
,
‖·‖
ℋ(
휔
푒
)
)→(ℋ(
Ω
)
,
‖·‖
ℋ(
Ω
)
)
decays as
휆
푒,푚
≤
퐶
exp
(−
푏푚
1
푑
+
1
)
in dimension
푑
, for some generic constant
퐶,푏
inde-
pendent of
푚
and
퐻
. Equivalently, if we write the left singular vectors as
푣
푒,푚
∈
퐻
1
(
Ω
)
,
which is local and supported in the neighboring elements of the edge
푒
, then there
exists some coefficient
푏
푒,푗
such that
(3.7)
‖
푄
퐸
퐻
푅
푒
푢
h
휔
푒
−
’
1
≤
푗
≤
푚
푏
푒,푗
푣
푒,푗
‖
ℋ(
Ω
)
≤
퐶
exp
(−
푏푚
1
푑
+
1
)‖
푢
h
휔
푒
‖
ℋ(
휔
푒
)
.
For more details, see Theorem 3.10 in [12]. Then, summing these local errors up, we
get
(3.8)
’
푒
∈ℰ
퐻
‖
푢
h
휔
푒
‖
2
ℋ(
휔
푒
)
≤
2
’
푒
∈ℰ
퐻
(‖
푢
|
휔
푒
‖
2
ℋ(
휔
푒
)
+‖
푢
b
휔
푒
‖
2
ℋ(
휔
푒
)
)
=
푂
(‖
푢
‖
2
ℋ(
Ω
)
+‖
푓
‖
2
퐿
2
(
Ω
)
)
,
where we used the fact that
‖
푢
b
휔
푒
‖
ℋ(
휔
푒
)
=
푂
(‖
푓
‖
퐿
2
(
휔
푒
)
)
by the elliptic estimate.
Combining the above estimates with edge coupling in Remark 3.2, we get the repre-
sentation
(3.9)
푢
=
푢
h
+
푢
b
=
’
푒
∈ℰ
퐻
’
1
≤
푗
≤
푚
푏
푒,푗
푣
푒,푗
+
’
푥
푖
∈풩
퐻
푢
(
푥
푖
)
휓
푖
+
푢
n
+
푂
(
exp
(−
푏푚
1
푑
+
1
)(‖
푢
‖
ℋ(
Ω
)
+‖
푓
‖
퐿
2
(
Ω
)
)
)
,
where
푢
n
:
=
푢
b
+
Õ
푒
∈ℰ
퐻
푄
퐸
퐻
푅
푒
푢
b
휔
푒
is a part that depends on
푓
locally.
Remark
3.3
.
We discuss several theoretical aspects and the implication of the above
representation.
•
The proof of the exponentially decaying singular values of
푄
퐸
퐻
푅
푒
is based on two
steps. The first step is the iterative argument of Caccioppoli’s inequality, first pro-
posed in [2] and then refined in [31]. It shows that the singular values of the restriction
operator on
푈
(
휔
푒
)
, which restricts a function from the original domain
휔
푒
to a subdo-
main
휔
∗
⊃
푒
, decay nearly exponentially fast. The second step is based on a stability
estimate of the operator
푄
퐸
퐻
푅
푒
acting on
푈
(
휔
∗
)
; see Lemma 3.10 in [11] or Lemma
6.1, 6.2 in [12].
•
We can understand that the oversampling technique is used to take advantage of
the low complexity property of the restriction operator. Historically, the idea of
oversampling was proposed in [22] to reduce the resonance error in MsFEM.
•
The remarkable thing about the representation in (3.9) is the exponentially decaying
error bound.
First, for elliptic equations with rough coefficients, the error bound implies that
these basis functions can capture the behavior of the solution, which is a hard task
for FEMs. Therefore, ExpMsFEM overcomes the difficulty of rough coefficients.
EXPONENTIALLY CONVERGENT MULTISCALE FINITE ELEMENT METHOD
7
Second, for the Helmholtz equation, the stability constant of the solution operator
can depend on
푘
; indeed, this is the main cause of the pollution effect [6]. Denote
the stability constant by
퐶
stab
(
푘
)
such that
‖
푢
‖
ℋ(
Ω
)
≤
퐶
stab
(
푘
)‖
푓
‖
퐿
2
(
Ω
)
. A prevalent
and reasonable assumption on the constant is that of polynomial growth, namely
퐶
stab
(
푘
)≤
퐶
(
1
+
푘
훾
)
for some constants
훾
and
퐶
; see, for example, [27]. In such case,
we can further bound the error by
exp
(−
푏푚
1
푑
+
1
)(‖
푢
‖
ℋ(
Ω
)
+‖
푓
‖
퐿
2
(
Ω
)
)≤
exp
(−
푏푚
1
푑
+
1
)(
퐶
(
1
+
푘
훾
)+
1
)‖
푓
‖
퐿
2
(
Ω
)
.
Therefore, once the number of basis functions per edge
푚
∼
log
푑
+
1
(
푘
)
(logarithmically
on
푘
only), the approximation error can be uniformly small for all
푘
. It implies
that the quantity
휂
(
푆
)
in (3.1) is small, which is important in determining the error
of Galerkin’s methods. In this sense, ExpMsFEM overcomes the difficulty of the
pollution effect by using basis functions whose number scales at most
log
푑
+
1
(
푘
)
.
•
The exponentially accurate representation in (3.9) will not be possible if we do not use
terms dependent on the right-hand side. Indeed, using basis functions independent
of
푓
, the optimal approximation error rate will be algebraic if the right-hand side is in
퐿
2
(
Ω
)
only, due to well-known results in approximation theory (the Kolmogorov
푛
-
width [43, 34]); see also the complexity analysis of the Green function of Helmholtz’s
equation [15]. From this perspective, we can understand that ExpMsFEM breaks the
Kolmogorov barrier by using
nonlinear model reduction
[42], i.e., the basis functions
can depend on the input of the model, here the right-hand side.
3.5.
The solver based on ExpMsFEM.
Now, we can use the representation in (3.9) to
solve the equation efficiently. First, we form
휓
푖
,푣
푒,푗
by computing the local extension
푄
퐸
퐻
̃
휓
푖
for each node and the top-
푚
left singular vectors
푣
푒,푗
,
1
≤
푗
≤
푚
of the local
operator
푄
퐸
퐻
푅
푒
for each
푒
; problems on different nodes and edges are independent
and parallelizable. These become our offline basis functions.
For any right-hand side
푓
, we compute the online part
푢
n
by solving local linear
equations involving
푓
. This step can be parallelized.
Then, we form an effective equation for
푢
−
푢
n
as
(3.10)
푎
(
푢
−
푢
n
,푣
)
=
(
푓 ,푣
)
Ω
−
푎
(
푢
n
,푣
)
,
for any
푣
∈ ℋ(
Ω
)
. We solve the equation for
푢
−
푢
n
using a Galerkin method. As an
example, using the Ritz-Galerkin method, we choose
푆
=
span
{
휓
푖
for
푥
푖
∈풩
퐻
, 푣
푒,푗
for
1
≤
푗
≤
푚, 푒
∈ℰ
퐻
}
,
and find a numerical solution
푢
푆
∈
푆
that satisfies
(3.11)
푎
(
푢
푆
,푣
)
=
(
푓 ,푣
)
Ω
−
푎
(
푢
n
,푣
)
,
for any
푣
∈
푆
. The final numerical solution is given by
푢
푆
+
푢
n
. We call
푢
n
the online
part and
푢
푆
the offline part since
푢
푆
lies in a space that is independent of
푓
.
Note that in the Galerkin method for solving
푢
푆
, the stiffness matrix only needs to be
assembled once and can be used for different
푓
afterward. We can understand (3.10)
as a reduced model of the original equation.
Remark
3.4
.
We discuss several theoretical aspects regarding the effectiveness of the
above method.
•
The accuracy of the numerical solution is due to the quasi-optimality property men-
tioned earlier in subsection 3.1: once
휂
(
푆
)
is small, the solution error is of the same
order compared to the optimal approximation using the basis functions, which is
exponentially small according to the representation (3.9).
8
YIFAN CHEN, THOMAS Y. HOU, YIXUAN WANG
•
When the solution is complex-valued, such as in the Helmholtz equations, we can
use both the Ritz and Petrov versions of the Galerkin methods; for the former, if
푆
≠
푆
, we need to replace
푆
by
푆
+
푆
; see discussions in [12].
•
One thing worth noting is that
‖
푢
n
‖
ℋ(
Ω
)
is of order
푂
(
퐻
)
, due to the standard elliptic
estimate [11, 12]. Therefore, if we aim for
푂
(
퐻
)
accuracy only, we can ignore this
part, and simply setting
푢
n
= 0 in the above algorithm will lead to a solution accurate
up to
푂
(
퐻
)
.
4.
Numerical Experiments
In this section, we present some numerical experiments to demonstrate the effective-
ness of ExpMsFEM. For all the experiments, we consider the domain
Ω =
[
0
,
1
]×[
0
,
1
]
and discretize it by a uniform two-level quadrilateral mesh; see a fraction of this mesh
in Figure 2, where we also show an edge
푒
and its oversampling domain
휔
푒
in solid
lines. The coarse and fine mesh sizes are denoted by
퐻
and
ℎ
, respectively.
푒
휔
푒
coarse mesh
fine mesh
Figure 2.
Two level mesh: a fraction
For a given equation, we compute the reference solution
푢
ref
using the classical FEM
on the fine mesh with a sufficiently small
ℎ
, which we choose to be
ℎ
=
1
/
1024
. By
a
posteriori
estimates, we can check that the fine mesh indeed resolves the corresponding
problems; thus, the associated fine mesh solutions could serve as accurate reference
solutions for all of our numerical examples. In our numerical computation, we solve
local problems that are required in the ExpMsFEM framework using the fine mesh. For
detailed implementation, we refer to [11, 12].
Remark
4.1 (Accuracy on the discrete level)
.
For simplicity of presentation, we do not
provide error analysis of ExpMsFEM on the fully discrete level, where the accuracy of
the local problems can depend on the resolution of the fine grid. For a detailed error
estimate on the fully discrete level in the context of partition of unity methods, see, for
example, [30, 29].
The accuracy of a numerical solution
푢
sol
is computed by comparing it with the
reference solution
푢
ref
on the fine mesh. The accuracy will be measured both in the
퐿
2
norm and energy norm:
(4.1)
푒
퐿
2
=
‖
푢
ref
−
푢
sol
‖
퐿
2
(
Ω
)
‖
푢
ref
‖
퐿
2
(
Ω
)
,
푒
ℋ
=
‖
푢
ref
−
푢
sol
‖
ℋ(
Ω
)
‖
푢
ref
‖
ℋ(
Ω
)
.
In subsection 4.1, we consider an elliptic equation where the coefficient
퐴
(
푥
)
is
periodic but contains multiple scales. This example demonstrates the exponential
accuracy of ExpMsFEM. In subsection 4.2, we consider an elliptic equation where
퐴
(
푥
)