1 Introduction

Spatial networks are a useful tool for constructing simplified discrete representations of complex geometric structures. Blood vessels may, for example, be modeled as connected one dimensional line segments forming a network of nodes and edges, see [15]. Paper consists of a web of wooden fibers that may be modeled as one dimensional beams forming a network, see [22]. Also porous materials, such as sandstone, may be represented by a pore network model, see [6]. Such simplifications reduce the complexity from a full three dimensional geometry to a discrete model for which computer simulations can be performed more easily. Nevertheless, highly heterogeneous materials and complicated geometries typically cause the underlying linear systems to be poorly conditioned.

This paper considers spatial network models that can be described by weighted graph Laplacians, arising from applications modeled by elliptic partial differential equations (PDEs). There already exists a variety of well-established iterative multilevel solvers for such problems. A prominent example is algebraic multigrid, see, e.g., [20, 25, 33] and the references therein. In a purely algebraic setting, the construction of multiple discretization levels can be challenging. For spatial networks with sufficient structure, it is possible to embed the network into a domain \(\Omega \subset {\mathbb {R}}^d\) and introduce scales by interpolating the network onto a family of finite element meshes, see, e.g., [17]. On coarser scales, the spatial network behaves like a continuous material and therefore inherits many advantageous properties from the continuous setting. This can be utilized for transferring successful algorithms for solving PDEs to spatial network models. We will focus on numerical homogenization, with the goal to compute an accurate effective coarse scale model of the full network.

Numerical homogenization is about the construction of problem-adapted, optimally approximating ansatz spaces possessing almost local computable bases. Near optimal numerical homogenization is achieved by the Localized Orthogonal Decomposition (LOD) [1, 3, 18, 23, 27, 28] and Gamblets [30, 31] which construct a fixed number of basis functions per mesh entity that decay exponentially relative to the mesh. For the computation of the basis, a localization of the basis functions to element patches is performed. The number of element layers in the element patches needs to be increased logarithmically with the desired accuracy. An alternative approach is taken by the AL-basis [16] and G(ms)FEM methods [5, 29] which solve local spectral problems and construct the global ansatz space using a partition of unity. Here, the support of the basis functions is fixed by the choice of partition of unity and the number of local eigenfunctions needs to be increased logarithmically with the desired accuracy.

Recently, the Super-localized Orthogonal Decomposition (SLOD) has been proposed in [19] performing a localization of the same space as the LOD, but utilizing a novel localization strategy allowing for super-exponentially decaying basis functions. This improved localization results in smaller local patch problems for the basis computation and a sparser coarse system matrix. The method proved its effectiveness also beyond elliptic multiscale problems for example in [2, 4, 14]; see also [13].

There are some contributions specifically targeting numerical homogenization of spatial network models, see [10, 12, 21, 22]. In particular, a spatial network version of the LOD has been developed in [10]. Therein a rigorous proof of the exponential decay of the LOD basis functions in the spatial network setting is provided utilizing the techniques developed in [24] and [17].

In this paper, we develop a SLOD for spatial network models. As model problem, we consider a weighted graph Laplacian K on the spatial network \({\mathcal {G}}=({\mathcal {N}},{\mathcal {E}})\), i.e., we seek the solution to the possibly poorly conditioned linear system

$$\begin{aligned} Ku=f, \end{aligned}$$
(1.1)

where u fulfills certain homogeneous Dirichlet boundary conditions and f is a given source term. For a depiction of a possible spatial network, see Fig. 1. Let the spatial network \({\mathcal {G}}\) be embedded into a domain \(\Omega \); for the ease of presentation, let \(\Omega = [0,1]^d\). We consider coarse finite element meshes \({\mathcal {T}}_H\) of \(\Omega \) and define problem-adapted ansatz spaces by applying the solution operator of (1.1) to \({\mathcal {T}}_H\)-piecewise constants. Such ansatz spaces have uniform approximation rates, independently of the smoothness of the material data. The SLOD then identifies local \({\mathcal {T}}_H\)-piecewise constant right-hand sides that yield rapidly decaying responses under the solution operator of the problem. These responses are then localized to element patches and used as basis functions of the SLOD.

We derive an a posteriori error bound for the SLOD error in terms of the size of the element patches and the coarse mesh size which holds under mild assumptions on the underlying network. In a series of numerical experiments, we illustrate the network assumptions, the super-locality of the SLOD basis and the method’s performance in presence of high-contrast data, in particular channeled data.

The outline of the paper is as follows. In Sect. 2, we introduce the spatial network model. A coarse scale discretization together with assumptions on the underlying network is presented in Sect. 3. Section 4 then constructs a prototypical problem-adapted ansatz space. We identify rapidly decaying basis functions and perform a localization of the basis in Sect. 5. Section 6 states the SLOD method along with an a posteriori error estimate and, finally, Sect. 7 presents numerical experiments.

2 A spatial network model

We consider a connected graph \({\mathcal {G}} = ({\mathcal {N}}, {\mathcal {E}})\) of nodes \({\mathcal {N}}\) and edges

$$\begin{aligned} {\mathcal {E}} = \{ \{x, y\} :\, \text {an edge connects nodes } x,y\in {\mathcal {N}}\} \end{aligned}$$

that is embedded into the unit hyper-cube \(\Omega = [0,1]^d\), \(d \in {\mathbb {N}}\). The presented methodology can also be generalized to polygonal and polyhedral domains, however, for the ease of presentation, we only consider hyper-cubes. We write \(x \sim y\) if two nodes \(x,y\in {\mathcal {N}}\) are connected by an edge in \({\mathcal {E}}\) and denote by \(|x-y|\) the Euclidean distance between the nodes. By \({\mathcal {N}}(\omega )\), we denote all nodes that are contained in the subset \(\omega \subset \Omega \). We impose homogeneous Dirichlet boundary conditions for nodes on the boundary segment \(\Gamma \subset \partial \Omega \). For a depiction of an example of a spatial network, see Fig. 1.

Fig. 1
figure 1

Illustration of spatial network with Dirichlet nodes in \(\Gamma \) marked black

Full size image

For the subsequent presentation, we introduce function spaces on the network. Let \({{\hat{V}}}\) denote the space of all real valued functions defined on the set of nodes \({\mathcal {N}}\) and denote by

$$\begin{aligned} V = \{ v \in {{\hat{V}}}\,:\,v(x) = 0,\; x \in \Gamma \} \end{aligned}$$

the subset of \({\hat{V}}\) satisfying homogeneous Dirichlet boundary conditions on the boundary segment \(\Gamma \). Furthermore, let \(V_\omega \) and \({{\hat{V}}}_\omega \) denote the spaces of functions in V and \({{\hat{V}}}\), respectively, that are supported in the subdomain \(\omega \).

2.1 Linear operators

We define a diagonal edge length weighted mass operator \(M:{{\hat{V}}} \rightarrow {{\hat{V}}}\) by the following sum of nodal contributions

$$\begin{aligned} M:=\sum _{x\in {\mathcal {N}}}M_x\quad \text { with }\quad M_x:{{\hat{V}}} \rightarrow {{\hat{V}}},\quad (M_x v, v) :=\frac{1}{2}\sum _{y \sim x} |x - y| v(x)^2, \end{aligned}$$
(2.1)

where we denote by \((\cdot ,\cdot )\) the Euclidean inner product of the vectors of nodal values of two elements in \({{\hat{V}}}\). Note that the nodal contributions \(M_x\) are uniquely defined by their associated quadratic forms. For subdomains \(\omega \subset \Omega \), we define a local version of M as

$$\begin{aligned} M_\omega :=\sum _{x \in {\mathcal {N}}(\omega )} M_x. \end{aligned}$$

The bilinear form \({( M \cdot \,,\,\cdot )}\) is an inner product on the space \({{\hat{V}}}\) with induced norm \({| \cdot |}_{M}^2 :={( M \cdot \,,\,\cdot )}\). By restriction to \(\omega \), we obtain the semi-norm \({| \cdot |}_{M,\omega }^2 :={( M_\omega \cdot \,,\,\cdot )}\) which can be interpreted as the mass of the network in subdomain \(\omega \).

Furthermore, we define a reciprocal edge length weighted graph Laplacian \(L:{{\hat{V}}} \rightarrow {{\hat{V}}}\) by

$$\begin{aligned} L:=\sum _{x\in {\mathcal {N}}}L_x\quad \text { with }\quad L_x:{{\hat{V}}} \rightarrow {{\hat{V}}},\quad (L_x v, v) :=\frac{1}{2}\sum _{y\sim x} \frac{(v(x) - v(y))^2}{|x - y|}, \end{aligned}$$
(2.2)

