\newsiamremark

remarkRemark \newsiamremarkhypothesisHypothesis \newsiamthmclaimClaim \headersLocalized Subspace Iteration Xiaofei Guan, Lijian Jiang, Yajun Wang, Zihao Yang

Localized subspace iteration methods for elliptic multiscale problems ^†^†thanks: Submitted to the editors DATE. \fundingThe work of the first author was supported by the National Science Foundation of China under grant 12271409, the Interdisciplinary Project in Ocean Research of Tongji University and the Fundamental Research Funds for the Central Universities. The work of the second author was supported by the National Science Foundation of China under grant 12271408.

Xiaofei Guan School of Mathematical Sciences, and Key Laboratory of Intelligent Computing and Applications (Ministry of Education), Tongji University, Shanghai 200092, China (). guanxf@tongji.edu.cn, ljjiang@tongji.edu.cn, 1910733@tongji.edu.cn Lijian Jiang²²footnotemark: 2 Yajun Wang²²footnotemark: 2 Zihao Yang School of Mathematics and Statistics, Northwestern Polytechnical University, Xi’an 710072, China (). yangzihao@nwpu.edu.cn

Abstract

This paper proposes localized subspace iteration (LSI) methods to construct generalized finite element basis functions for elliptic problems with multiscale coefficients. The key components of the proposed method consist of the localization of the original differential operator and the subspace iteration of the corresponding local spectral problems, where the localization is conducted by enforcing the local homogeneous Dirichlet condition and the partition of the unity functions. From a novel perspective, some multiscale methods can be regarded as one iteration step under approximating the eigenspace of the corresponding local spectral problems. Vice versa, new multiscale methods can be designed through subspaces of spectral problem algorithms. Then, we propose the efficient localized standard subspace iteration (LSSI) method and the localized Krylov subspace iteration (LKSI) method based on the standard subspace and Krylov subspace, respectively. Convergence analysis is carried out for the proposed method. Various numerical examples demonstrate the effectiveness of our methods. In addition, the proposed methods show significant superiority in treating long-channel cases over other well-known multiscale methods.

keywords:

Multiscale elliptic problems; Generalized Finite Element Method (GFEM); Localized subspace iteration (LSI); Krylov subspace; Spectral problems

{MSCcodes}

65N99, 65N30, 34E13

1 Introduction

Multiscale problems are important problems in both scientific computing and engineering applications. Examples include diffusion in fractured media [21] and deformation of composite materials with multiple nonseparated length scales [11], etc. However, the computation simulation by traditional numerical methods poses significant challenges because of the expensive computational cost for these problems [10]. This issue is caused by extremely fine computational grids that are required to resolve all relevant scales. To this end, numerical homogenization [31, 32] has been developed and replaced polynomial finite element ansatz functions with more general ansatz functions to overcome the global fine-scale computation. The essence of these methods lies in designing particular finite element ansatz functions that can efficiently capture the problem’s multiscale information.

Babuska, et al.’s pioneering work on the Generalized Finite Element Method (GFEM) [9] suggested that a specific form of basis function is desirable for one-dimensional problems with rough coefficients. Then the idea is extended to two-dimensional problems [5]. The multiscale finite element method (MsFEM) [22, 16] by Hou et al. is another significant milestone. They proposed to generate the harmonic extension of the standard Lagrangian finite element basis functions and to use these harmonic extensions afterward as ansatz functions in the Galerkin method. This method has broad applicability, but it can lead to so-called resonance errors [23, 20] due to the artificial boundary conditions of the local harmonic extension problems. After that, an important method to avoid resonance errors was suggested in [33, 1]. It involves using known global information to figure out better boundary conditions for local problems. If such global information is not available, more sophisticated constructions are necessary to derive suitable multiscale ansatz functions from local problems.

More recently, some new multiscale methods have emerged to tackle this challenge. One way is to use local spectral problems to figure out what information is redundant and what information is important. Then, you can use this information to build multiscale ansatz functions, such as SGFEM [6, 8, 29] and GMsFEM [14], etc. This significantly improves computation accuracy. Another method based on the orthogonal decomposition of the solution space was developed in [30, 19]. It is capable of converting arbitrary finite element (FE) basis functions into ansatz functions incorporating multiscale information. The exponential decay of these multiscale ansatz functions allows them to be localized, which is also known as localized orthogonal decomposition (LOD). By combining this technique with spectral problems, CEM-GMsFEM proposed in [13, 27] is able to improve POD and GMsFEM. There are also numerous methods for constructing some other multiscale basis functions and their variations [2, 24, 38, 37, 4, 17]. The reference list is incomplete. Subsequently, one crucial question is: What are the fundamental principles governing the construction of these basis functions? In our opinion, there are two fundamental principles for the construction of multiscale basis functions. The first is that the basis functions must have localized support, and this ensures that the stiffness matrix is sparse. Furthermore, solving localized problems is computationally efficient and desirable for parallel techniques. This will significantly decrease the CPU time for building the basis functions. Another fundamental principle is that local problems are connected to the inverse operator corresponding to the original problem. Several multiscale methods [30, 19, 13] use the inverse operator in advance, either locally or indirectly. By integrating these two principles, the local inverse operator is crucial in the design of multiscale basis functions.

An easy way to get multiscale basis functions from a standard set of finite element basis functions is to use the orthogonal decomposition technique in the LOD ([30, 19]). Employing the orthogonal decomposition technique several times will result in a more accurate basis function space sequence, and we refer to it as iterative orthogonal decomposition. As proved in [3], the trial and test function spaces obtained by the orthogonal decomposition technique are equivalent to the function spaces obtained by applying the inverse operator and the inverse conjugate operator to the initial basis function space, respectively. The iterative function space sequence that results from iterative orthogonal decomposition converges to the eigenfunction subspace of the inverse operator. At the same time, this iterative function space sequence is consistent with the subspace sequence $\{A^{k}X_{0}\}$ for solving the matrix eigenvalue problem $AX=\lambda X$ [34]. These motivate us to construct multiscale basis function spaces using a variety of subspaces derived from corresponding spectral problems.

The basic idea of the LSI starts with designing the local inverse operators in each local domain $\omega_{i}$ , where $\left\{\omega_{i}\right\}_{i}$ is an overlapping open cover of domain $\Omega$ . Moreover, we opt to enforce a local homogeneous Dirichlet boundary condition on the operator to accomplish localization. Then, novel multiscale methods can be developed by the subspace, which is initially utilized to deal with the spectral problem corresponding to the local inverse operator. Specifically, a localized standard subspace iteration (LSSI) method is proposed by utilizing standard subspaces [34] of local inverse operators. Furthermore, a localized Krylov subspace iteration (LKSI) method is also developed by utilizing the Krylov subspaces [34, 28] of local inverse operators. The proposed multiscale methods are applied to multiscale diffusion and elasticity equations, and they provide better stability than many other multiscale methods in dealing with long-channel problems.

The paper is organized as follows: In Section 2, an inherent relationship between the LOD and the subspace iteration is established. Then, in Section 3, two classic subspaces for solving spectral problems are introduced. Section 4 presents two multiscale methods based on the subspace iteration. In Section 5, the convergence analysis of proposed methods is carried out. In Section 6, a few numerical examples are provided to demonstrate the efficacy and accuracy of the proposed methods. Subsequently, conclusions are made in Section 7.

2 From LOD to spectral problem

2.1 Introduction to LOD

Localized Orthogonal Decomposition (LOD) [30, 19], a highly effective method for constructing multiscale basis functions, has greatly driven the development of many multiscale models. Consider the following elliptic problem:

(1)

\mathcal{L}u(x)=f(x)\quad\text{ in }\Omega,

where $\Omega\in\mathbb{R}^{d}$ is a polyhedral Lipschitz domain, $f\in L^{2}(\Omega)$ denotes a given source term, and $\mathcal{L}$ is a linear partial differential operator with some high-contrast or high-oscillation multiscale coefficients. For simplicity, suppose the above equation satisfies the homogeneous Dirichlet boundary conditions $u=0\text{ in }\partial\Omega$ . Let $V$ denote a Sobolev space that match the problem, then the variational form of (1) is given as

(2)

a(u,v)=\left(f,v\right)_{L^{2}(\Omega)}\text{ for all }v\in V,

where $a(u,v)$ is a bounded sesquilinear form. Let $\mathcal{T}_{H}$ denote a regular finite element mesh of the domain $\Omega$ into closed simplices. Let $\left\{l_{j}\right\}$ denote a set of finite element basis functions, such as the nodal basis functions (hat functions) of the Lagrange finite element space, and define $U_{H}=\text{span}\left\{l_{j}\right\}$ .

The LOD starts with a number of macroscopic quantities of interest, which extract the desired information from the exact solution [3]. These continuous linear functionals are denoted as $q_{j}\in\mathcal{V}^{*},j\in J$ , where $J$ is the finite index set with $N:=|J|$ . Without loss of generality, we assume that these functionals are linearly independent. A canonical choice of functional $q_{j}$ is

(3)

q_{j}:=\left(l_{j},\bullet\right)_{L^{2}(\Omega)}.

In fact, there are numerous alternative selections for the quantities of interest.

Remark 2.1.

For another set of quantities of interest $\tilde{q}_{j}\in\mathcal{V}^{*}$ , it is necessary to make the assumption $\tilde{q}_{j}\in\left[L^{2}(\Omega)\right]^{*}$ . According to Riesz representation theorem, there exists a unique $\tilde{l}_{j}\in\mathcal{V}$ such that

(4)

\tilde{q}_{j}(v)=\left(\tilde{l}_{j},v\right)_{L^{2}(\Omega)},\ \forall\ v\in V.

