# An improved high-order method for elliptic multiscale problems

Zhaonan Dong, Moritz Hauck, Roland Maier

► To cite this version:

Zhaonan Dong, Moritz Hauck, Roland Maier. An improved high-order method for elliptic multiscale problems. SIAM Journal on Numerical Analysis, 2023, 61 (4), pp.1918-1937. 10.1137/22M153392X. hal-03839664v2

HAL Id: hal-03839664

https://inria.hal.science/hal-03839664v2

Submitted on 31 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# AN IMPROVED HIGH-ORDER METHOD FOR ELLIPTIC MULTISCALE PROBLEMS

ZHAONAN DONG*, MORITZ HAUCK†, ROLAND MAIER‡

**ABSTRACT.** In this work, we propose a high-order multiscale method for an elliptic model problem with rough and possibly highly oscillatory coefficients. Convergence rates of higher order are obtained using the regularity of the right-hand side only. Hence, no restrictive assumptions on the coefficient, the domain, or the exact solution are required. In the spirit of the Localized Orthogonal Decomposition, the method constructs coarse problem-adapted ansatz spaces by solving auxiliary problems on local subdomains. More precisely, our approach is based on the strategy presented by Maier [SIAM J. Numer. Anal. 59(2), 2021]. The unique selling point of the proposed method is an improved localization strategy curing the effect of deteriorating errors with respect to the mesh size when the local subdomains are not large enough. We present a rigorous a priori error analysis and demonstrate the performance of the method in a series of numerical experiments.

**Keywords.** multiscale method, numerical homogenization, high-order method, localization

**AMS subject classifications.** 65N12, 65N30

## 1. INTRODUCTION

In this work, we consider the numerical solution of elliptic partial differential equations (PDEs) with possibly rough and highly varying coefficients. The numerical treatment of such problems with classical finite element approaches suffers from suboptimal approximation rates and pre-asymptotic effects if the discretization parameters do not resolve the oscillation scales. As globally resolving these scales can be computationally unfeasible, we aim to obtain reasonable errors already on coarse scales. Computational multiscale methods tackle this problem by constructing coarse problem-adapted ansatz spaces that contain local information on the underlying coefficient. Under minimal regularity assumptions on the coefficient, such methods achieve optimal orders of approximation while introducing only a moderate computational overhead (e.g., an increased support of the basis functions or an increased number of basis functions per mesh entity). Prominent first-order methods include, e.g., Generalized (Multiscale) Finite Element Methods [BO83, BCO94, BL11, EGH13], Adaptive Local Bases [GGS12, Wey16], the Localized Orthogonal Decomposition (LOD) method [MP14, HP13], or gamblets [Owh17]. For an overview of such methods, see also the recent textbooks [OS19, MP20] and the review article [AHP21].

Under additional smoothness assumptions, higher-order multiscale methods have been proposed based on, e.g., the Heterogeneous Multiscale Method [LMT12, AB12] or the Multiscale Finite Element Method [AB05, HZZ14]. Other strategies include Multiscale

<footnote>The work of Moritz Hauck is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 865751 - RandomMultiScales).</footnote>
Hybrid-Mixed Methods [HPV13, AHPV13] as well as the Multiscale Hybrid High-Order Method [CEL19]. However, for general $L^\infty$-coefficients, these methods only provide first-order convergence results in $H^1$. This issue is overcome in [Mai21] (see also [Mai20]), where a method is proposed that achieves high-order rates for arbitrarily rough coefficients. The idea is to exploit appropriate orthogonality properties in the spirit of the LOD which allows one to extract higher-order rates from the right-hand side only. Basis functions of the coarse problem-adapted ansatz space can be constructed as solutions to certain constrained energy minimization problems. The constraints enforce that the $L^2$-projection of a particular basis function into the space of piecewise polynomials coincides with a tensor-product Legendre polynomial of a given degree that is only supported on a single (quadrilateral) element. Since these basis functions are globally supported, a localization of the basis is required in order to obtain a practically feasible approach. In [Mai21], it is shown that the problem-adapted basis functions decay exponentially fast in modulus relative to the coarse mesh, which allows one to approximate these functions on localized subdomains. However, the approach theoretically and practically suffers from deteriorating errors as the mesh size is decreased when the local subdomains are not suitably increased.

In this work, we cure the above-mentioned error deterioration of the method in [Mai21] by proposing an improved localization technique. It slightly enlarges the support of the lowest-order basis functions, while keeping the other basis functions unchanged. Notably, this cures the issue of deteriorating errors for an arbitrary but fixed polynomial degree. In practice, this modification avoids an overestimation of the required size of the local subdomains and thereby reduces the computational costs and increases the reliability of the method.

This work is structured as follows. First, we present the elliptic model problem (Section 2) and preliminary notation (Section 3). A prototypical (non-localized) multiscale method is then presented in Section 4. In contrast to the definition in [Mai21], which is based on non-conforming polynomial spaces, our construction utilizes conforming bubble spaces. These spaces are designed such that the resulting prototypical multiscale method coincides with the one in [Mai21]. In Section 5, we propose an improved localization technique, which allows one to derive a practical version of the method (Section 6) which is able to remedy the issues of deteriorating errors for subdomains of fixed size. In Section 7, we finally show the performance of the novel strategy in a set of numerical experiments and, in particular, compare it to the original version.

## 2. MODEL PROBLEM

On a polygonal Lipschitz domain $D \subset \mathbb{R}^d$, $d \in \mathbb{N}$, we consider the prototypical second-order elliptic problem in weak form that seeks $u \in H_0^1(D)$ such that

$$ (2.1) \qquad a(u, v) := \int_D A \nabla u \cdot \nabla v \,dx = (f, v)_{L^2(D)} \quad \text{for all } v \in H_0^1(D) $$

with a right-hand side $f \in L^2(D)$ and a scalar-valued coefficient $A \in L^\infty(D)$ that fulfills $0 < \alpha \le A(x) \le \beta < \infty$ for almost all $x \in D$ and some positive constants $\alpha, \beta$. Apart from these bounds, we do not pose any further assumptions on $A$. However, we implicitly assume that the coefficient oscillates on multiple scales. It is well-known that in such a setting, classical finite element methods can perform arbitrarily badly [BO00], in particular if the
oscillations in the coefficient are not resolved. In the context of the model problem (2.1), this is also illustrated in [MP20, Ch. 2].

By the Lax–Milgram theorem, the solution $u$ exists and is uniquely defined. Further, the stability estimate

$$ (2.2) \qquad \|\nabla u\|_{L^2(D)} \le C_{\mathrm{PF}} \alpha^{-1} \|f\|_{L^2(D)} $$

holds, where the Poincaré–Friedrichs inequality with constant $C_{\mathrm{PF}} > 0$ is used to bound the $H^{-1}(D)$-norm of the right-hand side by its $L^2(D)$-norm.

### 3. PRELIMINARIES

Before we state a prototypical multiscale method in the following section, we first need to introduce coarse-scale spaces and interpolation operators that are required for the construction of the method.

3.1. **Non-conforming spaces.** Let $\mathcal{T}_H$ be a quasi-uniform and non-degenerate quadrilateral mesh with mesh size $H > 0$. We define for $\ell \in \mathbb{N}$ the element patch of order $\ell$ around $S \subset D$ by

$$ (3.1) \qquad \mathbf{N}^{\ell}(S) := \mathbf{N}^{1}(\mathbf{N}^{\ell-1}(S)), \quad \ell \ge 2, \qquad \mathbf{N}^{1}(S) := \bigcup \{T \in \mathcal{T}_H : \bar{S} \cap \bar{T} \neq \emptyset\} $$

and set $\mathbf{N}(S) := \mathbf{N}^{1}(S)$. For a fixed (but arbitrary) polynomial degree $p$, we denote with

$$ \mathcal{V}_H := \{v \in L^2(D) : v|_T, T \in \mathcal{T}_H, \text{ is a polynomial of coordinate degree } \le p\} $$

the non-conforming space (with respect to $H_0^1(D)$) consisting of element-wise defined polynomials. The restriction of $\mathcal{V}_H$ to a subdomain $S \subset D$ is given by

$$ \mathcal{V}_H(S) := \{v \in \mathcal{V}_H : \mathrm{supp}(v) \subset S\}. $$

As orthonormal basis $\mathcal{D}_T := \{\Lambda_{j,T}\}_{j=1}^N$ of the space $\mathcal{V}_H(T)$, $T \in \mathcal{T}_H$, we henceforth use the shifted tensor-product Legendre polynomials. Thus, $\mathcal{D} := \bigcup_{T \in \mathcal{T}_H} \mathcal{D}_T$ defines a local basis of $\mathcal{V}_H$. Denoting with $M := |\mathcal{T}_H|$ the number of elements in $\mathcal{T}_H$, this means that $\dim \mathcal{V}_H = MN$. Note that it holds $N = (p+1)^d$.