where the nodal contributions \(L_x\) are symmetric and are again uniquely defined by their associated quadratic forms. A local version of L is given by

$$\begin{aligned} L_\omega :=\sum _{x \in {\mathcal {N}}(\omega )} L_x. \end{aligned}$$
(2.3)

We remark that \((L_\omega v)(x)\) is nonzero for nodes x outside of \(\omega \) that are adjacent to nodes in \(\omega \). Since the graph \({\mathcal {G}}\) is assumed to be connected, the kernel of L consists only of the globally constant functions. Hence, \({| \cdot |}_{L}^2 = (L\cdot , \cdot )\) defines a semi-norm on \({{\hat{V}}}\). By restriction to \(\omega \), we obtain \({| \cdot |}_{L,\omega }^2 :=(L_\omega \cdot , \cdot )\).

Remark 2.1

(Weighting of M and L) The non-standard weightings in (2.1) and (2.2) are chosen such that, in one dimension, the operators L and M resemble the first order finite element stiffness matrix of the Poisson equation and the corresponding lumped mass matrix, respectively. This weighting enables an analysis which is similar in style to well-established PDE analysis, see, e.g., [22].

By combining the mass matrix and the graph Laplacian, we obtain another inner product on \({{\hat{V}}}\), namely \({( (L+M)\,\cdot \,,\,\cdot )}\). For its induced norm, we write \({| \cdot |}_{V}^2 :={| \cdot |}_{M}^2+{| \cdot |}_{L}^2\) and the restriction of the norm to a subdomain \(\omega \) is defined by \({| \cdot |}_{V,\omega }^2 :={| \cdot |}_{M,\omega }^2 + {| \cdot |}_{L,\omega }^2\).

2.2 Model problem

For demonstrating the extension of the SLOD to the spatial network setting, this paper considers a weighted graph Laplacian as model problem. It is defined as follows

$$\begin{aligned} K:=\sum _{x\in {\mathcal {N}}}K_x\quad \text { with }\quad K_x:{{\hat{V}}} \rightarrow {{\hat{V}}},\quad (K_x v, v) :=\frac{1}{2}\sum _{y\sim x}\gamma _{xy} \frac{(v(x) - v(y))^2}{|x - y|}\nonumber \\ \end{aligned}$$
(2.4)

with edge weights

$$\begin{aligned} 0<\alpha \le \gamma _{xy}\le \beta < \infty \end{aligned}$$
(2.5)

determining the edge’s conductivity. For subdomains \(\omega \), local versions \(K_\omega \) of can be defined analogously to (2.3).

Given a right-hand side \(f \in {{\hat{V}}}\), we seek \(u \in V\) such that, for all \(v \in V\),

$$\begin{aligned} (K u, v) = (f, v). \end{aligned}$$
(2.6)

In the spatial network setting it is natural to weight the right-hand side with the local edge length by introducing the function

$$\begin{aligned} {{\tilde{f}}}:=M^{-1}f\quad \Longrightarrow \quad (K u, v) = (M {{\tilde{f}}},v). \end{aligned}$$

This alternative weak formulation is useful because it simplifies the presentation of the proposed method and its analysis. Note that \({{\tilde{f}}}\) is well-defined since M is diagonal with positive entries. We will use both notations, f and \({{\tilde{f}}}\), throughout the paper. The following relationship between the M-norm of \({{\tilde{f}}}\) and the negative M-norm of f holds

$$\begin{aligned} |f|_{M^{-1}}:=\sup _{|v|_M=1}(f,v)=|{{\tilde{f}}}|_M. \end{aligned}$$

The unique solvability of problem (2.6) can be ensured under the minimal requirements that the network \({\mathcal {G}}\) is connected and that there exists at least one Dirichlet node in \(\Gamma \), i.e., \({\mathcal {N}}(\Gamma ) \ne \emptyset \). Subsequently, these assumptions are always assumed to hold. Using the bounds on the edge weights (2.5) we obtain that, for all \(v \in {{\hat{V}}}\),

$$\begin{aligned} \alpha {( Lv\,,\,v )}\le {( Kv\,,\,v )}\le \beta {( Lv\,,\,v )} \end{aligned}$$
(2.7)

which, in particular implies that on V, the energy norm \(|\cdot |^2 :={( K\cdot \,,\,\cdot )}\) is equivalent to the norm \({| \cdot |}_{L}\).

2.3 Global Friedrichs’ inequality

It is possible to derive a Friedrichs’ inequality in the spatial network setting. It immediately follows from the fact that the graph \({\mathcal {G}}\) is connected and that the value of one or more nodes at the boundary is fixed to zero.

Lemma 2.2

(Friedrichs’ inequality) There exists \(C_\textrm{fr}>0\), such that, for all \(v \in V\),

$$\begin{aligned} {| v |}_{M} \le C_{\textrm{fr}} {| v |}_{L}. \end{aligned}$$
(2.8)

Note that Lemma 2.2 does not provide an explicit characterization of \(C_\textrm{fr}\) in terms of the properties of the underlying network. However, one can establish a connection to the first eigenvalue of the well-studied normalized graph Laplacian for which such explicit characterizations exist, see, e.g., [8, 9].

2.4 Stability estimates

Let us also introduce a local variant of (2.6) which, for a subdomain \(\omega \subset \Omega \) and a given \(f\in V_\omega \), seeks \(u_\omega \in V_\omega \) such that, for all \(v \in V_\omega \),

$$\begin{aligned} {( K_\omega u_\omega \,,\,v )} = {( f\,,\,v )}. \end{aligned}$$
(2.9)

Existence and uniqueness of the solution \(u_\omega \) follow since \(K_\omega \), by construction, is also a weighted graph Laplacian and \(V_\omega \subset V\). Friedrichs’ inequality then allows us to show the stability of the model problem (2.6) and its local counterpart (2.9) with respect to \({| \cdot |}_{L}\) and \({| \cdot |}_{L,\omega }\) which are norms on V and \(V_\omega \), respectively.

Lemma 2.3

(Stable solvability) The solution to (2.6) is stable in the sense that

$$\begin{aligned} {| u |}_{L} \le C_\textrm{fr}\alpha ^{-1} {| {{\tilde{f}}} |}_{M} \end{aligned}$$

with \(C_\textrm{fr}\) from Lemma 2.2. This also holds for the local version (2.9), i.e.,

$$\begin{aligned} {| u_\omega |}_{L,\omega } \le C_\textrm{fr}\alpha ^{-1} {| {{\tilde{f}}} |}_{M,\omega }. \end{aligned}$$

Proof

Using the uniform lower bound (2.5), definition (2.6) and Lemma 2.2, we obtain

$$\begin{aligned} \alpha {| u |}_{L}^2 \le {( K u\,,\,u )} = (M {{\tilde{f}}},u) \le {| {{\tilde{f}}} |}_{M} {| u |}_{M} \le C_\textrm{fr} {| {{\tilde{f}}} |}_{M} {| u |}_{L}. \end{aligned}$$

Dividing by \({| u |}_{L}\), the global stability estimate follows.

The local stability estimate can be obtained similarly noting that \(u_\omega \), in particular, is also an element of V for which we can apply Lemma 2.2. Using that \(M = M_\omega \) on \(V_\omega \), \({| u_\omega |}_{M} = {| u_\omega |}_{M,\omega }\) and \({| u_\omega |}_{L} = {| u_\omega |}_{L,\omega }\), the local stability estimate can be concluded. \(\square \)

For later use, we define the solution operator \({\mathcal {K}}^{-1}:{{\hat{V}}} \rightarrow V,\) \({{\tilde{f}}} \mapsto u\) mapping a right-hand side \({{\tilde{f}}}\) to its corresponding solution u solving (2.6). On subsets \(\omega \subset \Omega \), we also define the local solution operator \({\mathcal {K}}_\omega ^{-1}:{{\hat{V}}}_\omega \rightarrow V_\omega \), \({{\tilde{f}}}\mapsto u_\omega \) mapping a local right-hand side to the local solution \(u_\omega \) satisfying (2.9).

3 Coarse scale discretization and network assumptions

3.1 Coarse mesh

The proposed method constructs its problem-adapted basis functions with respect to some coarse mesh \({\mathcal {T}}_H\). For simplicity, we restrict ourselves to Cartesian meshes which we define, unlike classical textbooks on finite element theory, as

$$\begin{aligned} {\mathcal {T}}_H:=\{S_H(x) \,:\,x = (x_1, \ldots , x_d) \in \Omega \text { and } H^{-1}x_i + 1/2 \in {\mathbb {Z}} \text { for } i = 1, \dots , d\} \end{aligned}$$