Hence, each quantity of interest satisfying $\tilde{q}_{j}\in\left[L^{2}(\Omega)\right]^{*}$ is uniquely associated with a corresponding basis function, and vice versa. It is worth noting that the $L^{2}$ inner product used here can be replaced by other well-defined inner products, such as weighted $L^{2}$ inner product defined by

(5)

\left(u,v\right)_{L^{2}(\Omega,\kappa)}:=\int_{\Omega}\kappa uv.

Given the macroscopic quantities of interest $q_{j}$ , we proceed to establish the kernel space

(6)

W:=\left\{v\in V\mid q_{j}(v)=0\text{ for all }j\in J\right\}=\bigcap_{j\in J}% \operatorname{ker}q_{j}.

This space is sometimes referred to as the fine-scale space [25], which encompasses fine-scale information that cannot be captured by the original basis function space $V_{H}$ . It inspired adaptive methods involving the addition of basis functions [26]. The LOD method’s core idea is to find a function space that is orthogonal to the given kernel space $W$ in the sense of sesquilinear form $a\left(\bullet,\bullet\right)$ , serving as the multiscale basis function space. To this end, we define two projections $\mathcal{C}:{V}\rightarrow{W}$ and $\mathcal{C}^{*}:{V}\rightarrow{W}$ , such that

(7)

a(\mathcal{C}v,w)=a(v,w)\text{ and }a(w,\mathcal{C}^{*}v)=a(w,v)\ \forall\ v% \in V,w\in W.

A natural conclusion is that $\mathcal{C}=\mathcal{C}^{*}$ if $a\left(\bullet,\bullet\right)$ is Hermitian. Using the kernel space $W$ and the projections $C$ and $C^{*}$ , the trial and test spaces are constructed by

(8)

\widetilde{U}_{H}:=(1-\mathcal{C})V\quad\text{ and }\quad\widetilde{V}_{H}:=% \left(1-\mathcal{C}^{*}\right)V.

Lemma 2.1.

Given the inf-sup condition, the spaces $\widetilde{U}_{H}$ and $\widetilde{V}_{H}$ possess a dimension of $N:=|J|$ and establish conforming decompositions of the overall space, thereby satisfying

(9)

V=\widetilde{U}_{H}\oplus W\quad\text{ and }\quad V=\widetilde{V}_{H}\oplus W.

Furthermore, we have the ‘orthogonality’ relations

(10)

a\left(\widetilde{U}_{H},W\right)=0\quad\text{ and }\quad a\left(W,\widetilde{% V}_{H}\right)=0.

Proof.

The proof of this lemma can be found in [3].

Lemma 2.2.

Let $\mathcal{L}^{*}$ be the conjugate operator of $\mathcal{L}$ , an alternative characterization of the trial and test spaces is provided by

(11)

\widetilde{U}_{H}=\operatorname{span}\left\{\mathcal{L}^{-1}q_{j}\mid j\in J% \right\}=\operatorname{span}\left\{\mathcal{L}^{-1}l_{j}\mid j\in J\right\}

and

(12)

\widetilde{V}_{H}=\operatorname{span}\left\{\mathcal{L^{*}}^{-1}q_{j}\mid j\in J% \right\}=\operatorname{span}\left\{\mathcal{L^{*}}^{-1}l_{j}\mid j\in J\right\}.

Proof.

The proof of equation (11) can be referenced in [3], and the proof of equation (12) follows a similar approach. Noticed that $a(u,v)=\left(\mathcal{L}u,v\right)_{L^{2}(\Omega)}=\left(u,\mathcal{L^{*}}v% \right)_{L^{2}(\Omega)}$ . Let $v=\mathcal{L^{*}}^{-1}q_{j}$ , then we have $a(w,v)=\left(w,\mathcal{L^{*}}v\right)_{L^{2}(\Omega)}=q_{j}(w)=0$ for all $w\in W$ , which implies that $\operatorname{span}\left\{\mathcal{L^{*}}^{-1}q_{j}\mid j\in J\right\}% \subseteq\widetilde{V}_{H}$ . Notice that both spaces possess dimension N, which leads to the conclusion of the lemma.

Given the trial and test spaces, the discrete variational problem can be expressed as follows: Find $\tilde{u}_{H}\in\widetilde{U}_{H}$ such that

(13)

a\left(\tilde{u}_{H},\tilde{v}_{H}\right)=\left(f,\tilde{v}_{H}\right)_{L^{2}(% \Omega)}\quad\text{ for all }\tilde{v}_{H}\in\widetilde{V}_{H}.

Lemma 2.3.

The bases of $\widetilde{U}_{H}$ and $\widetilde{V}_{H}$ can be obtained by two set of saddle point problems. For $j\in J$ , seek $\tilde{u}_{k}\in V$ and $\mu\in\mathbb{C}^{N}$ , such that

(14)		$\displaystyle a\left(\tilde{u}_{k},v\right)+\sum_{j\in J}\mu_{j}q_{j}(v)$	$\displaystyle=0$		$\displaystyle\text{ for all }v\in V,$
(14)		$\displaystyle q_{j}\left(\tilde{u}_{k}\right)$	$\displaystyle=\delta_{jk}$		$\displaystyle\text{ for all }j\in J.$

And for $j\in J$ , seek $\tilde{v}_{k}\in V$ and $\nu\in\mathbb{C}^{N}$ , such that

(15)		$\displaystyle a\left(v,\tilde{v}_{k}\right)+\sum_{j\in J}\nu_{j}q_{j}(v)$	$\displaystyle=0$		$\displaystyle\text{ for all }v\in V,$
(15)		$\displaystyle q_{j}\left(\tilde{v}_{k}\right)$	$\displaystyle=\delta_{jk}$		$\displaystyle\text{ for all }j\in J.$

Then the spaces $\widetilde{U}_{H}$ and $\widetilde{V}_{H}$ can be represented as follows

(16)

\widetilde{U}_{H}=\text{span}\left\{\tilde{u}_{k}\right\}\text{ and }% \widetilde{V}_{H}=\text{span}\left\{\tilde{v}_{k}\right\}.

Proof.

The proof of this lemma can be found in [3].

Indeed, these saddle point problems can be equivalently formulated as energy minimization problems subject to certain constraints. Exploiting this equivalence, the CEM-GMsFEM (Constrained Energy Minimization Generalized Multiscale Finite Element Method) [13] is introduced.

The aforementioned process, commonly referred to as orthogonal decomposition, allows us to construct the trial and test spaces $\widetilde{U}_{H}$ and $\widetilde{V}_{H}$ from a general basis function space $U_{H}$ , which is capable of accurately solving the original multiscale problem. If the basis functions $l_{j}$ exhibit localization, which means $\text{supp}(l_{j})$ represents a small part of the domain $\Omega$ , it can be demonstrated that the operators $\mathcal{C}$ and $\mathcal{C}^{*}$ possess exponential decay properties. This allows us to localize the projection operators and multiscale basis functions, which is referred to as localized orthogonal decomposition.

Before discussing localization, an interesting question arises: What results can be obtained by repeatedly applying orthogonal decomposition?

2.2 Iterative orthogonal decomposition

An iterative sequence of spaces can be constructed based on multiple iterations of orthogonal decomposition, referred to as iterative orthogonal decomposition. Firstly, we initialize the basis functions

(17)

u_{j}^{0}:=l_{j}\text{ and }v_{j}^{0}:=l_{j}\text{ for all }j\in J.

Define

(18)

\ U_{H}^{0}:=\text{span}\left\{u_{j}^{0}\right\}\text{ and }V_{H}^{0}:=\text{% span}\left\{v_{j}^{0}\right\}.

Then the quantities of interest are defined by

(19)

q_{j}^{n}:=\left(u_{j}^{n},\bullet\right)_{L^{2}(\Omega)}\text{ and }p_{j}^{n}% :=\left(v_{j}^{n},\bullet\right)_{L^{2}(\Omega)}\text{ for }j\in J,n=0,1,2,\cdots.

The kernel spaces are defined by

(20)		$\displaystyle W^{n}:=\left\{v\in V\mid q_{j}^{n}(v)=0\text{ for all }j\in J% \right\}=\bigcap_{j\in J}\operatorname{ker}q_{j}^{n}\text{ for }n=0,1,2,\cdots,$
(20)		$\displaystyle X^{n}:=\left\{v\in V\mid p_{j}^{n}(v)=0\text{ for all }j\in J% \right\}=\bigcap_{j\in J}\operatorname{ker}p_{j}^{n}\text{ for }n=0,1,2,\cdots.$

The projections $\mathcal{C}^{n}:{V}\rightarrow{W^{n}}$ and ${\mathcal{C}^{*}}^{n}:{V}\rightarrow{W^{n}}$ are defined by

(21)

a(\mathcal{C}^{n}v,w)=a(v,w)\text{ and }a(w,{\mathcal{C}^{*}}^{n}v)=a(w,v)\ \ % \forall\ v\in V,w\in W^{n},

for $n=0,1,2,\cdots$ . The trial and test spaces are constructed by

(22)

{U}_{H}^{n+1}:=(1-\mathcal{C}^{n})V\quad\text{ and }\quad{V}_{H}^{n+1}:=\left(% 1-{\mathcal{C}^{*}}^{n}\right)V\text{ for }n=0,1,2,\cdots.

Similarly, the bases of these two sequences of spaces can be constructed by solving two sets of saddle point problems. For $k\in J,n=0,1,2,\cdots$ , seek $u_{k}^{n+1}\in V$ and $\mu^{n+1}\in\mathbb{C}^{N}$ , such that

(23)		$\displaystyle a\left(u_{k}^{n+1},v\right)+\sum_{j\in J}\mu_{j}^{n+1}q_{j}^{n}(v)$	$\displaystyle=0$		$\displaystyle\text{ for all }v\in V,$
(23)		$\displaystyle q_{j}^{n}\left(u_{k}^{n+1}\right)$	$\displaystyle=\delta_{jk}$		$\displaystyle\text{ for all }j\in J.$