Let $\Pi_H: L^2(D) \to \mathcal{V}_H$ be the $L^2$-projection defined for any $v \in L^2(D)$ by the element-wise equation

$$ (\Pi_H v, v_H)_{L^2(T)} = (v, v_H)_{L^2(T)} \quad \text{for all } v_H \in \mathcal{V}_H(T), T \in \mathcal{T}_H. $$

From this definition, we directly obtain the following stability estimate,

$$ \|\Pi_H v\|_{L^2(T)} \le \|v\|_{L^2(T)}, \quad v \in L^2(T), T \in \mathcal{T}_H. $$

Additionally, there exists a constant $C_a > 0$ (independent of the mesh size) such that

$$ (3.2) \qquad \|(1 - \Pi_H)v\|_{L^2(T)} \le C_a H \|\nabla v\|_{L^2(T)}, \quad v \in H^1(T), T \in \mathcal{T}_H, $$

see, e.g., [HSS02]. We emphasize that the constant $C_a$ depends on the polynomial degree $p$. More precisely, $C_a \sim (p+1)^{-1}$.

For late use, we also define the broken Sobolev space $H^k(\mathcal{T}_H)$, $k \in \mathbb{N}$, by

$$ H^k(\mathcal{T}_H) := \{v \in L^2(D) : \forall T \in \mathcal{T}_H : v|_T \in H^k(T)\} $$

with the seminorm

$$ |\cdot|_{H^k(\mathcal{T}_H)} := \sum_{T \in \mathcal{T}_H} |\cdot|_{H^k(T)}. $$
3.2. **Conforming spaces.** While the locality of the non-conforming spaces is very beneficial for designing multiscale methods, it makes their theoretical investigation more cumbersome. Therefore, we replace the space $V_H$ by a conforming space $U_H \subset H_0^1(D)$ with the properties that $\Pi_H U_H = V_H$ and $\dim U_H = \dim V_H = MN$. The space $U_H$ is constructed based on a set of local basis bubble functions $\mathfrak{B}_T := \{b_{j,T}\}_{j=1}^N$ for each element $T \in \mathcal{T}_H$. These functions fulfill, for $1 \le j \le N$,

$$ (3.3) \qquad b_{j,T} \in H_0^1(T), \quad \Pi_H b_{j,T} = \Lambda_{j,T}. $$

Note that, throughout this paper, we do not distinguish between $H_0^1(S)$-functions on a subdomain $S \subset D$ and their $H_0^1(D)$-conforming extension by zero to the full domain $D$.

For the moment, an explicit characterization of the basis $\mathfrak{B} := \bigcup_{T \in \mathcal{T}_H} \mathfrak{B}_T$ of the space $U_H$ is not required. However, it is important that such functions actually exist, which is stated in the following lemma.

**Lemma 3.1 (Local bubble functions).** There exists a constant $C_b > 0$ which solely depends on the regularity of $\mathcal{T}_H$ and the polynomial degree $p$, such that for all $T \in \mathcal{T}_H$ and $j = 1, \dots, N$, there exists a function $b_{j,T} \in H_0^1(T)$ with $\Pi_H b_{j,T} = \Lambda_{j,T}$ and

$$ \left\| \nabla \sum_{j=1}^{N} c_j b_{j,T} \right\|_{L^2(T)} \le C_b H^{-1} \left\| \sum_{j=1}^{N} c_j \Lambda_{j,T} \right\|_{L^2(T)} $$

for any coefficients $c_j \in \mathbb{R}$. Regarding the scaling of $C_b$ with respect to $p$, we have $C_b \sim (p+1)^2$.

*Proof.* The proof is a special case of [Mai21, Cor. 3.6]. $\square$

<figure><img src="image_2.png" /></figure>

Based on the bases $\mathfrak{D}$ of $V_H$ and $\mathfrak{B}$ of $U_H$, we can construct a continuous mapping $\mathcal{B}_H : V_H \to U_H \subset H_0^1(D)$ that is uniquely defined by the property that

$$ \mathcal{B}_H \Lambda_{j,T} = b_{j,T} \quad \text{for all } T \in \mathcal{T}_H, 1 \le j \le N. $$

This bubble operator $\mathcal{B}_H$ can be extended to functions in $L^2(D)$ with the definition

$$ \mathcal{B}_H v := \mathcal{B}_H \Pi_H v \quad \text{for all } v \in L^2(D) $$

and it can be shown that it has the following representation,

$$ (3.4) \qquad \mathcal{B}_H v = \sum_{T \in \mathcal{T}_H} \sum_{j=1}^{N} \left( \int_T v \Lambda_{j,T} \, dx \right) b_{j,T}, \quad v \in L^2(D). $$

It is straightforward that $\ker \Pi_H = \ker \mathcal{B}_H$. Using Lemma 3.1, we directly have

$$ (3.5) \qquad \|\nabla \mathcal{B}_H v\|_{L^2(T)} \le C_{\text{bo}} H^{-1} \|v\|_{L^2(T)}, \quad v \in L^2(T), T \in \mathcal{T}_H. $$

The operator $\mathcal{B}_H$ will be important for the definition of a prototypical multiscale method in the following section and to link it to the method stated in [Mai21].

For later use, we introduce a quasi-interpolation operator $\mathcal{P}_H : L^2(D) \to H_0^1(D)$ which has the same kernel as the $L^2$-projection $\Pi_H$. The construction of $\mathcal{P}_H$ is based on another (low-order) quasi-interpolation operator defined as $\mathcal{I}_H := \mathcal{E}_H \circ \Pi_H^0$, where $\Pi_H^0$ denotes the $L^2$-projection onto piecewise constants and $\mathcal{E}_H$ is the averaging operator mapping piecewise constants to the space of continuous piecewise linear functions. For inner nodes $z$, we define

$$ (\mathcal{E}_H v)(z) := \frac{1}{\text{card}\{T \in \mathcal{T}_H : z \in \bar{T}\}} \sum_{T \in \mathcal{T}_H : z \in \bar{T}} v|_T(z) $$
and we set $(\mathcal{E}_{H}v)(z) := 0$ for nodes at the boundary. The operator $\mathcal{E}_{H}$ is well-known from the theory of domain decomposition methods, see [Osw93, Bre94]. Using the bubbles $b_{j,T}$, we then locally manipulate the operator $\mathcal{I}_{H}$ such that its kernel coincides with the kernel of $\Pi_{H}$. More precisely, for any $v \in L^{2}(D)$ we define

$$ (3.6) \qquad \mathcal{P}_{H}v := \mathcal{I}_{H}v + \mathcal{B}_{H}(v - \mathcal{I}_{H}v). $$

The operator $\mathcal{P}_{H}$ has been introduced in [AHP21, eq. (3.12)] employing very general functionals, so-called quantities of interest.

Remark (Effect of $\mathcal{P}_{H}$). In the definition of $\mathcal{P}_{H}$, the additional term involving $\mathcal{I}_{H}$ only has an effect on the zero-order bubble function (cf. Lemma 6.1 below for a proof). That is, the image of $\mathcal{P}_{H}$ is given by the span of the basis functions $\mathfrak{B}$ of $\mathcal{U}_{H}$, except for the zero-order bubble function, which is replaced by a modified function with a slightly increased support. The motivation for this modification is that $\mathcal{P}_{H}$ is better suited than the operator $\mathcal{B}_{H}$ to approximate (piecewise) constant functions. This property is crucial for an $H$-independent stability estimate as proved in Lemma 3.2 and will be important to obtain an improved version of the method in [Mai21].

The precise properties of $\mathcal{P}_{H}$ are investigated in the following lemma.

**Lemma 3.2** (Properties of $\mathcal{P}_{H}$). The operator $\mathcal{P}_{H}$ is a projection, i.e., it satisfies $\mathcal{P}_{H}^{2} = \mathcal{P}_{H}$. Moreover, its kernel coincides with the kernel of $\Pi_{H}$, i.e.,

$$ \ker \Pi_{H} = \ker \mathcal{P}_{H}. $$

Further, there exists a constant $C_{\mathcal{P}_{H}} > 0$, which solely depends on the regularity of $\mathcal{T}_{H}$ and the polynomial degree $p$, more precisely, $C_{\mathcal{P}_{H}} \sim (p+1)^{2}$, such that for all $v \in H_{0}^{1}(D)|_{\mathcal{N}(T)}$, $T \in \mathcal{T}_{H}$

$$ (3.7) \qquad \|\nabla \mathcal{P}_{H}v\|_{L^{2}(T)} + H^{-1}\|(1 - \mathcal{P}_{H})v\|_{L^{2}(T)} \le C_{\mathcal{P}_{H}} \|\nabla v\|_{L^{2}(\mathcal{N}(T))}. $$

*Proof.* We first prove that $\mathcal{P}_{H}$ is a projection. A technical calculation shows that $\mathcal{P}_{H}$ admits the representation

$$ (3.8) \qquad \mathcal{P}_{H}v = \sum_{T \in \mathcal{T}_{H}} \sum_{j=1}^{N} \left( \int_{T} v \Lambda_{j,T} \, dx \right) \mathcal{P}_{H}b_{j,T}. $$