with elements that are neither closed nor open but are rather defined as

$$\begin{aligned} S_H(x) = [x_1 - H/2, x_1 + H/2) \times \cdots \times [x_d - H/2, x_d + H/2). \end{aligned}$$

If \(x_i + H/2 = 1\), we replace \([x_i - H/2, x_i + H/2)\) by \([x_i - H/2, x_i + H/2]\). This definition ensures that the elements form a true partition of \(\Omega \) meaning that any point in \({\Omega }\) is contained in exactly one element (Fig. 2).

Fig. 2
figure 2

Coarse mesh \({\mathcal {T}}_{2^{-3}}\) with underlying spatial network

Full size image

We also introduce the concept of patches which is based on neighborhood relations of elements. The first order patch of a subset \(\omega \subset \Omega \) is defined as

$$\begin{aligned} {\textsf{N}}(\omega ) :=\{ x \in \Omega \,:\,\exists T \in {\mathcal {T}}_H \,:\,x \in T,\; {\overline{T}} \cap {\overline{\omega }} \ne \emptyset \}. \end{aligned}$$

By recursion, we can then define the \(\ell \)-th order patch as \({\textsf{N}}^\ell (\omega ) :={\textsf{N}}^{\ell -1}({\textsf{N}}(\omega ))\) with \({\textsf{N}}^1 ={\textsf{N}}\).

3.2 Network assumption

The proposed method relies on the assumption that the mesh \({\mathcal {T}}_H\) is coarse compared to the underlying network. On the length scale of the coarse mesh, the network then inherits many properties known from the continuous setting, e.g., an element-wise Poincaré inequality.

Hence, we restrict ourselves to meshes with mesh sizes larger than some parameter \(H_0>0\). More precisely, we consider a hierarchy of meshes

$$\begin{aligned} \{{\mathcal {T}}_H\,:\,H \in {\mathcal {H}}\} \end{aligned}$$
(3.1)

with \({\mathcal {H}}\) denoting a finite set of positive mesh size parameters with the smallest element being \(H_0\). The parameter \(H_0\) can be interpreted as the smallest mesh size for which desired properties from the continuous case still carry over. The following assumption provides an insight into the choice of \(H_0\).

Assumption 3.1

(Network connectivity) Assume that the smallest mesh size \(H_0\) of the hierarchy of meshes (3.1) is chosen such that, for any element \(T\in {\mathcal {T}}_H\), \(H \in {\mathcal {H}}\), there is a connected subgraph \(\bar{{\mathcal {G}}} = (\bar{{\mathcal {N}}}, \bar{{\mathcal {E}}})\) of \({\mathcal {G}}\) that contains

  • all edges with one or both endpoint in T,

  • only edges with endpoints contained within \({\textsf{N}}(T)\).

3.3 Element-wise Poincaré inequality

Under the previous assumption, it is possible to prove the following element-wise Poincaré inequality.

Lemma 3.2

(Element-wise Poincaré inequality) Let Assumption 3.1 be satisfied. Then, for all \(T \in {\mathcal {T}}_H\), \(H \in {\mathcal {H}}\), there exists a constant \(C_\textrm{po}>0\) such that, for all \(v \in {{\hat{V}}}\), there exists a constant c such that

$$\begin{aligned} |v - c|_{M,T} \le C_\textrm{po} |v|_{L,{\textsf{N}}(T)}. \end{aligned}$$
(3.2)

The constant c resembles the element-average in the continuous case.

Proof

The proof can be found in [17, Lemma 3.5]. For the sake of completeness, it is also presented here. Let \({{\bar{M}}}\), \({{\bar{L}}}\) be the operators defined on the subgraph \(\bar{{\mathcal {G}}}\) from Assumption 3.1. We consider the generalized eigenvalue problem \({{\bar{L}}} v = \lambda {{\bar{M}}} v\) posed in the space \({{\hat{V}}}(\bar{{\mathcal {N}}})\) with \(\lambda _1 \le \lambda _2 \le \dots \) denoting its eigenvalues. Due to Assumption 3.1, the subgraph \(\bar{{\mathcal {G}}}\) is connected and we have that \(\lambda _1 = 0\) (corresponding to constant eigenfunctions) and \(\lambda _2>0\). By the min–max theorem, the second eigenvalue admits the characterization

$$\begin{aligned} \lambda _2 = \inf _{\begin{array}{c} 0\ne v \in V\\ {( {{\bar{M}}} v\,,\,1 )} = 0 \end{array}} \frac{{( {{\bar{L}}} v\,,\,v )}}{{( {{\bar{M}}} v\,,\,v )}}. \end{aligned}$$

Denoting with c the \({{\bar{M}}}\)-orthogonal projection of v onto the constant functions, we obtain

$$\begin{aligned} |v-c|_{M,T}^2 \le (\bar{M}(v-c),v-c) \le \lambda _2^{-1} (\bar{L} v,v) \le =\lambda _2^{-1} {| v |}_{L,{\textsf{N}}(T)}^2 \end{aligned}$$

which is the desired inequality. \(\square \)

This lemma states, similarly as Lemma 2.2, only the existence of a constant \(C_\textrm{po}\) and not its qualitative behavior in dependence of H. In Sect. 7.1, a numerical experiment confirms that \(C_\textrm{po}\) is proportional to H provided that the considered meshes are sufficiently coarse, see Assumption 3.1. Theoretically, such a scaling can be proved under the assumption of a d dimensional isoperimetric inequality, which we will elaborate in the following. Let us denote by \({\bar{d}}_x\) the degree (number of connected edges) of the node x in the subgraph \(\bar{{\mathcal {G}}}\). Further for any node subset \(X \subset \bar{{\mathcal {N}}}\), we define

$$\begin{aligned} {\text {vol}}(X) :=\sum _{x \in X}{\bar{d}}_x. \end{aligned}$$

The set of subgraph edges with one endpoint in X and the other one in \(X^\prime \) is denoted by \(\bar{{\mathcal {E}}}(X,X^\prime )\subset \bar{{\mathcal {E}}}\).

Assumption 3.3

(Isoperimetric inequality) Assume that \(\max _{(x,y)\in {\mathcal {E}}} |x-y|\le H_0\), which can be ensured by adding additional nodes to the network. Further, we assume that there exist constants \(\nu _1, \nu _2>0\) such that for any \(T \in {\mathcal {T}}_H\), \(H \in {\mathcal {H}}\), there is a subgraph \(\bar{{\mathcal {G}}}\) (see Assumption 3.1) such that it holds

$$\begin{aligned} {\text {vol}}(\bar{{\mathcal {N}}})\le \nu _1 \left( \frac{H}{H_0}\right) ^d \end{aligned}$$

and the following isoperimetric inequality is satisfied

$$\begin{aligned} {\text {vol}}(X)^{(d-1)/d}\le \nu _2\,|\bar{{\mathcal {E}}}(X,\bar{{\mathcal {N}}}\backslash X)| \end{aligned}$$

for all \(X \subset \bar{{\mathcal {N}}}\) assuming that \({\text {vol}}(X)\le {\text {vol}}(\bar{{\mathcal {N}}}\backslash X)\).

Lemma 3.4

(Poincaré constant) Let Assumptions 3.1 and 3.3 be satisfied. Then there exists \(\mu > 0\) independent of H such that the constant \(C_\textrm{po}\) from Lemma 3.2 satisfies

$$\begin{aligned} C_\textrm{po} \le \mu H. \end{aligned}$$

Proof

The proof can be found in [17, Lemma 3.6]. For the sake of completeness, it is also presented here. Using the min–max characterization of the second eigenvalue and the definition of the weighted graph Laplacian with respect to \(\bar{{\mathcal {G}}}\), we obtain that

$$\begin{aligned} \lambda _2 = \inf _{\begin{array}{c} 0\ne v \in V\\ {( {{\bar{M}}} v\,,\,1 )} = 0 \end{array}} \frac{{( {{\bar{L}}} v\,,\,v )}}{{( {{\bar{M}}} v\,,\,v )}}\ge H_0^{-2}\inf _{\begin{array}{c} 0\ne v \in V\\ {( {{\bar{M}}} v\,,\,1 )} = 0 \end{array}} \frac{\sum _{x \in \bar{{\mathcal {N}}}} \sum _{y \sim x} (v(x)-v(y))^2}{\sum _{x \in \bar{{\mathcal {N}}}}{{\bar{d}}}_xv(x)^2} = H_0^{-2} {{\hat{\lambda }}}_2. \end{aligned}$$