And for $k\in J,n=0,1,2,\cdots$ , seek $v_{k}^{n+1}\in V$ and $\nu^{n}\in\mathbb{C}^{N}$ , such that

(24)		$\displaystyle a\left(v,v_{k}^{n+1}\right)+\sum_{j\in J}\nu_{j}^{n+1}p_{j}^{n}(v)$	$\displaystyle=0$		$\displaystyle\text{ for all }v\in V,$
(24)		$\displaystyle p_{j}^{n}\left(v_{k}^{n+1}\right)$	$\displaystyle=\delta_{jk}$		$\displaystyle\text{ for all }j\in J.$

Through iterative orthogonal decomposition, we obtain two sequences of spaces

(25)

U_{H}^{0}\rightarrow U_{H}^{1}\rightarrow U_{H}^{2}\rightarrow\cdots\text{ and% }V_{H}^{0}\rightarrow V_{H}^{1}\rightarrow V_{H}^{2}\rightarrow\cdots.

Naturally, we can select two spaces from the two sequences for a specific step to serve as the trial and test spaces. However, this method is ineffective as a finite element method for two reasons. First, the iterative process usually requires a lot of computational resources; second, the obtained functions lack the localizable property. A very important question is whether limits exist for these two space sequences. If limits do exist, what are their limits?

By Lemma 2.2, we can deduce

(26)

U_{H}^{n+1}=\operatorname{span}\left\{\mathcal{L}^{-1}u_{j}^{n}\mid j\in J% \right\}=\mathcal{L}^{-1}U_{H}^{n}

and

(27)

V_{H}^{n+1}=\operatorname{span}\left\{\mathcal{L^{*}}^{-1}v_{j}^{n}\mid j\in J% \right\}=\mathcal{L^{*}}^{-1}V_{H}^{n}

for $n=0,1,2,\cdots.$ Then we have

(28)

U_{H}^{n}=\mathcal{L}^{-n}U_{H}^{0}\text{ and }V_{H}^{n}=\mathcal{L^{*}}^{-n}V% _{H}^{0}.

The above formulas are naturally associated with the power and subspace iteration methods, which are commonly used to solve spectral problems. Before delving into the discussion of the spectral problem, it is necessary to make some assumptions:

Assumption 2.1.

The linear partial differential operator $\mathcal{L}$ is self-adjoint and positive definite, which means $\mathcal{L}=\mathcal{L^{*}}$ .

Assumption 2.2.

The problem (1) is well-posed, indicating the existence of a unique solution $u\in V$ corresponding to any given $f\in L^{2}(\Omega)$ .

Assumption 2.3.

The Sobolev space $V$ is compactly embedded in $L^{2}(\Omega)$ .

Assumption 2.2 states that the operator $\mathcal{L}^{-1}$ is a bounded operator from $L^{2}(\Omega)$ to $V$ . Together with Assumption 2.3, it can be proven that the operator $\mathcal{L}^{-1}$ is a compact operator defined on $L^{2}(\Omega)$ . These three assumptions make sure that the operator $\mathcal{L}^{-1}$ is both compact and self-adjoint. This allows us to use the spectral decomposition theorem, which is designed to work with compact self-adjoint operators [12].

Let $\left\{\left(\lambda_{i},\phi_{i}\right)|i=0,1,2,\cdots\right\}$ denote the set of all eigenpairs of the inverse operator $\mathcal{L}^{-1}$ in $V$ , where

(29)

\lambda_{0}\geq\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0.

Consequently, it can be deduced that the set $\left\{\phi_{i}\right\}$ forms a set of complete orthogonal bases of $V$ . We can make the assumption that $u_{j}^{0}\in V\text{ for all }j\in J$ . In the event that this assumption is not met, we redefine $u_{j}^{0}:=u_{j}^{1}\in V$ , ensuring that the iterative sequence remains unaltered.

Theorem 2.1.

Let $P_{i}$ is the spectral projector associated with the eigenvalues $\lambda_{1},\lambda_{2},\cdots,\lambda_{i}$ . Assume that $\operatorname{rank}\left(P_{i}\left[u_{1}^{0},u_{2}^{0},\ldots,u_{i}^{0}\right% ]\right)=i$ for $i=1,2,\ldots,N$ and $\lambda_{N}>\lambda_{N+1}$ . We can conclude that the sequence of spaces $\left\{U_{H}^{n}\right\}$ converges to the eigensubspace $V_{\text{eig}}$ spanned by the first $N$ eigenfunctions of the operator $\mathcal{L}^{-1}$ .

Proof.

Let $u_{j}^{0}=\sum_{i=1}^{\infty}\zeta_{j}^{i}\phi_{i}$ for $j\in J$ . Based on the assumption, there exists another set of functions $\left\{\bar{u}_{j}^{0}\right\}$ that can be expressed as linear combinations of $\left\{{u}_{j}^{0}\right\}$ , satisfying

(30)

\bar{u}_{j}^{0}=\sum_{i=j}^{\infty}\bar{\zeta}_{j}^{i}\phi_{i},\text{ for }j% \in J,

where $\bar{\zeta}_{j}^{j}\neq 0$ . Then we have

(31)

\mathcal{L}^{-n}\bar{u}_{j}^{0}=\sum_{i=j}^{\infty}\bar{\zeta}_{j}^{i}\lambda_% {i}^{n}\phi_{i},\text{ for }j\in J.

After normalization,

(32)

\hat{u}_{j}^{n}:=\frac{1}{\lambda_{j}^{n}}\mathcal{L}^{-n}\bar{u}_{j}^{0}=\bar% {\zeta}_{j}^{j}\phi_{j}+\sum_{i=j+1}^{\infty}\bar{\zeta}_{j}^{i}\left(\frac{% \lambda_{i}}{\lambda_{j}}\right)^{n}\phi_{i},\text{ for }j\in J.

It is easy to verify that $U_{H}^{n}=\operatorname{span}\left\{\hat{u}_{j}^{n}\right\}$ . If $\lambda_{j}>\lambda_{j+1}$ , we have $\lim_{n\to\infty}\hat{u}_{j}^{n}=\bar{\zeta}_{j}^{j}\phi_{j}$ . If $\lambda_{j}=\lambda_{j+1}=\cdots=\lambda_{j+k}>\lambda_{j+k+1}$ , we have $\lim_{n\to\infty}\hat{u}_{j+l}^{n}=\sum_{i=j+l}^{j+k}\bar{\zeta}_{j+l}^{i}\phi% _{i}$ for $l=0,1,\cdots,k$ . In conclusion,

(33)

\lim_{n\to\infty}U_{H}^{n}=\operatorname{span}\left\{\lim_{n\to\infty}\hat{u}_% {j}^{n}|j\in J\right\}=\operatorname{span}\left\{\phi_{j}|j\in J\right\}.

the proof is completed.

In fact, orthogonal decomposition can be viewed as a special case of the subspace iteration method. Start with $U_{H}^{n}=\left\{u_{j}^{n}\right\}$ , apply the operator $\mathcal{L}^{-1}$ , and we obtain $\left\{\mathcal{L}^{-1}u_{j}^{n}\right\}$ . In traditional subspace iteration methods, some orthogonalization techniques are used to obtain new basis functions, such as the Gram-Schmidt orthogonalization method. In orthogonal decomposition, based on the saddle problem (23), we seek

(34)

u_{k}^{n+1}=\sum_{j\in J}-\mu_{j}^{n+1}\mathcal{L}^{-1}u_{j}^{n},

such that $\left(u_{k}^{n+1},u_{j}^{n}\right)_{L^{2}(\Omega)}=\delta_{jk}$ . This can be considered as a special orthogonalization method, different from the traditional orthogonalization methods that satisfy $\left(u_{k}^{n+1},u_{j}^{n+1}\right)_{L^{2}(\Omega)}=\delta_{jk}$ .

Remark 2.2.

Orthogonalization techniques in the traditional subspace iteration methods are executed sequentially, facilitating the convergence of the basis functions towards their respective eigenfunctions. In the iterative orthogonal decomposition, the sequence of functions $\left\{u_{j}^{n}\right\}_{n=1}^{\infty}$ does not converge. It is worth noting that the sequence of functions can be decomposed into two alternating subsequences, each of which exhibits convergence, based on empirical observations from numerical experiments.

3 Two subspaces of spectral problem algorithms

In the preceding section, we provide a novel perspective that some multiscale methods can be be regarded as one iteration step under approximating the eigenspace of the corresponding local spectral problems. In fact, a multitude of multiscale methodologies exhibit a strong interconnection with spectral problems [15, 7, 13]. This impetus drives us to delve into broader possibilities within multiscale modeling, commencing with spectral problem algorithms.

We consider a spectral problem

(35)

Ax=\lambda x,

where $A\in\mathbb{C}^{m\times m}$ is hermitian. Let us assume that the eigenvalues are arranged in descending order, which means

(36)

\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{m}.

A classic mission is to find the leading $p$ eigenpairs $\left(\lambda_{i},x_{i}\right)$ for $i=1,2,\cdots,p$ . To achieve this, a natural choice is the subspace iteration method, a foundational and simple method.

3.1 Standard subspace iteration

Algorithm 1 shows the standard subspace iteration, which is capable of approximating the leading $p$ eigenpairs $\left(\lambda_{i},x_{i}\right),i=1,2,\cdots,p$ . The convergence speed for each corresponding eigenfunction $x_{i}$ can be expressed as $\left(\left|\frac{\lambda_{p+1}}{\lambda_{i}}\right|+\epsilon_{n}\right)^{n}$ , where $n$ denotes the number of iteration steps and $\epsilon_{n}$ tends to zero [35].

Algorithm 1 Standard Subspace Iteration