Using this representation, the projection property follows from the identity

$$ \int_{T} \mathcal{P}_{H}b_{i,K}\Lambda_{j,T} \, dx = \delta_{ij}\delta_{KT}, $$

where $\delta$ denotes the Kronecker delta function.

The equality of the kernels follows directly from (3.8) using the fact that the values $\int_{T} v \Lambda_{j,T} \, dx$ are also the coefficients of $\Pi_{H}v$ when expanded with respect to the basis function $\Lambda_{j,T}$.

For the proof of (3.7), we use the stability and approximation properties of $\mathcal{I}_{H}$, see, e.g., [EG17]. There exists a constant $C_{\mathcal{I}_{H}} > 0$ such that for all $v \in H_{0}^{1}(D)|_{\mathcal{N}(T)}$ and $T \in \mathcal{T}_{H}$,

$$ (3.9) \qquad \|\nabla \mathcal{I}_{H}v\|_{L^{2}(T)} + H^{-1}\|(1 - \mathcal{I}_{H})v\|_{L^{2}(T)} \le C_{\mathcal{I}_{H}} \|\nabla v\|_{L^{2}(\mathcal{N}(T))}. $$
Using the stability estimate of the bubble operator (3.5) and the properties of $\mathcal{I}_H$ stated in (3.9), we obtain

$$ \begin{aligned} \|\nabla \mathcal{P}_H v\|_{L^2(T)} &\le \|\nabla \mathcal{I}_H v\|_{L^2(T)} + C_{\text{bo}} H^{-1} \|v - \mathcal{I}_H v\|_{L^2(T)} \\ &\le C_{\mathcal{I}_H} (1 + C_{\text{bo}}) \|\nabla v\|_{L^2(\mathbb{N}(T))}, \end{aligned} $$

which proves the $H^1$-stability of $\mathcal{P}_H$ in (3.7). The scaling of $C_{\mathcal{P}_H}$ with respect to $p$ follows directly by the corresponding estimate for $C_{\text{bo}}$. The second ingredient, the $L^2$-approximation property in (3.7), follows with similar arguments and is omitted for the sake of brevity. $\square$

## 4. A PROTOTYPICAL MULTISCALE METHOD

The conformity of the space $\mathcal{U}_H$ is essential for the construction of a prototypical multiscale method, which directly relates to the classical LOD methodology as presented in [MP14, HP13, MP20]. We will adapt ideas from [HP22] and prove that the constructed $p$th-order prototypical multiscale method is equivalent to the one proposed in [Mai21]. The benefit of this alternative construction is that it enables an improved localization compared to [Mai21] that will be shown in the subsequent section based on the operator $\mathcal{P}_H$. To provide a better overview of the important spaces and operators, we have provided a summary of the most important definitions in Table 4.1 below.

In the spirit of the LOD, we define a so-called fine-scale space $\mathcal{W}$ of functions that are not captured by the $L^2$-projection $\Pi_H$,

$$ \mathcal{W} := \ker \Pi_H |_{H_0^1(D)}. $$

Note that the spaces $\mathcal{U}_H$ and $\mathcal{W}$ define an $L^2$-orthogonal decomposition of the space $H_0^1(D)$, that is

$$ H_0^1(D) = \mathcal{U}_H \oplus \mathcal{W}, \qquad H_0^1(D) \ni v = \mathcal{B}_H v + (1 - \mathcal{B}_H)v. $$

However, this decomposition is not suited for the approximation of the solution to (2.1) and we would like to modify the space $\mathcal{U}_H$ such that the decomposition is a-orthogonal instead and thus includes information on the problem at hand. Based on the space $\mathcal{W}$, we therefore introduce a so-called correction operator $\mathcal{C}: H_0^1(D) \to \mathcal{W}$. For any $v \in H_0^1(D)$, the correction $\mathcal{C}v \in \mathcal{W}$ is defined as the solution to

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

By construction, this defines an a-orthogonal decomposition of $H_0^1(D)$ into two subspaces. More precisely, we set $\tilde{U}_H := (1 - \mathcal{C})H_0^1(D)$ and observe that

$$ H_0^1(D) = \tilde{U}_H \oplus \mathcal{W}, \qquad H_0^1(D) \ni v = (1 - \mathcal{C})v + \mathcal{C}v. $$

Note that for the definition of $\tilde{U}_H$, it suffices to consider the coarse functions $\mathcal{U}_H$ in the domain of $(1 - \mathcal{C})$. That is, with the above operators, the property $\Pi_H(1 - \mathcal{B}_H)v = 0$ for all $v \in H_0^1(D)$, and the fact that $\mathcal{C}$ is the identity on $\mathcal{W}$, we deduce

$$ \begin{aligned} \tilde{U}_H &= (1 - \mathcal{C})H_0^1(D) = (1 - \mathcal{C})(1 - \mathcal{B}_H)H_0^1(D) + (1 - \mathcal{C})\mathcal{B}_H H_0^1(D) \\ &= (1 - \mathcal{C})\mathcal{B}_H H_0^1(D) = (1 - \mathcal{C})\mathcal{B}_H \mathcal{V}_H = (1 - \mathcal{C})\mathcal{U}_H. \end{aligned} $$
<table><caption>TABLE 4.1. Important spaces and operators</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>V<sub>H</sub></td><td>piecewise polynomial space with coordinate degree less or equal than p</td></tr><tr><td>U<sub>H</sub></td><td>H<sup>1</sup><sub>0</sub>(D)-conforming space with &Pi;<sub>H</sub>U<sub>H</sub> = V<sub>H</sub> and dimU<sub>H</sub> = dimV<sub>H</sub></td></tr><tr><td>&Pi;<sub>H</sub></td><td>L<sup>2</sup>-projection onto V<sub>H</sub></td></tr><tr><td>I<sub>H</sub></td><td>first-order conforming quasi-interpolation operator</td></tr><tr><td>B<sub>H</sub></td><td>bubble operator B<sub>H</sub>: V<sub>H</sub> &rarr; U<sub>H</sub> &subset; H<sup>1</sup><sub>0</sub>(D).</td></tr><tr><td>P<sub>H</sub></td><td>P<sub>H</sub>&bull; := I<sub>H</sub> &bull; + B<sub>H</sub>(&bull; &minus; I<sub>H</sub>&bull;) with ker P<sub>H</sub> = ker &Pi;<sub>H</sub></td></tr><tr><td>W</td><td>fine-scale space W := ker &Pi;<sub>H</sub>|<sub>H<sup>1</sup><sub>0</sub>(D)</sub></td></tr><tr><td>C</td><td>correction operator C: H<sup>1</sup><sub>0</sub>(D) &rarr; W</td></tr><tr><td>U&#771;<sub>H</sub></td><td>a-orthogonal complement of W in H<sup>1</sup><sub>0</sub>(D), U&#771;<sub>H</sub> = (1 &minus; C)U<sub>H</sub></td></tr></tbody></table>

Therefore, $(1 - C): U_H \to \tilde{U}_H$ is a bijective operator with inverse $B_H|_{\tilde{U}_H}$ and it holds $\dim U_H = \dim \tilde{U}_H$. Further,

$$ (4.3) \qquad \tilde{\mathfrak{B}} := \bigcup_{b \in \mathfrak{B}} (1 - C)b $$

defines a basis of $\tilde{U}_H$.

The newly defined space $\tilde{U}_H$ is a suitable candidate for approximating the solution $u$ to (2.1) and we define a Galerkin method based on the space $\tilde{U}_H$ as follows: find $\tilde{u}_H \in \tilde{U}_H$ such that

$$ (4.4) \qquad a(\tilde{u}_H, \tilde{v}_H) = (f, \tilde{v}_H)_{L^2(D)} \quad \text{for all } \tilde{v}_H \in \tilde{U}_H. $$

The well-posedness of (4.4) follows directly from the Lax-Milgram theorem with a stability estimate as in (2.2). The method stated in (4.4) is typically called the *prototypical method*, because it is defined based on the global correction operator $C$ that needs to be suitably localized and discretized in practice. This will be investigated in the following sections.

Before we turn to an optimal error estimate for the solution $\tilde{u}_H$ to (4.4), we show that the approximation space $\tilde{U}_H$ coincides with the multiscale space $RV_H$ from [Mai21, Sec. 3.1]. Therein, the operator $\mathcal{R}$ is defined by seeking for any $v_H \in V_H$ a pair $(Rv_H, \lambda_{v_H}) \in H_0^1(D) \times V_H$ solving

$$ (4.5) \qquad \begin{aligned} &a(Rv_H, v) + (\lambda_{v_H}, v)_{L^2(D)} = 0, \\ &(Rv_H, \mu_H)_{L^2(D)} = (v_H, \mu_H)_{L^2(D)} \end{aligned} $$