Here \(y \sim x\) means that \(x,y\in \bar{{\mathcal {N}}}\) are adjacent with respect to \(\bar{{\mathcal {G}}}\) and \({{\hat{\lambda }}}_2\) denotes the second eigenvalue of the normalized graph Laplacian of \(\bar{{\mathcal {G}}}\). Using the estimate \({{\hat{\lambda }}}_2 \ge C_{\nu _2}{\text {vol}}(\bar{{\mathcal {G}}})^{-2/d}\) from [9, Theorem 4], we obtain that

$$\begin{aligned} C_\textrm{po} = \lambda _2^{-1/2} \le H_0{{{\hat{\lambda }}}_2}^{-1/2} \le C_{\nu _2}^{-1/2}H_0{\text {vol}}(\bar{{\mathcal {G}}})^{1/d}\le C_{\nu _2}^{-1/2} \nu _1^{1/d}H. \end{aligned}$$

Thus, the result follows with the constant \(\mu :=C_{\nu _2}^{-1/2} \nu _1^{1/d}\). \(\square \)

If a graph has isoperimetric dimension d it essentially means that the graph mimics a d dimensional manifold. The graphs we consider in the numerical examples mimic two dimensional materials, i.e. \(d=2\).

4 Prototypical approximation

In this section, we construct prototypical problem-adapted ansatz spaces with uniform approximation properties, independent of the material data. The word prototypical emphasizes that, without modification, the constructed ansatz spaces are not feasible in a practical implementation.

A common technique for the construction of such problem-adapted ansatz spaces is the application of the inverse operator \({\mathcal {K}}^{-1}\) to some classical finite element spaces. Here, we adapt this technique to the spatial network setting. Let us consider the space of \({\mathcal {T}}_H\)-piecewise constants given by

$$\begin{aligned} {\mathbb {P}}^0({\mathcal {T}}_H)={\text {span}}\{{\textbf{1}}_T\,:\,T\in {\mathcal {T}}_H\} \subset {{\hat{V}}}, \end{aligned}$$

with the indicator function equaling one for \(x \in {\mathcal {N}}\) that are contained in T and zero otherwise.

4.1 \(L^2\)-type projection

Let us denote by \(\Pi _H:{{\hat{V}}} \rightarrow {\mathbb {P}}^0({\mathcal {T}}_H)\) the \((\cdot ,\cdot )_M\)-orthogonal projection onto \({\mathbb {P}}^0({\mathcal {T}}_H)\). It can be characterized by

$$\begin{aligned} \Pi _Hv = \sum _{T\in {\mathcal {T}}_H}\frac{{( M_T v\,,\,1 )}}{{| 1 |}_{M,T}^2}{\textbf{1}}_T. \end{aligned}$$
(4.1)

The following lemma states global stability and approximation estimates for \(\Pi _H\).

Lemma 4.1

(Properties of \(\Pi _H\)) The projection \(\Pi _H\) is stable, i.e., for all \(v \in {{\hat{V}}}\), it holds that

$$\begin{aligned} {| \Pi _Hv |}_{M}\le {| v |}_{M}. \end{aligned}$$
(4.2)

Further, if Assumptions 3.1 and 3.3 are satisfied, then there is a constant \(C_\Pi >0\) independent of H such that the following approximation estimate holds, for all \(v \in {{\hat{V}}}\),

$$\begin{aligned} {| v-\Pi _Hv |}_{M} \le C_\Pi H {| v |}_{L}. \end{aligned}$$
(4.3)

Proof

The stability estimate (4.2) is a standard property of orthogonal projections. For proving (4.3), we again split the norm into a sum of element contributions and employ (3.2) and Lemma 3.4 to obtain that

$$\begin{aligned} {| v-\Pi _Hv |}_{M}^2 = \sum _{T\in {\mathcal {T}}_H} {| v-\Pi _Hv |}_{M,T}^2\le \mu ^2 H^2 \sum _{T\in {\mathcal {T}}_H}{| v |}_{L,{\textsf{N}}(T)}^2 \le 3^d\mu ^2H^2{| v |}_{L}^2, \end{aligned}$$

where the constant \(3^d\) reflects the finite overlap of the patches \({\textsf{N}}(T)\). \(\square \)

4.2 Prototypical method

We define the prototypical problem-adapted ansatz space as

$$\begin{aligned} V_H :={\mathcal {K}}^{-1}{\mathbb {P}}^0({\mathcal {T}}_H). \end{aligned}$$
(4.4)

The corresponding Galerkin method seeks a discrete approximation \(u_H\in V_H\) such that, for all \(v_H\in V_H\),

$$\begin{aligned} {( K u_H\,,\,v_H )} = {( f\,,\,v_H )}. \end{aligned}$$
(4.5)

When using problem-adapted ansatz spaces as \(V_H\), the approximation problem of the solution u in \(V_H\) is transformed into an approximation problem of the right-hand side \({{\tilde{f}}}=M^{-1} f\) in \({\mathbb {P}}^0({\mathcal {T}}_H)\). This allows to show the desired uniform approximation properties without pre-asymptotic effects.

Lemma 4.2

(Uniform approximation) Let the network satisfy Assumptions 3.1 and 3.3. Then, for any \({{\tilde{f}}} \in {{\hat{V}}}\), the prototypical approximation \(u_H\) defined in (4.5) converges quadratically in H, i.e.,

$$\begin{aligned} {| u-u_H |}_{L} \le C_\Pi \alpha ^{-1}H {| {{\tilde{f}}}-\Pi _H{{\tilde{f}}} |}_{M} \le C_\Pi ^2 \alpha ^{-1}H^2{| {{\tilde{f}}} |}_{L}. \end{aligned}$$
(4.6)

Proof

This is the spatial network counterpart of [19, Lemma 3.2]. For the error estimate, we introduce \({{\bar{u}}}_H :={\mathcal {K}}^{-1}\Pi _H{{\tilde{f}}}\in V_H\) and employ Céa’s lemma for estimating the energy approximation error of \(u_H\) against the one for \({{\bar{u}}}_H\), i.e.,

$$\begin{aligned} {| u-u_H |}_{K} \le {| u-{{\bar{u}}}_H |}_{K}. \end{aligned}$$

Denoting \(e :=u-{{\bar{u}}}_H\), we further obtain using (2.7) and (4.3) that

$$\begin{aligned} {| e |}_{K}^2&= {( Ke\,,\,e )} = {( M({{\tilde{f}}}-\Pi _H{{\tilde{f}}})\,,\,e )} \le {( M({{\tilde{f}}}-\Pi _H{{\tilde{f}}})\,,\,e - \Pi _He )}\\&\le {| {{\tilde{f}}}-\Pi _H{{\tilde{f}}} |}_{M}{| e - \Pi _He |}_{M} \le C_\Pi \alpha ^{-1/2} H{| {{\tilde{f}}}-\Pi _H{{\tilde{f}}} |}_{M} |e|_K. \end{aligned}$$

After a division by \({| e |}_{K}\), inferring that

$$\begin{aligned} {| {{\tilde{f}}}-\Pi _H{{\tilde{f}}} |}_{M} \le C_\Pi H {| {{\tilde{f}}} |}_{L}, \end{aligned}$$

the assertion follows. \(\square \)

Remark 4.3

(\({\mathcal {T}}_H\)-piecewise right-hand sides) For \({{\tilde{f}}} \in {\mathbb {P}}^0({\mathcal {T}}_H)\), the prototypical method is exact, as for \({\mathcal {T}}_H\)-piecewise constant right-hand sides, it holds \({| {{\tilde{f}}}-\Pi _H{{\tilde{f}}} |}_{M} = 0\) in (4.6).

5 Localization

This section provides a localization strategy that turns the prototypical multiscale method introduced in the previous section into a practical method. Inspired by [19], we introduce a localization strategy for spatial network models that identifies local \({\mathcal {T}}_H\)-piecewise constant source terms with rapidly decaying responses under the solution operator \({\mathcal {K}}^{-1}\). This rapid decay enables a localization of the basis functions to element patches, which paves the way to an efficiently computable localized ansatz space.

In order to simplify the notation in the subsequent derivation, we fix an element \(T \in {\mathcal {T}}_H\) and its surrounding \(\ell \)-th order patch \(\omega :={\textsf{N}}^\ell (T)\). Furthermore, let \({\mathcal {T}}_{H,\omega }\) denote the submesh of \({\mathcal {T}}_H\) consisting of elements contained in \(\omega \).

5.1 Localization ansatz

For the construction of an almost local basis of the prototypical ansatz space \(V_H\) defined in (4.4), we assign one basis function to each element \(T \in {\mathcal {T}}_H\). The (in general global) basis function \(\psi \in V_H\) associated to element T is determined by the following ansatz