1. Start: Select an initial set of vectors

X_{0}=\left[x_{1}^{0},\ldots,x_{p}^{0}\right]

2. Iterate: Repeat until convergence is achieved, (a) Calculate

X_{k}:=AX_{k-1}

(b) Calculate

X_{k}=QR

the

QR

factorization of

X_{k}

, and set

X_{k}:=Q

The standard subspace iteration is the most fundamental algorithm for solving spectral problems. Compared to other methods, it may not offer a significant advantage. However, its significance lies in its ability to provide insights into multiscale method formulation, which we will discuss in the next section. Furthermore, applying specific projection or preprocessing strategies has the potential to accelerate computational processes.

3.2 Krylov subspace iteration

The Krylov subspace iteration, as a trivial extension of the standard subspace iteration, is one of the most important methods available for computing the eigenvalues and eigenvectors of large matrices, particularly in the Hermitian case. In comparison to the standard subspace iteration, the Krylov subspace iteration demonstrates enhanced efficiency, memory utilization, and flexibility.

For the same spectral problem (35), the Krylov subspace is defined by

(37)

\mathcal{K}_{r}(A,x):=\operatorname{span}\left\{x,Ax,A^{2}x,\ldots A^{r-1}x% \right\}.

If there is no possibility of ambiguity, $\mathcal{K}_{r}(A,x)$ is denoted as $\mathcal{K}_{r}$ . In contrast with the standard subspace iteration, the Krylov subspace iteration necessitates only a single matrix or operator operation at each iteration step, as opposed to multiple operations. A few well-known of these Krylov subspace methods are Arnoldi’s method and Lanczos’ method. This paper will introduce Arnoldi’s method as an illustrative example.

Arnoldi’s method is an orthogonal projection method onto $\mathcal{K}_{r}$ for large sparse matrices. Specifically, the procedure introduced by Arnoldi in 1951 starts by building an orthogonal basis of the Krylov subspace $\mathcal{K}_{r}$ . Subsequently, we approximate eigenpairs within the subspace $\mathcal{K}_{r}$ using orthogonal projection techniques, such as the Rayleigh-Ritz procedure. There are several distinct implementations of Arnoldi’s method, which are all mathematically equivalent, and Algorithm 2 is one of them.

Algorithm 2 Arnoldi’s Krylov subspace iteration

1. Start: Select an initial vector

x_{1}

2. Iterate: Compute for

j=1,2,\cdots,l-1

: (a)

h_{ij}=\left(Ax_{j},x_{i}\right),\quad i=1,2,\ldots,j

, (b)

y_{j}=Ax_{j}-\sum_{i=1}^{j}h_{ij}x_{i}

, (c)

h_{j+1,j}=\left\|y_{j}\right\|_{2},\quad\text{ if }h_{j+1,j}=0\text{ stop }

(d)

x_{j+1}=y_{j}/h_{j+1,j}

3. Orthogonal projection: Use the Rayleigh-Ritz procedure to obtain desired eigenpairs in

\mathcal{K}_{r}=\operatorname{span}\left\{x_{1},x_{2},\cdots,x_{l}\right\}

Indeed, there are numerous efficient methods for spectral problems, such as the Jacobi-Davidson method [36], and this article only highlights a few of the most common methods.

4 From spectral problem to LSI

Virtually every iteration method used for spectral problems has the potential to be extended to a multiscale modeling method. Specifically, we can extract a step or multiple steps from the iterative process to construct a subspace that approximates the eigenfunction subspace. Based on the standard subspace iteration method, we proposed the localized standard subspace iteration (LSSI) method. Furthermore, based on the Krylov subspace iteration, we proposed the localized Krylov subspace iteration (LKSI) method. Before introducing these methods, we start with a localization process, analogous to the technique adopted in the majority of multiscale methods.

Let $\left\{\omega_{i}\right\}_{i=1}^{N_{c}}$ be a open cover of domain $\Omega$ , such that $\Omega=\bigcup_{i=1}^{N_{c}}\omega_{i}$ , where $N_{c}$ is the number of subdomains. An elementary choice for the set $\left\{\omega_{i}\right\}$ is to extend each element $K_{i}$ of the finite element partition $\mathcal{T}_{H}$ by one or several layers. There is a set of partition of unity $\left\{\chi_{i}\right\}_{i=1}^{N_{c}}$ , such that

(38)

1=\sum_{i=1}^{N_{c}}\chi_{i}\text{ and }\operatorname{supp}(\chi_{i})=\omega_{% i}\text{ for }i=1,2,\cdots,N_{c}.

Let $V(\omega_{i}):=\{v\in V|\operatorname{supp}(v)\in\omega_{i}\}$ . The operator $\mathcal{L}$ restricted to the space $V(\omega_{i})$ is denoted as $\mathcal{L}_{i}$ . It is clear that $\mathcal{L}_{i}^{-1}$ , being defined on $L^{2}(\omega_{i})$ , is also compact and self-adjoint. Therefore, we can define local spectral problem

(39)

\mathcal{L}_{i}^{-1}\phi_{i}^{j}=\lambda_{i}^{j}\phi_{i}^{j},\quad\phi_{i}^{j}% \in V\left(\omega_{i}\right),

for $i=1,2,\cdots,N_{c},j=1,2,\cdots$ .

4.1 Localized standard subspace iteration (LSSI) method

For each subdomain $\omega_{i}$ , there is a set of elementary basis functions denoted as $\phi_{i}^{j,0}\in V\left(\omega_{i}\right)$ , where $j=1,2,\cdots,L_{i}$ and $L_{i}$ represents the number of basis functions in subdomain $\omega_{i}$ . Through the standard subspace iteration in each subdomain, a sequence of basis function sets can be acquired. To be specific, for $i=1,2,\cdots,N_{c},k=1,2,\cdots,L_{i},n=0,1,2,\cdots$ , we seek $\phi_{i}^{k,n+1}\in V\left(\omega_{i}\right)$ and $\mu_{i}^{n+1}\in\mathbb{C}^{L_{i}}$ , such that

(40)		$\displaystyle a\left(\phi_{i}^{k,n+1},v\right)+\sum_{j=1}^{L_{i}}\mu_{i}^{j,n+% 1}q_{i}^{j,n}(v)$	$\displaystyle=0$		$\displaystyle\text{ for all }v\in V\left(\omega_{i}\right),$
(40)		$\displaystyle q_{i}^{j,n}\left(\phi_{i}^{k,n+1}\right)$	$\displaystyle=\delta_{jk}$		$\displaystyle\text{ for all }j=1,2,\cdots,L_{i}.$

In the above equation, $q_{i}^{j,n}$ is a linear functional defined by

(41)

q_{i}^{j,n}:=\left(\phi_{i}^{j,n},\bullet\right)_{L^{2}(\Omega)}.

The local multiscale space $V_{\text{S}}^{i,n}$ is constructed by

(42)

V_{\text{S}}^{i,n}:=\text{span}\{\phi_{i}^{j,n}|j=1,2,\cdots L_{i}\},

$\text{ for }i=1,2,\cdots,N_{c},n=0,1,2,\cdots$ . The multiscale space $V_{\text{S}}^{n}$ is constructed by

(43)

V_{\text{S}}^{n}:=\bigoplus_{i=1}^{N_{c}}V_{\text{S}}^{i,n}.

Similar to Eq.(26), it can be deduced that

(44)

V_{\text{S}}^{i,n+1}=\operatorname{span}\left\{\mathcal{L}_{i}^{-1}\phi_{i}^{j% ,n}\mid j=1,2,\cdots,L_{i}\right\}=\mathcal{L}_{i}^{-1}V_{\text{S}}^{i,n}=% \mathcal{L}_{i}^{-(n+1)}V_{\text{S}}^{i,0}.

As a direct corollary of Theorem 2.1,

(45)

\lim_{n\to\infty}V_{\text{S}}^{i,n}=\operatorname{span}\left\{\phi_{i}^{j}|j=1% ,2,\cdots,L_{i}\right\}.

Once the multiscale space $V_{\text{S}}^{n}$ is obtained, the multiscale solution can be obtained using the Galerkin method. Due to its strong resemblance to the standard subspace iteration in spectral problem algorithms, this method is referred to as the localized standard subspace iteration (LSSI) method. The convergence rate of the LSSI depends on the separation of eigenvalues. Given the rapid decay characteristic of eigenvalues in multiscale spectral problems, it often only takes a few iterations to achieve highly satisfactory results.

4.2 Localized Krylov subspace iteration (LKSI) method

When employing the subspace iteration method to solve eigenvalue problems, it is common practice to exclusively use the result of the final iteration step for computation. In reality, the functions obtained at each iteration step can all be used for eigenfunction calculation. This fundamental principle reveals the core of the Krylov method, wherein it exhibits a distinct advantage in terms of computational efficiency and memory usage.

In each subdomain $\omega_{i}$ , suppose there is an initial basis function $\psi_{i}^{0}\in V\left(\omega_{i}\right)$ . For $i=1,2,\cdots,N_{c},n=0,1,\cdots$ , we seek $\psi_{i}^{n+1}$ and $\mu_{i}^{n+1}\in\mathbb{C}$ , such that

(46)		$\displaystyle a\left(\psi_{i}^{n+1},v\right)+\mu_{i}^{n+1}q_{i}^{n}(v)$	$\displaystyle=0$		$\displaystyle\text{ for all }v\in V\left(\omega_{i}\right),$
(46)		$\displaystyle q_{i}^{n}\left(\psi_{i}^{n+1}\right)$	$\displaystyle=1.$

Similarly, the definition of $q_{i}^{n}$ is as follows

(47)

q_{i}^{n}:=\left(\psi_{i}^{n},\bullet\right)_{L^{2}(\Omega)}.

The local multiscale space $V_{\text{K}}^{i,n}$ is constructed by