for all $v \in H_0^1(D)$ and all $\mu_H \in V_H$. Note that the unique solvability of (4.5) is shown in [Mai21].

The equivalence of spaces is stated in the following lemma.

**Lemma 4.1** (Equivalence of spaces). It holds $RV_H = \tilde{U}_H$. More precisely, for any $v_H \in V_H$ we have the identity

$$ Rv_H = (1 - C)B_H v_H. $$

*Proof.* Let $v_H \in V_H$. We show that $(1 - C)B_H v_H$ solves (4.5) together with an appropriate Lagrange multiplier $\lambda_{v_H}$. The assertion then follows by the uniqueness of (4.5).

By (4.1), we have that $a((1 - C)B_H v_H, w) = 0$ for all $w \in W$. Let now $\mu_H \in V_H$ and $v \in H_0^1(D)$. Due to the $L^2$-orthogonality of the spaces $V_H$ and $W$, we obtain

$$ (4.6) \qquad ((1 - C)B_H v_H, \mu_H)_{L^2(D)} = (B_H v_H, \mu_H)_{L^2(D)} = (v_H, \mu_H)_{L^2(D)}, $$
which shows that $(1 - C)B_H v_H$ solves the second equation of (4.5). We now turn to the first equality. Since $(1 - B_H)v \in W$, we have that

$$ (4.7) \qquad a((1 - C)B_H v_H, v) = a((1 - C)B_H v_H, B_H v). $$

To fulfill (4.5), we define a suitable Lagrange multiplier. Let $\lambda_{v_H} \in V_H$ be the solution to

$$ (4.8) \qquad (\lambda_{v_H}, q_H)_{L^2(D)} = -a((1 - C)B_H v_H, B_H q_H) \quad \text{for all } q_H \in V_H, $$

which is well-defined due to the coercivity of $(\cdot, \cdot)_{L^2(D)}$ with respect to $V_H \subset L^2(D)$ and the boundedness of $a$ and $B_H$. The combination of (4.7) and (4.8) leads to

$$ (4.9) \qquad a((1 - C)B_H v_H, v) + (\lambda_{v_H}, v)_{L^2(D)} = 0, $$

which concludes the proof.

<figure><img src="image_3.png" /></figure>

**4.1. Error estimates.** The equality of the space $RV_H$ and $\tilde{U}_H$ implies that the error estimates of the prototypical method (4.4) follow directly from the corresponding ones in [Mai21, Sec. 3.2]. In particular, we have the following theorem.

**Theorem 4.2 (Error of the prototypical method).** Assume that $f \in H^k(\mathcal{T}_H)$, $k \in \mathbb{N}_0$, and set $s := \min\{k, p+1\}$. Further, let $u \in H_0^1(D)$ and $\tilde{u}_H \in \tilde{U}_H$ be the solutions of (2.1) and (4.4), respectively. Then there exists a constant $C_e > 0$, which solely depends on the regularity of $\mathcal{T}_H$ and the polynomial degree $p$, such that

$$ (4.10) \qquad \|\nabla(u - \tilde{u}_H)\|_{L^2(D)} \le C_e H^{s+1} |f|_{H^s(\mathcal{T}_H)}, $$

$$ (4.11) \qquad \|u - \tilde{u}_H\|_{L^2(D)} \le C_e H^{s+2} |f|_{H^s(\mathcal{T}_H)} $$

with the notation $H^0(\mathcal{T}_H) := L^2(D)$ and $|\cdot|_{H^0(\mathcal{T}_H)} := \|\cdot\|_{L^2(D)}$. Note that the constant $C_e$ generally depends in a positive way on $p$ and $s$ and scales (at least) like $(p+1)^{-1}$.

*Proof.* The result is a direct consequence of [Mai21, Thm. 3.1].

<figure><img src="image_4.png" /></figure>

## 5. DECAY AND LOCALIZATION

As mentioned above, the prototypical method defined in (4.4) is not a practically feasible method since the corresponding basis functions are globally defined. In this section, we show that the Green's function of the correction decays exponentially fast, which justifies the localization of the basis functions. For enhanced stability properties, this paper employs a localization technique based on localized element corrections.

**5.1. Exponential decay.** For any $T \in \mathcal{T}_H$, the element correction operator $C_T: H_0^1(D) \to W$ is defined, for any $v \in H_0^1(D)$, by the solution $C_T v \in W$ to

$$ (5.1) \qquad a(C_T v, w) = a_T(v, w) \quad \text{for all } w \in W, $$

where we use the restricted bilinear form $a_T$ defined by

$$ a_T(v, w) := \int_T A \nabla v \cdot \nabla w \,dx, \quad v, w \in H_0^1(D). $$

We can then split the correction operator $C$ into its element contributions,

$$ (5.2) \qquad C = \sum_{T \in \mathcal{T}_H} C_T. $$

The following lemma proves that the moduli of the element corrections decay exponentially fast relative to the coarse mesh.
**Lemma 5.1** (Exponential decay of element corrections). There exists a constant $C_d > 0$, which is independent of $H, \ell, T$, such that for all $T \in \mathcal{T}_H$, $v \in H_0^1(D)$, $\ell \in \mathbb{N}$,

$$ \|\nabla C_T v\|_{L^2(D \setminus \mathbb{N}^\ell(T))} \le \exp(-C_d \ell) \|\nabla C_T v\|_{L^2(D)}. $$

Regarding the polynomial degree $p$, we have $C_d \sim (p+1)^{-1}$.

*Proof.* The proof is presented in Appendix A. □

We emphasize that the scaling of the constant $C_d$ in Lemma 5.1 is similar to the scaling in [Mai21, Thm. 4.1]. However, as pointed out in [Mai21, Rem. 4.3], the negative scaling with respect to $p$ seems to be pessimistic and is possibly not sharp. In fact, a positive scaling with respect to $p$ is observed in our numerical experiments, see Figure 7.3.

**5.2. Localized corrections.** The exponential decay of the element corrections $C_T$ motivates a localization to element patches around the element $T$. For a given oversampling parameter $\ell \in \mathbb{N}$, we therefore define a localized fine-scale space as

$$ W_T^\ell := W(\mathbb{N}^\ell(T)) := \{w \in W : \text{supp}(w) \subset \mathbb{N}^\ell(T)\} \subset W. $$

For any $T \in \mathcal{T}_H$, the localized element correction operator $C_T^\ell: H_0^1(D) \to W_T^\ell$ is defined, for any $v \in H_0^1(D)$, by the solution $C_T^\ell v \in W_T^\ell$ to

$$ (5.3) \qquad a(C_T^\ell v, w) = a_T(v, w) \quad \text{for all } w \in W_T^\ell. $$

Without localization, the correction operator $C$ equals the sum of the elements contributions $C_T$; see (5.2). Thus, we also define the localized correction operator $C^\ell$ as the sum of the localized element contributions,

$$ C^\ell := \sum_{T \in \mathcal{T}_H} C_T^\ell. $$

The next lemma shows exponential approximation properties of this localized correction operator.

**Lemma 5.2** (Localization error of correctors). There exists a constant $C_1 > 0$ independent of $H, \ell$ such that for all $v \in H_0^1(D)$, $\ell \in \mathbb{N}$

$$ \|\nabla(C - C^\ell)v\|_{L^2(D)} \le C_1 \ell^{d/2} \exp(-C_d \ell) \|\nabla v\|_{L^2(D)} $$

with the constant $C_d$ from Lemma 5.1. With respect to the polynomial degree, we have the scaling $C_1 \sim (p+1)$.

*Proof.* The proof is stated in Appendix A. □

## 6. PRACTICAL MULTISCALE METHOD

In this section, we present a practical (localized) multiscale method based on the considerations in the previous sections. As mentioned above, the (prototypical) basis $\tilde{\mathfrak{B}}$ defined in (4.3) consists of non-local but rapidly decaying basis functions; see Lemma 5.1. In order to avoid unfeasible global computations, we replace the prototypical basis by its (quasi-)local counterpart defined as

$$ (6.1) \qquad \tilde{\mathfrak{B}}^\ell := \bigcup_{b \in \mathfrak{B}} (1 - C^\ell) \mathcal{P}_H b. $$
Similarly, we also define a localized ansatz space for the practical multiscale method,

$$ \tilde{U}_H^\ell := (1 - C^\ell) P_H U_H. $$

This choice is different from the choice in [Mai21] as it employs a localization based on element corrections and additionally utilizes the operator $P_H$ in the basis definition. Henceforth, for $T \in \mathcal{T}_H$, we denote by $\Lambda_{1,T}$ and $b_{1,T}$ the constant Legendre basis function and its associated bubble function, respectively. In the following lemma, it is shown that $P_H$ only has an effect on the bubble functions corresponding to the polynomial degree 0, i.e., the set

$$ \mathfrak{B}^0 := \{b_{1,T} : T \in \mathcal{T}_H\} \subset \mathfrak{B}. $$

**Lemma 6.1 (Application of $P_H$).** For all $T \in \mathcal{T}_H$ and $j > 1$, it holds