$$\begin{aligned} \psi = {\mathcal {K}}^{-1} g\quad \text { with }\quad g :=\sum _{T \in {\mathcal {T}}_{H,\omega }} g_T{\textbf{1}}_T\in {\mathbb {P}}^0({\mathcal {T}}_{H,\omega }) \end{aligned}$$
(5.1)

with coefficients \((g_T)_{T\in {\mathcal {T}}_H}\) to be determined subsequently. A local counterpart of \(\psi \), which is supported on the patch \(\omega \), can be defined by applying the patch-local solution operator \({\mathcal {K}}_\omega ^{-1}\) instead of \({\mathcal {K}}^{-1}\), i.e.,

$$\begin{aligned} \varphi :={\mathcal {K}}_\omega ^{-1} g. \end{aligned}$$

In general, the locally supported function \(\varphi \) is a poor approximation of its global counterpart \(\psi \). However, we aim to choose the coefficients \((g_T)_{T\in {\mathcal {T}}_H}\) such that \(\psi \) is well approximated by its patch-local counterpart \(\varphi \).

5.2 Conormal derivatives

For this purpose, we transfer the concept of conormal derivatives to the spatial network setting. Let \({{\tilde{V}}}_\omega \) denote the subset of V consisting of functions that are supported on nodes in \(\omega \) or its neighboring nodes.

The following preliminary result is inevitable for the definition of conormal derivatives.

Lemma 5.1

(Inner product on \({{\tilde{V}}}_\omega \)) The bilinear form

$$\begin{aligned} {( (L_\omega + M_\omega )\, \cdot \,,\,\cdot )} \end{aligned}$$

is an inner product on the space \({{\tilde{V}}}_\omega \).

Proof

We only need to show the positivity of the inner product which we do by contradiction. Assume that there exists \(0 \ne v \in {{\tilde{V}}}_\omega \) such that \({( (L_\omega +M_\omega ) v\,,\,v )} = 0\). Let \(\tilde{{\mathcal {G}}}\) denote the subgraph of \({\mathcal {G}}\) within \(\omega \) extended by its neighboring nodes and the respective edges. We can decompose \(\tilde{{\mathcal {G}}}\) into a finite number of connected components \(\tilde{{\mathcal {G}}_{i}}\). For each connected component, we have that v equals a constant \(c_i\) (otherwise \((L_\omega v,v)>0\)). Denoting with \(M_i\) the mass matrix with respect to \(\tilde{{\mathcal {G}}_{i}}\), we have that \({( M_i c_i\,,\,c_i )} = c_i^2{( M_i 1\,,\,1 )}\) which is only zero if \(c_i = 0\). We remark that a connected component cannot have nodes solely outside of \(\omega \) since, by definition, these nodes are connected to nodes within \(\omega \). The assertion follows immediately. \(\square \)

In the spatial network setting, we define the conormal derivative of \(\varphi \), denoted by \(B_\omega \varphi \in {{\tilde{V}}}_\omega ^\prime \), as the functional that satisfies, for all \(v \in {{\tilde{V}}}_\omega \),

$$\begin{aligned} {( B_{\omega } \varphi \,,\,v )}\ :={( K \varphi \,,\,v )} - {( Mg\,,\,v )}. \end{aligned}$$
(5.2)

Further, we define its operator norm as

$$\begin{aligned} |B_{ \omega } \varphi |_{{{\tilde{V}}}_{ \omega }^\prime } :=\sup _{v \in {{\tilde{V}}}_{ \omega }} \frac{{( B_{ \omega } \varphi \,,\,v )}}{{\left| v \right| }_{V, \omega }}, \end{aligned}$$

where we recall that \({\left| v \right| }_{V, \omega }^2 = |v|_{M,\omega }^2+|v|_{L,\omega }^2\) for any \(v \in {{\tilde{V}}}_\omega \). The following lemma shows that the operator norm of \(B_\omega \varphi \) can be bounded in terms of \(\varphi \) and the corresponding right-hand side g.

Lemma 5.2

(Continuity of conormal derivative) The operator norm of \(B_\omega \varphi \) can be bounded as follows

$$\begin{aligned} |B_\omega \varphi |_{{{\tilde{V}}}_\omega ^\prime } \le \max \{2\beta ,1\}\sqrt{{| \varphi |}_{L,\omega }^2 + {| g |}_{M,\omega }^2}. \end{aligned}$$

Proof

We denote by \(\tilde{{\mathcal {G}}} = (\tilde{{\mathcal {N}}},\tilde{{\mathcal {E}}})\) the subgraph of \({\mathcal {G}}\) within \(\omega \) extended by its neighboring nodes and the respective edges. We define a semi-norm on \({{\tilde{V}}}_\omega \) by

$$\begin{aligned} |v|_{{{\tilde{L}}},\omega }^2 :=\sum _{x \in \tilde{{\mathcal {N}}}} \frac{1}{2}\sum _{\tilde{{\mathcal {N}}} \ni y\sim x} \frac{(v(x)-v(y))^2}{|x-y|} \end{aligned}$$

which is equivalent to \({| \cdot |}_{L,\omega }\). More precisely, for all \(v \in {{\tilde{V}}}_\omega \), it holds \({| v |}_{L,\omega } \le |v|_{{{\tilde{L}}},\omega }\) and

$$\begin{aligned} |v|_{{{\tilde{L}}},\omega }^2 \le 2 \sum _{x \in {\mathcal {N}}(\omega )} \frac{1}{2}\sum _{y\sim x} \frac{(v(x)-v(y))^2}{|x-y|} = 2 {| v |}_{L,\omega }^2. \end{aligned}$$

Using that \(g \in {{\hat{V}}}_\omega \), \(\varphi \in V_\omega \), and \(v \in {{\tilde{V}}}_\omega \), the norm equivalence, and the discrete Cauchy–Schwarz inequality, we obtain

$$\begin{aligned} {( B_\omega \varphi \,,\,v )}&= {( K \varphi \,,\,v )}-{( M g\,,\,\omega )} \le \beta |\varphi |_{{{\tilde{L}}},\omega }|v|_{{{\tilde{L}}},\omega } + {| g |}_{M,\omega }{| v |}_{M,\omega } \\&\le 2\beta {| \varphi |}_{L,\omega }{| v |}_{L,\omega } + {| g |}_{M,\omega }{| v |}_{M,\omega }\\&\le \max \{2\beta ,1\}{| v |}_{V,\omega }\sqrt{{| \varphi |}_{L,\omega }^2 + {| g |}_{M,\omega }^2} \end{aligned}$$

which is the boundedness of the conormal derivative. \(\square \)

The following result enables a computation of the \({{\tilde{V}}}_{ \omega }^\prime \)-norm of the conormal derivative and is key for the method’s practical implementation.

Lemma 5.3

(Computation of the conormal derivative’s norm) The operator norm of the conormal derivative \(B_{ \omega } \varphi \in {{\tilde{V}}}_{ \omega }^\prime \) can be computed as

$$\begin{aligned} |B_{ \omega } \varphi |_{{{\tilde{V}}}_{ \omega }^\prime } = {| \tau |}_{V,\omega }, \end{aligned}$$

where \(\tau \in {{\tilde{V}}}_\omega \) solves, for all \(v \in {{\tilde{V}}}_{ \omega }\),

$$\begin{aligned} {( (L_{ \omega } + M_\omega ) \tau \,,\,v )} = {( B_{ \omega } \varphi \,,\,v )}. \end{aligned}$$
(5.3)

Proof

Lemma 5.1 yields the unique existence of a solution \(\tau \) to problem (5.3). Furthermore, using (5.3), we can compute the operator norm of \(B_\omega \varphi \) as follows

$$\begin{aligned} |B_{ \omega } \varphi |_{{{\tilde{V}}}_{ \omega }^\prime } :=\sup _{v \in {{\tilde{V}}}_{ \omega }} \frac{{( B_{ \omega } \varphi \,,\,v )}}{{\left| v \right| }_{V, \omega }} = \sup _{v \in {{\tilde{V}}}_{ \omega }} \frac{{( L_{ \omega } \tau \,,\,v )} + {( M_{ \omega } \tau \,,\,v )}}{{| v |}_{V, \omega }} = {| \tau |}_{V, \omega }, \end{aligned}$$

where we used that the supremum is attained for \(v = \tau \). \(\square \)

5.3 Choice of local basis

It turns out that \(\varphi \) approximates \(\psi \) well provided that the parameters \((g_T)_{T \in {\mathcal {T}}_H}\) are chosen such that the conormal derivative of \(\varphi \) is small, since, for all \(v \in V\), it holds