(48)

V_{\text{K}}^{i,n}:=\operatorname{span}\{\psi_{i}^{k}|k=0,1,2,\cdots n\},

$\text{ for }i=1,2,\cdots,N_{c},n=0,1,2,\cdots$ . The multiscale space $V_{\text{K}}^{n}$ is constructed by

(49)

V_{\text{K}}^{n}:=\bigoplus_{i=1}^{N_{c}}V_{\text{K}}^{i,n}.

It is easy to deduce that

(50)

V_{\text{K}}^{i,n}=\operatorname{span}\left\{\mathcal{L}_{i}^{-k}\psi_{i}^{0}% \mid k=0,1,2,\cdots,n\right\}.

In fact, $V_{\text{K}}^{i,n}$ is the Krylov subspace of the operator $\mathcal{L}_{i}^{-1}$ with the initial function $\psi_{i}^{0}$ .

5 Convergence analysis

In this section, we will establish the convergence of our proposed methods. In order to demonstrate the convergence of the proposed methods, we will initially show the interpolation error of using the local eigenfunctions as basis functions.

5.1 Interpolation error

Define the local eigenfunction space $V_{\text{eig}}$ as follows:

(51)

V_{\text{eig}}:=\operatorname{span}\{\phi_{i}^{j}|i=1,2,\cdots,N_{c},j=1,2,% \cdots,L_{i}\},

where $\phi_{i}^{j}$ is the local eigenfunction difined by Eq. Eq. 39. Then for all $u\in V$ , define interpolation operator $\mathcal{I}_{\text{eig}}:V\rightarrow V_{\text{eig}}$ as follows:

(52)

\mathcal{I}_{\text{eig}}u:=\sum_{i=1}^{N_{c}}\sum_{j=1}^{L_{i}}\left\langle% \chi_{i}u,\phi_{i}^{j}\right\rangle\phi_{i}^{j}.

Theorem 5.1.

Suppose $u\in V$ , then we have an estimation for the interpolation error

(53)

\|u-\mathcal{I}_{\text{eig}}u\|_{a}\leq\sqrt{\lambda^{L+1}}\sum_{i=1}^{N_{c}}% \|\mathcal{L}\chi_{i}u\|_{L^{2}(\Omega)},

where $\lambda^{L+1}:=\mathop{\max}\limits_{i}\ \lambda_{i}^{L_{i}+1}$ , and the energy norm $\|\bullet\|_{a}$ is defined by $\|u\|_{a}^{2}=a(u,u)$ .

Proof.

Noticed that

(54)		$\displaystyle u-\mathcal{I}_{\text{eig}}u$	$\displaystyle=\sum_{i=1}^{N_{c}}\chi_{i}u-\sum_{i=1}^{N_{c}}\sum_{j=1}^{L_{i}}% \left\langle\chi_{i}u,\phi_{i}^{j}\right\rangle\phi_{i}^{j}$
(54)			$\displaystyle=\sum_{i=1}^{N_{c}}\sum_{j=L_{i}+1}^{\infty}\left\langle\chi_{i}u% ,\phi_{i}^{j}\right\rangle\phi_{i}^{j}.$

It is easy to verify that

(55)	$\displaystyle\\|\sum_{j=L_{i}+1}^{\infty}\left\langle\chi_{i}u,\phi_{i}^{j}% \right\rangle\phi_{i}^{j}\\|_{a}^{2}$	$\displaystyle=\sum_{j=L_{i}+1}^{\infty}\left\langle\chi_{i}u,\phi_{i}^{j}% \right\rangle^{2}\\|\phi_{i}^{j}\\|_{a}^{2}$
		$\displaystyle=\sum_{j=L_{i}+1}^{\infty}\left\langle\chi_{i}u,\phi_{i}^{j}% \right\rangle^{2}\frac{1}{\lambda_{i}^{j}}$
		$\displaystyle\leq\lambda_{i}^{L_{i}+1}\sum_{j=L_{i}+1}^{\infty}\left\langle% \chi_{i}u,\phi_{i}^{j}\right\rangle^{2}\left(\frac{1}{\lambda_{i}^{j}}\right)^% {2}$
		$\displaystyle\leq\lambda_{i}^{L_{i}+1}\sum_{j=1}^{\infty}\left\langle\chi_{i}u% ,\phi_{i}^{j}\right\rangle^{2}\left(\frac{1}{\lambda_{i}^{j}}\right)^{2}$
		$\displaystyle=\lambda_{i}^{L_{i}+1}\\|\mathcal{L}\chi_{i}u\\|_{L^{2}(\Omega)}^{2}.$

The last equation is based on the expansion of the operator $\mathcal{L}$ on $V(\omega_{i})$ ,

(56)

\mathcal{L}\chi_{i}u=\sum_{j=1}^{\infty}\frac{1}{\lambda_{i}^{j}}\left\langle% \chi_{i}u,\phi_{i}^{j}\right\rangle\phi_{i}^{j}.

Then we have

(57)	$\displaystyle\\|u-\mathcal{I}_{\text{eig}}u\\|_{a}$	$\displaystyle\leq\sum_{i=1}^{N_{c}}\\|\sum_{j=L_{i}+1}^{\infty}\left\langle\chi% _{i}u,\phi_{i}^{j}\right\rangle\phi_{i}^{j}\\|_{a}$
		$\displaystyle\leq\sum_{i=1}^{N_{c}}\sqrt{\lambda_{i}^{L_{i}+1}}\\|\mathcal{L}% \chi_{i}u\\|_{L^{2}(\Omega)}$
		$\displaystyle\leq\sqrt{\lambda^{L+1}}\sum_{i=1}^{N_{c}}\\|\mathcal{L}\chi_{i}u% \\|_{L^{2}(\Omega)}.$

5.2 Convergence result of the LSSI method

Before the convergence result of the LSSI is obtained, let us give a lemma without proof.

Lemma 5.1.

Suppose $\phi_{i}^{j}$ is the $j$ th eigenfunction of the local spectral problem Eq. 39, and the local multiscale space $V_{\text{S}}^{i,n}$ is defined by Eq. 42. There exists a $\hat{\phi}_{i}^{j,n}$ in $V_{\text{S}}^{i,n}$ , such that the following inequality is satisfied:

(58)

\|\hat{\phi}_{i}^{j,n}-\phi_{i}^{j}\|_{a}\leq C(\phi_{i}^{j},V_{\text{S}}^{i,0% })\left(\frac{\lambda_{i}^{L_{i}+1}}{\lambda_{i}^{j}}+\epsilon_{n}\right)^{n},

where $\epsilon_{n}$ tends to zero as n tends to infinity.

In fact, Let $P_{i}$ is the spectral projector associated with the invariant subspace associated with $\lambda_{i}^{1},\dots,\lambda_{i}^{L^{i}}$ . Then for each $\phi_{i}^{j}$ , there exists a unique $s\in V_{\text{S}}^{i,0}$ such that $P_{i}s=\phi_{i}^{j}$ . The constant $C(\phi_{i}^{j},V_{\text{S}}^{i,0})$ is defined by

(59)

C(\phi_{i}^{j},V_{\text{S}}^{i,0})=\|s-\phi_{i}^{j}\|_{a}.

For the sake of simplicity, $C(\phi_{i}^{j},V_{\text{S}}^{i,0})$ is replaced by a constant $C$ in the following. For details of the proof, we refer the reader to Yousef’s book ([35]; Theorem 5.2).

Based on the auxiliary function $\hat{\phi}_{i}^{j,n}$ in $V_{\text{S}}^{i,n}$ , we define Interpolation operator $\mathcal{I}_{S}:V\rightarrow V_{\text{S}}^{n}$ ,

(60)

\mathcal{I}_{\text{S}}u:=\sum_{i=1}^{N_{c}}\sum_{j=1}^{L_{i}}\left\langle\chi_% {i}u,\phi_{i}^{j}\right\rangle\hat{\phi}_{i}^{j}.

Theorem 5.2.

Suppose $u$ is the solution of Eq. 1, and $u_{\text{S}}$ is the finite element solution in the multiscale space $V_{\text{S}}^{n}$ . The following convergence result holds:

(61)

\|u-u_{S}\|_{a}\leq\sqrt{\lambda^{L+1}}\sum_{i=1}^{N_{c}}\|\mathcal{L}\chi_{i}% u\|_{L^{2}(\Omega)}+C\left(\lambda^{\frac{L+1}{L}}+\epsilon_{n}\right)^{n}\|u% \|_{L^{2}(\Omega)},

where $\lambda^{\frac{L+1}{L}}$ is defined by

(62)

\lambda^{\frac{L+1}{L}}:=\mathop{\max}\limits_{i}\frac{\lambda_{i}^{L_{i}+1}}{% \lambda_{i}^{L_{i}}}.

Proof.

By the céa Lemma and triangular inequality, we have

(63)		$\displaystyle\\|u-u_{S}\\|_{a}$	$\displaystyle=\mathop{\inf}\limits_{v\in V_{\text{S}}^{n}}\\|u-v\\|_{a}\leq\\|u-% \mathcal{I}_{\text{S}}u\\|_{a}$
(63)			$\displaystyle\leq\\|u-\mathcal{I}_{\text{eig}}u\\|_{a}+\\|\mathcal{I}_{\text{eig}% }u-\mathcal{I}_{\text{S}}u\\|_{a}.$

Based on the definitions of $\mathcal{I}_{\text{eig}}$ and $\mathcal{I}_{\text{S}}$ , we have

(64)

\displaystyle\mathcal{I}_{\text{eig}}u-\mathcal{I}_{\text{S}}u=\sum_{i=1}^{N_{% c}}\sum_{j=1}^{L_{i}}\left\langle\chi_{i}u,\phi_{i}^{j}\right\rangle\left(\phi% _{i}^{j}-\hat{\phi}_{i}^{j}\right).