$$ P_H b_{j,T} = b_{j,T}. $$

*Proof.* We consider a fixed element $T \in \mathcal{T}_H$. By definition, the bubble function $b_{j,T}$ satisfies the relation $\Pi_H b_{j,T} = \Lambda_{j,T}$. Using this, we obtain for all $j > 1$

$$ \Pi_H^0 b_{j,T} = (b_{j,T}, \Lambda_{1,T})_{L^2(T)} \Lambda_{1,T} = (\Lambda_{j,T}, \Lambda_{1,T})_{L^2(T)} \Lambda_{1,T} = 0. $$

Hence, we have $\mathcal{I}_H b_{j,T} = 0$, which implies $P_H b_{j,T} = b_{j,T}$ using (3.6). $\square$

<figure><img src="image_5.png" /></figure>

We show that, for any polynomial degree $p$, the construction in (6.1) ensures a more robust localization compared to [Mai21]. Notably, for this positive effect only the basis functions corresponding to elements in $\mathfrak{B}^0$ need to be modified. Their supports are enlarged by one additional layer of elements. This strategy has been used in [HP22] to cure instabilities for the case $p = 0$; see also [AHP21].

The final practical multiscale method approximates the solution $u$ of (2.1) by performing a Galerkin approach based on the space $\tilde{U}_H^\ell$. It seeks $\tilde{u}_H^\ell \in \tilde{U}_H^\ell$ such that

$$ (6.2) \qquad a(\tilde{u}_H^\ell, \tilde{v}_H^\ell) = (f, \tilde{v}_H^\ell)_{L^2(D)} \quad \text{for all } \tilde{v}_H^\ell \in \tilde{U}_H^\ell. $$

The next theorem proves an error estimate for the solution $\tilde{u}_H^\ell$.

**Theorem 6.2 (A priori error estimate).** Assume that $f \in H^k(\mathcal{T}_H)$, $k \in \mathbb{N}_0$ and set $s := \min\{k, p+1\}$. Further, let $u \in H_0^1(D)$ and $\tilde{u}_H^\ell \in \tilde{U}_H^\ell$ be the solutions to (2.1) and (6.2), respectively. Then, there exists a constant $C_{el} > 0$, which solely depends on the regularity of $\mathcal{T}_H$ and the polynomial degree $p$, such that

$$ (6.3) \qquad \|\nabla(u - \tilde{u}_H^\ell)\|_{L^2(D)} \le C_e H^{s+1} |f|_{H^s(\mathcal{T}_H)} + C_{el} \ell^{d/2} \exp(-C_d \ell) \|f\|_{L^2(D)}, $$

$$ (6.4) \qquad \|u - \tilde{u}_H^\ell\|_{L^2(D)} \le \beta (C_e H + C_{el} \ell^{d/2} \exp(-C_d \ell)) \|\nabla(u - \tilde{u}_H^\ell)\|_{L^2(D)} $$

with the constants $C_e$ from Theorem 4.2 and $C_d$ from Lemma 5.1. With respect to the polynomial degree, we have $C_{el} \sim (p+1)^3$.

*Proof.* By Céa's Lemma, we obtain for any $v \in \tilde{U}_H^\ell$

$$ \|\nabla(u - \tilde{u}_H^\ell)\|_{L^2(D)} \le \sqrt{\frac{\beta}{\alpha}} \|\nabla(u - v)\|_{L^2(D)}. $$

Adding and subtracting $\tilde{u}_H$ and employing the triangle inequality, we obtain

$$ (6.5) \qquad \|\nabla(u - v)\|_{L^2(D)} \le \|\nabla(u - \tilde{u}_H)\|_{L^2(D)} + \|\nabla(\tilde{u}_H - v)\|_{L^2(D)}. $$
The first term is the error of the prototypical method for which we use the bound from Theorem 4.2. To bound the second term, we require the equality $\tilde{u}_H = (1 - C)\mathcal{P}_H u$, which can be shown as follows. First, we observe that

$$ (1 - C)\mathcal{P}_H u = (1 - C)u - (1 - C)(1 - \mathcal{P}_H)u = (1 - C)u $$

since $(1 - \mathcal{P}_H)u \in W = \ker \Pi_H|_{H_0^1(D)}$. Further, we use that any $\tilde{v}_H \in \tilde{U}_H$ can be written as $\tilde{v}_H = (1 - C)b_H$ for some $b_H \in \mathcal{U}_H$ and calculate

$$ \begin{aligned} a((1 - C)u, \tilde{v}_H) &= a(u, \tilde{v}_H) - a(Cu, \tilde{v}_H) \\ &= (f, \tilde{v}_H)_{L^2(D)} - a(Cu, (1 - C)b_H) \\ &= (f, \tilde{v}_H)_{L^2(D)}, \end{aligned} $$

where we use the symmetry of $a$ and the definition of $C$ (cf. (4.1)) in the last step. Since $\tilde{u}_H$ is uniquely defined, we deduce that $\tilde{u}_H = (1 - C)u = (1 - C)\mathcal{P}_H u$.

We now go back to the second term in (6.5). Since $\mathcal{P}_H(u - B_H u) = 0$ (cf. the representation of $\mathcal{P}_H$ in (3.8)), $v := (1 - C^\ell)\mathcal{P}_H u = (1 - C^\ell)\mathcal{P}_H B_H u \in \tilde{U}_H^\ell$ is a valid approximation of $\tilde{u}_H$ in the space $\tilde{U}_H^\ell$. Using this and the above equality for $\tilde{u}_H$, we obtain

$$ \begin{aligned} \|\nabla(\tilde{u}_H - v)\|_{L^2(D)} &= \|\nabla(C - C^\ell)\mathcal{P}_H u\|_{L^2(D)} \\ &\le C_1 \ell^{d/2} \exp(-C_d \ell) \|\nabla \mathcal{P}_H u\|_{L^2(D)} \\ &\le C_{\mathcal{P}_H} C_1 \ell^{d/2} \exp(-C_d \ell) \|\nabla u\|_{L^2(D)} \\ &\le \alpha^{-1} C_{\text{PF}} C_{\mathcal{P}_H} C_1 \ell^{d/2} \exp(-C_d \ell) \|f\|_{L^2(D)} \end{aligned} $$

employing also Lemma 5.2 as well as the estimates (2.2) and (3.7). Combining the above error bounds yields the desired $H^1$-error estimate (6.3). Moreover, the scaling of the constant $C_{\text{el}}$ with respect to $p$ follows directly from the scaling of $C_{\mathcal{P}_H}$ and $C_1$.

To prove (6.4), we require the Aubin-Nitsche duality argument. We denote by $z_g \in H_0^1(D)$ the solution to (2.1) for the right-hand side $g := u - \tilde{u}_H^\ell \in L^2(D)$. Using the Galerkin orthogonality, we obtain for any $v \in \tilde{U}_H^\ell$

$$ \|u - \tilde{u}_H^\ell\|_{L^2(D)}^2 = (g, u - \tilde{u}_H^\ell)_{L^2(D)} = a(z_g, u - \tilde{u}_H^\ell) = a(z_g - v, u - \tilde{u}_H^\ell). $$

Choosing $v$ as the solution of the practical multiscale method (6.2) for the particular right-hand side $g$ and using (6.3) for $s = 0$, we get

$$ \begin{aligned} \|u - \tilde{u}_H^\ell\|_{L^2(D)}^2 &\le \beta \|\nabla(z_g - v)\|_{L^2(D)} \|\nabla(u - \tilde{u}_H^\ell)\|_{L^2(D)} \\ &\le \beta(C_e H + C_{\text{el}} \ell^{d/2} \exp(-C_d \ell)) \|g\|_{L^2(D)} \|\nabla(u - \tilde{u}_H^\ell)\|_{L^2(D)} \\ &= \beta(C_e H + C_{\text{el}} \ell^{d/2} \exp(-C_d \ell)) \|u - \tilde{u}_H^\ell\|_{L^2(D)} \|\nabla(u - \tilde{u}_H^\ell)\|_{L^2(D)}. \end{aligned} $$

Dividing by $\|u - \tilde{u}_H^\ell\|_{L^2(D)}$ on both sides, the desired $L^2$-norm estimate (6.4) follows. $\square$

Remark (Choice of $\ell$). The estimates in Theorem 6.2 suggest to choose $\ell \gtrsim (p+1)(s+1)|\log(H)|$ in order to retain the overall convergence rates as stated in Theorem 4.2. This is in line with the theoretical requirements in [Mai21, Cor. 4.6]. We emphasize, however, that this dependence on $p$ is rather pessimistic as also observed in numerical experiments.
<figure><img src="image_6.png" /><figcaption>FIGURE 6.1. Bubble functions $b_{1,T}$ for polynomial degrees $p = 0, 1$ (left) and polynomial degrees $p = 2, 3$ (right).</figcaption></figure>

## 7. NUMERICAL EXAMPLES

