A super-localized generalized finite element method

Abstract

This paper presents a novel multi-scale method for elliptic partial differential equations with arbitrarily rough coefficients. In the spirit of numerical homogenization, the method constructs problem-adapted ansatz spaces with uniform algebraic approximation rates. Localized basis functions with the same super-exponential localization properties as the recently proposed Super-Localized Orthogonal Decomposition enable an efficient implementation. The method’s basis stability is enforced using a partition of unity approach. A natural extension to higher order is presented, resulting in higher approximation rates and enhanced localization properties. We perform a rigorous a priori and a posteriori error analysis and confirm our theoretical findings in a series of numerical experiments. In particular, we demonstrate the method’s applicability for challenging high-contrast channeled coefficients.

1 Introduction

We consider the numerical solution of a second-order linear elliptic partial differential equation with a strongly heterogeneous coefficient. The coefficient may be non-periodic, with oscillations appearing on several non-separated scales. For such coefficients, classical finite element methods based on problem-independent polynomial ansatz spaces typically yield unsatisfactory approximations, cf. [1]. It is possible to overcome this issue by incorporating problem-specific information into the method’s ansatz space, which is commonly known under the term numerical homogenization and has been an active research field throughout the past decades. For an overview on numerical homogenization, we refer to the recent textbooks [2, 3] and the review article [4].

Under minimal assumptions on the coefficient, numerical homogenization is able to achieve optimal orders of approximation without any pre-asymptotic effects. However, this is not possible without a computational overhead. Compared to classical finite element methods, it is either necessary to consider basis functions with an enlarged support, or to increase the number of basis functions per mesh entity. As prominent examples, we mention the Multiscale Spectral Generalized Finite Element Method (MS-GFEM) [5,6,7], the Adaptive Local Basis (AL-Basis) [8, 9], the Localized Orthogonal Decomposition method (LOD) [10,11,12,13], Rough Polyharmonic Splines (RPS) [14], and gamblets [15].

The above-mentioned approaches can be distinguished into two classes. Methods like MS-GFEM and the AL-Basis first solve local spectral problems in the space of (locally) operator-harmonic functions. The respective ansatz spaces are then constructed by gluing together local eigenfunctions by means of a partition of unity [16, 17]. For such methods, the support of the basis functions is fixed by the choice of partition of unity. For convergence of optimal order, the number of local eigenfunctions taken into account needs to be increased logarithmically with the desired accuracy. In order to make these approaches computationally more efficient, one can use random sampling strategies as proposed, e.g., in [18].

The second class of methods includes the LOD, RPS, and gamblets. The idea is to construct problem-adapted ansatz spaces by applying the solution operator to specific classical finite element spaces with respect to some coarse mesh $\mathcal {T}_H$, which typically does not resolve the coefficient’s oscillations. Due to their connection to isogeometric analysis in the case of constant coefficients, such methods are sometimes referred to as spline-type approaches. While optimal approximation orders of these methods are achieved by design, the true challenge is to construct a local basis of the problem-dependent ansatz space. An almost optimal solution is provided by the LOD, which constructs a fixed number of basis functions per mesh entity that decay exponentially fast with respect to the coarse mesh. This rapid decay enables a localization of the basis functions to $\ell$-th order element patches with diameters of order $\ell H$. For convergence of optimal order, the oversampling parameter $\ell$ needs to be increased logarithmically with the desired accuracy.

Recently, the Super-Localized Orthogonal Decomposition (SLOD) has been proposed in [19] (see also [20,21,22,23,24]). The key contribution of the SLOD is a novel localization strategy that enables a significantly improved localization compared to the LOD. The SLOD constructs rapidly decaying basis functions, which in practice yields to super-exponentially decaying localization errors. Note that a rigorous proof of these super-exponential localization properties is still open, and only a justification based on a conjecture from spectral geometry is available, cf. [19]. Nevertheless, using LOD theory, one can prove the practically pessimistic result that the SLOD at least achieves the exponential localization properties of the LOD. The advantageous localization properties of SLOD result in smaller local patch problems for the basis computation and a sparser coarse system matrix. Numerical experiments indicate that the SLOD outperforms the LOD, achieving similar magnitude errors for significantly smaller oversampling parameters. Until now, for the best practical realization of the SLOD in [19], the stability of the SLOD basis functions cannot be guaranteed a priori. For high-contrast channeled coefficients or convection-dominated regimes, these basis stability issues may deteriorate the method’s approximation quality (see the numerical experiments in Sect. 7).

This paper proposes a novel multi-scale method that, on the one hand, preserves the unique localization properties of the SLOD and, on the other hand, resolves the aforementioned basis stability issues. This is achieved by combining the SLOD with a partition of unity approach. More precisely, locally on nodal patches, we apply the respective local solution operator to classical finite element source terms. Multiplying these spaces with the corresponding hat-functions yields local ansatz spaces with a low effective dimension. Consequently, low-dimensional optimally approximating spaces are constructed by solving local spectral problems. Compared to MS-GFEM methods, the proposed method has the major advantage that the local spectral problems are posed in a space spanned by a small number of deterministic snapshots. Hence, possible random sampling strategies can be avoided. Furthermore, due to their low dimension, the local spectral problems are easy to solve. The global problem-adapted ansatz space is then obtained by gluing together the local optimal approximation spaces using a partition of unity. Compared to, e.g., the LOD and SLOD, where the supports of the basis functions need to be increased with the desired accuracy, the proposed method has local basis functions and instead the number of basis functions per mesh entity is increased. We highlight that, from an application point of view, the proposed multi-scale method is conceptually simple and straightforward to implement. Further, by adapting the polynomial degree of the finite element source terms, one can easily construct higher-order versions of the method. Similarly as for the higher-order LOD [25, 26], one obtains higher-order convergence rates using the regularity of the source only.

We prove that the proposed method possesses the advantageous localization and convergence properties of the SLOD which can be quantified a posteriori. Building on the well-understood theoretical foundation of the LOD, we additionally perform a pessimistic a priori error analysis, proving that the proposed method at least recovers the convergence and localization properties of the LOD. For the method’s higher-order versions, we observe that solely increasing the polynomial degree significantly improves the localization properties. Numerical experiments even suggest that an (almost) local basis exists for sufficiently large polynomial degrees. Another noteworthy contribution of this work is the implementation of the proposed method as well as the SLOD in the Python-library gridlod [27]. This serves the principles of open access and reproducibility while enabling computations on large parallel clusters.

The outline of this paper is as follows. In Sect. 2, we introduce the prototypical elliptic model problem. Section 3 recalls preliminary results, which are then used in Sect. 4 for the definition of the novel multi-scale method. An a posteriori error analysis is presented in Sect. 5 followed by a pessimistic a priori error analysis in Sect. 6. Finally, Sect. 7 presents a series of numerical experiments which confirm our theoretical findings.

2 Model problem

We consider the prototypical second-order elliptic PDE $-\textrm{div}A\nabla u = f$ in weak form with homogeneous Dirichlet boundary conditions on a polygonal/polyhedral Lipschitz domain $\Omega \subset {\mathbb {R}}^d$, $d \in \{2,3\}$. Furthermore, without loss of generality, we assume that $\Omega$ is scaled such that its diameter is of order one. The coefficient function $A\in L^\infty (\Omega ,{\mathbb {R}}^{d\times d})$ may be matrix-valued and is assumed to be symmetric and positive definite almost everywhere. More specifically, we assume that there exist constants $0<\alpha \le \beta < \infty$ such that

$$\begin{aligned} \alpha |\eta \vert ^2 \le (A(x)\eta )\cdot \eta \le \beta |\eta \vert ^2,\qquad x \in \Omega ,\;\eta \in {\mathbb {R}}^d \end{aligned}$$

(2.1)

with $|\cdot \vert$ denoting the Euclidean norm. The weak formulation of the elliptic model problem uses the Sobolev space $\mathcal {V}:=H^1_0(\Omega )$ and the bilinear form $a:\mathcal {V}\times \mathcal {V}\rightarrow {\mathbb {R}}$, given by

$$\begin{aligned} a(u, v) :=\int _\Omega (A\nabla u)\cdot \nabla v\,\text {d}x. \end{aligned}$$

The symmetry and the condition (2.1) ensure that the above bilinear form is an inner product on $\mathcal {V}$. Its induced norm is the energy norm ${\left\Vert \cdot \right\Vert }_{a,\Omega } :=\sqrt{a(\cdot , \cdot )}$, which is equivalent to the canonical Sobolev norm on $\mathcal {V}$. The Lax–Milgram theorem ensures that, for all source terms $f\in L^2(\Omega )$, there exists a unique weak solution $u \in \mathcal {V}$ to the boundary value problem, satisfying

$$\begin{aligned} a(u, v) = (f,v)_{L^2(\Omega )},\qquad \text {for all }v \in \mathcal {V}. \end{aligned}$$

(2.2)

Note that the moderate restriction to source terms in $L^2(\Omega )$ (rather than the dual space $\mathcal {V}^\prime =H^{-1}(\Omega )$) will be essential for the uniform convergence of the numerical homogenization method. We note that the possibly rough coefficient generally prevents $H^2(\Omega )$-regularity of the solution, which would be required by classical finite elements. For a generalization to source terms with less regularity, we refer to [4]. Let us mention that the proposed method is not restricted to the class of elliptic PDEs with Dirichlet boundary conditions. Considering the extensions of the SLOD [20, 21], also an extension of the proposed method to (non-symmetric) coercive operators and even Helmholtz-type problems seems possible. In particular, more general boundary conditions of Neumann and Robin type can be taken into account.

Henceforth, we refer to ${\mathcal {A}}^{-1}:L^2(\Omega )\rightarrow \mathcal {V}$ as the solution operator that maps $f\in L^2(\Omega )$ to the unique solution $u\in \mathcal {V}$ of (2.2). Moreover, for a subdomain $\omega \subset \Omega$, we denote by $a_\omega (\cdot ,\cdot )$ and $\Vert \cdot \Vert _{a,\omega }$ the restriction of the bilinear form $a(\cdot ,\cdot )$ to $\omega$ and the restricted energy norm, respectively. The restricted solution operator subject to homogeneous Dirichlet boundary conditions on $\partial \omega$ is denoted by $\mathcal {A}^{-1}_{\omega }:L^2(\omega )\rightarrow H^1_0(\omega )$.

3 Preliminaries

Let $\mathcal {T}_H$ denote a quasi-uniform coarse mesh of $\Omega$ consisting of closed, simplicial, or quadrilateral shape-regular elements. The subscript H denotes the maximal element diameter, i.e., $H :=\max _{T\in \mathcal {T}_H}\textrm{diam}(T)$ and by $\mathcal {N}_{H}$, we denote the set of all (interior and boundary) vertices of $\mathcal {T}_H$. For the ease of presentation, we henceforth only consider quadrilateral meshes; the extension to triangular meshes is straightforward.

The proposed multi-scale method utilizes the concept of patches. For any $\ell \in {\mathbb {N}}_0$, we define the $\ell$-th order (element) patch of a union of elements $\omega \subset \Omega$ recursively by

$$\begin{aligned} \textsf{N}^0(\omega ) :=\omega ,\qquad \textsf{N}^{\ell + 1}(\omega ) :=\bigcup \left\{ T\in \mathcal {T}_H\,:\,T \cap \textsf{N}^\ell (\omega ) \ne \emptyset \right\} . \end{aligned}$$

3.1 Discontinuous finite element spaces