Using the 5.1,

(65)	$\displaystyle\\|\mathcal{I}_{\text{eig}}u-\mathcal{I}_{\text{S}}u\\|_{a}$	$\displaystyle\leq\sum_{i=1}^{N_{c}}\sum_{j=1}^{L_{i}}\left\langle\chi_{i}u,% \phi_{i}^{j}\right\rangle\\|\hat{\phi}_{i}^{j,n}-\phi_{i}^{j}\\|_{a}$
		$\displaystyle\leq C\sum_{i=1}^{N_{c}}\sum_{j=1}^{L_{i}}\left\langle\chi_{i}u,% \phi_{i}^{j}\right\rangle\left(\frac{\lambda_{i}^{L_{i}+1}}{\lambda_{i}^{j}}+% \epsilon_{n}\right)^{n}$
		$\displaystyle\leq C\sum_{i=1}^{N_{c}}\left(\frac{\lambda_{i}^{L_{i}+1}}{% \lambda_{i}^{L_{i}}}+\epsilon_{n}\right)^{n}\sum_{j=1}^{L_{i}}\left\langle\chi% _{i}u,\phi_{i}^{j}\right\rangle$
		$\displaystyle\leq C\left(\lambda^{\frac{L+1}{L}}+\epsilon_{n}\right)^{n}\sum_{% i=1}^{N_{c}}\\|\chi_{i}u\\|_{L^{2}(\Omega)}$
		$\displaystyle\leq C\left(\lambda^{\frac{L+1}{L}}+\epsilon_{n}\right)^{n}\\|u\\|_% {L^{2}(\Omega)}$

Combining Eqs. Eq. 53 and Eq. 65, the proof is completed.

5.3 Convergence result of the LKSI method

Similar to 5.1, we first give the following lemma without proof.

Lemma 5.2.

Suppose $\phi_{i}^{j}$ is the $j$ th eigenfunction of the local spectral problem Eq. 39, and the local multiscale space $V_{\text{K}}^{i,n}$ is defined by Eq. 48. There exists a $\bar{\phi}_{i}^{j,n}$ in $V_{\text{K}}^{i,n}$ , such that the following inequality is satisfied

(66)

\|\bar{\phi}_{i}^{j,n}-\phi_{i}^{j}\|_{a}\leq C\tan\theta\left(\phi_{i}^{j},% \psi_{i}^{0}\right)\frac{\alpha_{i}^{j}}{\left(1+4\gamma_{i}^{j}\right)^{n-j}},

where

(67)

\alpha_{i}^{1}=1,\ \alpha_{i}^{j}=\prod\limits_{k=1}^{j-1}\frac{\lambda_{i}^{k% }}{\lambda_{i}^{k}-\lambda_{i}^{j}}\text{ for }j>1,

and

(68)

\gamma_{i}^{j}=\frac{\lambda_{i}^{j}-\lambda_{i}^{j+1}}{\lambda_{i}^{j+1}}.

For the sake of simplicity, $C\tan\theta\left(\phi_{i}^{j},\psi_{i}^{0}\right)\frac{\alpha_{i}^{j}}{\left(1% +4\gamma_{i}^{j}\right)^{-j}}$ is replaced by a constant $C$ in the following. For details of the proof, we also refer the reader to Yousef’s book ([35]; Theorems 4.8 and 6.3).

Based on the auxiliary function $\bar{\phi}_{i}^{j,n}$ in $V_{\text{K}}^{i,n}$ , we define Interpolation operator $\mathcal{I}_{\text{K}}:V\rightarrow V_{\text{K}}^{n}$ ,

(69)

\mathcal{I}_{\text{K}}u:=\sum_{i=1}^{N_{c}}\sum_{j=1}^{L_{i}}\left\langle\chi_% {i}u,\phi_{i}^{j}\right\rangle\bar{\phi}_{i}^{j}.

Theorem 5.3.

Suppose $u$ is the accuracy solution of Eq. 1, $u_{\text{K}}$ is the finite element solution in the multiscale space $V_{\text{K}}^{n}$ . The following convergence result holds:

(70)

\|u-u_{K}\|_{a}\leq\sqrt{\lambda^{L+1}}\sum_{i=1}^{N_{c}}\|\mathcal{L}\chi_{i}% u\|_{L^{2}(\Omega)}+C\left(\frac{1}{1+4\Gamma}\right)^{n}\|u\|_{L^{2}(\Omega)},

where $\Gamma$ is defined by

(71)

\Gamma:=\mathop{\min}\limits_{i,j}\gamma_{i}^{j}.

The proof of this theorem is similar to the previous theorem, so we won’t go into details here.

6 Numerical experiments

In this section, we present several numerical examples to evaluate the performance of the proposed LSSI and LKSI methods. The methods proposed in this article are effective for symmetric positive-definite differential operators.

For each element $K_{i}$ of the coarse finite element mesh $\mathcal{T}_{H}$ , define the oversampling block $K_{i}^{m}$ as follows,

(72)			$\displaystyle K_{i}^{0}:=K_{i},$
(72)			$\displaystyle K_{i}^{m}:=\text{int}\left(\bigcup T\in\mathcal{T}_{H}\|T\cap% \overline{K}_{i}^{m-1}\neq\emptyset\right),m=1,2,3....$

In this work, the subdomain $\omega_{i}$ is chosen as $K_{i}^{m}$ , where $m$ is the number of oversampling layers. In addition, a fine finite element mesh $\mathcal{T}_{h}$ is used to obtain the reference solution and solve local problems. The errors in the following numerical examples are relative errors compared to the reference solution. Our numerical examples are performed on a desktop workstation with 16 GB of memory and a 3.4GHz Core i7 CPU.

6.1 Diffusion problem

We consider a diffusion problem

(73)

-\nabla\cdot\left(\kappa\nabla u\right)=f\text{ in }\Omega,

where $\Omega=[0,1]^{2}$ and $f=\sin(\pi x)\sin(\pi y)$ . $\kappa$ is a high-contrast permeability coefficient with multiscale characteristics, and is shown in Fig. 1. The fine mesh size is $h=1/100$ , and the coarse mesh size is $H=1/10$ . In the LOD and LSSI, we select $4$ bilinear functions on $K_{i}$ as the initial basis functions $\left\{\phi_{i}^{j,0}\right\}_{j=1}^{4}$ for each subdomain $K_{i}^{m}$ . In the LKSI, we select a piecewise constant on $K_{i}$ as the initial basis function $\psi_{i}^{0}$ for each subdomain $K_{i}^{m}$ . We use the notation ‘ $LSSI\text{-}n$ ’ to represent the LSSI and the ‘ $LKSI\text{-}n$ ’ to represent the LKSI with $n$ iteration steps.

Refer to caption — Figure 1: The high-contrast permeability coefficient $\kappa$ .

Fig. 2 displays solutions of several multiscale methods, including the LOD, LSSI and LKSI, where the number of oversampling layers is $m=4$ . Compared to the reference solution, all multiscale methods have captured the multiscale characteristics of the solution successfully and effectively. For further comparison, Table 1 lists the energy error, $L^{2}$ error, degree of freedom (DoF), CPU time and number of local problems (NoLP) in multiscale methods. With an equivalent degree of freedom, both the LSSI and LKSI exhibit exceptional accuracy, and as the number of iteration steps $n$ increases, the error of the LSSI decays. Table 2 lists the results we are focusing on for the LKSI with different numbers of iteration steps (NoIS) $n$ . As the number of iteration steps increases, there is a linear growth in the degrees of freedom, CPU time, and the number of local problems. Simultaneously, the energy error and $L^{2}$ error decrease at a decelerating rate.

Multiscale method	Energy error	$L^{2}$ error	DoF	CPU time (s)	NoLP
$LOD$	1.6780E-01	3.7891E-02	400	7.75	400
$LSSI\text{-}1$	1.8494E-02	1.0610E-03	400	4.46	400
$LSSI\text{-}2$	1.3449E-02	5.5521E-04	400	6.99	800
$LSSI\text{-}4$	1.2695E-02	5.3744E-04	400	10.97	1600
$LKSI\text{-}4$	1.1997E-02	5.3476E-04	400	9.88	400

Table 1: Comparison of different multiscale methods in terms of the energy error,

L^{2}

error, degree of freedom(DoF), CPU time and number of local problems(NoLP).

NoIS	Energy error	$L^{2}$ error	DoF	CPU time (s)	NoLP
$LKSI\text{-}1$	2.9165E-02	3.3226E-03	100	2.53	100
$LKSI\text{-}2$	1.6337E-02	1.1029E-03	200	5.16	200
$LKSI\text{-}3$	1.2730E-02	6.6714E-04	300	7.35	300
$LKSI\text{-}4$	1.1997E-02	5.3476E-04	400	9.88	400
$LKSI\text{-}5$	1.2298E-02	5.1922E-04	500	11.67	500

Table 2: The energy error,

L^{2}

error, degree of freedom(DoF), CPU time and number of local problems(NoLP) of the LKSI with different number of iteration steps (NoIS)

n

Among the many multiscale methods that employ oversampling techniques, the number of oversampling layers $m$ deserves significant consideration. When the coefficient $\kappa$ is high-contrast, the choice of $m$ in LOD is directly correlated with the contrast ${\kappa_{\max}}/{\kappa_{\min}}$ . This is because the exponential decay rate of global basis functions is associated with the contrast [3]. Fig. 3 shows the relative errors of multiscale methods versus the number of oversampling layers $m$ . In this numerical example, for the LSSI and LKSI, selecting $m=3$ is sufficient, whereas for the LOD, $m$ needs to be at least 4.