This section numerically investigates the proposed improved high-order method. For all numerical experiments, we use the domain $D = (0, 1)^2$ and consider a coarse Cartesian mesh $\mathcal{T}_H$ of $D$. Note that, henceforth, the mesh size denotes the element side-length instead of the diameter. The subsequent numerical experiments require an explicit construction of bubble functions. One possible choice is outlined in the following remark.

*Remark (Construction of bubbles).* For the construction of the bubble functions associated to an element $T \in \mathcal{T}_H$, we introduce the function $\theta \in H_0^1(T)$ obtained by multiplying all nodal $Q_1$ basis functions within $T$. The bubble functions $b_{j,T}$ are then obtained by the ansatz $b_{j,T} = \sum_{i=1}^N c_i \theta \Lambda_{i,T}$, where the coefficients are chosen such that

$$ \Pi_H b_{j,T} = \Lambda_{j,T}. $$

This condition translates to a $N \times N$ linear system, which is invertible under the condition that $\{\Pi_H(\theta \Lambda_{j,T}) : j = 1, \dots, N\}$ are linearly independent. The linear independence can be proved utilizing the positivity of $\theta$ in the interior of $T$. We emphasize that, by Lemma 6.1 and the localized counterpart of (4.5), only the zero-order bubble functions are explicitly required for the computation of the basis functions (6.1). The remaining basis functions can be computed as in [Mai21]. Zero-order bubble functions corresponding to different polynomial degrees which fulfill the above conditions are depicted in Fig. 6.1. We note that a symmetry argument shows that the bubble functions for consecutive polynomial degrees ($p = k$ and $p = k + 1, k = 0, 2, 4, \dots$) coincide.

Further, we need to discretize the infinite-dimensional localized corrections as defined in (5.3). For this, we employ the $Q_1$-finite element method on a fine mesh of the respective patches with mesh size $h = 2^{-9}$ obtained by uniform refinement of the coarse patch mesh.
<figure><img src="image_7.png" /><figcaption>FIGURE 7.1. Coefficients $A_1$, $A_2$ used in the numerical experiments.</figcaption></figure>

For details regarding the discretization of these problems, we refer to [Mai21, Sec. 4.3]. We emphasize that also the construction of the bubbles in Section 7 carries over to the fully discrete setting. All errors are computed with respect to the $Q_1$-finite element reference solution on the global Cartesian mesh with mesh size $h = 2^{-9}$. All considered errors are relative errors measured in the energy norm

$$ \| \cdot \|_a := \| A^{1/2} \nabla \cdot \|_{L^2(D)} $$

For the numerical experiments, we use the two scalar coefficients $A_1$, $A_2$ as depicted in Figure 7.1. Both coefficients are piecewise constant with respect to the Cartesian mesh $\mathcal{T}_{2-7}$. The coefficient $A_1$ is the sum of random contribution on several scales with $2^{-7}$ being the smallest scale. It varies within the interval $[1, 4]$. The coefficient $A_2$ features some channels of conductivity 10 which are placed on top of a noisy background. On each element of $\mathcal{T}_{2-7}$, the noisy background takes values generated by realizations of independent random variables which are uniformly distributed in the interval $[0.1, 1]$. As right-hand side, we consider the smooth function

$$ f(x_1, x_2) = (x_1 + \cos(3\pi x_1)) \cdot x_2^3. $$

For the first numerical experiment, we consider the diffusion coefficient $A_1$ and investigate the convergence of the proposed stabilized high-order method (referred to as sp-LOD) for polynomial degrees $p = 0, \dots, 3$. The results are compared to the localized high-order method in [Mai21] (referred to as p-LOD). For fixed oversampling parameters, Figure 7.2 plots the relative energy errors of both methods as a function of the coarse mesh size $H$. As reference, also a line indicating the expected rate of convergence is shown. Provided that the oversampling parameter $\ell$ is chosen large enough, one observes convergence of order $p+2$ for the methods of degree $p$. This numerically confirms the theoretical prediction from Theorem 4.2. Most importantly, Figure 7.2 confirms the improved stability properties of the sp-LOD compared to the p-LOD. For fixed oversampling parameters, the error of the sp-LOD first decreases and then stagnates as $H$ is decreased. In contrast, one observes an increasing error of the p-LOD for mesh sizes below a certain critical threshold depending on the oversampling parameter, which is in line with the theoretical results in [Mai21,
<figure><img src="image_8.png" /><figcaption>FIGURE 7.2. Relative errors with respect to $H$ for the coefficient $A_1$, different values of $\ell$, and polynomial degrees $p = 0, \dots, 3$ (top left to bottom right).</figcaption></figure>

Thm. 4.4]. Note that for both localization strategies, the stagnation (respectively deterioration) of the error can be avoided if the oversampling parameter is appropriately increased when $H$ is decreased (i.e., logarithmically in $1/H$; cf. also Theorem 6.2).

Next, we specifically investigate the localization errors. Therefore, we choose the right-hand side $f \equiv 1$ for which the first term on the right-hand side of error estimate (6.3) vanishes and, hence, only the localization error remains. For different choices of $\ell$ and $p$, we numerically investigate the decay of the localization error. In Fig. 7.3, one observes that the localization error exponentially decays as $\ell$ is increased which numerically confirms Theorem 6.2. Further, one observes that for larger values of $p$, the decay behavior is
<figure><img src="image_9.png" /><figcaption>FIGURE 7.3. Relative localization errors for the coefficient $A_1$ and $H = 2^{-4}$ for fixed $p$ and varying $\ell$ (left) and for fixed $\ell$ and varying $p$ (right).</figcaption></figure>

further improved. This can also be seen in Fig. 7.3, where the errors are plotted with respect to $p$ for different fixed values of $\ell$. This is in line with the numerical results in [Mai21] and indicates that the scaling of $C_d$ in Lemma 5.1 with respect to $p$ is probably not sharp.

Next, we depict selected prototypical basis functions. Therefore, we choose a fixed element $T$ of the coarse Cartesian mesh $\mathcal{T}_{2-4}$. For polynomial degrees $p = 0, \dots, 3$, Figure 7.4 shows plots of the discretized counterparts of the prototypical basis functions

$$ (1 - \mathcal{C})b_{1,T}, $$

where $b_{1,T} \in H_0^1(T)$ is the bubble function satisfying $\Pi_H b_{1,T} = \Lambda_{1,T}$ with $\Lambda_{1,T}$ denoting the $L^2$-normalized characteristic function of the element $T$ which is the first Legendre basis function. Note that the $p$-LOD and the $sp$-LOD have the same prototypical basis functions as both methods differ only with respect to the localization strategy. For all plotted basis functions in Figure 7.4, we use the same color coding indicating the absolute value of the respective basis functions. We employ a logarithmic color map which we restrict to the interval $[10^{-8}, 3.5]$. This means that all values below $10^{-8}$ are displayed in dark blue indiscriminate of their actual value. In Figure 7.4, one observes that, for larger polynomial degrees $p$, the basis functions tend to be more localized, i.e., a more rapid decay can be observed. This is in line with the results in [Mai21], where improved localization is observed for larger polynomial degrees.

In the second numerical experiment, we investigate the convergence of the $sp$-LOD for the more challenging diffusion coefficient $A_2$. For polynomial degrees $p = 0, \dots, 3$, Figure 7.5 plots the errors for oversampling parameters $\ell = 3, 5$ as a function of $H$. As reference, lines indicating the respective expected convergence rates are depicted. One observes that the errors plotted in Figure 7.5 are qualitatively similar to the ones in Figure 7.2. One obtains the expected optimal convergence rates of order $p+2$ provided that the oversampling parameter is chosen sufficiently large. Again, the method demonstrates
<figure><img src="image_10.png" /><figcaption>FIGURE 7.4. Ideal basis functions $(1 - C)b_{1,T}$ for the polynomial degrees $p = 0, \dots, 3$ (top left to bottom right).</figcaption></figure>

its improved stability properties, as the error first decreases and then stagnates for decreasing $H$ and does not show a deteriorating behavior. As mentioned above, the convergence rate of order $p + 2$ can be preserved if the oversampling parameter $\ell$ is logarithmically coupled to the mesh size $H$ as predicted by the theoretical result in Theorem 6.2.

## 8. CONCLUSIONS

In this work, we have presented an improved version of the high-order localized multi-scale method presented in [Mai21] which can be used to numerically solve elliptic multiscale problems with rough coefficients. The method only requires minimal assumptions on the coefficient and achieves higher-order approximation rates exploiting regularity properties
<figure><img src="image_11.png" /><figcaption>FIGURE 7.5. Convergence plot for the coefficient $A_2$ and selected values of $\ell$ and $p$.</figcaption></figure>

of the right-hand side only. The method is defined in the spirit of the Localized Orthogonal Decomposition methodology and is built upon appropriate high-order bubble functions. An improved localization strategy has been presented that cures the issue of deteriorating errors for smaller mesh sizes as theoretically and practically observed in [Mai21]. It slightly enlarges the support of the lowest-order basis functions and thereby notably achieves a stabilization for arbitrary polynomial degree. We have provided a complete priori error analysis and have demonstrated the improvements of the method in a set of numerical experiments.