For a fixed (but arbitrary) polynomial degree p, we denote with

$$\begin{aligned} \mathcal {P}(\mathcal {T}_H) :=\{v \in L^2(\Omega )\,:\,v\vert _T, T \in \mathcal {T}_H, \text { is polynomial of coordinate degree}\le p\} \end{aligned}$$

the non-conforming space (with respect to $\mathcal {V}$) consisting of element-wise defined polynomials. We define the restriction of $\mathcal {P}(\mathcal {T}_H)$ to a subdomain $\omega \subset \Omega$ by

$$\begin{aligned} \mathcal {P}(\omega ) :=\{v \in \mathcal {P}(\mathcal {T}_H)\,:\,{\text {supp}}(v) \subset \omega \}. \end{aligned}$$

One can characterize the space $\mathcal {P}(T)$, $T \in \mathcal {T}_H$ in terms of a suitable orthonormal basis $\{\Theta _{T,j}\,:\,j=1,\dots ,J\}$ with $J:=\dim (\mathcal {P}(T))$, e.g., shifted tensor-product Legendre polynomials. Hence, a local orthonormal basis of $\mathcal {P}(\mathcal {T}_H)$ is given by $\{\Theta _{T,j}\,:\,T \in \mathcal {T}_H, j = 1,\dots ,J\}$.

Let $\Pi _H:L^2(\Omega ) \rightarrow \mathcal {P}(\mathcal {T}_H)$ denote the $L^2$-orthogonal projection, which for each $v \in L^2(\Omega )$, is given by the element-wise equation

$$\begin{aligned} (\Pi _Hv, w)_{L^2(T)} = (v,w)_{L^2(T)}, \qquad \text {for all }w \in \mathcal {P}(T),\;T \in \mathcal {T}_H. \end{aligned}$$

The projection satisfies the following stability and approximation estimates

$$\begin{aligned}{} & {} \Vert \Pi _Hv \Vert _{L^2(T)}\le \Vert v \Vert _{L^2(T)},\qquad \qquad \qquad \quad \qquad v \in L^2(T), T\in \mathcal {T}_H, \end{aligned}$$

(3.1)

$$\begin{aligned}{} & {} \Vert (1-\Pi _H)v \Vert _{L^2(T)}\le C_\textrm{a}H \Vert \nabla v \Vert _{L^2(T)},\qquad v \in H^1(T), T \in \mathcal {T}_H\end{aligned}$$

(3.2)

with constant $C_\textrm{a}>0$ depending only on the regularity of the mesh $\mathcal {T}_H$ and the polynomial degree p, see, e.g., [28].

In addition, we define the broken Sobolev space $H^k(\mathcal {T}_H)$, $k \in {\mathbb {N}}$, by

$$\begin{aligned} H^k(\mathcal {T}_H) :=\big \{ v \in L^2(\Omega )\,:\,v\vert _T \in H^k(T),\;T \in \mathcal {T}_H\big \} \end{aligned}$$

with the seminorm

$$\begin{aligned} |\cdot \vert ^2_{H^k(\mathcal {T}_H)} :=\sum _{T \in \mathcal {T}_H} |\cdot \vert ^2_{H^k(T)}. \end{aligned}$$

3.2 Conforming companion spaces

We next define local conforming companions of the functions $\Theta _{T,j}$, so-called bubble functions, which we denote by $b_{T,j}$. For each element $T \in \mathcal {T}_H$, these functions fulfill, for $1 \le j \le J$,

$$\begin{aligned} b_{T,j}\in H^1_0(T), \qquad \Pi _Hb_{T,j}= \Theta _{T,j}. \end{aligned}$$

(3.3)

We do not require an explicit characterization. However, it is important that such functions actually exist. This is guaranteed by [25, Cor. 3.6] stating that, for all $\Theta _{T,j}$, there exist a corresponding bubble function $b_{T,j}$ such that

$$\begin{aligned} \Vert b_{T,j} \Vert _{L^2(T)} + H\Vert \nabla b_{T,j} \Vert _{L^2(T)} \le C_\textrm{b}\Vert \Theta _{T,j} \Vert _{L^2(T)},\qquad T\in \mathcal {T}_H,j=1,\dots ,J\nonumber \\ \end{aligned}$$

(3.4)

with constant $C_\textrm{b}>0$ depending solely on the mesh regularity of $\mathcal {T}_H$ and the polynomial degree p.

By means of the bubble functions, we can define the operator ${\mathcal {B}}_H$ mapping possibly non-conforming functions to $\mathcal {T}_H$-piecewise bubble functions with the same $L^2$-projection. For any function in $\mathcal {P}(\mathcal {T}_H)$, we uniquely define ${\mathcal {B}}_H$ by setting ${\mathcal {B}}_H\Theta _{T,j}:=b_{T,j}$ for all $T \in \mathcal {T}_H$, $j =1,\dots ,J$. We can extend the operator to $L^2(\Omega )$ by setting

$$\begin{aligned} {\mathcal {B}}_Hv :={\mathcal {B}}_H\Pi _Hv,\qquad \text {for all }v \in L^2(\Omega ). \end{aligned}$$

Clearly, the kernels of the operators $\Pi _H$ and ${\mathcal {B}}_H$ coincide and one can prove the local stability estimate

$$\begin{aligned} \Vert {\mathcal {B}}_Hv \Vert _{L^2(T)} + H\Vert \nabla {\mathcal {B}}_Hv \Vert _{L^2(T)} \le C_\textrm{bo}\Vert v \Vert _{L^2(T)},\qquad v \in L^2(T), T\in \mathcal {T}_H\qquad \end{aligned}$$

(3.5)

with another constant $C_\textrm{bo}>0$ depending solely on the mesh regularity of $\mathcal {T}_H$ and the polynomial degree p. By the definition of ${\mathcal {B}}_H$, we obtain, for all $v \in L^2(\Omega )$, $q \in \mathcal {P}(\mathcal {T}_H)$,

$$\begin{aligned} {( {\mathcal {B}}_Hv\,,\,q )}_{L^2(\Omega )} = {( {\mathcal {B}}_H\Pi _Hv\,,\,q )}_{L^2(\Omega )} = {( \Pi _Hv\,,\,q )}_{L^2(\Omega )} = {( q\,,\,v )}_{L^2(\Omega )}, \end{aligned}$$

(3.6)

i.e., ${\mathcal {B}}_H$ is the $L^2$-projection onto the space of bubble functions.

3.3 Partition of unity

The proposed multi-scale method is based on the framework of partition of unity methods, cf. [16, 17]. Although, there is great flexibility in the choice of such a partition, for simplicity, we restrict ourselves to the hat-functions $\{\Lambda _z\,:\,z \in \mathcal {N}_{H}\}$ corresponding to all (interior and boundary) nodes of $\mathcal {T}_H$. Recall that the hat-function $\Lambda _z$ associated with node $z \in \mathcal {N}_{H}$ is a continuous $\mathcal {T}_H$-piecewise bilinear function uniquely defined by setting its nodal values for all $y \in \mathcal {N}_{H}$ to $\Lambda _z(y) = \delta _{yz}$ with $\delta$ denoting the Kronecker symbol. By definition, the hat-functions have an $L^\infty$-norm of one and, due to the shape-regularity of $\mathcal {T}_H$, their gradients satisfy

$$\begin{aligned} \Vert \nabla \Lambda _z\Vert _{L^{\infty }(\Omega )} \le C_\Lambda H^{-1},\qquad z\in \mathcal {N}_{H} \end{aligned}$$

(3.7)

with constant $C_\Lambda >0$ depending solely on the mesh regularity of $\mathcal {T}_H$. Denoting $\omega _z:={\text {supp}}(\Lambda _z)$, the shape-regularity of $\mathcal {T}_H$ also implies that the supports $\{\omega _z\,:\,z \in \mathcal {N}_{H}\}$ have a finite overlap, i.e., the maximal number of overlapping supports

$$\begin{aligned} C_\textrm{ol}:=\max _{T\in \mathcal {T}_H}\#\{z \in \mathcal {N}_{H}\,:\,T \subset \omega _z\} \end{aligned}$$

(3.8)

is uniformly bounded. In what follows, we abbreviate the node patches around $z \in \mathcal {N}_{H}$ and the element patches around $T \in \mathcal {T}_H$ by $\omega _z^\ell :=\textsf{N}^\ell (\omega _z)$ and $\omega _T^\ell :=\textsf{N}^\ell (T)$, respectively.

4 Multi-scale method

This section introduces the proposed multi-scale method. The local ansatz spaces of the method are constructed by applying the local solution operator on an oversampling domain to piecewise polynomial source terms and by subsequent restriction to a subdomain. Due to the oversampling, the resulting local spaces have a low effective dimension, and thus, low-dimensional optimally approximating spaces are utilized. The ansatz space of the method is obtained by gluing together the low-dimensional local approximation spaces. Note that we consider a fixed polynomial degree p and do not track the dependence of constants on p. Explicitly tracking this dependence would make the analysis less clear and also add no value, as it relies on estimates that are pessimistic in p.

4.1 Local approximation spaces

For any $z \in \mathcal {N}_{H}$, we aim to approximate the restriction of the solution space $\mathcal {V}|_{\omega _z}$ using local approximation spaces. This is accomplished by choosing the local approximation space $\mathcal {V}_{H,z}^{\ell }|_{\omega _z}$ with

$$\begin{aligned} \mathcal {V}_{H,z}^{\ell }:=\textrm{span}\big \lbrace \mathcal {A}^{-1}_{\omega _z^\ell } \, q \,:\,q \in \mathcal {P}(\omega _z^\ell ) \big \rbrace \subset H^1_0(\omega _z^\ell ) \end{aligned}$$

(4.1)

being defined on the oversampling domain $\omega _z^\ell$.

Due to the oversampling, the restricted space contains many redundant functions. This holds, in particular, after the multiplication with the hat-function $\Lambda _z$ when gluing the local approximation spaces together. Hence, we investigate the optimal approximation of $\Lambda _z\mathcal {V}_{H,z}^{\ell }$ by n-dimensional subspaces $Q(n) \subset H^1_0(\omega _z)$. Given the subspace Q(n), its worst-case best approximation error is defined as

$$\begin{aligned} \sup _{v \, \in \, \mathcal {V}_{H,z}^{\ell }} \inf _{w \in Q(n)} \frac{\Vert \Lambda _zv -w \Vert _{a,\omega _z}}{\Vert v \Vert _{a,\omega _z^\ell }}. \end{aligned}$$

Typically, the minimal worst-case best approximation error is referred to as Kolmogorov n-width, cf. [29], and is defined as

$$\begin{aligned} d_n^z(H,\ell ):=\inf _{Q(n) \subset H^1_0(\omega _z)} \sup _{v \, \in \, \mathcal {V}_{H,z}^{\ell }} \inf _{w \in Q(n)} \frac{\Vert \Lambda _zv -w \Vert _{a,\omega _z}}{\Vert v \Vert _{a,\omega _z^\ell }}. \end{aligned}$$

(4.2)

Indeed, there exists a corresponding optimal local approximation space of dimension n, which we explicitly compute. For this, we solve the low-dimensional eigenvalue problem, which seeks eigenpairs $(v,\lambda ) \in \mathcal {V}_{H,z}^{\ell }\times {{\mathbb {R}}}$ such that