Table 3 and Table 4 present the energy errors and $L^{2}$ errors of multiscale methods in relation to the contrast of $\kappa$ , with a fixed number of oversampling layers $m=4$ . In the LOD, as the power exponent of contrast increases, there is a sharp rise in both the energy error and the $L^{2}$ error, reaching unacceptable levels. It is noteworthy that this issue can be mitigated by employing larger number of oversampling layers $m$ . However, the relative errors are very consistent within the frameworks of LSSI and LKSI, even when there are big changes in the power exponent of contrast. Under appropriate conditions, we can argue that the relative errors of our proposed LSSI and LKSI are independent of contrast.

Contrast	$LOD$	$LSSI\text{-}1$	$LSSI\text{-}2$	$LSSI\text{-}4$	$LKSI\text{-}4$
1E+02	8.9027E-02	2.3644E-02	1.9436E-02	1.1528E-02	1.1330E-02
1E+03	1.2928E-01	2.0283E-02	1.5518E-02	1.2251E-02	1.1951E-02
1E+04	1.6780E-01	1.8494E-02	1.3449E-02	1.2695E-02	1.1997E-02
1E+05	3.0026E-01	1.7854E-02	1.2740E-02	1.2738E-02	1.1955E-02
1E+06	6.4146E-01	1.7532E-02	1.2615E-02	1.2621E-02	1.1906E-02
1E+07	9.0074E-01	1.7396E-02	1.2557E-02	1.5008E-02	1.1891E-02

Table 3: Energy errors of multiscale methods versus the contrast of

\kappa

, where the number of oversampling layers

m=4

Contrast	$LOD$	$LSSI\text{-}1$	$LSSI\text{-}2$	$LSSI\text{-}4$	$LKSI\text{-}4$
1E+02	1.3417E-02	1.2870E-03	9.8800E-04	4.5969E-04	4.9329E-04
1E+03	2.6269E-02	1.1496E-03	6.7508E-04	4.9034E-04	5.3315E-04
1E+04	3.7892E-02	1.0610E-03	5.5521E-04	5.3744E-04	5.3476E-04
1E+05	9.3800E-02	1.0351E-03	5.1829E-04	5.5103E-04	5.2758E-04
1E+06	4.1972E-01	1.0314E-03	5.2722E-04	5.3654E-04	5.2040E-04
1E+07	8.2493E-01	1.0340E-03	5.2610E-04	6.3565E-04	5.1874E-04

Table 4:

L^{2}

errors of multiscale methods versus the contrast of

\kappa

, where the number of oversampling layers

m=4

Fig. 4 shows the relative errors of multiscale methods versus coarse mesh size $H$ , where the fine mesh size is $h=1/180$ and the number of oversampling layers is $m=\lceil 2\log(1/H)\rceil$ . For the LOD, the convergence rate of relative energy error is 1, and that of the relative $L^{2}$ error is 2, under an appropriate number of oversampling layers $m$ , just as concluded in [30]. Despite the fact that the convergence rates of relative errors in our proposed methods are lower than those of the LOD, the values themselves are significantly lower.

6.2 Discussion of the channel length

In multiscale problems with high-contrast coefficients, the focal challenge revolves around dealing with contrast, exemplified by the investigation of methods that remain independent of contrast variations. Interestingly, the geometric features of the large coefficient regions can also affect the results of multiscale methods. The influence of geometric features is very complex, and this article only focuses on channel length in simple cases.

We also consider problem Eq. 73, and Fig. 5 shows multiscale coefficients $\kappa$ with different channel lengths (from $2H$ to $10H$ ). The fine mesh size is $h=1/200$ , the coarse mesh size is $H=1/20$ , the number of oversampling layers is $m=5$ and the contrast is $10^{4}$ . Other problem settings are the same as in the previous numerical example. Fig. 6 shows the relative errors of multiscale methods versus the channel length. When the channel length is no higher than 5, the relative errors of the LOD are significantly better than those of our proposed methods. But when the channel length is greater than 5, the errors of the LOD increase sharply, while those of our proposed method increase slightly. This shows that our proposed method has stronger stability for long-channel cases.

6.3 Elasticity problem

The methods we introduced are not limited to diffusion problems; in fact, they hold for general positive definite operators. We consider an elasticity problem

(74)

-\nabla\cdot\mathbf{\sigma}\big{(}\mathbf{u}\big{)}=\mathbf{f},\text{ in }{% \Omega},\\

where $\Omega=[0,1]^{2}$ and $\mathbf{f}=[\sin(\pi x)\sin(\pi y),1]$ . The stress-strain relationship is given by

(75)

\mathbf{\sigma}\big{(}\mathbf{u}\big{)}=2{\mu}\epsilon(u)+{\lambda}\nabla\cdot% \mathbf{u}I,

where $\lambda>0$ and $\mu>0$ are the Lam $\acute{\text{e}}$ constants. The strain tensor $\epsilon(\mathbf{u})=(\epsilon_{ij}(\mathbf{u}))_{1\leq i,j\leq 2}$ is defined by

(76)

\epsilon(\mathbf{u})=\frac{1}{2}\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}% \right).

The Lam $\acute{\text{e}}$ constants $\lambda$ and $\mu$ are the same as $\kappa$ in Section 6.1. The fine mesh size is $h=1/100$ , and the coarse mesh size is $H=1/10$ . In the LOD and LSSI, we select $8$ bilinear functions on $K_{i}$ as the initial basis functions $\left\{\phi_{i}^{j,0}\right\}_{j=1}^{8}$ for each subdomain $K_{i}^{m}$ . In the LKSI, we select piecewise constant on $K_{i}$ as the initial basis function $\psi_{i}^{0}$ for each subdomain $K_{i}^{m}$ .

Multiscale method	Energy error	$L^{2}$ error	DoF	CPU time (s)	NoLP
$LOD$	1.8414E-01	4.5567E-02	800	49.62	800
$LSSI\text{-}1$	1.7168E-02	1.3079E-03	800	18.30	800
$LSSI\text{-}2$	1.0592E-02	8.1541E-04	800	27.10	1600
$LSSI\text{-}8$	1.7368E-02	1.5965E-03	800	78.94	6400
$LKSI\text{-}6$	9.7085E-03	9.7971E-04	600	42.54	600

Table 5: Comparison of different multiscale methods in terms of the energy error,

L^{2}

error, degree of freedom(DoF), CPU time and number of local problems(NoLP).

NoIS	Energy error	$L^{2}$ error	DoF	CPU time	NoLP
$LSSI\text{-}1$	1.7168E-02	1.3079E-03	800	18.30	800
$LSSI\text{-}2$	1.0592E-02	8.1541E-04	800	27.10	1600
$LSSI\text{-}3$	1.1081E-02	7.9146E-04	800	34.63	2400
$LSSI\text{-}4$	1.2494E-02	9.3155E-04	800	43.46	3200
$LSSI\text{-}5$	1.3972E-02	1.0856E-03	800	51.97	4000
$LSSI\text{-}6$	1.5448E-02	1.2735E-03	800	61.36	4800
$LSSI\text{-}7$	1.6625E-02	1.4493E-03	800	70.01	5600
$LSSI\text{-}8$	1.7308E-02	1.5965E-03	800	78.94	6400

Table 6: The energy error,

L^{2}

error, degree of freedom(DoF), CPU time and number of local problems(NoLP) of the LSSI with different number of iteration steps (NoIS)

n

NoIS	Energy error	$L^{2}$ error	DoF	CPU time	NoLP
$LKSI\text{-}1$	6.0448E-02	2.4621E-02	100	7.59	100
$LKSI\text{-}2$	2.7115E-02	3.9845E-03	200	13.96	200
$LKSI\text{-}3$	1.5114E-02	1.7718E-03	300	20.44	300
$LKSI\text{-}4$	1.2045E-02	1.2045E-02	400	26.7	400
$LKSI\text{-}5$	1.0695E-02	1.1438E-03	500	33.47	500
$LKSI\text{-}6$	9.7085E-03	9.7971E-04	600	42.54	600

Table 7: The energy error,

L^{2}

error, degree of freedom(DoF), CPU time and number of local problems(NoLP) of the LKSI with different number of iteration steps (NoIS)

n

Fig. 7 displays solutions of several multiscale methods, where the number of oversampling layers is $m=4$ . Table 5 lists the energy error, $L^{2}$ error, degree of freedom (DoF), CPU time and number of local problems (NoLP) in various multiscale methods. The results obtained are almost consistent with those obtained in the previous numerical examples. Table 6 and Table 7 list the results of the LSSI and LKSI that we focus on with different number of iteration steps (NoIS) $n$ . For the LSSI, the increase in the number of iteration steps means that the obtained basis functions are closer to the local eigenfunctions, which does not necessarily lead to smaller errors. Typically, achieving excellent results requires only two to three steps. For LKSI, the space obtained at the current iteration step is included in the space obtained at the next iteration step. Therefore, as the number of iteration steps increases, the errors will decrease, and the rate of decrease will slow down.

7 Conclusions

In this paper, we proposed two local subspace iteration methods for elliptic multiscale problems. Localization and the inverse operator are fundamental components of several multiscale methods. Multiple implementations of orthogonal decomposition establish the relationship between the LOD and spectral problem algorithms. Orthogonal decomposition can be regarded as an iteration step in an algorithm designed to solve spectral problems. We presented two compelling examples (the LSSI and LKSI) to illustrate our novel perspective: new multiscale methods can be designed through spectral problem algorithms. Numerical examples demonstrated that the proposed methods exhibit exceptional efficiency and applicability in long-channel diffusion fields, which are challenging for most of previous multiscale methods .

This study focused on the implementation of localization by enforcing the local homogeneous Dirichlet boundary condition. More investigation of localization would be a worthwhile research in the future. For example, in the Generalized Multiscale Finite Element Method (GMsFEM), homogeneous Neumann boundary conditions were used in the localization process. Furthermore, for asymmetric and non-positive definite operators, multiscale methods typically required a special design [18]. Future research will investigate asymmetric and non-positive-definite problems using the proposed multiscale methods.