## REFERENCES

[AB05] G. Allaire and R. Brizzi. A multiscale finite element method for numerical homogenization. *Multiscale Model. Simul.*, 4(3):790–812, 2005.

[AB12] A. Abdulle and Y. Bai. Reduced basis finite element heterogeneous multiscale method for high-order discretizations of elliptic homogenization problems. *J. Comput. Phys.*, 231(21):7014–7036, 2012.

[AHP21] R. Altmann, P. Henning, and D. Peterseim. Numerical homogenization beyond scale separation. *Acta Numer.*, 30:1–86, 2021.

[AHPV13] R. Araya, C. Harder, D. Paredes, and F. Valentin. Multiscale hybrid-mixed method. *SIAM J. Numer. Anal.*, 51(6):3505–3531, 2013.

[BCO94] 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 J. Numer. Anal.*, 31(4):945–981, 1994.

[BL11] I. Babuška and R. Lipton. Optimal local approximation spaces for generalized finite element methods with application to multiscale problems. *Multiscale Model. Simul.*, 9(1):373–406, 2011.

[BO83] I. Babuška and J. E. Osborn. Generalized finite element methods: their performance and their relation to mixed methods. *SIAM J. Numer. Anal.*, 20(3):510–536, 1983.
[BO00] I. Babuška and J. E. Osborn. Can a finite element method perform arbitrarily badly? *Math. Comp.*, 69(230):443–462, 2000.

[Bre94] S. C. Brenner. Two-level additive Schwarz preconditioners for nonconforming finite elements. *Contemp. Math.*, 180:9–14, 1994.

[CEL19] M. Cicuttin, A. Ern, and S. Lemaire. A hybrid high-order method for highly oscillatory elliptic problems. *Comput. Methods Appl. Math.*, 19(4):723–748, 2019.

[EG17] A. Ern and J.-L. Guermond. Finite element quasi-interpolation and best approximation. *ESAIM Math. Model. Numer. Anal.*, 51(4):1367–1385, 2017.

[EGH13] Y. Efendiev, J. Galvis, and T. Y. Hou. Generalized multiscale finite element methods (GMsFEM). *J. Comput. Phys.*, 251:116–135, 2013.

[GGS12] L. Grasedyck, I. Greff, and S. Sauter. The AL basis for the solution of elliptic problems in heterogeneous media. *Multiscale Model. Simul.*, 10(1):245–258, 2012.

[HP13] P. Henning and D. Peterseim. Oversampling for the multiscale finite element method. *Multiscale Model. Simul.*, 11(4):1149–1175, 2013.

[HP22] M. Hauck and D. Peterseim. Multi-resolution localized orthogonal decomposition for Helmholtz problems. *Multiscale Model. Simul.*, 20(2):657–684, 2022.

[HPV13] C. Harder, D. Paredes, and F. Valentin. A family of multiscale hybrid-mixed finite element methods for the Darcy equation with rough coefficients. *J. Comput. Phys.*, 245:107–130, 2013.

[HSS02] P. Houston, C. Schwab, and E. Süli. Discontinuous *hp*-finite element methods for advection-diffusion-reaction problems. *SIAM J. Numer. Anal.*, 39(6):2133–2163, 2002.

[HZZ14] J. S. Hesthaven, S. Zhang, and X. Zhu. High-order multiscale finite element method for elliptic problems. *Multiscale Model. Simul.*, 12(2):650–666, 2014.

[LMT12] R. Li, P. Ming, and F. Tang. An efficient high order heterogeneous multiscale method for elliptic problems. *Multiscale Model. Simul.*, 10(1):259–283, 2012.

[Mai20] R. Maier. *Computational Multiscale Methods in Unstructured Heterogeneous Media*. PhD thesis, University of Augsburg, 2020.

[Mai21] R. Maier. A high-order approach to elliptic multiscale problems with general unstructured coefficients. *SIAM J. Numer. Anal.*, 59(2):1067–1089, 2021.

[MP14] A. Målqvist and D. Peterseim. Localization of elliptic multiscale problems. *Math. Comp.*, 83(290):2583–2603, 2014.

[MP20] A. Målqvist and D. Peterseim. *Numerical homogenization by localized orthogonal decomposition, volume 5 of SIAM Spotlights*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2020.

[OS19] H. Owhadi and C. Scovel. *Operator-adapted wavelets, fast solvers, and numerical homogenization, volume 35 of Cambridge Monographs on Applied and Computational Mathematics*. Cambridge University Press, Cambridge, 2019.

[Osw93] P. Oswald. On a BPX-preconditioner for P1 elements. *Computing*, 51(2):125–133, 1993.

[Owh17] H. Owhadi. Multigrid with rough coefficients and multiresolution operator decomposition from hierarchical information games. *SIAM Rev.*, 59(1):99–149, 2017.

[Wey16] M. Weymuth. *Adaptive local basis for elliptic problems with $L^\infty$-coefficients*. PhD thesis, University of Zurich, 2016.
## APPENDIX A. PROOFS OF LEMMAS 5.1 AND 5.2

*Proof of Lemma 5.1.* This proof uses classical LOD techniques and is similar to the proofs of [MP20, Thm. 4.1] and [HP22, Lem. 5.2]. For shorter notation, let us denote $\varphi := C_T v$. We define the finite element cut-off function $\eta \in W^{1,\infty}(D, [0, 1])$ such that

$$ 
\begin{array}{ll}
\eta \equiv 0 & \text{in } \mathbb{N}^{\ell-1}(T), \\
\eta \equiv 1 & \text{in } D \setminus \mathbb{N}^{\ell}(T), \\
0 \le \eta \le 1 & \text{in } R := \mathbb{N}^{\ell}(T) \setminus \mathbb{N}^{\ell-1}(T).
\end{array}
$$

with $\|\nabla \eta\|_{L^\infty(D)} \le C_\eta H^{-1}$ for some $C_\eta > 0$. We estimate

$$ (A.1) \qquad \begin{aligned} \alpha \|\nabla \varphi\|_{L^2(D \setminus \mathbb{N}^\ell(T))}^2 &\le (A\nabla \varphi, \eta\nabla \varphi) = a(\varphi, \eta\varphi) - (A\nabla \varphi, \varphi\nabla \eta) \\ &\le |a(\varphi, (1 - \mathcal{B}_H)(\eta\varphi))| + |a(\varphi, \mathcal{B}_H(\eta\varphi))| + |(A\nabla \varphi, \varphi\nabla \eta)| \end{aligned} $$

For the first term we utilize definition (5.1) and the fact that $w := (1 - \mathcal{B}_H)(\eta\varphi) \in W$ and $w \equiv 0$ in $T$. For $\ell \ge 1$, this yields

$$ a(\varphi, (1 - \mathcal{B}_H)(\eta\varphi)) = a_T(\varphi, w) = 0. $$

Using that $\mathcal{B}_H(\eta\varphi) \equiv 0$ in $D \setminus \mathbb{N}^\ell(T)$ and employing the estimates (3.5) and (3.2), we obtain for the second term in (A.1)

$$ 
\begin{aligned} |a(\varphi, \mathcal{B}_H(\eta\varphi))| &\le C_{\text{bo}}\beta H^{-1} \|\nabla\varphi\|_{L^2(R)} \|\eta\varphi\|_{L^2(R)} \\ &\le C_{\text{bo}}\beta H^{-1} \|\nabla\varphi\|_{L^2(R)} \|(1 - \Pi_H)\varphi\|_{L^2(R)} \le C_{\text{bo}}C_a\beta \|\nabla\varphi\|_{L^2(R)}^2. \end{aligned} $$

For the third term in (A.1), we obtain with similar arguments

$$ |(A\nabla\varphi, \varphi\nabla\eta)| \le C_\eta \beta H^{-1} \|\nabla\varphi\|_{L^2(R)} \|(1 - \Pi)\varphi\|_{L^2(R)} \le C_\eta C_a \beta \|\nabla\varphi\|_{L^2(R)}^2. $$

Combining the previous estimates yields

$$ 
\begin{aligned} \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^\ell(T))}^2 &\le C \|\nabla\varphi\|_{L^2(R)}^2 = C \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))} - C \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^\ell(T))}^2 \\ \Leftrightarrow (1 + C) \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^\ell(T))}^2 &\le C \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))}^2 \\ \Leftrightarrow \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^\ell(T))} &\le \sqrt{\frac{C}{1+C}} \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))}, \end{aligned} $$

where $C = C_a(C_{\text{bo}} + C_\eta)\beta\alpha^{-1}$. Defining $C_d := \frac{1}{2} |\log(\frac{C}{1+C})| > 0$ and iterating the last inequality, we finally obtain

$$ \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^\ell(T))} \le \exp(-C_d\ell) \|\nabla\varphi\|_{L^2(D)}. $$