$$\begin{aligned} a_{\omega _z}(\Lambda _zv, \Lambda _zw) = \lambda \, a_{\omega _z^\ell }(v,w),\qquad \text {for all }w \in \mathcal {V}_{H,z}^{\ell }. \end{aligned}$$

(4.3)

We denote the eigenfunctions by $\{v_i\,:\,i = 1,\dots ,N:=\dim (\mathcal {V}_{H,z}^{\ell })\}$ assuming an ordering such that the corresponding eigenvalues satisfy $\lambda _1\ge \lambda _2\ge \dots \ge \lambda _N \ge 0$. Consequently, denoting

$$\begin{aligned} \mathcal {V}_{H,z}^{\ell ,n}:=\textrm{span}\{v_i\,:\,i=1,\dots ,n\}, \end{aligned}$$

the optimal local approximation space of dimension n is given by $\Lambda _z\mathcal {V}_{H,z}^{\ell ,n}$.

4.2 Global approximation space

A global approximation space is obtained by gluing together the above local approximation spaces using the partition of unity, i.e.,

$$\begin{aligned} \mathcal {V}_{H}^{\ell ,n}:=\sum _{z \in \mathcal {N}_{H}} \Lambda _z\mathcal {V}_{H,z}^{\ell ,n}. \end{aligned}$$

We measure the overall error when approximating $\Lambda _z\mathcal {V}_{H,z}^{\ell }$ by spaces of dimension n by

$$\begin{aligned} d_n(H,\ell ) :=\max _{z\in \mathcal {N}_{H}} d_n^z(H,\ell ). \end{aligned}$$

(4.4)

Finally, the proposed method seeks $u_{H}^{\ell ,n}\in \mathcal {V}_{H}^{\ell ,n}$ such that

$$\begin{aligned} a(u_{H}^{\ell ,n},v) = {( f\,,\,v )}_{L^2(\Omega )},\qquad \text {for all }v \in \mathcal {V}_{H}^{\ell ,n}. \end{aligned}$$

(4.5)

Note that the mesh size H determines the accuracy of the approximation. The oversampling parameter $\ell$ specifies the size of the local patch problems and determines the method’s localization error, whereas the number of local functions is given by n and needs to be chosen sufficiently large. For a precise choice of the parameters $H,\ell$, and n, we refer to Remarks 5.6 and 6.3. The following two sections are devoted to the theoretical analysis of the proposed multi-scale method. In Sect. 5, we present an a posteriori error analysis, while in Sect. 6, we present a priori error bounds.

5 A posteriori error analysis

Subsequently, we derive an a posteriori error analysis of the proposed method by establishing a connection to the SLOD introduced in [19]. The SLOD is conceptually related to the proposed method, as it also constructs its basis functions by applying the local solution operator to $\mathcal {T}_H$-piecewise polynomial source terms.

5.1 Higher-order SLOD

For the a posteriori error analysis, we first briefly introduce a higher-order variant of the SLOD. Note that this variant only serves theoretical purposes and is not investigated numerically. Let us fix an arbitrary element $T\in \mathcal {T}_H$ and oversampling parameter $\ell \in {\mathbb {N}}$. Henceforth, we drop all fixed indices and denote the $\ell$-th order patch around T just by $\omega :=\omega _T^\ell$. Furthermore, we make the meaningful assumption that no patch $\omega$ coincides with the entire domain $\Omega$.

For its prototypical (global) basis functions $\{\varphi _j \,:\,j =1,\dots ,J \}$ associated to the element T, the SLOD uses the following ansatz

$$\begin{aligned} \varphi _j :={\mathcal {A}}^{-1}g_j \end{aligned}$$

with source terms $g_j\in \mathcal {P}(\omega )$ to be determined subsequently. We obtain a localized approximation $\psi _j \in H^1_0(\omega )$ of the basis function $\varphi _j$ by computing its Galerkin projection onto the local subspace $H^1_0(\omega )$, i.e., $\psi _j \in H^1_0(\omega )$ satisfies

$$\begin{aligned} a_\omega (\psi _j,v) = {( g_j\,,\,v )}_{L^2(\omega )},\qquad \text {for all }v\in H^1_0(\omega ). \end{aligned}$$

(5.1)

For the choice of $g_j$, we recall some notation and results on traces of $H^1(\omega )$-functions (see, e.g., [30] for details). Denoting $U :=\mathcal {V}|_\omega \subset H^1(\omega )$, we introduce the trace operator on $\omega$ restricted to U as

$$\begin{aligned} {\text {tr}}:U \rightarrow X :=\textrm{range}\,{\text {tr}}\subset H^{1/2}(\partial \omega ). \end{aligned}$$

An example of a continuous right-inverse of ${\text {tr}}$ is the ${\mathcal {A}}$-harmonic extension, henceforth denoted by ${\text {tr}}^{-1}$. Given $w \in X$, it satisfies ${\text {tr}} {\text {tr}}^{-1}w= w$ and

$$\begin{aligned} a_\omega ({\text {tr}}^{-1}w,v)= 0,\qquad \text {for all }v \in H^1_0( \omega ). \end{aligned}$$

(5.2)

Using definitions (5.1) and (5.2), and that $(v - {\text {tr}}^{-1}{\text {tr}} v) \in H^1_0(\omega )$, it holds

$$\begin{aligned} a(\psi _j,v) = a_\omega (\psi _j,v) = a_\omega (\psi _j,v-{\text {tr}}^{-1} {\text {tr}} v) = (g_j,v- {\text {tr}}^{-1} {\text {tr}} v)_{L^2(\omega )}. \end{aligned}$$

This result yields, together with the definition of $\varphi _j$ and the local support of $g_j$, the following key observation

$$\begin{aligned} a(\varphi _j-\psi _j,v) = (g_j,v)_{L^2(\omega )}-a_\omega (\psi _j,v)= (g_j, {\text {tr}}^{-1} {\text {tr}}\, v)_{L^2(\omega )},\qquad v \in H^1_0(\Omega ).\nonumber \\ \end{aligned}$$

(5.3)

Hence, we can rephrase the smallness of the localization error as the (almost) $L^2(\omega )$-orthogonality of $g_j$ to the space

$$\begin{aligned} Y:={\text {tr}}^{-1}X\subset U \end{aligned}$$

(5.4)

of ${\mathcal {A}}$-harmonic functions on $\omega$ (which satisfy the homogeneous Dirichlet boundary condition on $\partial \Omega \cap \partial \omega$). Since the restricted $L^2$-projection $\Pi _H|_Y$ has a finite rank of dimension less or equal to $K :=J\cdot \#(\mathcal {T}_{H}\cap \omega )$, there exists a singular value decomposition (SVD) such that

$$\begin{aligned} \Pi _H|_Y\,v = \sum _{k = 1}^{K} \sigma _k (v,w_k)_{H^1(\omega )}\,q_k \end{aligned}$$

(5.5)

where $\sigma _1\ge \dots \ge \sigma _{K}\ge 0$ denote the singular values, $\{q_1,\ldots ,q_K\}$ the $L^2(\omega )$-orthonormal left singular vectors, and $\{w_1,\ldots ,w_K\}$ the $H^1(\omega )$-orthonormal right singular vectors.

We choose the source terms $g_j$ as the left singular vectors corresponding to the J smallest singular values, i.e.,

$$\begin{aligned} g_j :=q_{K-J+j},\qquad j =1,\dots ,J. \end{aligned}$$

(5.6)

This yields

$$\begin{aligned} \sup _{v\in Y:\Vert v\Vert _{H^1(\omega )} = 1} {{( g_j\,,\,v )}_{L^2(\omega )}} \le \sigma _{K-J+1},\qquad j =1,\dots ,J. \end{aligned}$$

which follows directly from the properties of the SVD. We define the quantity $\sigma$ measuring the (quasi-)orthogonality between the $g_j$ and Y as

$$\begin{aligned} \sigma (H,\ell ) :=\max _{T\in \mathcal {T}_H}\sigma _{T}(H,\ell ),\qquad \sigma _{T}(H,\ell ):=\sigma _{K-J+1}. \end{aligned}$$

(5.7)

The quantity $\sigma$ is crucial for the error analysis of the SLOD as it determines the localization error. Note that the dependence of $\sigma$ on the (fixed) polynomial degree p is not made explicit in (5.7). The following remark deals with the decay of $\sigma$ with respect to the oversampling parameter $\ell$ and, for the sake of curiosity, also the polynomial degree p.

Remark 5.1

(Decay of $\sigma$) In [19], it has been numerically observed and conjectured that, for $p = 0$, the quantity $\sigma$ decays super-exponentially as $\ell$ is increased, cf. Fig. 1 (left). A similar decay in $\ell$ can also be observed for $p>0$. In accordance with [19, 20], we state the following conjecture: there exists $C_\sigma >0$ depending algebraically on $H, \ell$ and $C>0$ independent of $H,\ell$ such that

$$\begin{aligned} \sigma (H,\ell ) \le C_\sigma (H,\ell ) \exp \big (-C \ell ^{\frac{d}{d-1}}\big ). \end{aligned}$$

(5.8)

Conversely, for fixed $\ell$, a rapid decay of $\sigma$ in p can be observed, cf. Fig. 1 (right). The low level of magnitude of the last singular values may suggest that the respective source terms correspond to fully local basis functions. However, due to the low levels or singular values, this is difficult to verify numerically.

Fig. 1

Singular values $\sigma _k$ of operator $\Pi _H|_Y$ defined in (5.5) for an interior patch, for different pairs of $\ell ,$ p. The singular values $\sigma _{K-J+1}$ relevant for (5.7) are marked by dashed horizontal lines

Choosing for each patch the basis functions corresponding to the J source terms (5.6), we obtain the SLOD ansatz space

$$\begin{aligned} \mathcal {V}_{H}^{\ell ,\,\textrm{SLOD}}:=\textrm{span}\{\psi _{T,j}^{\ell }\,:\,T\in \mathcal {T}_H, j = 1,\dots ,J\}. \end{aligned}$$

(5.9)

The Galerkin SLOD solution $u_{H}^{\ell ,\,\textrm{SLOD}}\in \mathcal {V}_{H}^{\ell ,\,\textrm{SLOD}}$ then satisfies

$$\begin{aligned} a(u_{H}^{\ell ,\,\textrm{SLOD}},w) = {( f\,,\,w )}_{L^2(\Omega )},\qquad \text {for all }w \in \mathcal {V}_{H}^{\ell ,\,\textrm{SLOD}}. \end{aligned}$$

(5.10)

Note that a reasonable SLOD approximation requires a stable choice of the basis functions in (5.9). However, for large oversampling parameters and patches intersecting the boundary, the choice (5.6) may be insufficient, and a special treatment is required. In [19, App. B], such an algorithm is proposed, curing possible stability and uniqueness issues in practice. Since we still cannot guarantee stability in an a priori manner, we assume that the source terms corresponding to the basis function of  (5.9) form a Riesz basis of $\mathcal {P}(\mathcal {T}_H)$, which is formulated in the following.

Assumption 5.2

(Riesz stability) The set

$$\begin{aligned} \{g_{T,j}^{\ell }\,:\,T\in \mathcal {T}_H, j = 1,\dots ,J\} \end{aligned}$$

is a Riesz basis of $\mathcal {P}(\mathcal {T}_H)$, i.e., there is $C_{\textrm{r}}(H,\ell )>0$ depending polynomially on H and $\ell$ such that, for all $(c_{T,j})_{T \in \mathcal {T}_H,j=1,\dots ,J}$,