References

[1] J. E. Aarnes, Y. Efendiev, and L. Jiang, Mixed multiscale finite element methods using limited global information, Multiscale Modeling & Simulation, 7 (2008), pp. 655–676.
[2] G. Allaire, Homogenization and two-scale convergence, SIAM Journal on Mathematical Analysis, 23 (1992), pp. 1482–1518.
[3] R. Altmann, P. Henning, and D. Peterseim, Numerical homogenization beyond scale separation, Acta Numerica, 30 (2021), pp. 1–86.
[4] T. Arbogast, G. Pencheva, M. F. Wheeler, and I. Yotov, A multiscale mortar mixed finite element method, Multiscale Modeling & Simulation, 6 (2007), pp. 319–346.
[5] I. Babuška, G. Caloz, and J. E. Osborn, Special finite element methods for a class of second order elliptic problems with rough coefficients, SIAM Journal on Numerical Analysis, 31 (1994), pp. 945–981.
[6] I. Babuska and R. Lipton, Optimal local approximation spaces for generalized finite element methods with application to multiscale problems, Multiscale Modeling & Simulation, 9 (2011), pp. 373–406.
[7] I. Babuska and R. Lipton, Optimal local approximation spaces for generalized finite element methods with application to multiscale problems, Multiscale Modeling & Simulation, 9 (2011), pp. 373–406.
[8] I. Babuška, R. Lipton, P. Sinz, and M. Stuebner, Multiscale-spectral gfem and optimal oversampling, Computer Methods in Applied Mechanics and Engineering, 364 (2020), p. 112960.
[9] I. Babuška and J. E. Osborn, Generalized finite element methods: their performance and their relation to mixed methods, SIAM Journal on Numerical Analysis, 20 (1983), pp. 510–536.
[10] I. Babuška and J. E. Osborn, Can a finite element method perform arbitrarily badly?, Mathematics of computation, 69 (2000), pp. 443–462.
[11] M. A. Biot, Mechanics of deformation and acoustic propagation in porous media, Journal of applied physics, 33 (1962), pp. 1482–1498.
[12] H. Brézis, Functional analysis, Sobolev spaces and partial differential equations, New York: Springer, 2011.
[13] E. T. Chung, Y. Efendiev, and W. T. Leung, Constraint energy minimizing generalized multiscale finite element method, Computer Methods in Applied Mechanics and Engineering, 339 (2018), pp. 298–319.
[14] Y. Efendiev, J. Galvis, and T. Y. Hou, Generalized multiscale finite element methods (gmsfem), Journal of computational physics, 251 (2013), pp. 116–135.
[15] Y. Efendiev, J. Galvis, and T. Y. Hou, Generalized multiscale finite element methods (gmsfem), Journal of computational physics, 251 (2013), pp. 116–135.
[16] Y. Efendiev and T. Y. Hou, Multiscale finite element methods: theory and applications, vol. 4, Springer Science & Business Media, 2009.
[17] X. Guan, L. Jiang, and Y. Wang, Multiscale model reduction for stochastic elasticity problems using ensemble variable-separated method, Journal of Computational and Applied Mathematics, 421 (2023), p. 114895.
[18] X. Guan, L. Jiang, and Y. Wang, Regularized coupling multiscale method for thermomechanical coupled problems, Journal of Computational Physics, 499 (2024), p. 112737.
[19] P. Henning and A. Målqvist, Localized orthogonal decomposition techniques for boundary value problems, SIAM Journal on Scientific Computing, 36 (2014), pp. A1609–A1634.
[20] P. Henning and D. Peterseim, Oversampling for the multiscale finite element method, Multiscale Modeling & Simulation, 11 (2013), pp. 1149–1175.
[21] U. Hornung and R. E. Showalter, Diffusion models for fractured media, Journal of mathematical analysis and applications, 147 (1990), pp. 69–80.
[22] T. Y. Hou and X.-H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, Journal of computational physics, 134 (1997), pp. 169–189.
[23] T. Y. Hou, X.-H. Wu, and Y. Zhang, Removing the cell resonance error in the multiscale finite element method via a petrov-galerkin formulation, Communications in Mathematical Sciences, 2 (2004), pp. 185–205.
[24] T. J. Hughes, G. R. Feijóo, L. Mazzei, and J.-B. Quincy, The variational multiscale method—a paradigm for computational mechanics, Computer methods in applied mechanics and engineering, 166 (1998), pp. 3–24.
[25] T. J. Hughes, G. R. Feijóo, L. Mazzei, and J.-B. Quincy, The variational multiscale method—a paradigm for computational mechanics, Computer methods in applied mechanics and engineering, 166 (1998), pp. 3–24.
[26] M. G. Larson and A. Målqvist, Adaptive variational multiscale methods based on a posteriori error estimation: energy norm estimates for elliptic problems, Computer methods in applied mechanics and engineering, 196 (2007), pp. 2313–2324.
[27] M. Li, E. Chung, and L. Jiang, A constraint energy minimizing generalized multiscale finite element method for parabolic equations, Multiscale Modeling & Simulation, 17 (2019), pp. 996–1018.
[28] J. Liesen and Z. Strakos, Krylov subspace methods: principles and analysis, Numerical Mathematics and Scie, 2013.
[29] C. Ma, R. Scheichl, and T. Dodwell, Novel design and analysis of generalized finite element methods based on locally optimal spectral approximations, SIAM Journal on Numerical Analysis, 60 (2022), pp. 244–273.
[30] A. Målqvist and D. Peterseim, Localization of elliptic multiscale problems, Mathematics of Computation, 83 (2014), pp. 2583–2603.
[31] A. Målqvist and D. Peterseim, Numerical homogenization by localized orthogonal decomposition, SIAM, 2020.
[32] H. Owhadi and C. Scovel, Operator-Adapted Wavelets, Fast Solvers, and Numerical Homogenization: From a Game Theoretic Approach to Numerical Approximation and Algorithm Design, vol. 35, Cambridge University Press, 2019.
[33] H. Owhadi and L. Zhang, Metric-based upscaling, Communications on Pure and Applied Mathematics, 60 (2007), pp. 675–723.
[34] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.
[35] Y. Saad, Numerical methods for large eigenvalue problems: revised edition, SIAM, 2011.
[36] G. L. Sleijpen and H. A. Van der Vorst, A jacobi–davidson iteration method for linear eigenvalue problems, SIAM review, 42 (2000), pp. 267–293.
[37] E. Weinan and B. Engquist, The heterognous multiscale methods, Communications in Mathematical Sciences, 1 (2003), pp. 87–132.
[38] X.-H. Wu, Y. Efendiev, and T. Y. Hou, Analysis of upscaling absolute permeability, Discrete and Continuous Dynamical Systems Series B, 2 (2002), pp. 185–204.

(57)	$\displaystyle\\|u-\mathcal{I}_{\text{eig}}u\\|_{a}$	$\displaystyle\leq\sum_{i=1}^{N_{c}}\\|\sum_{j=L_{i}+1}^{\infty}\left\langle\chi% _{i}u,\phi_{i}^{j}\right\rangle\phi_{i}^{j}\\|_{a}$
		$\displaystyle\leq\sum_{i=1}^{N_{c}}\sqrt{\lambda_{i}^{L_{i}+1}}\\|\mathcal{L}% \chi_{i}u\\|_{L^{2}(\Omega)}$
		$\displaystyle\leq\sqrt{\lambda^{L+1}}\sum_{i=1}^{N_{c}}\\|\mathcal{L}\chi_{i}u% \\|_{L^{2}(\Omega)}.$

(63)		$\displaystyle\\|u-u_{S}\\|_{a}$	$\displaystyle=\mathop{\inf}\limits_{v\in V_{\text{S}}^{n}}\\|u-v\\|_{a}\leq\\|u-% \mathcal{I}_{\text{S}}u\\|_{a}$
(63)			$\displaystyle\leq\\|u-\mathcal{I}_{\text{eig}}u\\|_{a}+\\|\mathcal{I}_{\text{eig}% }u-\mathcal{I}_{\text{S}}u\\|_{a}.$

(65)	$\displaystyle\\|\mathcal{I}_{\text{eig}}u-\mathcal{I}_{\text{S}}u\\|_{a}$	$\displaystyle\leq\sum_{i=1}^{N_{c}}\sum_{j=1}^{L_{i}}\left\langle\chi_{i}u,% \phi_{i}^{j}\right\rangle\\|\hat{\phi}_{i}^{j,n}-\phi_{i}^{j}\\|_{a}$
		$\displaystyle\leq C\sum_{i=1}^{N_{c}}\sum_{j=1}^{L_{i}}\left\langle\chi_{i}u,% \phi_{i}^{j}\right\rangle\left(\frac{\lambda_{i}^{L_{i}+1}}{\lambda_{i}^{j}}+% \epsilon_{n}\right)^{n}$
		$\displaystyle\leq C\sum_{i=1}^{N_{c}}\left(\frac{\lambda_{i}^{L_{i}+1}}{% \lambda_{i}^{L_{i}}}+\epsilon_{n}\right)^{n}\sum_{j=1}^{L_{i}}\left\langle\chi% _{i}u,\phi_{i}^{j}\right\rangle$
		$\displaystyle\leq C\left(\lambda^{\frac{L+1}{L}}+\epsilon_{n}\right)^{n}\sum_{% i=1}^{N_{c}}\\|\chi_{i}u\\|_{L^{2}(\Omega)}$
		$\displaystyle\leq C\left(\lambda^{\frac{L+1}{L}}+\epsilon_{n}\right)^{n}\\|u\\|_% {L^{2}(\Omega)}$