Regarding the scaling of $C_d$ with respect to $p$, we point out that the scaling of $C$ above is mainly determined by the factor $C_a C_{\text{bo}}$ and, thus, $|\log(\frac{C}{1+C})|$ scales at most like $(p+1)^{-1}$. $\square$

For the proof of Lemma 5.2, we infer the following intermediate result.

**Lemma A.1 (Localization error of element corrections).** There exists $C_t > 0$ independent of $H, \ell, T$ such that for all $T \in \mathcal{T}_H$, $v \in H_0^1(D)$, and $\ell \in \mathbb{N}$

$$ \|\nabla(C_T^\ell v - C_T v)\|_{L^2(D)} \le C_t \exp(-C_d \ell) \|\nabla C_T v\|_{L^2(D)}, $$

where $C_d$ is the constant from Lemma 5.1 and $C_t \le C(p+1)$.
Proof. This proof is again similar to the proofs of [MP20, Cor. 4.2] and [HP22, Lem. A.1]. We abbreviate $\varphi := C_T v$, $\varphi^\ell := C_T^\ell v$ and consider an arbitrary $w \in W_T^\ell$. Using the Galerkin orthogonality, we obtain

$$ \begin{aligned} \alpha \|\nabla(\varphi^\ell - \varphi)\|_{L^2(D)}^2 &\le a(\varphi^\ell - \varphi, \varphi^\ell - \varphi) = a(\varphi^\ell - \varphi, w - \varphi) \\ &\le \beta \|\nabla(\varphi^\ell - \varphi)\|_{L^2(D)} \|\nabla(w - \varphi)\|_{L^2(D)}. \end{aligned} $$

We define a particular function $w := (1 - B_H)(\eta\varphi)$ with the finite element cut-off function $\eta \in W^{1,\infty}(D, [0, 1])$ such that

$$ \begin{array}{ll} \eta = 1 & \text{in } \mathbb{N}^{\ell-1}(T), \\ \eta = 0 & \text{in } D \setminus \mathbb{N}^{\ell}(T), \\ 0 \le \eta \le 1 & \text{in } R := \mathbb{N}^{\ell}(T) \setminus \mathbb{N}^{\ell-1}(T) \end{array} $$

and $\|\nabla\eta\|_{L^\infty(D)} \le C_\eta H^{-1}$. Since $\varphi \in W$ and $w \in W_T^\ell$, we obtain using the estimates (3.5) and (3.2)

$$ \begin{aligned} \|\nabla(\varphi^\ell - \varphi)\|_{L^2(D)} &\le \beta\alpha^{-1}\|\nabla(w - \varphi)\|_{L^2(D)} \\ &= \beta\alpha^{-1}\|\nabla((1 - B_H)(1 - \eta)\varphi)\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))} \\ &\le (1 + C_a(C_\eta + C_{bo}))\beta\alpha^{-1}\|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))} \\ &\le C_t \exp(-C_d\ell)\|\nabla\varphi\|_{L^2(D)}, \end{aligned} $$

where $C_t := (1 + C_a(C_\eta + C_{bo}))\beta\alpha^{-1} \exp(C_d)$. In the previous estimate, we have used that

$$ \begin{aligned} &\|\nabla((1 - B_H)(1 - \eta)\varphi)\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))} \\ &\le \|\nabla((1 - \eta)\varphi)\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))} + C_{bo}H^{-1}\|(1 - \eta)\varphi\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))} \\ &(A.2) \qquad \le \|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))} + (C_\eta + C_{bo})H^{-1}\|\varphi\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))} \\ &\le (1 + C_a(C_\eta + C_{bo}))\|\nabla\varphi\|_{L^2(D \setminus \mathbb{N}^{\ell-1}(T))}, \end{aligned} $$

using once again the estimates (3.5) and (3.2). As in the proof of Lemma 5.1, we observe that the scaling of $C_t$ with respect to $p$ is determined by the scaling of $C_a C_{bo}$, which proves the stated scaling of $C_t$. $\square$

Proof of Lemma 5.2. This proof is similar to the proof of [MP20, Thm. 4.3] and [HP22, Lem. 5.4]. We denote $z := C^\ell v - C v$ and, for $T \in \mathcal{T}_H$, $z_T := C_T^\ell v - C_T v$. It holds

$$ \alpha \|\nabla z\|_{L^2(D)}^2 \le a(z, z) = \sum_{T \in \mathcal{T}_H} a(z_T, z). $$

Henceforth, we fix an element $T \in \mathcal{T}_H$ and define the finite element the cut-off function $\eta \in W^{1,\infty}(D, [0, 1])$ such that

$$ \begin{array}{ll} \eta = 0 & \text{in } \mathbb{N}^{\ell}(T), \\ \eta = 1 & \text{in } D \setminus \mathbb{N}^{\ell+1}(T), \\ 0 \le \eta \le 1 & \text{in } R := \mathbb{N}^{\ell+1}(T) \setminus \mathbb{N}^{\ell}(T) \end{array} $$

and $\|\nabla\eta\|_{L^\infty(D)} \le C_\eta H^{-1}$. Using that $\text{supp}((1 - B_H)(\eta z)) \subset D \setminus \mathbb{N}^\ell(T)$ and $(1 - B_H)(\eta z) \in W$ and employing definition (5.1), we get

$$ a(z_T, (1 - B_H)(\eta z)) = -a(C_T v, (1 - B_H)(\eta z)) = 0. $$
Hence, using that $z = (1 - \mathcal{B}_H)z$ since $z \in W$, we arrive at

$$a(z_T, z) = a(z_T, z - (1 - \mathcal{B}_H)(\eta z)) = a(z_T, (1 - \mathcal{B}_H)((1 - \eta)z)).$$

Proceeding similarly as in estimate (A.2), we obtain

$$|a(z_T, (1 - \mathcal{B}_H)((1 - \eta)z))| \le C\beta \|\nabla z_T\|_{L^2(D)} \|\nabla z\|_{L^2(\mathbb{N}^{\ell+1}(T))}.$$

with $C = (1 + C_a(C_\eta + C_{bo}))$. Note that the element corrections satisfy the estimate $\|\nabla C_T v\|_{L^2(D)} \le \beta \alpha^{-1} \|\nabla v\|_{L^2(T)}$, which follows from

$$\alpha \|\nabla C_T v\|_{L^2(D)}^2 \le a(C_T v, C_T v) = a_T(v, C_T v) \le \beta \|\nabla v\|_{L^2(T)} \|\nabla C_T v\|_{L^2(D)}$$

using the definition (5.1). Making use of the above inequalities, Lemma A.1, the discrete Cauchy-Schwarz inequality and the finite finite overlap of the patches, we obtain after summation over all elements

$$\begin{aligned} \alpha \|\nabla z\|_{L^2(D)}^2 &\le \sum_{T \in \mathcal{T}_H} a(z_T, z) \le C\beta \sum_{T \in \mathcal{T}_H} \|\nabla z_T\|_{L^2(D)} \|\nabla z\|_{L^2(\mathbb{N}^{\ell+1}(T))} \\ &\le CC_t \beta \exp(-C_d \ell) \sum_{T \in \mathcal{T}_H} \|\nabla C_T v\|_{L^2(D)} \|z\|_{V(\mathbb{N}^{\ell+1}(T))} \\ &\le CC_t \beta^2 \alpha^{-1} \exp(-C_d \ell) \sum_{T \in \mathcal{T}_H} \|\nabla v\|_{L^2(T)} \|\nabla z\|_{L^2(\mathbb{N}^{\ell+1}(T))} \\ &\le CC_t C_{ol} \beta^2 \alpha^{-1} \ell^{d/2} \exp(-C_d \ell) \|\nabla v\|_{L^2(D)} \|\nabla z\|_{L^2(D)}, \end{aligned}$$

with a constant $C_{ol} > 0$, which only depends on the regularity of the mesh $\mathcal{T}_H$. The assertion now follows with the constant $C_1 = CC_t C_{ol} \beta^2 \alpha^{-2}$, which scales like $C_t$ with respect to the polynomial degree $p$. $\square$

<footnote>* INRIA, 2 RUE SIMONE IFF, 75589 PARIS, FRANCE & CERMICS, ECOLE DES PONTS, 77455 MARNE-LA-VALLEE, FRANCE</footnote>

<footnote>Email address: zhaonan.dong@inria.fr</footnote>

<footnote>† INSTITUTE OF MATHEMATICS, UNIVERSITY OF AUGSBURG, UNIVERSITÄTSSTR. 12A, 86159 AUGSBURG, GERMANY</footnote>

<footnote>Email address: moritz.hauck@uni-a.de</footnote>

<footnote>‡ INSTITUTE OF MATHEMATICS, FRIEDRICH SCHILLER UNIVERSITY JENA, ERNST-ABBE-PLATZ 2, 07743 JENA, GERMANY</footnote>

<footnote>Email address: roland.maier@uni-jena.de</footnote>