$$\begin{aligned} C^{-1}_{\textrm{r}}(H,\ell )\sum _{\begin{array}{c} T\in {\mathcal {T}}_H\\ j=1,\dots ,J \end{array}} c_{T,j}^2 \le \bigg \Vert \sum _{\begin{array}{c} T\in {\mathcal {T}}_H\\ j=1,\dots ,J \end{array}} c_{T,j} g_{T,j}^{\ell }\biggl \Vert _{L^2(\Omega )}^2 \;\le C_{\textrm{r}}(H,\ell ) \sum _{\begin{array}{c} T\in {\mathcal {T}}_H\\ j=1,\dots ,J \end{array}} c_{T,j}^2.\nonumber \\ \end{aligned}$$

(5.11)

5.2 A posteriori error bound using SLOD

We provide an a posteriori error analysis of the proposed method based on SLOD techniques. Conceptually, it is similar to the one for the SLOD, cf. [19, Thm. 6.1], but it additionally includes the local optimal approximation error $d_n$ defined in (4.4).

Theorem 5.3

(A posteriori error bound) Let Assumption 5.2 be satisfied and let u and $u_{H}^{\ell ,n}$ denote the solutions to (2.2) and (4.5), respectively. Then, there exists a constant $C>0$ independent of $H,\ell ,$ and n, such that, for any $f \in H^k(\mathcal {T}_H)$, $k \in {{\mathbb {N}}}_0$,

$$\begin{aligned} \Vert u-&u_{H}^{\ell ,n} \Vert _{a,\Omega }\\&\le C\Big ( H^{s+1}|f\vert _{H^s(\mathcal {T}_H)} + \ell ^{d+1} C_\textrm{r}^{1/2}(H,\ell ) \big (\sigma (H,\ell ) + Hd_n(H,\ell )\big )\Vert f \Vert _{L^2(\Omega )}\Big ) \end{aligned}$$

with $s :=\min \{k,p+1\}$ and the notation $H^0(\mathcal {T}_H) :=L^2(\Omega )$ and $|\cdot \vert _{H^0(\mathcal {T}_H)}:=\Vert \cdot \Vert _{L^2(\Omega )}$.

Proof

The application of Céa’s Lemma yields, for arbitrary $v \in \mathcal {V}_{H}^{\ell ,n}$,

$$\begin{aligned} \Vert u-u_{H}^{\ell ,n} \Vert _{a,\Omega } \le \Vert u-v \Vert _{a,\Omega }. \end{aligned}$$

(5.12)

Let the solution to the (higher-order) collocation variant of the SLOD, cf. [19, Rem. 5.1], with oversampling parameter $m:=\lfloor \ell /2\rfloor$ be denoted by $u_{H}^{m,\,\textrm{SLOD}}$. Given a source term $f \in L^2(\Omega )$, its solution is obtained by the following linear combination of localized basis functions $\psi _{T,j}^{m}$:

$$\begin{aligned} u_{H}^{m,\,\textrm{SLOD}}= \sum _{\begin{array}{c} T\in {\mathcal {T}}_H\\ j=1,\dots ,J \end{array}} c_{T,j} \psi _{T,j}^{m}\end{aligned}$$

with coefficients $c_{T,j}$ that are uniquely defined by

$$\begin{aligned} \sum _{\begin{array}{c} T\in {\mathcal {T}}_H\\ j=1,\dots ,J \end{array}}c_{T,j} g_{T,j}^{m}= \Pi _H f. \end{aligned}$$

(5.13)

Adding and subtracting $u_{H}^{m,\,\textrm{SLOD}}$ in (5.12) and employing the triangle inequality yields

$$\begin{aligned} \Vert u-u_{H}^{\ell ,n} \Vert _{a,\Omega } \le \Vert u-u_{H}^{m,\,\textrm{SLOD}} \Vert _{a,\Omega } + \Vert u_{H}^{m,\,\textrm{SLOD}}-v \Vert _{a,\Omega }. \end{aligned}$$

(5.14)

The first term is the error of the collocation variant of the SLOD. It can be bounded using a higher-order version of [19, Thm. 6.1] stating the existence of a constant $C_\textrm{slod}>0$ independent of H and $m$ such that

$$\begin{aligned} \Vert u-u_{H}^{m,\,\textrm{SLOD}} \Vert _{a,\Omega }\le C_\textrm{slod}\big ( H^{s+1} |f\vert _{H^s(\mathcal {T}_H)}+ C_{\textrm{r}}^{1/2}(H,m){m}^{d/2}\sigma (H,m)\Vert f\Vert _{L^2(\Omega )}\big ). \end{aligned}$$

For the second term in (5.14), we choose $v \in \mathcal {V}_{H}^{\ell ,n}$ as sum of functions $v_z^n \in \mathcal {V}_{H,z}^{\ell ,n}$ to be specified later, i.e.,

$$\begin{aligned} v = \sum _{z \in \mathcal {N}_{H}} \Lambda _z v_z^n. \end{aligned}$$

Using the partition of unity property of the hat-functions $\sum _{z \in \mathcal {N}_{H}}\Lambda _z \equiv 1$, we obtain

$$\begin{aligned} \Vert u_{H}^{m,\,\textrm{SLOD}}-v \Vert _{a,\Omega }^2&= \Big \Vert \sum _{z \in {\mathcal {N}}_H} \Lambda _z(u_{H}^{m,\,\textrm{SLOD}}- v_z^n)\Big \Vert _{a,\Omega }^2 \\ {}&\le C_\textrm{ol}\sum _{z \in {\mathcal {N}}_H}\Vert \Lambda _z(u_{H}^{m,\,\textrm{SLOD}}- v_z^n) \Vert _{a,\omega _z}^2 \end{aligned}$$

with $C_\textrm{ol}$ defined in (3.8) denoting the maximal number of overlapping $\omega _z$. For any $z \in \mathcal {N}_{H}$, we can locally on $\omega _z$ replace $u_{H}^{m,\,\textrm{SLOD}}$ by $u_{H,z}^{m,\,\textrm{SLOD}}$ defined by

$$\begin{aligned} u_{H,z}^{m,\,\textrm{SLOD}}:=\sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}}c_{T,j}\psi _{T,j}^{m}\in H^1_0(\omega _z^\ell ) \end{aligned}$$

with the basis functions $\psi _{T,j}^{m}:={\mathcal {A}}^{-1}_{\omega ^m_T} g_{T,j}^{m}$. As approximation to $u_{H,z}^{m,\,\textrm{SLOD}}$, we use

$$\begin{aligned} v_z :=\sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}} c_{T,j} {\widehat{\psi }}_{T,j}^{m}\in \mathcal {V}_{H,z}^{\ell }\end{aligned}$$

with ${\widehat{\psi }}_{T,j}^{m}:={\mathcal {A}}^{-1}_{\omega _z^\ell } g_{T,j}^{m}\in \mathcal {V}_{H,z}^{\ell }$ being approximations of $\psi _{T,j}^{m}$. We choose $v_z^n \in \mathcal {V}_{H,z}^{\ell ,n}$ as the not necessarily unique element minimizing $\Vert \Lambda _z v_z -\Lambda _z v_z^n \Vert _{a,\omega _z^\ell }$. Abbreviating

$$\begin{aligned} e_z^1 :=u_{H,z}^{m,\,\textrm{SLOD}}- v_z,\qquad e_z^2:=v_z - v_z^n, \end{aligned}$$

and performing the above-mentioned local replacement, we obtain

$$\begin{aligned} \Vert \Lambda _z(u_{H}^{m,\,\textrm{SLOD}}- v_z^n) \Vert _{a,\omega _z}^2 \le 2\big ( \Vert \Lambda _z e_z^1 \Vert _{a,\omega _z}^2+\Vert \Lambda _z e_z^2 \Vert _{a,\omega _z}^2\big ), \end{aligned}$$

(5.15)

where we add and subtract $v_z$ and employ the triangle inequality. Using the product rule and the bound (3.7), we get for the first term

$$\begin{aligned} \begin{aligned}&\Vert \Lambda _z e_z^1 \Vert _{a,\omega _z}^2 \le 2\beta \big (C_\Lambda ^2 H^{-2} \Vert e_z^1 \Vert _{L^2(\omega _z)}^2 + \Vert \nabla e_z^1 \Vert _{L^2(\omega _z)}^2\big ). \end{aligned} \end{aligned}$$

(5.16)

Noting that by $e_z^1 \in H^1_0(\omega _z^\ell )$, we obtain for the first term in (5.16)

$$\begin{aligned} \Vert e_z^1 \Vert _{L^2(\omega _z)}^2 \le \Vert e_z^1 \Vert _{L^2(\omega _z^{\ell })}^2 \le C_\textrm{p}^2 \pi ^{-2}\ell ^2 H^2 \Vert \nabla e_z^1 \Vert _{L^2(\omega _z^\ell )}^2 \end{aligned}$$

(5.17)

using Friedrichs’ inequality on $\omega _z^\ell$ with ${\text {diam}}(\omega _z^\ell ) \le C_\textrm{p}\ell H$, $C_\textrm{p}>0$. For the second term in (5.16), we infer the trivial estimate

$$\begin{aligned} \Vert \nabla e_z^1 \Vert _{L^2(\omega _z)}^2 \le \Vert \nabla e_z^1 \Vert _{L^2(\omega _z^\ell )}^2, \end{aligned}$$

which implies that, in order to bound (5.16), it suffices to estimate $\Vert \nabla e_z^1 \Vert _{L^2(\omega _z^\ell )}$. Using the continuity estimate

$$\begin{aligned} \Vert {\text {tr}}_{\omega ^m_T}^{-1}{\text {tr}}_{\omega ^m_T} v\Vert _{H^1(\omega _T^m)} \le C_\textrm{tr} \Vert v\Vert _{H^1(\omega _T^m)},\qquad v \in \mathcal {V}|_{\omega _T^m} \end{aligned}$$

with a constant $C_\textrm{tr}>0$ independent of H, $\ell$ from the proof of [19, Thm. 6.1], one can show that

$$\begin{aligned} {( g_{T,j}^{m}\,,\,{\text {tr}}_{\omega ^m_T}^{-1}{\text {tr}}_{\omega ^m_T} e_z^1 )}_{L^2(\omega ^m_T)} \le C_\textrm{tr} \sigma (H,m) \Vert e_z^1\Vert _{H^1(\omega ^m_T)}. \end{aligned}$$

(5.18)

By (5.3) and (5.18), as well as the discrete Cauchy–Schwarz inequality, the finite overlap of the patches $\omega ^m_T$ and Friedrichs’ inequality, we get

$$\begin{aligned} \alpha \Vert \nabla e_z^1 \Vert _{L^2(\omega _z^\ell )}^2&\le \sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}} c_{T,j} a(\psi _{T,j}^{m}-{\widehat{\psi }}_{T,j}^{m},e_z^1) \\ {}&= -\sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}} c_{T,j} {( g_{T,j}^{m}\,,\,{\text {tr}}_{\omega ^m_T}^{-1}{\text {tr}}_{\omega ^m_T} e_z^1 )}_{L^2(\omega ^m_T)}\\&\le C_\textrm{tr}\sigma (H,m) \sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}} |c_{T,j}\vert \Vert e_z^1\Vert _{H^1(\omega ^m_T)}\\&\le C C_\textrm{tr} (C_\textrm{p}\ell H \pi ^{-1} \hspace{-2.22214pt}+ 1)\sigma (H,m) {m}^{d/2}J^{1/2}\sqrt{\sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}} c_{T,j}^2}\; \Vert \nabla e_z^1 \Vert _{L^2(\omega _z^\ell )} \end{aligned}$$