$$\begin{aligned} {( K(\varphi -\psi )\,,\,v )} = {( K \varphi \,,\,v )}-{( M g\,,\,v )} = {( B_{ \omega } \varphi \,,\,v )}. \end{aligned}$$

Hence, denoting by \({\mathcal {R}}:{\mathbb {P}}^0({\mathcal {T}}_{H,\omega }) \rightarrow {{\tilde{V}}}_\omega \), \(g\mapsto \tau \) the linear mapping that maps the right-hand side g to its corresponding function \(\tau \) defined in (5.3), we choose g as an \({| \cdot |}_{M,\omega }\)-normalized minimizer of the following quadratic constraint minimization problem

$$\begin{aligned} g \in \mathop {\textrm{argmin}}\limits _{q \in {\mathbb {P}}^0({\mathcal {T}}_{H,\omega })}\frac{{| {\mathcal {R}}q |}_{V,\omega }^2}{{| q |}_{M,\omega }^2}. \end{aligned}$$
(5.4)

Due to Lemma 5.3, this is equivalent to the minimization of the conormal derivative. For a depiction of selected basis functions and their corresponding right-hand sides, we refer to Fig. 3.

Fig. 3
figure 3

Basis functions \(\varphi \) for \(\ell = 1\) and several coarse mesh sizes H (top) and corresponding right-hand sides g (bottom)

Full size image

Remark 5.4

(Numerical implementation) Numerically, instead of solving (5.4), one can equivalently solve for the eigenvector corresponding to the smallest eigenvalue of the low dimensional generalized eigenvalue problem

$$\begin{aligned} Ax = \lambda Cx \end{aligned}$$
(5.5)