with constant $C>0$ reflecting the overlap of the patches $\omega _T^m$. Using (5.17), this yields an estimate for (5.16) and consequently bounds the first term in (5.15).

For the second expression in (5.15), using the definition of the Kolmogorov n-width (4.2), Friedrichs’ inequality on the patch $\omega _z^\ell$, the discrete Cauchy–Schwarz inequality, and that the $g_{T,j}^{m}$ are $L^2$-normalized, we obtain

$$\begin{aligned} \Vert \Lambda _ze_z^2 \Vert _{a,\omega _z}&\le d_n(H,\ell )\Vert v_z \Vert _{a,\omega _z^\ell } = d_n(H,\ell )\bigg \Vert {\mathcal {A}}_{\omega _z^\ell }^{-1}\sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}} c_{T,j} g_{T,j}^{m}\bigg \Vert _{a,\omega _z^\ell }\\&\le C_\textrm{p}d_n(H,\ell )\alpha ^{-1/2} \pi ^{-1} \ell H \bigg \Vert \sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}}c_{T,j} g_{T,j}^{m}\bigg \Vert _{L^2(\omega _z^\ell )}\\&\le C C_\textrm{p}d_n(H,\ell )\alpha ^{-1/2} \pi ^{-1} \ell H m^{d/2}J^{1/2} \sqrt{\sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}}c_{T,j}^2} \end{aligned}$$

with constant $C>0$ appearing in the bound $C^2 m^dJ$ of the number of terms in the above sum. Consequently, we conclude the estimate for the second term in (5.15).

The assertion can be finalized by combining all estimates utilizing

$$\begin{aligned} \sum _{z \in {\mathcal {N}}_H} \sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}}c_{T,j}^2 \le C m^d \sum _{\begin{array}{c} T\in \mathcal {T}_H\\ j=1,\dots ,J \end{array}} c_{T,j}^2 \le C C_\textrm{r}(H,\ell ) m^d {\left\Vert f \right\Vert }_{L^2(\Omega )}^2 \end{aligned}$$

with constant $C>0$ reflecting the overlap of the patches $\omega _z^m$. Here, we used Assumption 5.2 and (3.1) and (5.13). Finally, for the sake of readability, we substitute $m= \lfloor \tfrac{\ell }{2}\rfloor$ by $\ell$, which may introduce additional constants that change the decay rate of $\sigma$ by some constant factor. $\square$

5.3 Local approximation error bound using SLOD

The error estimate from Theorem 5.3 incorporates the spectral approximation error $d_n$ defined in (4.4). Subsequently, we derive a bound for $d_n$ utilizing the SLOD for constructing bases of the local spaces $\mathcal {V}_{H,z}^{\ell }$ defined in (4.1). For this purpose, we fix a node $z \in \mathcal {N}_{H}$ and treat $\omega _z^ \ell$ as the whole domain. For the oversampling parameters $m = 1,\dots ,\ell -1$, we denote the basis functions with a tilde to emphasize that, in general, they do not coincide with their global counterparts $\psi _{T,j}^{m}$, i.e.,

$$\begin{aligned} \{{{\tilde{\varphi }}}_{T,j}^{m}:={\mathcal {A}}_{\omega _z^\ell }^{-1} {{\tilde{g}}}_{T,j}^m\,:\,T \subset \omega _z^\ell ,\,j = 1,\dots ,J\} \end{aligned}$$

(5.19)

with source terms ${{\tilde{g}}}_{T,j}^m\in \mathcal {P}({{\tilde{\omega }}}_T^m)$, where $\tilde{\omega }_T^m :=\omega _T^m \cap \omega _z^\ell$. We denote the corresponding localized basis by

$$\begin{aligned} \{{{\tilde{\psi }}}_{T,j}^{m}:={\mathcal {A}}_{{\tilde{\omega }}^m_T}^{-1} {{\tilde{g}}}_{T,j}^m\,:\,T \subset \omega _z^\ell ,\,j = 1,\dots ,J\}. \end{aligned}$$

(5.20)

Similarly, as in the full domain setting, we need to measure, for all $T \subset \omega _z^\ell$, the (quasi-) orthogonality of the source terms $\{{{\tilde{g}}}_{T,j}^m\,:\,j = 1,\dots ,J\}$ on the corresponding space of ${\mathcal {A}}$-harmonic functions

$$\begin{aligned} {\text {tr}}^{-1}_{{\tilde{\omega }}^m_T}{\text {tr}}_{{\tilde{\omega }}^m_T}\big (H^1_0(\omega _z^\ell )|_{{\tilde{\omega }}^m_T}\big ). \end{aligned}$$

(5.21)

Similarly to (5.5), for a node $z \in \mathcal {N}_{H}$ and element $T \subset \omega _z^\ell$, we denote the singular values of $\Pi _H$ restricted to (5.21) by ${{\tilde{\sigma }}}_1\ge \tilde{\sigma }_2\ge \dots \ge {{\tilde{\sigma }}}_K\ge 0$. Analogously to (5.7), we define

$$\begin{aligned} {{\tilde{\sigma }}}(H,\ell ,m) :=\max _{z\in \mathcal {N}_{H}}\max _{T\in \mathcal {T}_H} {{\tilde{\sigma }}}_{z,T},\qquad {{\tilde{\sigma }}}_{z,T} :=\tilde{\sigma }_{K-J+1}. \end{aligned}$$

The quantity ${{\tilde{\sigma }}}$ is strongly related to its counterpart $\sigma$ from (5.7) and, in numerical experiments, exhibits the same qualitative behavior as described in Remark 5.1. A local variant of Assumption 5.2 is required to ensure the stability of the local basis.

Assumption 5.4

(Local Riesz stability) For all patches $\omega _z^\ell$, the set

$$\begin{aligned} \{{{\tilde{g}}}_{T,j}^m\,:\,T\subset \omega _z^\ell , j = 1,\dots ,J\} \end{aligned}$$

is a Riesz basis of $\mathcal {P}(\omega _z^\ell )$, i.e., there is ${{\tilde{C}}}_{\textrm{r}}(H,\ell ,m)>0$ depending polynomially on $H,\ell ,$ and $m$ such that, for all $z \in \mathcal {N}_{H}$ and all $(c_{T,j})_{T\subset \omega _z^\ell ,j=1,\dots ,J}$,

$$\begin{aligned} {{\tilde{C}}}^{-1}_{\textrm{r}}(H,\ell ,m)\sum _{\begin{array}{c} T\subset \omega _z^\ell \\ j=1,\dots ,J \end{array}} c_{T,j}^2 \le \bigg \Vert \sum _{\begin{array}{c} T\subset \omega _z^\ell \\ j=1,\dots ,J \end{array}} c_{T,j} {{\tilde{g}}}_{T,j}^m\biggl \Vert _{L^2(\omega _z^\ell )}^2 \;\le {{\tilde{C}}}_{\textrm{r}}(H,\ell ,m) \sum _{\begin{array}{c} T\subset \omega _z^\ell \\ j=1,\dots ,J \end{array}} c_{T,j}^2.\qquad \end{aligned}$$

(5.22)

Theorem 5.5

(Bound on $d_n$) Let Assumption 5.4 be satisfied. Then, there exists a constant $C>0$ independent of $H,\ell ,$ and m such that, for $m= 1,\dots ,\ell -1$

$$\begin{aligned} d_n(H,\ell ) \le C\ell \, m^{d/2} H^{-1} {{\tilde{C}}}^{1/2}_\textrm{r}(H,\ell ,m) {{\tilde{\sigma }}}(H,\ell ,m), \end{aligned}$$

where $n \approx m^d$.

Proof

Let us consider a fixed node $z \in \mathcal {N}_{H}$ and oversampling parameter $\ell$. As approximation space Q(n) of dimension $n \approx m^d$, we choose

$$\begin{aligned} Q(n) :=\textrm{span}\{\Lambda _z {{\tilde{\psi }}}_{T,j}^{m}\,:\,T \subset \omega _z^m, j = 1,\dots ,J\} \end{aligned}$$

with basis functions defined in (5.20). For the approximation of $v_z \in \mathcal {V}_{H,z}^{\ell }$, we choose the element $w_z \in Q(n)$ as

$$\begin{aligned} w_z = \Lambda _z \sum _{\begin{array}{c} T\subset \omega _z^m\\ j=1,\dots ,J \end{array}} c_{T,j}{{\tilde{\psi }}}_{T,j}^{m}, \end{aligned}$$

where the $c_{T,j}$ are the coefficients of the expansion of $v_z$ in terms of the basis functions ${{\tilde{\varphi }}}_{T,j}^{m}$ defined in (5.19). Note that, by Assumption 5.4, the coefficients $c_{T,j}$ are uniquely determined. Thus, we can estimate the spectral approximation error (4.2) using $v_z$ and $w_z$ as

$$\begin{aligned} d_n^z(H,\ell )&= \inf _{Q(n) \subset H^1_0(\omega _z)} \sup _{v_z \, \in \, \mathcal {V}_{H,z}^{\ell }} \inf _{w_z \in Q(n)} \frac{\Vert \Lambda _zv_z -w_z \Vert _{a,\omega _z}}{\Vert v_z \Vert _{a,\omega _z^\ell }} \\ {}&\le \sup _{v_z \, \in \, \mathcal {V}_{H,z}^{\ell }} \frac{\Vert \Lambda _zv_z -w_z \Vert _{a,\omega _z}}{\Vert v_z \Vert _{a,\omega _z^\ell }}. \end{aligned}$$

Denoting

$$\begin{aligned} e_z :=\sum _{\begin{array}{c} T\subset \omega _z^m\\ j=1,\dots ,J \end{array}} c_{T,j}\big ({{\tilde{\varphi }}}_{T,j}^{m}-{{\tilde{\psi }}}_{T,j}^{m}\big ) \in H^1_0(\omega _z^\ell ), \end{aligned}$$

we obtain for the numerator, using the product rule, the triangle inequality, and (3.7):

$$\begin{aligned} \Vert \Lambda _zv_z -w_z \Vert _{a,\omega _z}&\le \beta ^{1/2} {\left\Vert \nabla (\Lambda _z e_z) \right\Vert }_{L^2(\omega _z^\ell )} \\ {}&\le \beta ^{1/2}\big (C_\Lambda H^{-1} {\left\Vert e_z \right\Vert }_{L^2(\omega _z^\ell )} + {\left\Vert \nabla e_z \right\Vert }_{L^2(\omega _z^\ell )}\big ). \end{aligned}$$

We apply Friedrichs’ inequality on the patch $\omega _z^\ell$ using that ${\text {diam}}(\omega _z^\ell )\le C_\textrm{p}\ell H$, with a constant $C_\textrm{p}>0$. Hence, we can bound the first term against the second term, i.e.,

$$\begin{aligned} {\left\Vert e_z \right\Vert }_{L^2(\omega _z^\ell )} \le C_\textrm{p}H \ell \pi ^{-1} {\left\Vert \nabla e_z \right\Vert }_{L^2(\omega _z^\ell )}. \end{aligned}$$

We adapt estimate (5.18) to the local setting introducing a constant ${{\tilde{C}}}_\textrm{tr}>0$. Using the finite overlap of the patches ${\tilde{\omega }}^m_T$, the discrete Cauchy–Schwarz inequality and Friedrichs’ inequality on $\omega _z^\ell$, we obtain

$$\begin{aligned} \alpha \Vert \nabla e_z \Vert _{L^2(\omega _z^\ell )}^2&\le \sum _{\begin{array}{c} T\subset \omega _z^m\\ j=1,\dots ,J \end{array}} c_{T,j} a({{\tilde{\varphi }}}_{T,j}^{m}-{{\tilde{\psi }}}_{T,j}^{m},e_z) \\ {}&\le {{\tilde{C}}}_\textrm{tr} {{\tilde{\sigma }}}(H,\ell ,m) \sum _{\begin{array}{c} T\subset \omega _z^{m}\\ j=1,\dots ,J \end{array}}c_{T,j} \Vert e_z\Vert _{H^1({\tilde{\omega }}^m_T)}\\&\le C {{\tilde{C}}}_\textrm{tr} (C_\textrm{p}\pi ^{-1}H\ell + 1) m^{d/2} {{\tilde{\sigma }}}(H,\ell ,m) \sqrt{\sum _{\begin{array}{c} T\subset \omega _z^m\\ j=1,\dots ,J \end{array}} c_{T,j}^2}\Vert \nabla e_z\Vert _{L^2(\omega _z^\ell )}, \end{aligned}$$

where $C>0$ reflects the overlap of the patches ${\tilde{\omega }}^m_T$.

Adding the remaining coefficients $c_{T,j}$ from the expansion of $v_z$ in terms of the ${{\tilde{\varphi }}}_{T,j}^{m}$ and using Assumption 5.4, we get

$$\begin{aligned} {{\tilde{C}}}_\textrm{r}^{-1}(H,\ell ,m)\sum _{\begin{array}{c} T\subset \omega _z^\ell \\ j=1,\dots ,J \end{array}}c_{T,j}^2&\le \bigg \Vert \sum _{\begin{array}{c} T\subset \omega _z^\ell \\ j=1,\dots ,J \end{array}} c_{T,j} {{\tilde{g}}}_{T,j}^m\biggl \Vert _{L^2(\omega _z^\ell )}^2 \\ {}&\le C_\textrm{bo}^2 H^{-2} \bigg \Vert \sum _{\begin{array}{c} T\subset \omega _z^\ell \\ j=1,\dots ,J \end{array}} c_{T,j} {{\tilde{g}}}_{T,j}^m\biggl \Vert _{H^{-1}(\omega _z^\ell )}^2\\ {}&\le C_\textrm{bo}^2 \beta H^{-2} {\left\Vert v_z \right\Vert }_{a,\omega _z^\ell }^2. \end{aligned}$$

Here, we also employed that, by (3.5) and (3.6), we have, for all $q \in \mathcal {P}(\omega _z^\ell )$,

$$\begin{aligned} \Vert q \Vert _{L^2(\omega _z^\ell )}&= \sup _{v \, \in \, H^1_0(\omega _z^\ell )} \frac{{( q\,,\,v )}_{L^2(\omega _z^\ell )}}{\Vert v \Vert _{L^2(\omega _z^\ell )}} \le C_\textrm{bo}H^{-1}\sup _{v \, \in \, H^1_0(\omega _z^\ell )} \frac{{( q\,,\,{\mathcal {B}}_Hv )}_{L^2(\omega _z^\ell )}}{\Vert \nabla {\mathcal {B}}_Hv \Vert _{L^2(\omega _z^\ell )}}\\&\le C_\textrm{bo}H^{-1} \Vert q\Vert _{H^{-1}(\omega _z^\ell )}. \end{aligned}$$

Combining the estimate yields the assertion. $\square$

Remark 5.6

(Choice of parameters) This remark specifies how to choose the oversampling parameter $\ell$ and the number of local functions n in order to preserve the optimal order of convergence of $s+1$ in Theorem 5.3. For $\ell$, the super-exponential decay (5.8) implies that it needs to be chosen of order $|\log H\vert ^{(d-1)/d}$. Using Theorem 5.5 and that ${{\tilde{\sigma }}}$ has similar decay properties as $\sigma$, we obtain that n needs to be chosen of order $|\log H\vert ^{d-1}$. Note that these choices require the validity of (5.8) and Assumption 5.2 and 5.4. For a (pessimistic) estimate which is valid without any assumptions, we refer to Remark 6.3 below.

6 (Pessimistic) a priori error analysis

This section presents an a priori error analysis of the proposed method, which is based on the LOD framework, cf. [4, 10, 11]. Note that the exponential localization properties of the LOD cannot match the practically observed super-exponential localization properties of the SLOD, cf. Remark 5.1. Nevertheless, the LOD construction has the decisive advantage that the basis stability is guaranteed by construction and that the exponential localization can be rigorously proved. This enables an a priori analysis without assumptions on the stability of the SLOD basis and without conjectures on the decay of singular values, cf. Assumptions 5.2 and 5.4 and Remark 5.1.

6.1 Higher-order LOD

We briefly introduce a higher-order version of the LOD similar to the constructions in [25, 26]. The LOD constructs its problem-adapted basis functions by adding fine-scale information to coarse-scale finite element functions. We define the space of fine-scale functions as

$$\begin{aligned} \mathcal {W}:=\textrm{kernel}\,\Pi _H. \end{aligned}$$

(6.1)

The step of adding fine-scale information is called correction and utilizes the so-called correction operator $\mathcal {C}:\mathcal {V}\rightarrow \mathcal {W}$ defined as the a-orthogonal projection onto $\mathcal {W}$, i.e.,

$$\begin{aligned} a(\mathcal {C}v,w) = a(v,w),\qquad \text {for all }w\in \mathcal {W}. \end{aligned}$$

We split up the correction operator into a sum of element correction operators, i.e., $\mathcal {C}= \sum _{T \in \mathcal {T}_H} \mathcal {C}_T$ with element correction operators $\mathcal {C}_T:\mathcal {V}\rightarrow \mathcal {W}$ defined by

$$\begin{aligned} a(\mathcal {C}_Tv,w) = a_T(v,w),\qquad \text {for all }w\in \mathcal {W}. \end{aligned}$$

For any $v \in \mathcal {V}$, the correction $\mathcal {C}_Tv$ decays exponentially fast away from its associated element $T \in \mathcal {T}_H$, cf. [26, Lemma 5.1], which motivates a localization. For this purpose, we substitute the global space $\mathcal {W}$ by localized counterparts $\mathcal {W}_T^\ell :=\mathcal {W}(\omega _T^\ell )$, where, for a subdomain $\omega \subset \Omega$, we use the definition

$$\begin{aligned} \mathcal {W}(\omega ) :=\{w \in \mathcal {W}\,:\,{\text {supp}}(w)\subset \omega \}. \end{aligned}$$

(6.2)

The localized element correction operator $\mathcal {C}_T^\ell :\mathcal {V}\rightarrow \mathcal {W}_T^{\ell }$ is then defined by

$$\begin{aligned} a(\mathcal {C}_T^\ell v,w) = a_T(v,w),\qquad \text {for all }w\in \mathcal {W}_T^\ell . \end{aligned}$$

Similarly as for the correction operator $\mathcal {C}$ which can be decomposed into a sum of element correction operators, we define the localized correction operator by

$$\begin{aligned} \mathcal {C}^\ell :=\sum _{T \in \mathcal {T}_H} \mathcal {C}_T^\ell . \end{aligned}$$

The ansatz space of the LOD is then obtained by adding (localized) corrections to the bubble functions $\{b_{T,j}\,:\,T\in \mathcal {T}_H, j = 1,\dots ,J\}$ defined in (3.3), i.e.,

$$\begin{aligned} \mathcal {V}_{H}^{\ell ,\,\textrm{LOD}}:=\textrm{span}\{(1-\mathcal {C}^\ell ) b_{T,j}\,:\,T \in \mathcal {T}_H, j=1,\dots ,J\}. \end{aligned}$$

The Galerkin LOD approximation $u_{H}^{\ell ,\,\textrm{LOD}}\in \mathcal {V}_{H}^{\ell ,\,\textrm{LOD}}$ then satisfies

$$\begin{aligned} a(u_{H}^{\ell ,\,\textrm{LOD}},w) = {( f\,,\,w )}_{L^2(\Omega )},\qquad \text {for all }w \in \mathcal {V}_{H}^{\ell ,\,\textrm{LOD}}. \end{aligned}$$

(6.3)

6.2 A priori error bound using LOD

Using LOD techniques, we can prove the following a priori error estimate for the proposed multi-scale method. Numerical experiments show that this estimate is tentatively pessimistic. Nevertheless, compared to Theorem 5.3, it has the crucial advantage that it does not rely on additional assumptions or conjectures.

Theorem 6.1

(A priori error bound)  Let u and $u_{H}^{\ell ,n}$ denote the solutions of (2.2) and (4.5), respectively. There exist constants $C,C_\textrm{d}>0$ independent of $H,\ell ,$ and n such that, for any $f \in H^k(\mathcal {T}_H)$, $k \in {{\mathbb {N}}}_0$

$$\begin{aligned} \Vert u&-u_{H}^{\ell ,n} \Vert _{a,\Omega }\\ {}&\le C\Big ( H^{s+1} |f\vert _{H^s(\mathcal {T}_H)} + \ell ^{d/2}H^{-1} \big (\ell ^{d/2}\exp (-C_\textrm{d}\ell ) + d_n(H,\ell )\big )\Vert f \Vert _{L^2(\Omega )}\Big ). \end{aligned}$$

with $s :=\min \{k,p+1\}$.

Proof

We apply Céa’s Lemma, which yields, for arbitrary $v \in \mathcal {V}_{H}^{\ell ,n}$,

$$\begin{aligned} \Vert u-u_{H}^{\ell ,n} \Vert _{a,\Omega } \le \Vert u-v \Vert _{a,\Omega }. \end{aligned}$$

(6.4)

For the oversampling parameter $m:=\lfloor \ell /2\rfloor$, we define the approximation

$$\begin{aligned} u_{H}^{m,\,\textrm{LOD}}:=(1-\mathcal {C}^m){\mathcal {B}}_Hu \end{aligned}$$

which is not the Galerkin LOD solution (6.3) but has the same approximation properties, cf. proof of [25, Thm. 4.4]. Adding and subtracting $u_{H}^{m,\,\textrm{LOD}}$ in (6.4) and using the triangle inequality yields

$$\begin{aligned} \Vert u-u_{H}^{\ell ,n} \Vert _{a,\Omega } \le \Vert u-u_{H}^{m,\,\textrm{LOD}} \Vert _{a,\Omega } + \Vert u_{H}^{m,\,\textrm{LOD}}-v \Vert _{a,\Omega }. \end{aligned}$$

(6.5)

Using the above-mentioned approximation properties of $u_{H}^{m,\,\textrm{LOD}}$, we obtain the following estimate for the first term in (6.5)

$$\begin{aligned} \Vert u-u_{H}^{m,\,\textrm{LOD}} \Vert _{a,\Omega }\le C_\textrm{lod}\big ( H^{s+1} |f\vert _{H^s(\mathcal {T}_H)} + m^{d/2} H^{-1}\exp (-C_\textrm{d}m)\Vert f\Vert _{L^2(\Omega )}\big ). \end{aligned}$$

with constants $C_\textrm{lod},C_\textrm{d}>0$ independent of $H,\ell ,$ and m.