with matrices \(A,C \in {\mathbb {R}}^{N\times N}\), \(N :=\# {\mathcal {T}}_{H,\omega }\), defined as

$$\begin{aligned} A_{ij} = {( (L_\omega + M_\omega ) {\mathcal {R}}{\textbf{1}}_{T_j}\,,\,{\mathcal {R}}{\textbf{1}}_{T_i} )},\quad C_{ij} :={( M_\omega {\textbf{1}}_{T_j}\,,\,{\textbf{1}}_{T_i} )}, \end{aligned}$$

where \(\{T_i\,:\,i = 1,\dots N\}\) is some numbering of the elements in \({\mathcal {T}}_{H,\omega }\).

As the notation in (5.4) already indicates, the solution to the minimization problem is, in general, non-unique. Indeed, for large oversampling parameters \(\ell \) and for patches \(\omega \) close to the boundary \(\partial \Omega \), there might by several (linearly independent) g with similar in size Rayleigh quotients. In such cases, the appropriate choice of g can be difficult, cf. Section 6.1 for an implementation resolving this issue in practice.

5.4 Localization error

We define the quantity \(\sigma _T\) as the norm of the conormal derivative of the selected basis function

$$\begin{aligned} \sigma _T = \sigma _T(H,\ell ) = {| {\mathcal {R}}g |}_{V,\omega }. \end{aligned}$$
(5.6)

The quantity \(\sigma _T\) determines the local localization error, i.e., the error when approximating \(\psi \) by its local counterpart \(\varphi \).

6 Super-localized approximation

Using the localized basis function introduced in the previous section, this section transforms the prototypical method (4.5) into a practically feasible method. For a fixed oversampling parameter, we define the localized ansatz space as

$$\begin{aligned} V_{H,\ell } :=\textrm{span}\{\varphi _{T,\ell }\,:\,T\in {\mathcal {T}}_H\}\subset V. \end{aligned}$$

Note that, in the case of ambiguity, we write \(\varphi _{T,\ell }\), \(\psi _{T,\ell }\), and \(g_{T,\ell }\) for \(\varphi \), \(\psi \), and g in order to emphasize their dependence on \(T,\ell \).

The SLOD determines the Galerkin approximation of the solution u to (2.6) in the space \(V_{H,\ell }\), i.e., it seeks \(u_{H,\ell }\in V_{H,\ell }\) such that, for all \(v_{H,\ell }\in V_{H,\ell }\),

$$\begin{aligned} {( K u_{H,\ell }\,,\,v_{H,\ell } )} = {( f\,,\,v_{H,\ell } )}. \end{aligned}$$
(6.1)

6.1 Riesz basis property of right-hand sides

The choice of right-hand sides (5.4), in general, does not guarantee their linear independence. For the analysis, we assume that the right-hand sides \(\{g_{T,\ell }\,:\,T\in {\mathcal {T}}_H\}\) span \({\mathbb {P}}^0({\mathcal {T}}_H)\) in a stable way. Subsequently, we also present an implementation strategy that ensures the basis stability in practice.

Assumption 6.1

(Riesz stability) The set \(\{g_{T,\ell }\,:\,T\in {\mathcal {T}}_H\}\) is a Riesz basis of \({\mathbb {P}}^0({\mathcal {T}}_H)\), i.e., there is \(C_{\textrm{r}}(H,\ell )>0\) depending polynomially on \(H, \ell \) such that, for all \((c_T)_{T \in {\mathcal {T}}_H}\),

$$\begin{aligned} C^{-1}_{\textrm{r}}(H,\ell )\sum _{T\in {\mathcal {T}}_H} c_T^2 \le \Big \vert \sum _{T\in {\mathcal {T}}_H} c_Tg_{T,\ell }\Bigl \vert _{M}^2 \;\le C_{\textrm{r}}(H,\ell ) \sum _{T\in {\mathcal {T}}_H} c_T^2. \end{aligned}$$
(6.2)

The Riesz basis property is closely related to the eigenvalues of the matrix containing inner products of the right-hand sides \(g_{T,\ell }\) as the following remark shows.

Remark 6.2

(Riesz constant) The constant \(C_\textrm{r}\) is determined by the smallest and largest eigenvalue of the matrix \(G \in {\mathbb {R}}^{m\times m}\), \(m:=\#{\mathcal {T}}_H\) defined as

$$\begin{aligned} G_{ij} = {( Mg_{T_j,\ell }\,,\,g_{T_i,\ell } )}, \end{aligned}$$

where \(\{T_i\,:\,i = 1,\dots ,m\}\) is some numbering of the elements in \({\mathcal {T}}_H\). Denoting its eigenvalues by \(\lambda _1\le \lambda _2 \le \dots , \le \lambda _m\), the constant in the lower (resp. upper) bound of (6.2) is then the smallest (resp. largest) eigenvalue of G. Thus, we can set

$$\begin{aligned} C_\textrm{r} = \max \{\lambda _m,\lambda _1^{-1}\}. \end{aligned}$$

The stability issue of the basis of right-hand sides similarly appears in the continuous setting for the SLOD from [19]. Therein, an algorithm is proposed that selects the right-hand sides in a stable and efficient manner. The algorithm is based on the observation that stability issues only occur for patches close to the boundary \(\partial \Omega \). By grouping certain troubled patches and solving one minimization problem of the type (5.4) for such a group of patches, the stability of the right-hand sides within this group of patches can be ensured by construction. As this procedure only needs to be applied close to the global boundary for certain groups of patches, the algorithm only introduces minimal communication between the patches. For more information and a detailed illustrative description of the algorithm, see [19, Appendix B].

6.2 A-posteriori error estimate

We derive an error estimate for the SLOD solution (6.1) which is explicit in the quantity

$$\begin{aligned} \sigma = \sigma (H,\ell ) :=\max _{T\in {\mathcal {T}}_H} \sigma _T(H,\ell ) \end{aligned}$$
(6.3)

with \(\sigma _T\) defined in (5.6). The quantity \(\sigma \) determines the size of the method’s localization error. Remark 6.4 presents a summary of existing decay results for \(\sigma \) in the continuous setting and the spatial network setting.

Theorem 6.3

(Uniform localized approximation) Let the network satisfy Assumptions 3.1 and 3.3. Further, let \(\{g_{T,\ell }\,:\,T\in {\mathcal {T}}_H\}\) be stable in the sense of Assumption 6.1. Then, for any \({{\tilde{f}}} \in {{\hat{V}}}\), the SLOD approximation defined in (6.1) converges quadratically in H plus a localization error, i.e., there exists \(C,C^\prime >0\) such that, for all \({{\tilde{f}}} \in {{\hat{V}}}\),

$$\begin{aligned} \begin{aligned} {| u-{u_{H,\ell }} |}_{L}&\le C\big (H {| {{\tilde{f}}}-\Pi _H {{\tilde{f}}} |}_{M} + C_{\textrm{r}}^{1/2}(H,\ell )\ell ^{d/2} \sigma (H,\ell ){| {{\tilde{f}}} |}_{M}\big )\\&\le C^\prime \big (H^2 {| {{\tilde{f}}} |}_{L} + C_{\textrm{r}}^{1/2}(H,\ell )\ell ^{d/2} \sigma (H,\ell ){| {{\tilde{f}}} |}_{M}\big ) \end{aligned} \end{aligned}$$
(6.4)

with \(C_{\textrm{r}}\) from Assumption 6.1.

Proof

Using Céa’s Lemma, we can estimate the energy error of \(u_{H,\ell }\) approximating u by the respective energy error for \(v_{H,\ell }\), for any \(v_{H,\ell } \in V_{H,\ell }\). This yields

$$\begin{aligned} {| u-u_{H,\ell } |}_{K}\le {| u-v_{H,\ell } |}_{K}. \end{aligned}$$

Next, we add and subtract \({{\bar{u}}}_H :={\mathcal {K}}^{-1}\Pi _H{{\tilde{f}}}\) and obtain with the triangle inequality

$$\begin{aligned} {| u-u_{H,\ell } |}_{K}\le {| u-{{\bar{u}}}_H |}_{K}+{| {\bar{u}}_H-v_{H,\ell } |}_{K}. \end{aligned}$$

The first term has already been estimated in the proof of Lemma 4.2. Therein, it was shown that

$$\begin{aligned} {| u-{{\bar{u}}}_H |}_{K} \le C_\Pi ^2 \alpha ^{-1/2}H^2{| {{\tilde{f}}} |}_{L}. \end{aligned}$$

For the second term, we use that the prototypical method (4.5) is exact for piecewise constant right-hand sides, i.e., in particular for \(\Pi _H{{\tilde{f}}}\). This enables us to represent \({{\bar{u}}}_H\) using the functions \(\psi _{T,\ell } = {\mathcal {K}}^{-1} g_{T,\ell }\) from (5.1) as

$$\begin{aligned} {{\bar{u}}}_H = \sum _{T\in {\mathcal {T}}_H} c_T\, \psi _{T,\ell }, \end{aligned}$$

where \((c_T)_{T \in {\mathcal {T}}_H}\) are the coefficients of the expansion of \(\Pi _H{{\tilde{f}}}\) in terms of the right-hand sides \(g_{T,\ell }\). This is possible as the \(g_{T,\ell }\) are linearly independent which is, in particular, guaranteed by Assumption 6.1. For the particular choice

$$\begin{aligned} v_{H,\ell } :=\sum _{T\in {\mathcal {T}}_H} c_T\, \varphi _{T,\ell }\in V_{H,\ell }, \end{aligned}$$

we obtain defining \(e:={{\bar{u}}}_H - v_{H,\ell }\)

$$\begin{aligned} {| e |}_{K}^2 = \sum _{T\in {\mathcal {T}}_H}c_T\, {( K(\psi _{T,\ell }-\varphi _{T,\ell })\,,\,e )}. \end{aligned}$$
(6.5)

Writing \({\textsf{N}}^\ell (T)\) instead of \(\omega \), we obtain using Lemma 5.3 and the definitions of \(\psi _{T,\ell }\), the conormal derivative \(B_{{\textsf{N}}^\ell (T)}\) in (5.2), and \(\sigma _T\) from (5.6) that

$$\begin{aligned} \begin{aligned} {( K(\psi _{T,\ell }-\varphi _{T,\ell })\,,\,e )}&= {( M g_{T,\ell }\,,\,e )} - {( K \varphi _{T,\ell }\,,\,e )}\\&= -{( B_{{\textsf{N}}^\ell (T)} \varphi _{T,\ell }\,,\,e )} = - {( B_{{\textsf{N}}^\ell (T)}\varphi _{T,\ell }\,,\,{{\tilde{e}}} )}\\&\le \sigma _T(H,\ell ){| {{\tilde{e}}} |}_{V,{\textsf{N}}^\ell (T)} = \sigma _T(H,\ell ){| e |}_{V,{\textsf{N}}^\ell (T)}. \end{aligned} \end{aligned}$$
(6.6)

Here, we extended definition (5.2) to test functions in V and use the notation \({\tilde{e}}\) for the restriction of \(e \in V\) to the space \({{\tilde{V}}}_{{\textsf{N}}^\ell (T)}\). In the last step, we utilized \({| {{\tilde{e}}}\, |}_{V,{\textsf{N}}^\ell (T)} = {| e |}_{V,{\textsf{N}}^\ell (T)}\) which holds as the norm only considers nodes in the support of \({{\tilde{e}}}\).

The combination of (6.5), (6.6), Assumption 6.1, the discrete Cauchy–Schwarz inequality, the finite overlap of the patches, and (4.2) yields

$$\begin{aligned} \begin{aligned} {| e |}_{K}^2&= \sum _{T\in {\mathcal {T}}_H}c_T\, {( K(\psi _{T,\ell }-\varphi _{T,\ell })\,,\,e )}\le \sigma (H,\ell )\sum _{T\in {\mathcal {T}}_H}c_T{| e |}_{V,{\textsf{N}}^\ell (T)} \\&\le \sigma (H,\ell )\sqrt{\sum _{T\in {\mathcal {T}}_H}c_T^2}\sqrt{\sum _{T\in {\mathcal {T}}_H} {| e |}_{V,{\textsf{N}}^\ell (T)}^2}\\&\le \sigma (H,\ell )\,C_{\textrm{r}}^{1/2}(H,\ell ){| \Pi _H{{\tilde{f}}} |}_{M}\,C_\textrm{ol}\ell ^{d/2}{| e |}_{V} \\&\le C_\textrm{ol} C_{\textrm{r}}^{1/2}(H,\ell )\ell ^{d/2}\sigma (H, \ell )\alpha ^{-1/2}\sqrt{1+C_\textrm{fr}^2}{| e |}_{K} {| {{\tilde{f}}} |}_{M} \end{aligned} \end{aligned}$$

with constant \(C_\textrm{ol}>0\) reflecting the overlap of the patches \({\textsf{N}}^\ell (T)\). In the last step, we applied Friedrichs’ inequality (2.8) on the full network. Putting together the previous estimates and using (2.5), the assertion can be concluded. \(\square \)

Remark 6.4

(Decay of localization error) In the continuous setting, [19] conjectures that \(\sigma \) decays super-exponentially in \(\ell \). More precisely, there exists \(C_\sigma >0\) depending polynomially on \(H,\ell \) and \(C_\textrm{d}>0\) independent of \(H,\ell \) such that, for all \(\ell \),

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

In [19, Section 7], this result is justified theoretically using a conjecture from spectral geometry. Numerical experiments in [19, Section 8] confirm the super-exponential decay numerically. Furthermore, [19, Lemma 6.4] provides a rigorous proof that \(\sigma \) decays at least exponentially in \(\ell \).

In the spatial network setting, we also observe a rapid decay of the quantity \(\sigma \) as \(\ell \) is increased. Qualitatively, the decay behavior is is similar to the one in the continuous setting, see the numerical experiments in Sect. 7.2 and Sect. 7.3. Using techniques from [10], one can show, similarly as in the continuous setting, that \(\sigma \) decays at least exponentially.

Remark 6.5

(\({\mathcal {T}}_H\)-piecewise right-hand sides) For \({{\tilde{f}}} \in {\mathbb {P}}^0({\mathcal {T}}_H)\), only the localization error in (6.4) is present, as for \({\mathcal {T}}_H\)-piecewise constant right-hand sides, it holds \({| {{\tilde{f}}}-\Pi _H{{\tilde{f}}} |}_{M} = 0\).

Remark 6.6

(Relation to AMG coarse spaces) To improve the convergence of the outer iteration, many AMG methods use problem-adapted coarse spaces. The basis functions of such coarse spaces are computed, e.g., by solving local eigenproblems (see, e.g., [7, 11]) or by solving an energy minimization problem subject to a partition of unity constraint (see, e.g., [26, 32]). Using the coarse space in a one-shot method, i.e. computing the Galerkin approximation in the coarse space, typically does not yield satisfactory approximation orders. In contrast, the proposed method achieves optimal convergence orders introducing a computational overhead (\(\ell \) must be increased with the desired accuracy, see Theorem 6.3).

7 Numerical experiments

For the subsequent numerical experiments, we utilize a spatial network constructed as follows. First, we sample 20,000 lines of length 0.05 which are uniformly rotated and with midpoints uniformly distributed in \([-0.025,1.025]^2\). Next, we remove all line segments outside of the unit square so that all lines are contained in the domain \(\Omega = [0,1]^2\). We then define the network nodes as the line segments’ endpoints and intersections. Two nodes are connected by an edge if the nodes share a line segment. We only consider the largest connected component of the network in order to ensure that the network is connected. We also remove all hanging nodes (nodes of degree one) in the interior of the domain along with the respective edges. All nodes at the boundary \(\partial \Omega \) are Dirichlet nodes. The total number of nodes is around 80, 000.

7.1 Poincaré constant

This numerical experiment is to justify Lemma 3.4 numerically. Given \(T \in {\mathcal {T}}_H\), \(H \in {\mathcal {H}}\), we construct the subgraph \(\bar{{\mathcal {G}}} = (\bar{{\mathcal {N}}},\bar{{\mathcal {E}}})\) by means of a breadth-first search algorithm. Then, we compute the second eigenvalue \(\lambda _2\) of the generalized eigenvalue problem \({{\bar{L}}} v = \lambda {{\bar{M}}} v\) posed in the space \({{\hat{V}}}(\bar{{\mathcal {N}}})\) with \({{\bar{L}}}\), \({{\bar{M}}}\) being defined on \(\bar{{\mathcal {G}}}\), see also the proof of Lemma 3.2.

In Fig. 4, the reciprocal square root of the numerically computed eigenvalues \(\lambda _2\) corresponding to elements \(T \in {\mathcal {T}}_H\) is depicted for different mesh sizes \(H \in {\mathcal {H}}\). The dotted black line interpolates, for all mesh sizes H, the averaged eigenvalue for the respective H. Note that the scaling of the y-axis is chosen such that the black line should appear linear if the desired scaling from Lemma 3.4 holds.

Fig. 4
figure 4

Eigenvalues \(\lambda _2\) for different mesh sizes H

Full size image

Figure 4 confirms Lemma 3.4 numerically. We note that for smaller mesh sizes, there is a considerable variation in the eigenvalues. This happens when we reach the critical mesh size, which is \(H_0\approx 2^{-5}\) in this example.

7.2 Decay of \(\sigma \)

In this experiment, we numerically investigate the decay of \(\sigma \) defined in (6.3) which is the maximum of the local quantities \(\sigma _T\) from (5.6). For illustration purposes, we pick an element \(T \in {\mathcal {T}}_{2^{-4}}\) whose fourth order patch has no intersection with the global boundary \(\partial \Omega \). Figure 5 then depicts the square root of the eigenvalues of the eigenvalue problem (5.5) for the patches \({\textsf{N}}^\ell (T)\) with \(\ell = 1,\dots ,4\). By definition, the square root of the smallest eigenvalue coincides with \(\sigma _T\). The values of \(\sigma _T\) are marked using dashed horizontal lines.

Fig. 5
figure 5

Eigenvalues of the eigenvalue problem (5.5) for patches of orders \(\ell \). The dashed lines marks the values of \(\sigma _T\) for the respective patches

Full size image

In Fig. 5, one observes a rapid decay of \(\sigma _T\) as \(\ell \) is increased. Note that for \(\ell \ge 3\), the difference in magnitude of the smallest and largest eigenvalue of (5.5) is at least of order \(10^{14}\) which means that the respective matrices have a large condition number. This affects the accuracy of the eigenvalue solver and explain the flatting of the decay for large oversampling parameters as observed in Figs. 5 and 6.

7.3 Super-localization

For this numerical experiment, we consider a weighted graph Laplacian (2.4) with edge weights \(\gamma _{xy}\) independently and uniformly distributed in the interval [0.01, 1]. We choose the right-hand side \({{\tilde{f}}} \equiv 1\) or equivalently \(f=M1\). Since \({{\tilde{f}}}\in {\mathbb {P}}^0({\mathcal {T}}_H)\), the error is bounded solely by the localization error, see Remark 6.5. In Fig. 6 (left), the localization errors for the SLOD are shown. The localization errors are plotted for several coarse grids \({\mathcal {T}}_H\) with respect to \(\ell \). We additionally depict the localization errors when using the LOD for localizing the same prototypical ansatz space \(V_H\) from (4.4). As reference, we indicate lines showing the expected rates of decay of the localization errors for the SLOD and LOD which is super-exponential decay for the SLOD, cf. (6.7), and exponential decay for the LOD. In Fig. 6 (right), we again depict the localization errors for the SLOD but this time together with the values of its estimator

$$\begin{aligned} \textrm{est}(H,\ell ) :=C_\textrm{r}^{1/2}(H,\ell )\ell ^{d/2}\sigma (H,\ell ) \end{aligned}$$
(7.1)

appearing in the error estimate in Theorem 6.3. Note that the scaling of the x-axis in Fig. 6 is chosen such that a super-exponentially decaying curve appears to be a straight line.

Fig. 6
figure 6

Plot of the relative localization errors as a function of the oversampling parameter \(\ell \) of the SLOD and the LOD (left). Depiction of the relative localization errors of the SLOD with the estimator (7.1) (right)

Full size image

Figure 6 (left) numerically confirms the super-exponential decay rates of the localization errors as known for the continuous setting [19]. Note that, similarly as for the decay of \(\sigma \) in Sect. 7.2, we can also observe a flattening of the decay for \(\ell \ge 3\). This might again be explained by the high condition number of the matrices in (5.5), see Sect. 7.2. The localization error of the LOD, depicted in Fig. 6 (left), decays exponentially, see, e.g., [1, 10]. Much larger values of \(\ell \) are necessary in order to reach the accuracy level of the SLOD.

Figure 6 (right) shows that the error estimator is quite well aligned with the localization error and thereby underlines the sharpness of the estimator. For \(\ell = 4\), we observe that also the estimator is slightly affected by the aforementioned conditioning issue.

7.4 Optimal convergence

For demonstrating the convergence of the SLOD for spatial networks, we consider the edge weights from the previous numerical example and the right-hand side

$$\begin{aligned} {{\tilde{f}}}(x_1,x_2) = \sin (x_1)\sin (x_2). \end{aligned}$$

Figure 7 shows the convergence plots for the SLOD (left) and the LOD (right) for different oversampling parameters as the coarse mesh \({\mathcal {T}}_H\) is refined. Note that we only consider combinations of H, \(\ell \) for which there is no patch \({\textsf{N}}^\ell (T)\) that coincides with the whole domain \(\Omega \). As reference, lines of slope 2 are depicted.

Fig. 7
figure 7

Plot of the errors \({| u - u_{H,\ell } |}_{L}/{| u |}_{L}\) as a function of the mesh size H of the SLOD (left) and the LOD (right)

Full size image

Figure 7 confirms the method’s convergence properties as stated in Theorem 6.3, i.e., provided that the oversampling parameter \(\ell \) is chosen sufficiently large, optimal convergence of order 2 can be observed. Note that in Fig. 7 (left), the yellow line is partially below the purple line and thus is difficult to see. For the LOD, the optimal order convergence is difficult to see. The reason is that, for all considered localization parameters, the localization error still dominates the optimal order error.

7.5 High-contrast example

In this numerical experiment, we use edge weights \(\gamma _{xy}\) that are independently and uniformly distributed in [0.01, 1] and add several channels of high conductivity with edge weights of \(10^4\). Hence, the contrast in this numerical example is of order \(10^6\). For an illustration of the setup, we refer to Fig. 8 (left), where the high-conductivity channels are marked in green. The high conductivity effectively extends the homogeneous Dirichlet boundary conditions also to the channels. We consider this experiment as particularly challenging as the channels are not aligned with the considered Cartesian meshes \({\mathcal {T}}_H\) and thus, artificial boundary-like conditions are imposed on the SLOD basis functions. As right-hand side, we choose \({{\tilde{f}}} \equiv 1\), i.e., only the localization error is present, cf. Remark 6.5.

Fig. 8
figure 8

Spatial network with high conductivity channels indicated in green (left) and SLOD solution (right)

Full size image

For this experiment, the SLOD is again able to achieve very accurate approximations for much smaller oversampling parameters than the LOD. For example, for the mesh \({\mathcal {T}}_{2^{-4}}\) and \(\ell = 3\), the SLOD achieves a relative L-norm error of \(2.75\times 10^{-3}\), while the LOD, for the same discretization parameters, only achieves an error of \(1.73 \times 10^{-1}\). For reaching a similar accuracy as the SLOD, for the LOD, we would need to choose \(\ell \) so large that many patch problems are already global problems.

Note that in Fig. 8, for the ease of illustration, we only depicted a subnetwork and the correspondingly restricted solution. Plotting the full network is infeasible as the high density of the lines would make the lines nearly impossible to distinguish.