For the second term in (6.5), we choose $v \in \mathcal {V}_{H}^{\ell ,n}$ as sum of functions $v_z^n \in \mathcal {V}_{H,z}^{\ell ,n}$ to be specified later, i.e.,

$$\begin{aligned} v = \sum _{z \in \mathcal {N}_{H}} \Lambda _z v_z^n. \end{aligned}$$

Using the partition of unity property of the hat-functions $\sum _{z \in \mathcal {N}_{H}}\Lambda _z \equiv 1$, we obtain for the second term in (6.5),

$$\begin{aligned} \Vert u_{H}^{m,\,\textrm{LOD}}-v \Vert _{a,\Omega }^2&= \Big \Vert \sum _{z \in {\mathcal {N}}_H} \Lambda _z(u_{H}^{m,\,\textrm{LOD}}- v_z^n)\Big \Vert _{a,\Omega }^2 \\ {}&\le C_\textrm{ol}\sum _{z \in {\mathcal {N}}_H}\Vert \Lambda _z(u_{H}^{m,\,\textrm{LOD}}- v_z^n) \Vert _{a,\omega _z}^2 \end{aligned}$$

with $C_\textrm{ol}$ defined in (3.8). For any $z \in \mathcal {N}_{H}$, we can, locally on $\omega _z$, substitute $u_{H}^{m,\,\textrm{LOD}}$ by

$$\begin{aligned} u_{H,z}^{m,\,\textrm{LOD}}:=(1-\mathcal {C}^m)({\mathcal {B}}_Hu|_{\omega _z^m}). \end{aligned}$$

Hence, we define $v_z^n\in \mathcal {V}_{H,z}^{\ell ,n}$ as the (not necessarily unique) elements minimizing the expression $\Vert \Lambda _z v_z -\Lambda _z v_z^n \Vert _{a,\omega _z^\ell }$, where

$$\begin{aligned} v_z :=(1-\tilde{\mathcal {C}}^\ell _z)({\mathcal {B}}_Hu|_{\omega _z^m}) \end{aligned}$$

is an approximation to $u_{H,z}^{m,\,\textrm{LOD}}$. We denote $\mathcal {W}_z^\ell :=\mathcal {W}(\omega _z^\ell )$ and define the above used correction operator $\tilde{\mathcal {C}}^\ell _z:H^1_0(\Omega )\rightarrow \mathcal {W}_z^\ell$ by

$$\begin{aligned} a_{\omega _z^\ell }( \tilde{\mathcal {C}}^\ell _zv,w) = a_{\omega _z^\ell }(v,w)\qquad \text {for all }w \in \mathcal {W}(\omega _z^\ell ) . \end{aligned}$$

(6.6)

Note that it holds $v_z \in \mathcal {V}_{H,z}^{\ell }$ which is a non-trivial observation, cf. [4, Rem 3.7 ii].

Abbreviating

$$\begin{aligned} e_z^1 :=u_{H,z}^{m,\,\textrm{LOD}}- v_z,\qquad e_z^2:=v_z - v_z^n, \end{aligned}$$

we obtain after performing the above-mentioned local substitution

$$\begin{aligned} \Vert \Lambda _z(u_{H}^{m,\,\textrm{LOD}}- v_z^n) \Vert _{a,\omega _z}^2 \le 2\big ( \Vert \Lambda _z e_z^1 \Vert _{a,\omega _z}^2+\Vert \Lambda _z e_z^2 \Vert _{a,\omega _z}^2\big ). \end{aligned}$$

(6.7)

Following (5.16), in order to bound the first term in (6.7), it suffices to bound $\Vert e_z^1 \Vert _{L^2(\omega _z)}$ and $\Vert \nabla e_z^1 \Vert _{L^2(\omega _z)}$. It holds

$$\begin{aligned} e_z^1 = (\mathcal {C}^m- \tilde{\mathcal {C}}^\ell _z) ({\mathcal {B}}_Hu|_{\omega _z^m}) \in \mathcal {W}_z^\ell \end{aligned}$$

which implies that $e_z^1$ has vanishing element averages. Thus, by Poincare’s inequality

$$\begin{aligned} \Vert e_z^1 \Vert _{L^2(\omega _z)}^2 \le 4\pi ^{-2}H^2 \Vert \nabla e_z^1 \Vert _{L^2(\omega _z)}^2 \end{aligned}$$

it is sufficient to estimate $\Vert \nabla e_z^2 \Vert _{L^2(\omega _z)}$ in order to obtain a bound for the first term in (6.7). Given a function supported in ${\omega _z^m}$ (e.g., ${\mathcal {B}}_Hu|_{\omega _z^m}$) the correction operator $\mathcal {C}^m$ coincides with the localization of $\tilde{\mathcal {C}}^\ell _z$ to $m$-th order patches. Hence, we can apply the localization error estimate from the proof of [25, Thm. 4.4], here in the oversampling parameter $m$ and (3.5) to obtain

$$\begin{aligned} \Vert \nabla e_z^1 \Vert _{L^2(\omega _z)}^2&\le \Vert \nabla (\mathcal {C}^m-\tilde{\mathcal {C}}^\ell _z)({\mathcal {B}}_Hu|_{\omega _z^m}) \Vert _{L^2(\omega _z^\ell )}^2 \\&\le C_\textrm{lo}^2 m^{d}\exp (-C_\textrm{d}m)^2 \Vert \nabla ({\mathcal {B}}_Hu|_{\omega _z^m}) \Vert _{L^2(\omega _z^\ell )}^2\\&\le C_\textrm{lo}^2 C_\textrm{bo}^2 H^{-2} m^{d} \exp (-C_\textrm{d}m)^2 \Vert u \Vert _{L^2(\omega _z^m)}^2. \end{aligned}$$

with constants $C_\textrm{lo},C_\textrm{d}>0$ independent of $H,\ell ,$ and m.

For the second term in (6.7), we obtain by the definition of the Kolmogorov n-width, the stability of $(1-\tilde{\mathcal {C}}^\ell _z)$ and (3.5) that

$$\begin{aligned} \Vert \Lambda _z e_z^2 \Vert _{a,\omega _z}^2&\le d_n^2(H,\ell ,n)\Vert v_z \Vert _{a,\omega _z^\ell }^2 \le \beta d_n^2(H,\ell ,n)\Vert \nabla ({\mathcal {B}}_Hu|_{\omega _z^m}) \Vert _{L^2(\omega _z^\ell )}^2\\*&\le C_\textrm{bo}^2 H^{-2}\beta ^2 d_n^2(H,\ell ,n) \Vert u \Vert _{L^2(\omega _z^m)}^2. \end{aligned}$$

Combining the previous estimates, using Friedrichs’ inequality on $\Omega$ (recall that $\Omega$ is scaled to unit size), we get

$$\begin{aligned} \Vert u \Vert _{L^2(\Omega )} \le \pi ^{-1}\Vert \nabla u \Vert _{L^2(\Omega )} \le \pi ^{-2}\alpha ^{-1}\Vert f \Vert _{L^2(\Omega )}. \end{aligned}$$

Finally, we substitute $m = \lfloor \tfrac{\ell }{2}\rfloor$ by $\ell$, which introduces additional constants and changes the exponential decay rate $C_\textrm{d}$ by a factor of two. $\square$

6.3 Local approximation error bound using LOD

Subsequently, we derive an a priori bound for $d_n$, which is fully explicit in H and $\ell$. This is the analog to Theorem 5.5, which does not rely on additional assumptions or conjectures. Numerical experiments show that this estimate is tentatively pessimistic.

Theorem 6.2

(Bound on $d_n$) There exists $C,C_\textrm{d}>0$ independent of $H,\ell ,$ and m such that, for $m= 1,\dots ,\ell -1$

$$\begin{aligned} d_n(H,\ell ) \le C\ell \, m^{d/2} \exp (-C_\textrm{d}m), \end{aligned}$$

where $n \approx m^d$.

Proof

Let us consider a fixed node $z \in \mathcal {N}_{H}$ and oversampling parameter $\ell$. By [4, Rem. 3.7 ii], we can write any $v \in \mathcal {V}_{H,z}^{\ell }$ as

$$\begin{aligned} v = (1-\tilde{\mathcal {C}}^\ell _z){\mathcal {B}}_Hv \end{aligned}$$

with the correction operator $\tilde{\mathcal {C}}^\ell _z$ defined in (6.6). We define the patch-local localized correction operator by ${\tilde{\mathcal {C}}}^m:=\sum _{T \subset \omega _z^\ell } {\tilde{\mathcal {C}}}_T^m$, where, denoting ${{\tilde{\mathcal {W}}}}_T^m :=\mathcal {W}( \omega _T^m\cap \omega _z^\ell )$, the element correctors ${\tilde{\mathcal {C}}}_T^m:H^1_0(\omega _z^\ell )\rightarrow {{\tilde{\mathcal {W}}}}_T^m$ are defined by

$$\begin{aligned} a_{\omega _z^\ell }({{\tilde{\mathcal {C}}}}_T^mv,w) = a_T(v,w),\qquad \text {for all }w \in {{\tilde{\mathcal {W}}}}_T^m. \end{aligned}$$

As approximation space Q(n) of dimension $n \approx m^d$, we choose

$$\begin{aligned} Q(n) :=\textrm{span}\{\Lambda _z (1-{{\tilde{\mathcal {C}}}}^m)b_{T,j}\,:\,T \subset \omega _z^m, j = 1,\dots ,J\} \end{aligned}$$

and as approximation $w_z \in Q(n)$ of an element $v_z\in \mathcal {V}_{H,z}^{\ell }$, we use

$$\begin{aligned} w_z = \Lambda _z (1-{{\tilde{\mathcal {C}}}}^m)({\mathcal {B}}_Hv_z|_{\omega _z^m}). \end{aligned}$$

Using the approximation space Q(n) and the above defined choice of $w_z$, we can bound the Kolmogorov n-width as follows

$$\begin{aligned} d_n^z(H,\ell )&= \inf _{Q(n) \subset H^1_0(\omega _z)} \sup _{v_z \, \in \, \mathcal {V}_{H,z}^{\ell }} \inf _{w_z \in Q(n)} \frac{\Vert \Lambda _z v_z -w_z \Vert _{a,\omega _z}}{\Vert v_z \Vert _{a,\omega _z^\ell }} \\ {}&\le \sup _{v_z \, \in \, \mathcal {V}_{H,z}^{\ell }} \frac{\Vert \Lambda _z v_z -w_z \Vert _{a,\omega _z}}{\Vert v_z \Vert _{a,\omega _z^\ell }}. \end{aligned}$$

Abbreviating

$$\begin{aligned} e_z :=({{\tilde{\mathcal {C}}}}^m-\tilde{\mathcal {C}}^\ell _z) {\mathcal {B}}_Hv_z, \end{aligned}$$

we can estimate the numerator using (3.7) and

$$\begin{aligned} \Lambda _z(1-{{\tilde{\mathcal {C}}}}^m)({\mathcal {B}}_Hv_z|_{\omega _z^m}) = \Lambda _z(1-{{\tilde{\mathcal {C}}}}^m){\mathcal {B}}_Hv_z, \end{aligned}$$

as

$$\begin{aligned} \Vert \Lambda _zv_z -w_z \Vert _{a,\omega _z}&\le \beta ^{1/2}\Vert \nabla \Lambda _z e_z \Vert _{L^2(\omega _z)} \\ {}&\le \beta ^{1/2}\big (C_\Lambda H^{-1} \Vert e_z \Vert _{L^2(\omega _z)} + \Vert \nabla e_z \Vert _{L^2(\omega _z)}\big ). \end{aligned}$$

It holds that $e_z \in \mathcal {W}_z^\ell$, which implies that $e_z^1$ has vanishing element averages. Thus, by Poincare’s inequality

$$\begin{aligned} \Vert e_z \Vert _{L^2(\omega _z)} \le 2\pi ^{-1}H\Vert \nabla e_z \Vert _{L^2(\omega _z)}. \end{aligned}$$

Applying the localization error estimate from the proof of [25, Tehorem 4.4] to show that ${{\tilde{\mathcal {C}}}}^m$ approximates $\tilde{\mathcal {C}}^\ell _z$ exponentially and using (3.5), and Friedrichs’ inequality on the patch $\omega _z^\ell$ with ${\text {diam}}(\omega _z^\ell ) \le C_\textrm{p}\ell H$, $C_\textrm{p}>0$, we obtain

$$\begin{aligned} \Vert \Lambda _zv_z&-w_z \Vert _{a,\omega _z} \\ {}&\le \beta ^{1/2} (2C_\Lambda \pi ^{-1}+1)\Vert \nabla e_z \Vert _{L^2(\omega _z^\ell )}\\&\le C_\textrm{lo}m^{d/2} \beta ^{1/2}(2C_\Lambda \pi ^{-1}+1)\exp (-C_\textrm{d}m) \Vert \nabla {\mathcal {B}}_Hv \Vert _{L^2(\omega _z^\ell )}\\&\le C_\textrm{bo}C_\textrm{lo}C_\textrm{p}\ell \, m^{d/2} \beta ^{1/2}(2C_\Lambda \pi ^{-1}+1)\alpha ^{-1/2}\pi ^{-1}\exp (-C_\textrm{d}m) \Vert v_z \Vert _{a,\omega _z^\ell }, \end{aligned}$$

with constants $C_\textrm{lo},C_\textrm{d}>0$ independent of $H,\ell$, and m. The assertion follows immediately. $\square$

Remark 6.3

(Choice of parameters) This remark specifies how to choose the oversampling parameter $\ell$ and the number of local functions n in order to guarantee the optimal order of convergence of $s+1$ in Theorem 6.1. By Theorem 6.1, $\ell$ needs to be chosen of order $|\log H\vert$. Using Theorem 6.2, we obtain that n needs to be chosen of order $|\log H\vert ^{d}$. According to the experiments, these choices are pessimistic, cf. Remark 5.6.

7 Implementation and numerical experiments

In this section, we numerically investigate the proposed multi-scale method regarding the localization error, optimal convergence properties, high-contrast channeled coefficients, and higher-order polynomials using suitable benchmark problems. As a comparison, we use the SLOD from [19], which we consider as state-of-the-art. We refer to [19, Sec. 8] for a comparison of the SLOD to other multi-scale methods such as the LOD. For underlining the origin of the proposed method and its super-localization properties, cf. Theorems 5.3 and 5.5, we subsequently refer to it as Super-Localized Generalized Finite Element Method (SL-GFEM).

7.1 Implementation

For the practical implementation of the SL-GFEM, we need to perform a fine-scale discretization, i.e., we substitute the infinite-dimensional function space $\mathcal {V}$ by the finite element space $\mathcal {V}_h :=\mathcal {V}\cap \mathcal {P}_1(\mathcal {T}_h)$. Here, $\mathcal {T}_h$ denotes a fine-scale mesh obtained by uniform refinement of $\mathcal {T}_H$, where the number of refinements should be chosen such that the resulting mesh resolves all oscillations of A and f. For solving the patch problems (4.1) and (4.3), one considers patch-local subspaces of $\mathcal {V}_h$.

The SL-GFEM is straightforward to implement, as only very few technical details need to be addressed. The local spaces $\mathcal {V}_{H,z,h}^\ell$ (discrete counterparts of $\mathcal {V}_{H,z}^{\ell }$) can be computed in parallel. Their computation only requires the local stiffness and mass matrices on the respective patches.

In contrast to, e.g., MS-GFEM methods [5,6,7], the SL-GFEM solves local eigenvalue problems which are posed in the space spanned by a small number of deterministic snapshots. This results in a lower dimension of the eigenvalue problems (4.3) and, hence, makes them easier to solve numerically. After a multiplication with the respective hat-functions (partition of unity functions), we store the eigenfunctions corresponding to the n largest eigenvalues of the eigenvalue problems (4.3). These functions are then used as ansatz functions for computing the global approximation (4.5). Compared to the SLOD, by construction, no stability issues in the choice of basis can occur for the SL-GFEM and thus, no special treatment of boundary patches is required, cf. [19, App. B].

For the implementation, we use gridlod [27], which is a Python-library initially designed for the implementation of LOD-related methods. Although we do not require particular LOD-functionality from gridlod, it is convenient to use its flexible data structures for patches and its local discretization tools. Similarly, as in [31], our implementation can solve all local patch problems in parallel on an HPC cluster. As a comparison, we also implemented the SLOD from [19] in gridlod. All experiments are fully reproducible, and the corresponding source code can be found in [32].¹

7.2 Numerical experiments

We consider the domain $\Omega = (0,1)^2$ equipped with coarse Cartesian meshes $\mathcal {T}_H$ and a fine Cartesian mesh $\mathcal {T}_h$ obtained by uniform refinement of $\mathcal {T}_H$. Note that, for ease of presentation, H and h henceforth denote the elements side lengths instead of their diameters. For all numerical experiments, we use $h=2^{-10}$, which results in about one million degrees of freedom for the fine mesh. Note that our implementation also works for higher spatial dimensions d. For the numerical experiments, we consider two scalar diffusion coefficients A (realization of random field with short correlation length and high contrast channeled coefficient) and two source terms f (constant and non-polynomial). The precise definitions can be found in the respective experiments. Each configuration serves its own purpose for numerically investigating the SL-GFEM. For measuring the approximation quality, we use the relative energy error, i.e.,

$$\begin{aligned} e_{a}^{\text {rel}}({\tilde{u}}) :=\frac{\Vert u_h - {\tilde{u}} \Vert _{a,\Omega }}{\Vert u_h \Vert _{a,\Omega }}, \end{aligned}$$

where $u_h \in \mathcal {V}_h$ denotes the first order finite element approximation of (2.2) which we use as reference solution. Further, ${\tilde{u}}$ is a placeholder for the SL-GFEM approximation (4.5) or the SLOD approximation (5.10).

7.2.1 Super-exponential localization

First, we investigate the localization properties of the SL-GFEM given several choices of the local approximation space size n. For the choice $f \equiv 1$, the optimal order term in Theorems 5.3 and 6.1 disappears and only the localization error and the approximation error $d_n$ are present. As coefficient A, we consider a realization of the random field taking piecewise constant values on $\mathcal {T}_{2^{-8}}$, which are independent and identically distributed in the interval [1, 100]. This results in a maximum contrast of 100. We consider the fixed coarse mesh $\mathcal {T}_{2^{-5}}$ and the polynomial degree $p=0$.

Fig. 2

Localization errors of the SL-GFEM for multiple choices of n and of the SLOD for a fixed coarse mesh

For several sizes of local approximation spaces n, Fig. 2 depicts the relative energy errors of the SL-GFEM and the SLOD as a function of the oversampling parameter $\ell$. Clearly, the parameter n strongly impacts the approximation error of the SL-GFEM. One observes a large difference in the approximation quality, for instance, for $n=10$ and $n=15$. Conversely, choosing $n=20$ does not yield a significantly better approximation than $n=15$. This effect is related to the large jumps between the plateaus in Fig. 1 and is also visible Theorems 5.5 and 6.2. In conclusion, for a sufficiently large n, the SL-GFEM shows a rapid decay of the localization error confirming Theorems 5.3 and 6.1 numerically. For the SLOD, the super-localization property [19, Sec. 7] is visible in Fig. 2.

7.2.2 Optimal convergence

For investigating the convergence with respect to the coarse mesh size H, we use the same coefficient as in Sect. 7.2.1 but consider the non-polynomial source term

$$\begin{aligned} f_1(x_1,x_2) :=(x_1 + \cos (3 \pi x_1)) \cdot x_2^3\in H^1(\Omega ). \end{aligned}$$

Figure 3 depicts the errors of the SL-GFEM and SLOD for multiple choices of n and $\ell$ as a function of H. As a reference, a line with slope two indicates the expected order of convergence.

Fig. 3

Convergence plot of the SL-GFEM and SLOD for multiple choices of n and $\ell$

For n and $\ell$ sufficiently large, one observes that the SL-GFEM converges with an order of two which numerically confirms Theorems 5.3 and 6.1. Notably, the errors of the SL-GFEM are smaller by nearly one order of magnitude than the errors of the SLOD. This effect only appears for non-trivial coefficients A, i.e., the effect is most probably related to the contrast of the coefficient. The contrast dependence is investigated more closely in the following subsection.

7.2.3 High-contrast channeled coefficient

One of the major challenges for multi-scale methods is their sensitivity to high-contrast channeled coefficients. In this numerical experiment, we consider the coefficient $A_{\kappa }$ constructed by adding four channels of conductivity $\kappa$ in the coefficient from Sects. 7.2.1 and 7.2.2. Some of the channels touch the boundary, while others stop before. The number $\kappa$ is the maximum contrast of $A_\kappa$. Figure 4 depicts the coefficient for $\kappa = 10^5, 10^7$.

Fig. 4

Coefficient $A_{\kappa }$ for $\kappa =10^4, 10^7$ (left and right)

For this numerical experiment, we choose the same setup as in Sect. 7.2.1. Figure 5 depicts the localization errors for the above choices of $\kappa$.

Fig. 5

Localization errors of the SL-GFEM and the SLOD for the high-contrast channeled coefficient $A_\kappa$ for $\kappa = 10^4, 10^7$ (left and right)

The SL-GFEM appears to be largely unaffected by large values of $\kappa$. For the SLOD, in contrast, the best practical realization known until now [19] yields a basis with deteriorating stability as $\kappa$ is increased (growing constants in Assumptions 5.2 and 5.4). This explains the worse performance of the SLOD for $\kappa = 10^7$ compared to $\kappa = 10^4$. Notably, when compared to Sect. 7.2.1, the SL-GFEM does not need more local functions n to attain a good approximation quality, which suggests that the choice of n is not affected by the contrast.

7.2.4 Higher-order polynomials

One key benefit of the proposed method is its flexibility with regard to the choice of polynomial degree, i.e., the construction of higher-order methods is straightforward. While the previous numerical experiments have investigated the performance of the SL-GFEM for $p = 0$, this experiment also considers higher polynomial degrees. Using the setup from Sect. 7.2.2, Fig. 6 depicts the errors of the SL-GFEM for $p = 0,1,2$ as a function of H together with lines indicating the respective expected orders of convergence.

Fig. 6

Convergence plot of the higher-order SL-GFEM for multiple choices of n, $\ell$, and p

For n and $\ell$ sufficiently large, it can be observed that the method of degree p converges with an order of $p+2$ (recall that f is sufficiently smooth). This numerically confirms Theorems 5.3 and 6.1. Note that the choice of n needs to be adapted to p, which is related to the larger plateaus in Fig. 1. We observed that n needs to be increased linearly as p is increased. It is left to future research to find a (possibly adaptive) choice of n such that pessimistic choices (that may result in unnecessarily many basis functions) can be avoided.