Regularity and convergence analysis in Sobolev and H\"older spaces for generalized Whittle-Mat\'ern fields

TL;DR

This paper establishes optimal convergence rates for Galerkin approximations of generalized Whittle-Matérn Gaussian fields in Sobolev and Hölder spaces, based on minimal regularity assumptions, and validates results through numerical experiments.

Contribution

It provides the first rigorous convergence analysis for Galerkin methods approximating generalized Whittle-Matérn fields under minimal regularity conditions.

Findings

Optimal convergence rates in Sobolev and Hölder spaces.

Error estimates for covariance function approximations.

Numerical validation of theoretical convergence rates.

Abstract

We analyze several Galerkin approximations of a Gaussian random field $\mathcal{Z}\colon\mathcal{D}\times\Omega\to\mathbb{R}$ indexed by a Euclidean domain $\mathcal{D}\subset\mathbb{R}^d$ whose covariance structure is determined by a negative fractional power $L^{-2\beta}$ of a second-order elliptic differential operator $L:= -\nabla\cdot(A\nabla) + \kappa^2$. Under minimal assumptions on the domain $\mathcal{D}$, the coefficients $A\colon\mathcal{D}\to\mathbb{R}^{d\times d}$, $\kappa\colon\mathcal{D}\to\mathbb{R}$, and the fractional exponent $\beta>0$, we prove convergence in $L_q(\Omega; H^\sigma(\mathcal{D}))$ and in $L_q(\Omega; C^\delta(\overline{\mathcal{D}}))$ at (essentially) optimal rates for (i) spectral Galerkin methods and (ii) finite element approximations. Specifically, our analysis is solely based on $H^{1+\alpha}(\mathcal{D})$-regularity of the differential operator $L$, where $0<\alpha\leq 1$. For this setting, we furthermore provide rigorous estimates for the error in the covariance function of these approximations in $L_{\infty}(\mathcal{D}\times\mathcal{D})$ and in the mixed Sobolev space $H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})$, showing convergence which is more than twice as fast compared to the corresponding $L_q(\Omega; H^\sigma(\mathcal{D}))$-rate. For the well-known example of such Gaussian random fields, the original Whittle-Mat\'ern class, where $L=-\Delta + \kappa^2$ and $\kappa \equiv \operatorname{const.}$, we perform several numerical experiments which validate our theoretical results.

Figures

Tables

Table 1. Table 1. Expected rates of convergence, cf. Theorems  6.18 and  6.23

	$L_2$	$L_{\infty }$	$H_0^1$
${𝒵}_{h,k}^{\beta }$	$\min \{2\beta -1/2,p+1\}$	$\min \{2\beta -1/2,p+1\}$	$\min \{2\beta -3/2,p\}$
${\varrho }_{h,k}^{\beta }$	$\min \{4\beta -1/2,p+1\}$	$\min \{4\beta -1,p+1\}$	$\min \{4\beta -5/2,p\}$

Table 2. Table 2. Observed (resp. theoretical) rates of convergence for the errors of the field/covariance shown in Figures  1 – 2

		$\beta$ for field error studies
	$p$	$0.5$	$0.8$	$1.1$	$1.4$	$1.7$
$L_2$	$1$	0.54 (0.5)	1.10 (1.1)	1.67 (1.7)	1.94 (2)	1.96 (2)
	$2$	0.56 (0.5)	1.10 (1.1)	1.68 (1.7)	2.27 (2.3)	2.85 (2.9)
$L_{\infty }$	$1$	0.55 (0.5)	1.05 (1.1)	1.60 (1.7)	1.93 (2)	1.99 (2)
	$2$	0.68 (0.5)	1.14 (1.1)	1.67 (1.7)	2.25 (2.3)	2.79 (2.9)
$H_0^1$	$1$	–	0.22 (0.1)	0.70 (0.7)	1.00 (1)	1.05 (1)
	$2$	–	0.27 (0.1)	0.73 (0.7)	1.30 (1.3)	1.87 (1.9)
		$\beta$ for covariance error studies
	$p$	$0.5$	$0.6$	$0.7$	$0.8$	$0.9$	$1$
$L_2$	$1$	1.53 (1.5)	1.85 (1.9)	1.98 (2)	2.00 (2)	2.00 (2)	2.00 (2)
	$2$	1.57 (1.5)	1.94 (1.9)	2.32 (2.3)	2.69 (2.7)	2.94 (3)	3.00 (3)
$L_{\infty }$	$1$	1.07 (1)	1.41 (1.4)	1.72 (1.8)	1.91 (2)	1.98 (2)	1.99 (2)
	$2$	1.23 (1)	1.52 (1.4)	1.86 (1.8)	2.23 (2.2)	2.61 (2.6)	2.99 (3)

Equations

$$L^{\beta }\mathcal{Z}=\mathrm{d}\mathcal{W},\text{in }\mathcal{D},L:=-\Delta +{\kappa }^2,$$

$$Lu=-\nabla \cdot (A\nabla u)+{\kappa }^{2}u,u\in \mathcal{D}(L)\subseteq L_2(\mathcal{D}).$$

$$\displaystyle\left(\mathbb{E}\Bigl{[}\bigl{\|}\mathcal{Z}^{\beta}-\mathcal{Z}^{\beta}_{N}\bigr{\|}_{H^{\sigma}(\mathcal{D})}^{q}\Bigr{]}\right)^{\nicefrac{{1}}{{q}}}$$

$$\displaystyle\left(\mathbb{E}\Bigl{[}\bigl{\|}\mathcal{Z}^{\beta}-\mathcal{Z}^{\beta}_{N}\bigr{\|}_{C^{\delta}(\overline{\mathcal{D}})}^{q}\Bigr{]}\right)^{\nicefrac{{1}}{{q}}}$$

$$\displaystyle\bigl{\|}\varrho^{\beta}-\varrho^{\beta}_{N}\bigr{\|}_{H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})}$$

$$\displaystyle\sup_{x,y\in\overline{\mathcal{D}}}\bigl{|}\varrho^{\beta}(x,y)-\varrho^{\beta}_{N}(x,y)\bigr{|}$$

$$\beta >0,0\le \sigma \le \min \{1,2\beta \},-1\le \delta \le 1+\alpha ,\delta \ne \text{nicefrac}12\text{with}2\beta +\delta -\sigma >0,$$

$$\displaystyle\bigl{\|}L^{-\beta}g-\widetilde{L}_{N}^{-\beta}g\bigr{\|}_{H^{\sigma}(\mathcal{D})}$$

$$\Vert f{\Vert }_{L_q(\mathcal{D};E)}$$

$$\Vert X{\Vert }_{L_q(\Omega ;E)}$$

$$\Vert f{\Vert }_{L_{\infty }(\mathcal{D};E)}:=\mathrm{ess}\underset{x\in \mathcal{D}}{\sup }\Vert f(x){\Vert }_E,f\in L_{\infty }(\mathcal{D};E).$$

$$|\cdot {|}_{C^{\gamma }(\overline{\mathcal{D}};E)},\Vert \cdot {\Vert }_{C^{\gamma }(\overline{\mathcal{D}};E)}:C(\overline{\mathcal{D}};E)\to [0,\infty ]$$

$$\bigl{(}C(\overline{\mathcal{D}};E),\|\,\cdot\,\|_{C(\overline{\mathcal{D}};E)}\bigr{)},\qquad\|f\|_{C(\overline{\mathcal{D}};E)}:=\sup_{x\in\overline{\mathcal{D}}}\|f(x)\|_{E},$$

$$|f{|}_{C^{\gamma }(\overline{\mathcal{D}};E)}$$

$$\Vert f{\Vert }_{C^{\gamma }(\overline{\mathcal{D}};E)}$$

$$C^{\gamma }(\overline{\mathcal{D}};E)$$

$$\exists a_0>0:\mathrm{ess}\underset{x\in \mathcal{D}}{\inf }{\xi }^{\top }A(x)\xi \ge a_0|\xi {|}^2\forall \xi \in {\mathbb{R}}^d;$$

$$\exists a_{\mathrm{Lip}}>0:|A_{ij}(x)-A_{ij}(y)|\le a_{\mathrm{Lip}}|x-y|\forall x,y\in \overline{\mathcal{D}},$$

$$\left|\int_{\mathcal{D}}\bigl{[}(A(x)\nabla u(x),\nabla v(x))_{\mathbb{R}^{d}}+\kappa^{2}(x)u(x)v(x)\bigr{]}\,\mathrm{d}x\right|\leq C\|v\|_{L_{2}(\mathcal{D})}\quad\forall v\in H^{1}_{0}(\mathcal{D}),$$

$$\displaystyle\int_{\mathcal{D}}\bigl{[}(A(x)\nabla u(x),\nabla v(x))_{\mathbb{R}^{d}}+\kappa^{2}(x)u(x)v(x)\bigr{]}\,\mathrm{d}x=(Lu,v)_{L_{2}(\mathcal{D})}.$$

$${\lambda }_j{\eqsim }_{(A,\kappa ,\mathcal{D})}{j}^{\text{nicefrac}2d},j\in \mathbb{N}.$$

$$\mathscr{D}(L^{\sigma}):=\Biggl{\{}\psi\in L_{2}(\mathcal{D}):\sum_{j\in\mathbb{N}}\lambda_{j}^{2\sigma}(\psi,e_{j})_{L_{2}(\mathcal{D})}^{2}<\infty\Biggr{\}}$$

$$L^{\sigma }\psi :=\sum _{j\in \mathbb{N}}{\lambda }_j^{\sigma }(\psi ,e_j{)}_{L_2(\mathcal{D})}{e}_j,\psi \in \mathcal{D}(L^{\sigma }).$$

$$\displaystyle\bigl{(}\dot{H}^{\sigma}_{L},(\,\cdot\,,\,\cdot\,)_{\sigma}\bigr{)},\qquad\dot{H}^{\sigma}_{L}:=\mathscr{D}\bigl{(}L^{\nicefrac{{\sigma}}{{2}}}\bigr{)}\subset L_{2}(\mathcal{D}),$$

$$(\phi,\psi)_{\sigma}:=\bigl{(}L^{\nicefrac{{\sigma}}{{2}}}\phi,L^{\nicefrac{{\sigma}}{{2}}}\psi\bigr{)}_{L_{2}(\mathcal{D})}=\sum_{j\in\mathbb{N}}\lambda_{j}^{\sigma}(\phi,e_{j})_{L_{2}(\mathcal{D})}(\psi,e_{j})_{L_{2}(\mathcal{D})},$$

$$\bigl{(}\dot{H}^{\sigma}_{L},\|\,\cdot\,\|_{\sigma}\bigr{)}\cong\bigl{(}\left[L_{2}(\mathcal{D}),H^{1}_{0}(\mathcal{D})\right]_{\sigma},\|\,\cdot\,\|_{[L_{2}(\mathcal{D}),H^{1}_{0}(\mathcal{D})]_{\sigma}}\bigr{)},\qquad 0\leq\sigma\leq 1,$$

$$\bigl{(}\dot{H}^{\sigma}_{L},\|\,\cdot\,\|_{\sigma}\bigr{)}\hookrightarrow\bigl{(}H^{\sigma}(\mathcal{D}),\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}\bigr{)},\qquad 0\leq\sigma\leq 1,$$

$$\bigl{(}\dot{H}^{\sigma}_{L},\|\,\cdot\,\|_{\sigma}\bigr{)}\cong\bigl{(}H^{\sigma}(\mathcal{D})\cap H^{1}_{0}(\mathcal{D}),\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}\bigr{)},\qquad 1\leq\sigma\leq 2.$$

$$f\in L_2(\mathcal{D})\Rightarrow u:=L^{-1}f\in H^2(\mathcal{D})\cap H_0^1(\mathcal{D}).$$

$$(\mathcal{Z},\psi {)}_{L_2(\mathcal{D})}=\mathcal{W}(T^{\ast }\psi )\mathbb{P}\text{-a.s.}\forall \psi \in L_2(\mathcal{D}).$$

$$\mathcal{Z}(x)=\mathcal{W}(T^{\ast }{\text{dh}}_x)\mathbb{P}\text{-a.s.}\forall x\in \overline{\mathcal{D}}.$$

$$\left(\mathbb{E}\Bigl{[}\|\mathcal{Z}\|_{C^{\delta}(\overline{\mathcal{D}})}^{q}\Bigr{]}\right)^{\nicefrac{{1}}{{q}}}\lesssim_{(q,\gamma,\delta,\mathcal{D})}\|T\|_{\mathcal{L}(L_{2}(\mathcal{D});C^{\gamma}(\overline{\mathcal{D}}))}.$$

$$\displaystyle\left(\mathbb{E}\bigl{[}|\mathcal{Z}_{0}(x)-\mathcal{Z}_{0}(y)|^{2}\bigr{]}\right)^{\nicefrac{{1}}{{2}}}$$

Full text

For the well-known example of such Gaussian random fields, the original Whittle–Matérn class, where $L=-\Delta+\kappa^{2}$ and $\kappa\equiv\operatorname{const.}$, we perform several numerical experiments which validate our theoretical results.

Acknowledgment. The authors thank Christoph Schwab and Mark Veraar for helpful and valuable comments.

1. Introduction

1.1. Motivation and background

By virtue of their practicality owing to the full characterization by their mean and covariance structure, Gaussian random fields (GRFs for short) are popular models for many applications in spatial statistics and uncertainty quantification, see, e.g., [4, 7, 26, 33, 35]. As a result, several methodologies in these disciplines require the efficient simulation of GRFs at unstructured locations in various possibly non-convex Euclidean domains, and this topic has been intensively discussed in both areas, spatial statistics and computational mathematics, see, e.g., [2, 8, 9, 13, 17, 19, 24, 30]. In particular, sampling from non-stationary GRFs, for which methods based on circulant embedding are inapplicable, has become a central topic of current research, see, e.g., [2, 9, 17].

In order to capture both stationary and non-stationary GRFs, a new class of random fields has been introduced in [26], which is based on the following observation made by P. Whittle [40]: A GRF $\mathcal{Z}$ on $\mathcal{D}:=\mathbb{R}^{d}$ with covariance function of Matérn type solves the fractional-order stochastic partial differential equation (SPDE for short)

[TABLE]

where $\Delta$ denotes the Laplacian, $\mathrm{d}\mathcal{W}$ is white noise on $\mathbb{R}^{d}$, and $\kappa>0$, $\beta>\nicefrac{{d}}{{4}}$ are constants which determine the practical correlation length and the smoothness of the field. In [26] this relation has been exploited to formulate generalizations of Matérn fields, the generalized Whittle–Matérn fields, by considering the SPDE (1.1) for non-stationary differential operators $L$ (e.g., by allowing for a spatially varying coefficient $\kappa\colon\mathcal{D}\to\mathbb{R}$) on bounded domains $\mathcal{D}\subset\mathbb{R}^{d}$, $d\in\{1,2,3\}$. Note that the covariance structure of a GRF is uniquely determined by its covariance operator, in this case given by the negative fractional-order differential operator $L^{-2\beta}$. Furthermore, for the case $2\beta\in\mathbb{N}$, approximations based on a finite element discretization have been proposed in [26].

Subsequently, a computational approach which allows for arbitrary fractional exponents $\beta>\nicefrac{{d}}{{4}}$ has been suggested in [2, 3]. To this end, a sinc quadrature combined with a Galerkin discretization of the differential operator $L$ is applied to the Balakrishnan integral representation of the fractional-order inverse $L^{-\beta}$.

In this work, we investigate Sobolev and Hölder regularity of generalized Whittle–Matérn fields and we perform a rigorous error analysis in these norms for several Galerkin approximations, including the sinc-Galerkin approximations of [2, 3]. Specifically, we consider a GRF $\mathcal{Z}^{\beta}\colon\mathcal{D}\times\Omega\to\mathbb{R}$, indexed by a Euclidean domain $\mathcal{D}\subset\mathbb{R}^{d}$, whose covariance operator is given by the negative fractional power $L^{-2\beta}$ of a second-order elliptic differential operator $L\colon\mathscr{D}(L)\subseteq L_{2}(\mathcal{D})\rightarrow L_{2}(\mathcal{D})$ in divergence form with Dirichlet boundary conditions, formally given by

[TABLE]

Here, we solely assume that $\mathcal{D}\subset\mathbb{R}^{d}$ has a Lipschitz boundary, $\kappa\in L_{\infty}(\mathcal{D})$, and that A\in L_{\infty}\bigl{(}\mathcal{D};\mathbb{R}^{d\times d}\bigr{)} is symmetric and uniformly positive definite.

For a sequence \bigl{(}\mathcal{Z}_{N}^{\beta}\bigr{)}_{N\in\mathbb{N}} of Galerkin approximations for $\mathcal{Z}^{\beta}$ (namely, spectral Galerkin approximations in Section 5 and sinc-Galerkin approximations in Section 6) defined with respect to family $(V_{N})_{N\in\mathbb{N}}$ of subspaces $V_{N}\subset H^{1}_{0}(\mathcal{D})$ of finite dimension $\dim(V_{N})=N<\infty$, we prove convergence at (essentially) optimal rates. More precisely, under minimal regularity conditions on the operator $L$ in (1.2) and for $0\leq\sigma<2\beta-\nicefrac{{d}}{{2}}$, $\delta\in(0,\sigma)$, within a suitable parameter range we show that for all $\varepsilon,q>0$ there exists a constant $C>0$ such that, for all $N\in\mathbb{N}$,

[TABLE]

Here, $\varrho^{\beta},\varrho_{N}^{\beta}\colon\mathcal{D}\times\mathcal{D}\to\mathbb{R}$ denote the covariance functions of the Whittle–Matérn field $\mathcal{Z}^{\beta}$ and of the Galerkin approximation $\mathcal{Z}_{N}^{\beta}$, respectively. For details, see Corollaries 5.1–5.3 for spectral Galerkin approximations, and Theorems 6.18, 6.23 for the sinc-Galerkin approach. “Suitable parameter range” refers to the observations that

(i) if a finite element method of polynomial degree $p\in\mathbb{N}$ is used to define the sinc-Galerkin approximation or

(ii) if $L$ in (1.2) is $H^{1+\alpha}(\mathcal{D})$-regular for $0<\alpha\leq 1$ maximal (see Definition 6.20),

then the convergence rates of the sinc-Galerkin approximation cannot exceed $p+1-\sigma$ or $\min\{1+\alpha-\sigma,2\alpha\}$, where $0\leq\sigma\leq 1$.

We point out that due to the low regularity of white noise, $\mathrm{d}\mathcal{W}\in H^{-\nicefrac{{d}}{{2}}-\varepsilon}(\mathcal{D})$, which holds $\mathbb{P}$-almost surely and in $L_{q}(\Omega)$ (cf. [2, Prop. 2.3]) the convergence results (1.3)–(1.6) are (essentially, up to $\varepsilon>0$) optimal and they are also reflected in our numerical experiments, see Section 7 and the discussion in Section 8. Note furthermore that the convergence rates in (1.4), (1.6) of the field with respect to $L_{q}(\Omega;C^{\delta}(\overline{\mathcal{D}}))$ and of the covariance function in the $C(\overline{\mathcal{D}\times\mathcal{D}})$-norm, which we obtain via a Kolmogorov–Chentsov argument, are by $\nicefrac{{d}}{{2}}$ better than applying the results (1.3), (1.5) combined with the Sobolev embeddings $H^{\delta+\nicefrac{{d}}{{2}}}(\mathcal{D})\hookrightarrow C^{\delta}(\overline{\mathcal{D}})$ and $H^{\varepsilon+\nicefrac{{d}}{{2}},\,\varepsilon+\nicefrac{{d}}{{2}}}(\mathcal{D}\times\mathcal{D})\hookrightarrow C(\overline{\mathcal{D}\times\mathcal{D}})$, respectively. We remark that strong convergence of the sinc-Galerkin approximation with respect to the $L_{2}(\Omega;L_{2}(\mathcal{D}))$-norm, i.e., (1.3) for $\sigma=0$, at the rate $2\beta-\nicefrac{{d}}{{2}}$ has already been proven in [2, Thm. 2.10]. However, the assumptions made in [2, Ass. 2.6 and Eq. (2.19)] require the differential operator $L$ to be at least $H^{2}(\mathcal{D})$-regular. Thus, our results do not only generalize the analysis of [2] for the strong error to different norms, but also to less regular differential operators. This is of relevance for several practical applications, since the spatial domain, where the GRF is simulated, may be non-convex or the coefficient $A$ may have jumps. For this reason, in Subsection 6.3.2 we work under the assumption that $L$ is $H^{1+\alpha}(\mathcal{D})$-regular for some $0<\alpha\leq 1$ (for instance, $\alpha<\nicefrac{{\pi}}{{\omega}}$ if $\mathcal{D}$ is a non-convex domain with largest interior angle $\omega>\pi$).

As an interim result while deriving the error bounds (1.3)–(1.6) for the sinc-Galerkin approximation, we prove a non-trivial extension of one of the main results in [5]. Namely, we show that for all

[TABLE]

and for all $\varepsilon>0$, there exists a constant $C>0$ such that, for $N\in\mathbb{N}$ and $g\in H^{\delta}(\mathcal{D})$,

[TABLE]

Here, $\widetilde{L}_{N}^{-1}\colon H^{-1}(\mathcal{D})\to V_{N}$ denotes the approximation of the data-to-solution map $L^{-1}\colon H^{-1}(\mathcal{D})\to H^{1}_{0}(\mathcal{D})$ with respect to the Galerkin space $V_{N}\subset H^{1}_{0}(\mathcal{D})$. For details see Theorem 6.6, Remark 6.7 and Lemmata 6.21–6.22. This error estimate was proven in [5, Thm. 4.3 & Rem. 4.1] only for $\beta\in(0,1)$, $\sigma=0$, and $\delta\geq 0$, see also the comparison in Remark 6.8.

1.2. Outline

After specifying the mathematical setting as well as our notation in Subsections 1.3–1.4, we rigorously define the second-order elliptic differential operator $L$ from (1.2) under minimal assumptions on the coefficients $A,\kappa$ and the domain $\mathcal{D}\subset\mathbb{R}^{d}$ in Section 2; thereby collecting several auxiliary results for this type of operators. Section 3 is devoted to the regularity analysis of a GRF colored by a linear operator $T$ which is bounded on $L_{2}(\mathcal{D})$. These results are subsequently applied in Section 4 to the class of generalized Whittle–Matérn fields, where $T:=L^{-\beta}$ with $L$ defined as in Section 2 and $\beta>\nicefrac{{d}}{{4}}$. In Section 5 we derive the convergence results (1.3)–(1.6) for spectral Galerkin approximations where the finite-dimensional subspace $V_{N}$ is generated by the eigenvectors of the operator $L$ corresponding to the $N$ smallest eigenvalues. We then investigate sinc-Galerkin approximations in Section 6, where we first let $V_{N}$ be an abstract Galerkin space satisfying certain approximation properties, see Subsections 6.1–6.2. Subsequently, in Subsection 6.3 we show that these properties are indeed satisfied if the Galerkin spaces originate from a quasi-uniform family of finite element discretizations of polynomial degree $p\in\mathbb{N}$, and we discuss the convergence behavior for two cases in detail:

(i) the coefficients $A,\kappa$ and the domain $\mathcal{D}$ in (1.2) are smooth, and

(ii) $A,\kappa,\mathcal{D}$are such that the differential operator $L$ in (1.2) is only $H^{1+\alpha}(\mathcal{D})$-regular for some $0<\alpha\leq 1$.

In Section 7 we perform several numerical experiments for the model example (1.1), $d=1$, and sinc-Galerkin discretizations generated with a finite element method of polynomial degree $p\in\{1,2\}$. In Section 8 we reflect on our outcomes.

1.3. Setting

Throughout this article, we let $(\Omega,\mathcal{F},\mathbb{P})$ be a complete probability space with expectation operator $\mathbb{E}$, and $\mathcal{D}$ be a bounded, connected and open subset of $\mathbb{R}^{d}$, $d\in\mathbb{N}$, with closure $\overline{\mathcal{D}}$.

In addition, we let $\mathcal{W}\colon L_{2}(\mathcal{D})\rightarrow L_{2}(\Omega)$ be an $L_{2}(\mathcal{D})$-isonormal Gaussian process in the sense of [29, Def. 1.1.1].

1.4. Notation

For $B\subseteq\mathbb{R}^{d}$, $\mathcal{B}(B)$ denotes the Borel $\sigma$-algebra on $B$ (i.e., the $\sigma$-algebra generated by the sets that are relatively open in $B$). For two $\sigma$-algebras $\mathcal{F}$ and $\mathcal{G}$, $\mathcal{F}\otimes\mathcal{G}$ is the $\sigma$-algebra generated by $\mathcal{F}\times\mathcal{G}$.

If $(E,\|\,\cdot\,\|_{E})$ is a Banach space, then $(E^{*},\|\,\cdot\,\|_{E^{*}})$ denotes its dual, $\langle\,\cdot\,,\,\cdot\,\rangle_{E^{*}\times E}$ the duality pairing on $E^{*}\times E$, $\operatorname{Id}_{E}$ the identity on $E$, and $\mathcal{L}(E;F)$ the space of bounded linear operators from $(E,\|\,\cdot\,\|_{E})$ to another Banach space $(F,\|\,\cdot\,\|_{F})$. For $T\in\mathcal{L}(E;F)$ we write $T^{*}\in\mathcal{L}(F^{*};E^{*})$ for the adjoint of $T$. If $V$ is a vector space such that $E,F\subseteq V$ and if, in addition, $\operatorname{Id}_{V}|_{E}\in\mathcal{L}(E;F)$, then we write $(E,\|\,\cdot\,\|_{E})\hookrightarrow(F,\|\,\cdot\,\|_{F})$. Moreover, the notation $(E,\|\,\cdot\,\|_{E})\cong(F,\|\,\cdot\,\|_{F})$ indicates that $(E,\|\,\cdot\,\|_{E})\hookrightarrow(F,\|\,\cdot\,\|_{F})\hookrightarrow(E,\|\,\cdot\,\|_{E})$.

If not specified otherwise, $(\,\cdot\,,\,\cdot\,)_{H}$ is the inner product on a Hilbert space $H$ and $\mathcal{L}_{2}(H;U)\subseteq\mathcal{L}(H;U)$ denotes the Hilbert space of Hilbert–Schmidt operators between two Hilbert spaces $H$ and $U$. The adjoint of $T\in\mathcal{L}(H;U)$ is identified with $T^{*}\in\mathcal{L}(U;H)$ (via the Riesz maps on $H$ and on $U$). We write $\mathcal{L}(E)$ and $\mathcal{L}_{2}(H)$ whenever $E=F$ and $H=U$. The domain of a possibly unbounded operator $L$ is denoted by $\mathscr{D}(L)$.

For $1\leq q<\infty$, $L_{q}(\mathcal{D};E)$ is the space of (equivalence classes of) $E$-valued, Bochner measurable, $q$-integrable functions on $\mathcal{D}$ and $L_{q}(\Omega;E)$ denotes the space of (equivalence classes of) $E$-valued random variables with finite $q$-th moment, i.e.,

[TABLE]

The space $L_{\infty}(\mathcal{D};E)$ consists of all equivalence classes of $E$-valued, Bochner measurable functions which are essentially bounded on $\mathcal{D}$, i.e.,

[TABLE]

For $\gamma\in(0,1)$, we furthermore define the mappings

[TABLE]

on the Banach space

[TABLE]

of continuous functions from $\overline{\mathcal{D}}$ to $(E,\|\,\cdot\,\|_{E})$ via

[TABLE]

Note that the norm $\|\,\cdot\,\|_{C^{\gamma}(\overline{\mathcal{D}};E)}$ renders the subspace

[TABLE]

of $\gamma$-Hölder continuous functions a Banach space. Whenever the functions or random variables are real-valued, we omit the image space and write $C(\overline{\mathcal{D}})$, $C^{\gamma}(\overline{\mathcal{D}})$, $L_{q}(\mathcal{D})$, and $L_{q}(\Omega)$, respectively. For $\sigma>0$, the (integer- or fractional-order) Sobolev space is denoted by $H^{\sigma}(\mathcal{D})$ (see [12, Sec. 2], see also [41, Sec. 1.11.4/5]), and $H^{1}_{0}(\mathcal{D})\subset H^{1}(\mathcal{D})$ is the closure of the space $C^{\infty}_{\operatorname{c}}(\mathcal{D})$ of compactly supported smooth functions in \bigl{(}H^{1}(\mathcal{D}),\|\,\cdot\,\|_{H^{1}(\mathcal{D})}\bigr{)}.

We mark equations which hold almost everywhere or $\mathbb{P}$-almost surely with a.e. and $\mathbb{P}$-a.s., respectively. For two random variables $X,Y$, we write $X\overset{d}{=}Y$ whenever $X$ and $Y$ have the same probability distribution. The Dirac measure at $x\in\overline{\mathcal{D}}$ is denoted by $\text{\dh}_{x}$. Given a parameter set $\mathcal{P}$ and mappings $A,B\colon\mathcal{P}\to\mathbb{R}$, we let $A(p)\lesssim B(p)$ denote the relation that there exists a constant $C>0$, independent of $p\in\mathcal{P}$, such that $A(p)\leq CB(p)$ for all $p\in\mathcal{P}$. For a further parameter set $\mathcal{Q}$ and mappings $A,B\colon\mathcal{P}\times\mathcal{Q}\rightarrow\mathbb{R}$, we write $A(p,q)\lesssim_{q}B(p,q)$ if, for all $q\in\mathcal{Q}$, there exists a constant $C_{q}>0$, independent of $p\in\mathcal{P}$, such that $A(p,q)\leq C_{q}B(p,q)$ for all $p\in\mathcal{P}$ and $q\in\mathcal{Q}$. Finally, $A(p)\eqsim B(p)$ indicates that both relations, $A(p)\lesssim B(p)$ and $B(p)\lesssim A(p)$, hold simultaneously; and similarly for $A(p,q)\eqsim_{q}B(p,q)$.

2. Auxiliary results on second-order elliptic differential operators

As outlined in Subsection 1.1, the overall objective of this article is to study (generalized) Whittle–Matérn fields and Galerkin approximations for them. Here, we call a Gaussian random field a generalized Whittle–Matérn field if its covariance operator is given by a negative fractional power of a second-order elliptic differential operator. The purpose of this section is to present preliminary results on second-order differential operators which will be of importance for the regularity and error analysis of these fields.

Firstly, we specify the class of differential operators that we consider. We start by formulating assumptions on the coefficients of the operator.

Assumption 2.1** (on the coefficients $A$ and $\kappa$).**

Throughout this article we assume:

1. I.

A\in L_{\infty}\bigl{(}\mathcal{D};\mathbb{R}^{d\times d}\bigr{)} is symmetric and uniformly positive definite, i.e.,

[TABLE] 2. II.

$\kappa\in L_{\infty}(\mathcal{D})$.

Where explicitly specified, we require in addition:

1. III.

$A\colon\overline{\mathcal{D}}\to\mathbb{R}^{d\times d}$ is Lipschitz continuous on the closure $\overline{\mathcal{D}}$, i.e.,

[TABLE]

for all $i,j\in\{1,\ldots,d\}$.

Under Assumptions 2.1.I–II we let $L\colon\mathscr{D}(L)\subset L_{2}(\mathcal{D})\rightarrow L_{2}(\mathcal{D})$ denote the maximal accretive operator on $L_{2}(\mathcal{D})$ associated with $A$ and $\kappa^{2}$ with largest domain $\mathscr{D}(L)\subset H^{1}_{0}(\mathcal{D})$. By this we mean that $\mathscr{D}(L)$ consists of precisely those $u\in H^{1}_{0}(\mathcal{D})$ for which there exists a constant $C\geq 0$ such that

[TABLE]

and, for $u\in\mathscr{D}(L)$, $Lu$ is the unique element of $L_{2}(\mathcal{D})$ which, for all $v\in H^{1}_{0}(\mathcal{D})$, satisfies

[TABLE]

It is well-known that the operator $L\colon\mathscr{D}(L)\to L_{2}(\mathcal{D})$ defined via (2.2) is densely defined and self-adjoint (e.g., [31, Prop. 1.22 and Prop. 1.24]). Furthermore, by the Lax–Milgram lemma, its inverse exists and extends to a bounded linear operator $L^{-1}\colon H^{1}_{0}(\mathcal{D})^{*}\to H^{1}_{0}(\mathcal{D})$ (e.g., [31, Lem. 1.3]). By the Kondrachov compactness theorem $L^{-1}\colon L_{2}(\mathcal{D})\to L_{2}(\mathcal{D})$ is compact (e.g., [18, Thm. 7.22]).

For this reason, the spectrum of $L$ consists of a system of only positive eigenvalues $(\lambda_{j})_{j\in\mathbb{N}}$ with no accumulation point, whence we can assume them to be in nondecreasing order. The following asymptotic spectral behavior, known as Weyl’s law (see, e.g., [11, Thm. 6.3.1]), will be exploited several times in our analysis.

Lemma 2.2**.**

Let $L$ be the second-order differential operator in (2.2), defined with respect to the bounded open domain $\mathcal{D}\subset\mathbb{R}^{d}$, and with coefficients $A$ and $\kappa$ fulfilling Assumptions 2.1.I–II. Then, the eigenvalues of $L$ (in nondecreasing order) satisfy

[TABLE]

We let $\mathcal{E}:=\{e_{j}\}_{j\in\mathbb{N}}$ denote a system of eigenvectors of the operator $L$ in (2.2) which corresponds to the eigenvalues $(\lambda_{j})_{j\in\mathbb{N}}$ and which is orthonormal in $L_{2}(\mathcal{D})$. Note that, for $\sigma>0$, the fractional power operator $L^{\sigma}\colon\mathscr{D}(L^{\sigma})\subset L_{2}(\mathcal{D})\rightarrow L_{2}(\mathcal{D})$ is well-defined. Indeed, on the domain

[TABLE]

the action of $L^{\sigma}$ is given via the spectral representation

[TABLE]

The subspace

[TABLE]

is itself a Hilbert space with respect to the inner product

[TABLE]

and the corresponding induced norm $\|\,\cdot\,\|_{\sigma}$. In what follows, we let $\dot{H}^{0}_{L}:=L_{2}(\mathcal{D})$ and, for $\sigma>0$, $\dot{H}^{-\sigma}_{L}$ denotes the dual space $(\dot{H}^{\sigma}_{L})^{*}$ after identification via the inner product $(\,\cdot\,,\,\cdot\,)_{L_{2}(\mathcal{D})}$ on $L_{2}(\mathcal{D})$ which is continuously extended to a duality pairing.

In order to derive regularity and convergence results with respect to the Sobolev space $H^{\sigma}(\mathcal{D})$ and the space $C^{\gamma}(\overline{\mathcal{D}})$ of $\gamma$-Hölder continuous functions in (1.9), we wish to relate the various norms involved by well-known results from interpolation theory and by Sobolev embeddings. To this end, we need to consider various assumptions on the spatial domain $\mathcal{D}$, specified below.

Assumption 2.3** (on the domain $\mathcal{D}$).**

Throughout this article, we assume that

1. I.

$\mathcal{D}$ has a Lipschitz continuous boundary $\partial\mathcal{D}$.

Where explicitly specified, we additionally suppose one or both of the following:

1. II.

$\mathcal{D}$ is convex; 2. III.

$\mathcal{D}$ is a polytope.

Note that II. implies I. (see, e.g., [21, Cor. 1.2.2.3]).

In the following lemma we specify the relationship between the spaces $\dot{H}^{\sigma}_{L}$ in (2.4) and the Sobolev space $H^{\sigma}(\mathcal{D})$, under two sets of assumptions on the spatial domain $\mathcal{D}$ and on the coefficients $A$, $\kappa$ of the differential operator $L$ in (2.2). We recall that $[E,F]_{\sigma}$ denotes the complex interpolation space between $(E,\|\,\cdot\,\|_{E})$ and $(F,\|\,\cdot\,\|_{F})$ with parameter $\sigma\in[0,1]$, see, e.g., [27, Ch. 2].

Lemma 2.4**.**

Let Assumptions 2.1.I–II and 2.3.I be satisfied. Then

[TABLE]

holds for the space \bigl{(}\dot{H}^{\sigma}_{L},\|\,\cdot\,\|_{\sigma}\bigr{)} from (2.4). Furthermore,

[TABLE]

and the norms $\|\,\cdot\,\|_{\sigma}$, $\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}$ are equivalent on $\dot{H}^{\sigma}_{L}$ for $0\leq\sigma\leq 1$ and $\sigma\neq\nicefrac{{1}}{{2}}$. If, in addition, Assumptions 2.1.III and 2.3.II hold, then

[TABLE]

Proof.

First, note that [41, Cor. 2.4] implies (2.5). If $(E,\|\,\cdot\,\|_{E})$, $(F,\|\,\cdot\,\|_{F})$, and $(G,\|\,\cdot\,\|_{G})$ are Banach spaces such that $(F,\|\,\cdot\,\|_{F})\hookrightarrow(G,\|\,\cdot\,\|_{G})$, then by definition of complex interpolation we have \bigl{(}[E,F]_{\sigma},\|\,\cdot\,\|_{[E,F]_{\sigma}}\bigr{)}\hookrightarrow\bigl{(}[E,G]_{\sigma},\|\,\cdot\,\|_{[E,G]_{\sigma}}\bigr{)}. This observation in connection with [41, Thm. 1.35] (which collects several results from [39]) shows (2.6). Equivalence of $\|\,\cdot\,\|_{\sigma}$, $\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}$ on $\dot{H}^{\sigma}_{L}$ for $0\leq\sigma\leq 1$, $\sigma\neq\nicefrac{{1}}{{2}}$, is proven in [20, Thm. 8.1]. By combining (2.5) for $\sigma=1$, [27, Thm. 4.36] and [22, Lem. A2] (recalling Assumption 2.3.II) we find that (2.7) for $\sigma\in(1,2)$ follows once (2.7) is established for the case $\sigma=2$.

It thus remains to prove (2.7) for $\sigma=2$. To this end, we first observe that, for a vanishing coefficient $\kappa\equiv 0$ of the operator $L$ in (2.2), we have, e.g., by [21, Thm. 3.2.1.2] the regularity result

[TABLE]

If $\kappa\not\equiv 0$, then $u\in H^{1}_{0}(\mathcal{D})$ satisfies the equality $-\nabla\cdot(A\nabla u)=f-\kappa^{2}u$ in the weak sense so that [21, Thm. 3.2.1.2] applied to $\widetilde{f}:=f-\kappa^{2}u\in L_{2}(\mathcal{D})$ again yields (2.8). By the closed graph theorem, \bigl{(}\dot{H}^{2}_{L},\|\,\cdot\,\|_{2}\bigr{)}\hookrightarrow\bigl{(}H^{2}(\mathcal{D})\cap H^{1}_{0}(\mathcal{D}),\|\,\cdot\,\|_{H^{2}(\mathcal{D})}\bigr{)}.

We now establish the reverse embedding. By Assumption 2.1.III and, e.g., [16, Thm. 4 in Ch. 5.8] (note that the assumptions on the boundary posed therein can be circumvented by exploiting an extension argument as, e.g., in [36, Sec. VI.2.3 Thm. 3], see also the remark below [16, Thm. 4 in Ch. 5.8]), $A_{ij}$ is differentiable a.e. in $\mathcal{D}$ with essentially bounded weak derivatives $\partial_{x_{k}}A_{ij}\in L_{\infty}(\mathcal{D})$, $1\leq i,j,k\leq d$. Thus (by first approximating $A_{ij}$ in $H^{1}(\mathcal{D})$ with a sequence in $C^{\infty}(\mathcal{D})$ to obtain that $A_{ij}\partial_{x_{j}}u$ is weakly differentiable with $\partial_{x_{k}}(A_{ij}\partial_{x_{j}}u)=\partial_{x_{k}}A_{ij}\partial_{x_{j}}u+A\partial_{x_{k}x_{j}}u$), we conclude that $A\nabla u\in H^{1}(\mathcal{D})^{d}$ whenever $u\in H^{2}(\mathcal{D})\cap H^{1}_{0}(\mathcal{D})$. Finally, by the closed graph theorem \bigl{(}H^{2}(\mathcal{D})\cap H^{1}_{0}(\mathcal{D}),\|\,\cdot\,\|_{H^{2}(\mathcal{D})}\bigr{)}\hookrightarrow\bigl{(}\dot{H}^{2}_{L},\|\,\cdot\,\|_{2}\bigr{)} follows. ∎

3. General results on Gaussian Random Fields (GRFs)

In this section we address different notions of regularity (Hölder and Sobolev) for Gaussian random fields (GRFs) and their covariance functions. We first recall the definition of a GRF and specify then what we mean by a colored GRF. As usually, we work in the setting formulated in Subsection 1.3.

Definition 3.1**.**

Let $B\subseteq\mathbb{R}^{d}$. A family of $\mathcal{F}$-measurable $\mathbb{R}$-valued random variables $(\mathcal{Z}(x))_{x\in B}$ is called a random field (indexed by $B$). It is called Gaussian if the random vector $(\mathcal{Z}(x_{1}),\ldots,\mathcal{Z}(x_{n}))^{\top}$ is Gaussian for all finite sets $\{x_{1},\ldots,x_{n}\}\subset B$. It is called continuous if the mapping $x\mapsto\mathcal{Z}(x)(\omega)$ is continuous for all $\omega\in\Omega$.

Definition 3.2**.**

Let $T\in\mathcal{L}(L_{2}(\mathcal{D}))$. We call $\mathcal{Z}\colon\mathcal{D}\times\Omega\rightarrow\mathbb{R}$ a Gaussian random field (GRF) colored by $T$ if it is a GRF, a $\mathcal{B}(\mathcal{D})\otimes\mathcal{F}$-measurable mapping, and

[TABLE]

Remark 3.3*.*

It is well-known and easily verified (see also Proposition 3.7) that there exists a square-integrable GRF $\mathcal{Z}$ colored by $T$ if and only if $T\in\mathcal{L}_{2}(L_{2}(\mathcal{D}))$, and in this case \mathbb{E}\bigl{[}\|\mathcal{Z}\|_{L_{2}(\mathcal{D})}^{2}\bigr{]}=\operatorname{tr}(TT^{*}), where $\operatorname{tr}(\,\cdot\,)$ is the trace on $L_{2}(\mathcal{D})$.

3.1. Hölder regularity of GRFs

We now provide an abstract result on the construction and Hölder regularity of a GRF assuming that the color and, thus, the covariance structure of the field is given.

Proposition 3.4**.**

Assume that $T\in\mathcal{L}(L_{2}(\mathcal{D});C^{\gamma}(\overline{\mathcal{D}}))$ for some $\gamma\in(0,1)$. Then there exists a continuous GRF $\mathcal{Z}$ colored by $T$ such that

[TABLE]

Furthermore, for $q\in(0,\infty)$ and $\delta\in(0,\gamma)$, we have

[TABLE]

Proof.

We first define the random field $\mathcal{Z}_{0}\colon\overline{\mathcal{D}}\times\Omega\to\mathbb{R}$ by $\mathcal{Z}_{0}(x):=\mathcal{W}(T^{*}\text{\dh}_{x})$ for all $x\in\overline{\mathcal{D}}$. By the properties of an isonormal Gaussian process we find, for $x,y\in\overline{\mathcal{D}}$,

[TABLE]

Since $\mathcal{Z}_{0}(x)-\mathcal{Z}_{0}(y)=\mathcal{W}(T^{*}(\text{\dh}_{x}-\text{\dh}_{y}))$ is a real-valued Gaussian random variable, we can apply the Khintchine inequalities (see, e.g., [25, Thm. 4.7 and p. 103]) and conclude with (3.4) that, for all $q\in(0,\infty)$, the estimate

[TABLE]

holds, with a constant $C_{q}>0$ depending only on $q$.

Thus, by the Kolmogorov–Chentsov continuity theorem (e.g., [34, Thm. I.2.1], combined with an extension argument as discussed in the proof of [28, Thm. 2.1], see also [36, Thm. VI.2.3]), there exists a continuous random field $\mathcal{Z}\colon\overline{\mathcal{D}}\times\Omega\to\mathbb{R}$ such that $\mathcal{Z}(x)=\mathcal{Z}_{0}(x)$ $\mathbb{P}$-a.s. for all $x\in\overline{\mathcal{D}}$, and furthermore, for every $\delta\in(0,\gamma)$ and every finite $q>(\gamma-\delta)^{-1}$, we can find a constant $C_{q,\gamma,\delta,\mathcal{D}}>0$, depending only on $q$, $\gamma$, $\delta$, as well as the dimension and the diameter of $\mathcal{D}\subset\mathbb{R}^{d}$, such that

[TABLE]

Next, again by the Khintchine inequalities, we have, for all $x\in\overline{\mathcal{D}}$ and all $q\in(0,\infty)$,

[TABLE]

From (1.7)–(1.8) we deduce, for every $\delta\in(0,1)$ and all $f\in C^{\delta}(\overline{\mathcal{D}})$, the relation

[TABLE]

We combine this observation with (3.5), (3.6), and (3.7) to derive, for all $\delta\in(0,\gamma)$ and all finite $q>(\gamma-\delta)^{-1}$, the bound

[TABLE]

Note that Hölder’s inequality and (3.8) ensure that (3.3) holds for every $\delta\in(0,\gamma)$ and every $q\in(0,\infty)$. Furthermore, for every $\psi\in L_{2}(\mathcal{D})$, one readily verifies the identity \mathbb{E}\bigl{[}|(\mathcal{Z},\psi)_{L_{2}(\mathcal{D})}-\mathcal{W}(T^{*}\psi)|^{2}\bigr{]}=0, i.e., $\mathcal{Z}$ is colored by $T$. ∎

If Assumption 2.3.I is fulfilled, the Sobolev embedding theorem (see, e.g., [12, Thm. 5.4 and Thm. 8.2]) is applicable and we obtain $\gamma$-Hölder continuity (1.9) for elements in the fractional-order Sobolev space $H^{\gamma+\nicefrac{{d}}{{2}}}(\mathcal{D})$ for every $\gamma\in(0,1)$. This continuous embedding, $H^{\gamma+\nicefrac{{d}}{{2}}}(\mathcal{D})\hookrightarrow C^{\gamma}(\overline{\mathcal{D}})$, combined with Proposition 3.4 leads to the following result.

Corollary 3.5**.**

Let Assumption 2.3.I, $\gamma\in(0,1)$, and T\in\mathcal{L}\bigl{(}L_{2}(\mathcal{D});H^{\gamma+\nicefrac{{d}}{{2}}}(\mathcal{D})\bigr{)} be satisfied. Then there exists a continuous GRF $\mathcal{Z}\colon\overline{\mathcal{D}}\times\Omega\rightarrow\mathbb{R}$ colored by $T$, cf. (3.1), such that $\mathcal{Z}(x)=\mathcal{W}(T^{*}\text{\dh}_{x})$ $\mathbb{P}$-a.s. for all $x\in\overline{\mathcal{D}}$. Moreover, the stability estimate

[TABLE]

for the $q$-th moment of $\mathcal{Z}$ with respect to the $\delta$-Hölder norm (1.8) holds for every $\delta\in(0,\gamma)$ and $q\in(0,\infty)$.

We close the subsection with a brief discussion on

(i) the continuity of covariance functions of colored GRFs, and

(ii) the $L_{\infty}(\mathcal{D}\times\mathcal{D})$-distance between two covariance functions of GRFs colored by different operators.

By definition, the covariance function $\varrho\in L_{2}(\mathcal{D}\times\mathcal{D})$ of a square-integrable random field $\mathcal{Z}\in L_{2}(\mathcal{D}\times\Omega)$ satisfies

[TABLE]

We obtain the one-to-one correspondence

[TABLE]

with the covariance operator $\mathcal{C}\colon L_{2}(\mathcal{D})\to L_{2}(\mathcal{D})$ of the field $\mathcal{Z}$, which is defined via

[TABLE]

From this definition it is evident that a GRF $\mathcal{Z}$ colored by $T$ (note that $\mathbb{E}[\mathcal{Z}]=0$ by construction, see Definition 3.2) has the covariance operator $\mathcal{C}=TT^{*}$. In the next lemma, this relation is exploited to characterize continuity of the covariance function $\varrho$ in terms of the color $T$ of the GRF $\mathcal{Z}$.

Proposition 3.6**.**

Let $\mathcal{Z},\widetilde{\mathcal{Z}}$ be GRFs colored by $T$ and $\widetilde{T}$, respectively, see (3.1), with covariance functions denoted by $\varrho$ and $\widetilde{\varrho}$, cf. (3.10). Then,

1. (i)

$\varrho$* has a continuous representative on $\overline{\mathcal{D}\times\mathcal{D}}$ (again denoted by $\varrho$) if and only if $T\in\mathcal{L}(L_{2}(\mathcal{D});C(\overline{\mathcal{D}}))$. In this case,*

[TABLE] 2. (ii)

if $T,\widetilde{T}\in\mathcal{L}(L_{2}(\mathcal{D});C(\overline{\mathcal{D}}))$, then $\varrho,\widetilde{\varrho}\in C(\overline{\mathcal{D}\times\mathcal{D}})$ satisfy

[TABLE]

Proof.

By (3.10), the covariance function $\varrho$ of a GRF $\mathcal{Z}$ colored by $T$ is given by

[TABLE]

First, let $T\in\mathcal{L}(L_{2}(\mathcal{D});C(\overline{\mathcal{D}}))$. Then, we have $T^{*}\in\mathcal{L}(C(\overline{\mathcal{D}})^{*};L_{2}(\mathcal{D})^{*})$ and continuity of $\varrho\colon\overline{\mathcal{D}\times\mathcal{D}}\to\mathbb{R}$ follows from (3.13). Assume now that $\varrho\in C(\overline{\mathcal{D}\times\mathcal{D}})$. Then, again by (3.13), we obtain $\|T^{*}\text{\dh}_{x}\|_{L_{2}(\mathcal{D})^{*}}^{2}=\varrho(x,x)<\infty$ for all $x\in\overline{\mathcal{D}}$ and

[TABLE]

holds for all $\phi\in L_{2}(\mathcal{D})$ with $\|\phi\|_{L_{2}(\mathcal{D})}\leq 1$. Thus, $T\in\mathcal{L}(L_{2}(\mathcal{D});C(\overline{\mathcal{D}}))$ if $\varrho$ is continuous. Furthermore, by identifying $L_{2}(\mathcal{D})^{*}\cong L_{2}(\mathcal{D})$ via the Riesz map, the covariance operator $\mathcal{C}$ of $\mathcal{Z}$ satisfies $\mathcal{C}=TT^{*}\in\mathcal{L}(C(\overline{\mathcal{D}})^{*};C(\overline{\mathcal{D}}))$, and we can deduce (3.11) from (3.13) since, for all $x,y\in\overline{\mathcal{D}}$,

[TABLE]

Finally, the estimate (3.12) can be shown similarly since, for all $x,y\in\overline{\mathcal{D}}$,

[TABLE]

3.2. Sobolev regularity of GRFs and their covariances

After having characterized

1. (i)

the Hölder regularity (in $L_{q}(\Omega)$-sense) of a GRF $\mathcal{Z}$, and 2. (ii)

continuity of the covariance function $\varrho$ in (3.10),

in terms of the color of $\mathcal{Z}$, we now proceed with this discussion for Sobolev spaces. Specifically, we investigate the regularity of $\mathcal{Z}$ in $L_{q}(\Omega;H^{\sigma}(\mathcal{D}))$ and of the covariance function $\varrho$ with respect to the norm on the mixed Sobolev space

[TABLE]

Here, $\mathop{\hat{\otimes}}$ denotes the tensor product of Hilbert spaces. Thus, the inner product on $H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})$ inducing the norm $\|\,\cdot\,\|_{H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})}$ is uniquely defined via

[TABLE]

To this end, in the following proposition we first quantify the $\dot{H}^{\sigma}_{L}$-regularity (in $L_{q}(\Omega)$-sense) of a colored GRF in terms of its color, cf. (2.4) and Definition 3.2. In addition, we specify the regularity of the covariance function (3.10) in the Hilbert tensor product space

[TABLE]

cf. (3.14). Finally, we characterize the distance between two GRFs which are colored by different operators with respect to these norms.

Proposition 3.7**.**

Let $\mathcal{Z}\colon\mathcal{D}\times\Omega\to\mathbb{R}$ be a GRF colored by $T\in\mathcal{L}(L_{2}(\mathcal{D}))$, cf. (3.1). Then $\mathcal{Z}$ is square-integrable, i.e., $\mathcal{Z}\in L_{2}(\mathcal{D}\times\Omega)$, if and only if its covariance operator $\mathcal{C}=TT^{*}$ has a finite trace on $L_{2}(\mathcal{D})$. More generally, for $\sigma\geq 0$ and $q\in(0,\infty)$, we have

[TABLE]

Here, $\operatorname{tr}(\,\cdot\,)$ is the trace on $L_{2}(\mathcal{D})$, $L$ is the differential operator in (2.2) with coefficients $A,\kappa$ satisfying Assumptions 2.1.I–II, $\varrho$ is the covariance function of $\mathcal{Z}$, see (3.10), and $\mathcal{L}_{2}^{\theta;\sigma}$ is a short notation for the Hilbert–Schmidt space \mathcal{L}_{2}\bigl{(}\dot{H}^{\theta}_{L};\dot{H}^{\sigma}_{L}\bigr{)}.

If $\widetilde{\mathcal{Z}}\in L_{2}(\mathcal{D}\times\Omega)$ is another GRF colored by $\widetilde{T}\in\mathcal{L}(L_{2}(\mathcal{D}))$, with covariance function $\widetilde{\varrho}$ and covariance operator $\widetilde{\mathcal{C}}=\widetilde{T}\widetilde{T}^{*}$, we have, for $\sigma\geq 0$ and $q\in(0,\infty)$,

[TABLE]

Proof.

Assume first that $\mathcal{Z}\in L_{2}(\mathcal{D}\times\Omega)$. Since $\mathcal{Z}$ has mean zero and since it is colored by $T\in\mathcal{L}(L_{2}(\mathcal{D}))$, we obtain $\mathcal{C}=TT^{*}$, i.e.,

[TABLE]

By choosing $\phi=\psi:=\lambda_{j}^{\nicefrac{{\sigma}}{{2}}}e_{j}$, summing these equalities over $j\in\mathbb{N}$, and exchanging the order of summation and expectation, we obtain the identity

[TABLE]

and the first part of the proposition as well as (3.16) are proven. The estimate (3.17) follows from (3.16) by the Kahane–Khintchine inequalities (see, e.g., [25, Thm. 4.7 and p. 103]), since $\mathcal{Z}$ is an $\dot{H}^{\sigma}_{L}$-valued zero-mean Gaussian random variable.

Assume now that $\widetilde{\mathcal{Z}}\in L_{2}(\mathcal{D}\times\Omega)$ is another GRF colored by $\widetilde{T}\in\mathcal{L}(L_{2}(\mathcal{D}))$. Then we obtain (3.19) from (3.17), since $\mathcal{Z}-\widetilde{\mathcal{Z}}$ is again a GRF, colored by $T-\widetilde{T}$, see (3.1) and Definition 3.2. Furthermore, we find

[TABLE]

This proves (3.20) and (3.18) follows from this result for $\widetilde{\mathcal{Z}}\equiv 0$. ∎

Remark 3.8*.*

Note that if Assumptions 2.1.I–II, 2.3.I, and $0\leq\sigma\leq 1$ (or Assumptions 2.1.I–III, 2.3.II, and $0\leq\sigma\leq 2$) are satisfied and $\sigma\neq\nicefrac{{1}}{{2}}$, it follows from Lemma 2.4 that all assertions of Proposition 3.7 remain true if we replace the equalities with equivalences and the norms $\|\,\cdot\,\|_{\sigma}$, $\|\,\cdot\,\|_{\sigma,\sigma}$ (cf. the spaces in (2.4), (3.15)) with the Sobolev norm $\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}$ and with the norm $\|\,\cdot\,\|_{H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})}$ on the mixed Sobolev space (3.14), respectively. Furthermore, by (2.6) Proposition 3.7 provides upper bounds for these quantities if $\sigma=\nicefrac{{1}}{{2}}$.

4. Regularity of Whittle–Matérn fields

In this section we focus on the regularity of (generalized) Whittle–Matérn fields, i.e., of GRFs colored (cf. Definition 3.2) by a negative fractional power of the differential operator $L$ as provided in (2.2). Specifically, we consider

[TABLE]

for

[TABLE]

We emphasize the dependence of the covariance structure of $\mathcal{Z}^{\beta}$ on the fractional exponent $\beta>0$ by the index and write $\varrho^{\beta}$ for the covariance function (3.10) of $\mathcal{Z}^{\beta}$.

The first aim of this section is to apply Proposition 3.7 for specifying the regularity of $\mathcal{Z}^{\beta}$ in (4.1) and of its covariance function $\varrho^{\beta}$ with respect to the spaces $\dot{H}^{\sigma}_{L}$ and $\dot{H}^{\sigma,\sigma}_{L}$ in (2.4), (3.15). As already pointed out in Remark 3.8, provided that the assumptions of Lemma 2.4 are satisfied, this implies regularity in the Sobolev space $H^{\sigma}(\mathcal{D})$ and in the mixed Sobolev space $H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})$ in (3.14), respectively.

Besides this regularity result with respect to the spaces $\dot{H}^{\sigma}_{L}$ and $H^{\sigma}(\mathcal{D})$, we obtain a stability estimate with respect to the Hölder norm from Corollary 3.5 and continuity of the covariance function from Proposition 3.6. Although we believe that, at least in some specific cases, these results are well-known, for the sake of completeness, we derive them here in our general framework.

Lemma 4.1**.**

Let Assumptions 2.1.I–II be fulfilled, $\beta,q\in(0,\infty)$, $\sigma\geq 0$, and $\mathcal{Z}^{\beta}$ be the Whittle–Matérn field in (4.1), with covariance function $\varrho^{\beta}$. Then,

1. (i)

\mathbb{E}\Bigl{[}\bigl{\|}\mathcal{Z}^{\beta}\bigr{\|}_{\sigma}^{q}\Bigr{]}<\infty* if and only if $2\beta>\sigma+\nicefrac{{d}}{{2}}$, and* 2. (ii)

\bigl{\|}\varrho^{\beta}\bigr{\|}_{\sigma,\sigma}<\infty* if and only if $2\beta>\sigma+\nicefrac{{d}}{{4}}$.*

If, in addition, Assumption 2.3.I and $0\leq\sigma\leq 1$ (or Assumptions 2.1.I–III, 2.3.II, and $0\leq\sigma\leq 2$) hold, then the assertions (i)–(ii) remain true if we formulate them with respect to the Sobolev norms $\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}$, $\|\,\cdot\,\|_{H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})}$.

Proof.

By Proposition 3.7 we have, for any $\beta,q\in(0,\infty)$ and $\sigma\geq 0$,

[TABLE]

Combining the spectral behavior (2.3) of $L$ from Lemma 2.2 with (4.3)/(4.4) proves (i)/(ii) for $\|\,\cdot\,\|_{\sigma}$, $\|\,\cdot\,\|_{\sigma,\sigma}$. If the assumptions stated in the second part of the lemma are satisfied, then applying Lemma 2.4 completes the proof. ∎

Lemma 4.2**.**

Suppose that

1. (i)

Assumptions 2.1.I–II are satisfied, $0<2\gamma\leq 1$, and $d=1$, or 2. (ii)

Assumptions 2.1.I–III and 2.3.II are fulfilled, $d\in\{1,2,3\}$ and $\gamma\in(0,1)$ are such that $\gamma\leq 2-\nicefrac{{d}}{{2}}$.

In either of these cases and if $2\beta\geq\gamma+\nicefrac{{d}}{{2}}$, there exists a continuous Whittle–Matérn field $\mathcal{Z}^{\beta}\colon\overline{\mathcal{D}}\times\Omega\rightarrow\mathbb{R}$ satisfying (4.1) such that $\mathcal{Z}^{\beta}(x)=\mathcal{W}(L^{-\beta}\text{\dh}_{x})$ $\mathbb{P}$-a.s. for all $x\in\overline{\mathcal{D}}$, and, for every $\delta\in(0,\gamma)$ and $q\in(0,\infty)$, the bound

[TABLE]

for the $q$-th moment of $\mathcal{Z}^{\beta}$ with respect to the $\delta$-Hölder norm, cf. (1.8), holds.

Proof.

Note that by definition of $\dot{H}^{\sigma}_{L}$, see (2.4), for any $\beta>0$, the operator

[TABLE]

is an isometric isomorphism. For this reason, $L^{-\beta}\colon L_{2}(\mathcal{D})\to\dot{H}^{\gamma+\nicefrac{{d}}{{2}}}_{L}$ is bounded provided that $2\beta\geq\gamma+\nicefrac{{d}}{{2}}$. For $d$ and $\gamma$ as specified in (i)/(ii) above, we have \bigl{(}\dot{H}^{\gamma+\nicefrac{{d}}{{2}}}_{L},\|\,\cdot\,\|_{\gamma+\nicefrac{{d}}{{2}}}\bigr{)}\hookrightarrow\bigl{(}H^{\gamma+\nicefrac{{d}}{{2}}}(\mathcal{D}),\|\,\cdot\,\|_{H^{\gamma+\nicefrac{{d}}{{2}}}(\mathcal{D})}\bigr{)} by the relations (2.6)–(2.7) from Lemma 2.4 and we conclude that L^{-\beta}\in\mathcal{L}\bigl{(}L_{2}(\mathcal{D});H^{\gamma+\nicefrac{{d}}{{2}}}(\mathcal{D})\bigr{)}. The proof is then completed by applying Corollary 3.5 in both cases (i)/(ii). ∎

Lemma 4.3**.**

Let Assumptions 2.1.I–II be satisfied and $\beta>\nicefrac{{d}}{{4}}$. Suppose furthermore that a system of $L_{2}(\mathcal{D})$-orthonormal eigenvectors $\mathcal{E}=\{e_{j}\}_{j\in\mathbb{N}}$ corresponding to the eigenvalues $0<\lambda_{1}\leq\lambda_{2}\leq\ldots$ of $L$ in (2.2) is uniformly bounded in $C(\overline{\mathcal{D}})$, i.e.,

[TABLE]

Then the covariance function, cf. (3.10), of the Whittle–Matérn field $\mathcal{Z}^{\beta}$ in (4.1) has a continuous representative $\varrho^{\beta}\colon\overline{\mathcal{D}\times\mathcal{D}}\to\mathbb{R}$ and

[TABLE]

where $\operatorname{tr}(\,\cdot\,)$ denotes the trace on $L_{2}(\mathcal{D})$.

Proof.

By Proposition 3.6(i) we have to show boundedness of $L^{-\beta}\colon L_{2}(\mathcal{D})\to C(\overline{\mathcal{D}})$ to infer that $\varrho^{\beta}\in C(\overline{\mathcal{D}\times\mathcal{D}})$, with

[TABLE]

For $\psi\in L_{2}(\mathcal{D})$, the spectral representation $L^{-\beta}\psi=\sum_{j\in\mathbb{N}}\lambda_{j}^{-\beta}(\psi,e_{j})_{L_{2}(\mathcal{D})}\,e_{j}$ shows that, for all $x\in\overline{\mathcal{D}}$,

[TABLE]

is finite, provided that $\beta>\nicefrac{{d}}{{4}}$. Here, we have used the Cauchy–Schwarz inequality and the spectral behavior (2.3) from Lemma 2.2 in the last estimate. Similarly,

[TABLE]

for all $\varphi\in C(\overline{\mathcal{D}})^{*}$. Combining (4.7) and (4.8) completes the proof. ∎

Remark 4.4*.*

Note that if $\gamma\in(0,1)$ and $d\in\{1,2,3\}$ are such that Assumption (i) or (ii) of Lemma 4.2 is satisfied, then the Sobolev embedding $H^{\delta+\nicefrac{{d}}{{2}}}(\mathcal{D})\hookrightarrow C^{\delta}(\overline{\mathcal{D}})$ and Lemma 2.4 are applicable for any $0<\delta\leq\gamma$. Thus, if $2\beta\geq\delta+\nicefrac{{d}}{{2}}$, we find

[TABLE]

i.e., $L^{-\beta}\colon L_{2}(\mathcal{D})\to C(\overline{\mathcal{D}})$ is bounded. Thus, by Proposition 3.6(i) the covariance function $\varrho^{\beta}\colon\overline{\mathcal{D}\times\mathcal{D}}\to\mathbb{R}$ of the Whittle–Matérn field $\mathcal{Z}^{\beta}$ in (4.1) is a continuous kernel and $\dot{H}^{2\beta}_{L}$ is the corresponding reproducing kernel Hilbert space, cf. [37].

5. Spectral Galerkin approximations

In this section we investigate convergence of spectral Galerkin approximations for the Whittle–Matérn field $\mathcal{Z}^{\beta}$ in (4.1). Recall that the covariance structure of the GRF $\mathcal{Z}^{\beta}$ is uniquely determined via its color (3.1) given by the negative fractional power $L^{-\beta}$ of the second-order differential operator $L$ in (2.2) which is defined with respect to the bounded spatial domain $\mathcal{D}\subset\mathbb{R}^{d}$.

For $N\in\mathbb{N}$, the spectral Galerkin approximation $\mathcal{Z}^{\beta}_{N}$ of $\mathcal{Z}^{\beta}$ is ($\mathbb{P}$-a.s.) defined by

[TABLE]

i.e., it is a GRF colored by the finite-rank operator

[TABLE]

mapping to the finite-dimensional subspace $V_{N}:=\operatorname{span}\{e_{1},\ldots,e_{N}\}$ generated by the first $N$ eigenvectors of $L$ corresponding to the eigenvalues $0<\lambda_{1}\leq\ldots\leq\lambda_{N}$.

The following three corollaries, which provide explicit convergence rates of these approximations and their covariance functions with respect to the truncation parameter $N$, are consequences of the Propositions 3.4, 3.6 and 3.7. We first formulate the results in the Sobolev norms.

Corollary 5.1**.**

Suppose Assumptions 2.1.I–II and that $d\in\mathbb{N}$, $\sigma\geq 0$, $\beta,q\in(0,\infty)$. Let $\mathcal{Z}^{\beta}$ be the Whittle–Matérn field in (4.1) and, for $N\in\mathbb{N}$, let $\mathcal{Z}^{\beta}_{N}$ be the spectral Galerkin approximation in (5.1). If $2\beta-\sigma>\nicefrac{{d}}{{2}}$, the following bounds hold,

[TABLE]

where $\varrho^{\beta},\varrho^{\beta}_{N}$ denote the covariance functions of $\mathcal{Z}^{\beta}$ and $\mathcal{Z}^{\beta}_{N}$, respectively, cf. (3.10).

If, in addition, Assumption 2.3.I and $0\leq\sigma\leq 1$ (or Assumptions 2.1.I–III, 2.3.II, and $0\leq\sigma\leq 2$) are satisfied, then the assertions (5.3)–(5.4) remain true if we formulate them with respect to the Sobolev norms $\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}$, $\|\,\cdot\,\|_{H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})}$.

Proof.

The estimates (5.3)/(5.4) follow from (3.19)/(3.20) of Proposition 3.7 with $\mathcal{Z}:=\mathcal{Z}^{\beta}$, $\widetilde{\mathcal{Z}}:=\mathcal{Z}^{\beta}_{N}$, $T:=L^{-\beta}$, and $\widetilde{T}:=L_{N}^{-\beta}$ by exploiting the spectral behavior (2.3) from Lemma 2.2. Finally, applying Lemma 2.4 proves the last claim of this proposition. ∎

By Proposition 3.6 we furthermore obtain the following convergence result as $N\to\infty$ for the covariance function $\varrho^{\beta}_{N}$ in the $L_{\infty}(\mathcal{D}\times\mathcal{D})$-norm.

Corollary 5.2**.**

Suppose Assumptions 2.1.I–II and that the system $\mathcal{E}=\{e_{j}\}_{j\in\mathbb{N}}$ of $L_{2}(\mathcal{D})$-orthonormal eigenvectors of the operator $L$ in (2.2) is uniformly bounded in $C(\overline{\mathcal{D}})$ as in (4.6). Then, for $\beta>\nicefrac{{d}}{{4}}$, the covariance functions of $\mathcal{Z}^{\beta}$ in (4.1) and of $\mathcal{Z}^{\beta}_{N}$ in (5.1) have continuous representatives $\varrho^{\beta},\varrho^{\beta}_{N}\colon\overline{\mathcal{D}\times\mathcal{D}}\to\mathbb{R}$, and

[TABLE]

Proof.

By Lemma 4.3, $\varrho^{\beta}$ and $\varrho^{\beta}_{N}$ have continuous representatives. In addition, the estimate (3.12) from Proposition 3.6 proves (5.5) since, for all $x\in\overline{\mathcal{D}}$, $\varphi\in C(\overline{\mathcal{D}})^{*}$,

[TABLE]

Finally, for $\beta>\nicefrac{{d}}{{4}}$, the spectral behavior (2.3) of $L$ from Lemma 2.2 yields

[TABLE]

If Assumption (i) or (ii) of Lemma 4.2 is satisfied, we obtain not only Sobolev regularity of the GRF $\mathcal{Z}^{\beta}$ in ($L_{q}(\Omega)$-sense), but also Hölder continuity. The next proposition shows that in this case the sequence of spectral Galerkin approximations \bigl{(}\mathcal{Z}^{\beta}_{N}\bigr{)}_{N\in\mathbb{N}} converges also with respect to these norms.

Corollary 5.3**.**

Suppose that $d\in\{1,2,3\}$, $\gamma\in(0,1)$ satisfy Assumption (i) or (ii) of Lemma 4.2, and let $L,L_{N}^{-\beta}$ be the operators in (2.2) and (5.2). Then, for every $N\in\mathbb{N}$ and $2\beta\geq\gamma+\nicefrac{{d}}{{2}}$, there exist continuous random fields $\mathcal{Z}^{\beta},\mathcal{Z}^{\beta}_{N}\colon\overline{\mathcal{D}}\times\Omega\rightarrow\mathbb{R}$, colored by $L^{-\beta}$ and $L_{N}^{-\beta}$, respectively, such that

[TABLE]

for every $\delta\in(0,\gamma)$ and $q\in(0,\infty)$.

Proof.

By Lemma 4.2 there exist continuous random fields $\mathcal{Z}^{\beta},\mathcal{Z}^{\beta}_{N}\colon\overline{\mathcal{D}}\times\Omega\rightarrow\mathbb{R}$ colored by $L^{-\beta}$ and $L_{N}^{-\beta}$, respectively. Their difference $\mathcal{Z}^{\beta}-\mathcal{Z}^{\beta}_{N}$ is then a continuous random field colored by $T_{N}:=L^{-\beta}-L_{N}^{-\beta}=(L-L_{N})^{-\beta}$ and we obtain the convergence result in (5.6) from the stability estimate (4.5) of Lemma 4.2 applied to $\mathcal{Z}^{\beta}-\mathcal{Z}^{\beta}_{N}$, since, for every $\psi\in L_{2}(\mathcal{D})$,

[TABLE]

Here, we have used the spectral behavior (2.3) from Lemma 2.2 for $\lambda_{N}$. ∎

6. General Galerkin approximations

After having derived error estimates for spectral Galerkin approximations in the previous subsection, we now consider a family of general Galerkin approximations for the Whittle–Matérn field $\mathcal{Z}^{\beta}$ in (4.1) which, for the case $\beta\in(0,1)$, has been proposed in [2, 3]. Recall that the random field $\mathcal{Z}^{\beta}$ is indexed by the bounded spatial domain $\mathcal{D}\subset\mathbb{R}^{d}$.

6.1. Sinc-Galerkin approximations

The approximations proposed in [2, 3] are based on a Galerkin method for the spatial discretization $L_{h}$ of $L$ and a sinc quadrature for an integral representation of the resulting discrete fractional inverse $L_{h}^{-\beta}$. We recall that approach in this subsection, and formulate all assumptions and auxiliary results which are needed for the subsequent error analysis in Subsection 6.2.

6.1.1. Galerkin discretization

We assume that we are given a family $(V_{h})_{h>0}$ of subspaces of $H^{1}_{0}(\mathcal{D})$, with dimension $N_{h}:=\dim(V_{h})<\infty$. We let $\Pi_{h}\colon L_{2}(\mathcal{D})\to V_{h}$ denote the $L_{2}(\mathcal{D})$-orthogonal projection onto $V_{h}$. Since $V_{h}\subset H^{1}_{0}(\mathcal{D})=\dot{H}^{1}_{L}$, $\Pi_{h}$ can be uniquely extended to a bounded linear operator $\Pi_{h}\colon\dot{H}^{-1}_{L}\to V_{h}$. In addition, we let $L_{h}\colon V_{h}\to V_{h}$ be the Galerkin discretization of the differential operator $L$ in (2.2) with respect to $V_{h}$, i.e.,

[TABLE]

We arrange the eigenvalues of $L_{h}$ in nondecreasing order,

[TABLE]

and let $\{e_{j,h}\}_{j=1}^{N_{h}}$ be a set of corresponding eigenvectors which are orthonormal in $L_{2}(\mathcal{D})$. The operator $R_{h}\colon H^{1}_{0}(\mathcal{D})=\dot{H}^{1}_{L}\to V_{h}$ is the Rayleigh–Ritz projection, i.e., $R_{h}:=L_{h}^{-1}\Pi_{h}L$ and, for all $\psi\in\dot{H}^{1}_{L}$,

[TABLE]

All further assumptions on the finite-dimensional subspaces $(V_{h})_{h>0}$ are summarized below and explicitly referred to, when needed in our error analysis.

Assumption 6.1** (on the Galerkin discretization).**

1. I.

There exist $\theta_{1}>\theta_{0}>0$ and a linear operator $\mathcal{I}_{h}\colon H^{\theta_{1}}(\mathcal{D})\to V_{h}$ such that, for all $\theta_{0}<\theta\leq\theta_{1}$, $\mathcal{I}_{h}\colon H^{\theta}(\mathcal{D})\to V_{h}$ is a continuous extension, and

[TABLE]

holds for $0\leq\sigma\leq\min\{1,\theta\}$ and sufficiently small $h>0$. 2. II.

For all $h>0$ sufficiently small and all $0\leq\sigma\leq 1$ the following inverse inequality holds:

[TABLE] 3. III.

$\dim(V_{h})=N_{h}\eqsim_{\mathcal{D}}h^{-d}$ for sufficiently small $h>0$. 4. IV.

There exist $r,s_{0},t,C_{0},C_{\lambda}>0$ such that for all $h>0$ sufficiently small and for all $j\in\{1,\ldots,N_{h}\}$ the following error estimates hold:

[TABLE]

where $\{(\lambda_{j},e_{j})\}_{j\in\mathbb{N}}$ are the eigenpairs of the operator $L$ in (2.2).

We refer to Subsection 6.3 for explicit examples of finite element spaces $(V_{h})_{h>0}$, which satisfy these assumptions.

Remark 6.2*.*

It is a consequence of the min-max principle that the first inequality in (6.5), $\lambda_{j}\leq\lambda_{j,h}$, is satisfied for all conforming Galerkin spaces $V_{h}\subset\dot{H}^{1}_{L}$.

In Theorem 6.6 below, we bound the deterministic Galerkin error in the fractional case, i.e., we consider the error between $L^{-\beta}g$ and $L_{h}^{-\beta}\Pi_{h}g$. This theorem is one of our main results and it will be a crucial ingredient when analyzing general Galerkin approximations of the Whittle–Matérn field $\mathcal{Z}^{\beta}$ from (4.1) in Subsection 6.2. For its derivation, we need the following two lemmata.

Lemma 6.3**.**

Let $L$ be as in (2.2) and, for $h>0$, let $L_{h},R_{h}$ be as in (6.1) and (6.2). Suppose Assumptions 2.1.I–II and 2.3.I. Let $0<\alpha\leq 1$ be such that

[TABLE]

where $\dot{H}^{1+\delta}_{L}$ is defined as in in (2.4). Furthermore, let Assumption 6.1.I be satisfied with parameters $\theta_{0}\in(0,1)$ and $\theta_{1}\geq 1+\alpha$. Then, for every $0\leq\eta\leq\vartheta\leq\alpha$,

[TABLE]

for sufficiently small $h>0$.

Proof.

Since $R_{h}u\in V_{h}$ is the best approximation of $u\in\dot{H}^{1}_{L}$ with respect to $\|\,\cdot\,\|_{1}$, we find by Assumption 6.1.I and the assumed equivalence (6.7) that, for $e:=u-R_{h}u$ and any $0\leq\vartheta\leq\alpha$,

[TABLE]

i.e., (6.8) for $\eta=0$ follows. Furthermore, if we let $\psi:=L^{-\vartheta}e\in\dot{H}^{1+2\vartheta}_{L}$, the estimate above and the orthogonality of $e$ to $V_{h}$ in $\dot{H}^{1}_{L}$, combined with (2.6), Assumption 6.1.I and (6.7) yield

[TABLE]

which proves (6.8) for $\eta=\vartheta$ since $\|\psi\|_{1+\vartheta}=\|e\|_{1-\vartheta}$. For $\eta\in(0,\vartheta)$, the result (6.8) holds by interpolation. Now let $g\in\dot{H}^{\vartheta-1}_{L}$ be given. Then, (6.9) follows from (6.8) for $u:=L^{-1}g\in\dot{H}^{1+\vartheta}_{L}$, since $\|u\|_{1+\vartheta}=\|g\|_{\vartheta-1}$. ∎

Remark 6.4*.*

Note that if Assumptions 2.1.I–III and 2.3.II are satisfied, i.e., if the coefficient $A$ of the operator $L$ in (2.2) is Lipschitz continuous and the domain $\mathcal{D}$ is convex, then the equivalence (6.7) with $\alpha=1$ is part of Lemma 2.4, see (2.7).

Lemma 6.5**.**

Suppose Assumptions 2.1.I–II and 2.3.I. Let $L$ be as in (2.2) and, for $h>0$, let $L_{h}$ be as in (6.1). Then, for each $0\leq\gamma\leq\nicefrac{{1}}{{2}}$, we have

[TABLE]

Furthermore, if the $L_{2}(\mathcal{D})$-orthogonal projection $\Pi_{h}$ is $H^{1}(\mathcal{D})$-stable, i.e., if there exists a constant $C_{\Pi}>0$ such that

[TABLE]

for all sufficiently small $h>0$, then, for such $h>0$ and all $0\leq\gamma\leq\nicefrac{{1}}{{2}}$,

[TABLE]

If additionally Assumption 6.1.II is satisfied and if $0<\alpha\leq 1$ is as in (6.7), then (6.12) holds for $0\leq\gamma\leq\nicefrac{{(1+\alpha)}}{{2}}$.

Proof.

For $g\in L_{2}(\mathcal{D})=\dot{H}^{0}_{L}$, we find by the definition (6.1) of $L_{h}$ that

[TABLE]

Thus, (6.10) holds for $\gamma\in\{0,\nicefrac{{1}}{{2}}\}$. In other words, the canonical embedding $I_{h}$ of $V_{h}$ into $L_{2}(\mathcal{D})$ is a continuous mapping from $\dot{H}^{2\gamma}_{h}$ to $\dot{H}^{2\gamma}_{L}$, for $\gamma\in\{0,\nicefrac{{1}}{{2}}\}$, where $\dot{H}^{2\gamma}_{h}$ denotes the space $V_{h}$ equipped with the norm $\|\,\cdot\,\|_{\dot{H}^{2\gamma}_{h}}:=\|L_{h}^{\gamma}\,\cdot\,\|_{L_{2}(\mathcal{D})}$. Thus,

[TABLE]

follows by interpolation for all $0\leq\gamma\leq\nicefrac{{1}}{{2}}$, which completes the proof of (6.10).

If $\Pi_{h}$ is $H^{1}(\mathcal{D})$-stable, by Lemma 2.4 we have $\|\Pi_{h}\|_{\mathcal{L}\left(\dot{H}^{1}_{L}\right)}\lesssim_{(A,\kappa,\mathcal{D})}C_{\Pi}$, and

[TABLE]

follows, i.e., (6.12) holds for $\gamma\in\{0,\nicefrac{{1}}{{2}}\}$. By interpreting this result as continuity of $\Pi_{h}$ as a mapping from $\dot{H}^{2\gamma}_{L}$ to $\dot{H}^{2\gamma}_{h}$, again by interpolation, we obtain (6.12) for all $0\leq\gamma\leq\nicefrac{{1}}{{2}}$. Finally, if $\gamma=\nicefrac{{(1+\vartheta)}}{{2}}$ for some $0<\vartheta\leq\alpha$, we use the identity

[TABLE]

where $R_{h}=L_{h}^{-1}\Pi_{h}L$ is the Rayleigh–Ritz projection (6.2). Since $0<\vartheta\leq\alpha\leq 1$, we obtain for the first term by (6.10) that

[TABLE]

To estimate the second term, we write $E_{h}^{R}:=\operatorname{Id}_{\dot{H}^{1+\vartheta}_{L}}-R_{h}$. Then,

[TABLE]

Here, $\|L_{h}^{\nicefrac{{\vartheta}}{{2}}}\Pi_{h}L^{-\nicefrac{{\vartheta}}{{2}}}\|_{\mathcal{L}(L_{2}(\mathcal{D}))}\lesssim_{(\gamma,A,\kappa,\mathcal{D})}1$, since $0<\vartheta=2\gamma-1\leq 1$, and we can use Assumption 6.1.II, (6.11), and (6.8) to conclude for $\vartheta\neq\nicefrac{{1}}{{2}}$ as follows,

[TABLE]

If $\gamma=\nicefrac{{3}}{{4}}$ and, thus, $\vartheta=\nicefrac{{1}}{{2}}$, a slight modification completes the proof of (6.12) for all $\nicefrac{{1}}{{2}}<\gamma\leq\nicefrac{{(1+\alpha)}}{{2}}$. ∎

Theorem 6.6**.**

Let $L$ be as in (2.2) and, for $h>0$, let $L_{h}$ be as in (6.1). Suppose Assumptions 2.1.I–II, 2.3.I, 6.1.II and that $\Pi_{h}$ is $H^{1}(\mathcal{D})$-stable, see (6.11). Let Assumption 6.1.I be satisfied with parameters $\theta_{0}\in(0,1)$ and $\theta_{1}\geq 1+\alpha$, where $0<\alpha\leq 1$ is as in (6.7). Assume further that $\beta>0$, $0\leq\sigma\leq 1$, and $-1\leq\delta\leq 1+\alpha$ are such that $2\beta+\delta-\sigma>0$ and $2\beta-\sigma\geq 0$. Then, for all $g\in\dot{H}^{\delta}_{L}$, we have

[TABLE]

for arbitrary $\varepsilon>0$ and all $h>0$ sufficiently small.

Remark 6.7* (Sobolev bounds).*

By (2.6) of Lemma 2.4 and under the assumption given by (6.7), the result (6.13) implies an error bound with respect to the Sobolev norms, for all $0\leq\sigma\leq 1$ and $-1\leq\delta\leq 1+\alpha$, $\delta\neq\nicefrac{{1}}{{2}}$. Namely, for all $g\in H^{\delta}(\mathcal{D})$,

[TABLE]

for any $\varepsilon>0$ and all $h>0$ sufficiently small.

Remark 6.8* (Comparison with [5]).*

For the specific case $\beta\in(0,1)$, $\sigma=0$, and $\delta\geq 0$ the error in (6.13) has already been investigated in [5], where $(V_{h})_{h>0}$ are chosen as finite element spaces with continuous piecewise linear basis functions, defined with respect to a quasi-uniform family of triangulations $(\mathcal{T}_{h})_{h>0}$ of $\mathcal{D}$. If $g\in\dot{H}^{\delta}_{L}$, $\delta\geq 0$ and $\alpha<\beta$, the results of [5, Thm. 4.3] show convergence at the rate $2\alpha$, in accordance with (6.13). For $\alpha\geq\beta$ and $g\in\dot{H}^{\delta}_{L}$, by [5, Thm. 4.3 & Rem. 4.1]

[TABLE]

i.e., compared to (6.13), one obtains a log-term $\ln(1/h)$ instead of $h^{-\varepsilon}$ in the first case. We point out that the purpose of Theorem 6.6 was to allow for all $\beta>0$ and, in addition, for the wider range of parameters: $0\leq\sigma\leq 1$ and $-1\leq\delta\leq 1+\alpha$.

Remark 6.9* ($p$-FEM).*

Due to the term $2\alpha$ and $0<\alpha\leq 1$, (6.13) will be sharp for finite elements of first order, but not for finite elements of polynomial degree $p\geq 2$ when $\beta>1$ and the problem is “smooth” such that (6.7) holds with $\alpha>1$.

Proof of Theorem 6.6.

We first prove (6.13) for $0\leq\delta\leq 1+\alpha$. To this end, let $\beta>0$ and $0\leq\sigma\leq\min\{2\beta,1\}$ satisfying $2\beta+\delta>\sigma$ be given. Without loss of generality we may assume that $\varepsilon\in\left(0,2\beta+\delta-\sigma-\alpha\mathds{1}_{\{2\beta+\delta-\sigma-\alpha>0\}}\right)$. We write $I:=\operatorname{Id}_{L_{2}(\mathcal{D})}$ and split

[TABLE]

In order to estimate term (A), we first note that by Assumption 6.1.I, with $\theta=1+\alpha$, and by (6.7) we have, for $h>0$ sufficiently small,

[TABLE]

since $\Pi_{h}g\in V_{h}$ is the $L_{2}(\mathcal{D})$-best approximation of $g\in H^{\theta}(\mathcal{D})$. Furthermore, $\|I-\Pi_{h}\|_{\mathcal{L}(L_{2}(\mathcal{D}))}=1$, and by interpolation

[TABLE]

By exploiting the identity

[TABLE]

which holds for all $\phi,\psi\in L_{2}(\mathcal{D})$, we thus obtain, for all $h>0$ sufficiently small,

[TABLE]

where we set $\theta:=\min\{2\beta-\sigma,1+\alpha\}$ and, hence, $0\leq\theta,\delta\leq 1+\alpha$.

For deriving a bound for (B), we first note that by (6.12) of Lemma 6.5

[TABLE]

Next, we define the contour

[TABLE]

where $\omega\in(0,\pi)$ and $r:=\nicefrac{{\lambda_{1}}}{{2}}$. By, e.g., [32, Ch. 2.6, Eq. (6.3)] we have

[TABLE]

From the limit $\omega\to\pi$, we then obtain the representation

[TABLE]

which applied to $L^{-\beta}$ and $L_{h}^{-\beta}$ (recall that $\lambda_{1,h}\geq\lambda_{1}$, see Remark 6.2) implies that

[TABLE]

We exploit this integral representation as well as the identity

[TABLE]

which holds for any $z\in\mathbb{C}$, and bound term (B) as follows

[TABLE]

where $\mu:=\nicefrac{{(1+\eta+\sigma)}}{{2}}$, $\nu:=\nicefrac{{(1+\vartheta-\delta)}}{{2}}$ and $0\leq\eta\leq\vartheta\leq\alpha$ are chosen as follows

[TABLE]

By (6.12) and (6.9), we find for the term outside of the integral,

[TABLE]

for $h>0$ sufficiently small, where these three cases can be summarized as in (6.13), since $2\beta+\delta-\sigma-\varepsilon<\alpha\leq 1+\alpha-\sigma$ for all $0\leq\sigma\leq 1$ if $2\beta+\delta-\sigma\leq\alpha$ and $2\alpha<2\beta+\delta-\sigma-\varepsilon$ for $\varepsilon>0$ sufficiently small if $2\beta+\delta-\sigma>2\alpha$. It remains to show that the two integrals in (6.14) converge, uniformly in $h$. To this end, we first note that $0\leq\mu\leq 1$ and, thus, for any $t>0$,

[TABLE]

By the same argument we find that \bigl{\|}L_{h}^{\nu}(tI+L_{h})^{-1}\Pi_{h}\bigr{\|}_{\mathcal{L}(L_{2}(\mathcal{D}))}\leq t^{\nu-1}, for $t>0$, since also $0\leq\nu\leq 1$. Thus, we can bound the first integral arising in (6.14) by

[TABLE]

Here, we have used that $r=\nicefrac{{\lambda_{1}}}{{2}}$, $\mu+\nu-2-\beta=-1+\nicefrac{{(\eta+\vartheta+\sigma-\delta-2\beta)}}{{2}}\leq-1-\nicefrac{{\varepsilon}}{{2}}<-1$ if $2\beta+\delta-\sigma\leq 2\alpha$, and $\mu+\nu-2-\beta\leq-1-(\beta+\nicefrac{{\delta}}{{2}}-\nicefrac{{\sigma}}{{2}}-\alpha)<-1$ if $2\beta+\delta-\sigma>2\alpha$. To estimate the second integral in (6.14), we note that, for any $z\in\mathbb{C}$ with $|z|=\nicefrac{{\lambda_{1}}}{{2}}$,

[TABLE]

since $(x+y)^{\mu}\leq x^{\mu}+y^{\mu}$ if $0\leq\mu\leq 1$ and $x,y\geq 0$. Similarly, for $0\leq\nu\leq 1$,

[TABLE]

With these observations, we finally can bound the second integral in (6.14),

[TABLE]

which completes the proof of (6.13) for the case that $0\leq\delta\leq 1+\alpha$.

Assume now that $\delta=-\widetilde{\sigma}$ for some $0<\widetilde{\sigma}\leq 1$. Then, for $g\in\dot{H}^{\delta}_{L}$,

[TABLE]

After rewriting,

[TABLE]

we may exploit (6.13), which has already been proven for $0\leq\delta\leq 1+\alpha$, as follows,

[TABLE]

since $\widetilde{\delta}:=2\beta-\widetilde{\sigma}-\sigma=2\beta+\delta-\sigma>0$ by assumption. Furthermore, by (6.10) of Lemma 6.5 we have \bigl{\|}L_{h}^{-\nicefrac{{\widetilde{\sigma}}}{{2}}}\Pi_{h}L^{\nicefrac{{\widetilde{\sigma}}}{{2}}}\bigr{\|}_{\mathcal{L}(L_{2}(\mathcal{D}))}\lesssim_{\widetilde{\sigma}}1. We conclude that

[TABLE]

for the whole range of parameters $\sigma,\delta$ as stated in the theorem. ∎

6.1.2. Sinc quadrature and fully discrete scheme

After the Galerkin discretization (in space), we need a second component to approximate the generalized Whittle–Matérn field $\mathcal{Z}^{\beta}$ in (4.1). Namely, we have to numerically realize a fractional inverse of the Galerkin operator $L_{h}$ in (6.1). To this end, as proposed in [2], we introduce, for $\beta\in(0,1)$ and $k>0$, the sinc quadrature approximation of $L_{h}^{-\beta}$ from [5],

[TABLE]

where K^{-}:=\bigl{\lceil}\tfrac{\pi^{2}}{4\beta k^{2}}\bigr{\rceil}, K^{+}:=\bigl{\lceil}\frac{\pi^{2}}{4(1-\beta)k^{2}}\bigr{\rceil}. We also formally define this operator for the case $\beta=0$ by setting $Q^{0}_{h,k}:=\operatorname{Id}_{V_{h}}$.

For a general $\beta=n_{\beta}+\beta_{\star}>0$ as in (4.2), we then consider the approximations $\widetilde{\mathcal{Z}}^{\beta}_{h,k},\mathcal{Z}^{\beta}_{h,k}\colon\mathcal{D}\times\Omega\to\mathbb{R}$ of the Whittle–Matérn field $\mathcal{Z}^{\beta}$ in (4.1) which are ($\mathbb{P}$-a.s.) defined by

[TABLE]

i.e., $\mathcal{Z}^{\beta}_{h,k}$, $\widetilde{\mathcal{Z}}^{\beta}_{h,k}$ are GRFs colored by $Q^{\beta_{\star}}_{h,k}L_{h}^{-n_{\beta}}\Pi_{h}$ and $Q^{\beta_{\star}}_{h,k}L_{h}^{-n_{\beta}}\widetilde{\Pi}_{h}$, respectively, cf. Definition 3.2. Here, the finite-rank operator $\widetilde{\Pi}_{h}$ is given by

[TABLE]

For $\beta\in(0,1)$, the construction (6.17) of $\widetilde{\mathcal{Z}}^{\beta}_{h,k}$ gives the same approximation as considered in [2, 3]. Note furthermore that, in contrast to $\Pi_{h}$, the operator $\widetilde{\Pi}_{h}$ in (6.18) is neither a projection nor self-adjoint, and its definition depends on the particular choice of the eigenbases $\{e_{j}\}_{j\in\mathbb{N}}\subset L_{2}(\mathcal{D})$ and $\{e_{j,h}\}_{j=1}^{N_{h}}\subset V_{h}$. The reason for why we consider both approximations $\mathcal{Z}^{\beta}_{h,k}$, $\widetilde{\mathcal{Z}}^{\beta}_{h,k}$ will become apparent in the error analysis of Subsection 6.2. Although, in general, they do not coincide in $L_{q}(\Omega;L_{2}(\mathcal{D}))$-sense, i.e.,

[TABLE]

they have the same Gaussian distribution as shown in the following lemma.

Lemma 6.10**.**

Suppose Assumptions 2.1.I–II. Let $\Pi_{h}$ denote the $L_{2}(\mathcal{D})$-orthogonal projection onto $V_{h}$, and $\widetilde{\Pi}_{h}$ be the operator in (6.18). Then, if $T_{h}\in\mathcal{L}(V_{h})$,

[TABLE]

holds for all $\phi,\psi\in L_{2}(\mathcal{D})$. In particular, $\mathcal{Z}^{\beta}_{h,k}\overset{d}{=}\widetilde{\mathcal{Z}}^{\beta}_{h,k}$ as $L_{2}(\mathcal{D})$-valued random variables, where $\mathcal{Z}^{\beta}_{h,k}$ and $\widetilde{\mathcal{Z}}^{\beta}_{h,k}$ are as defined in (6.16)–(6.17).

Proof.

Note that $(T_{h}\Pi_{h})^{*}$ (resp. $(T_{h}\widetilde{\Pi}_{h})^{*}$) denotes the adjoint of $T_{h}\Pi_{h}$ (resp. of $T_{h}\widetilde{\Pi}_{h}$) when interpreted as an operator in $\mathcal{L}(L_{2}(\mathcal{D}))$. This means, we are identifying $T_{h}\Pi_{h}$ with $I_{h}T_{h}\Pi_{h}$ (resp. $T_{h}\widetilde{\Pi}_{h}$ with $I_{h}T_{h}\widetilde{\Pi}_{h}$), where $I_{h}$ denotes the canonical embedding of $V_{h}$ into $L_{2}(\mathcal{D})$. Since $I_{h}^{*}=\Pi_{h}$ we thus find that $(T_{h}\Pi_{h})^{*}=T_{h}^{*}\Pi_{h}$ and $(T_{h}\widetilde{\Pi}_{h})^{*}=\widetilde{\Pi}_{h}^{*}T_{h}^{*}\Pi_{h}$, which combined with $\widetilde{\Pi}_{h}\widetilde{\Pi}_{h}^{*}=\operatorname{Id}_{V_{h}}$ proves (6.19).

By definition of $\mathcal{Z}_{h,k}^{\beta}$, $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ in (6.16)–(6.17), for $M\in\mathbb{N}$ and $\psi_{1},\ldots,\psi_{M}\in L_{2}(\mathcal{D})$, the random vectors $\mathbf{z}$, $\widetilde{\mathbf{z}}$ with entries z_{j}=\bigl{(}\mathcal{Z}_{h,k}^{\beta},\psi_{j}\bigr{)}_{L_{2}(\mathcal{D})} and \widetilde{z}_{j}=\bigl{(}\widetilde{\mathcal{Z}}_{h,k}^{\beta},\psi_{j}\bigr{)}_{L_{2}(\mathcal{D})}, $1\leq j\leq M$, are multivariate Gaussian distributed. Furthermore, both vanish in expectation and their covariance matrices, $\mathbf{C}:=\operatorname{Cov}(\mathbf{z})$ and \widetilde{\mathbf{C}}:=\operatorname{Cov}\bigl{(}\widetilde{\mathbf{z}}\bigr{)}, coincide due to (6.19) applied to $T_{h}:=Q^{\beta_{\star}}_{h,k}L_{h}^{-n_{\beta}}$. This shows that $\mathcal{Z}_{h,k}^{\beta}\overset{d}{=}\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ as $L_{2}(\mathcal{D})$-valued random variables. ∎

Remark 6.11* (Simulation in practice).*

To simulate samples of the in (6.16)–(6.17) abstractly defined ($\mathbb{P}$-a.s.) $V_{h}$-valued Gaussian random variables $\mathcal{Z}_{h,k}^{\beta}$ and $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ in practice, in both cases, one first has to generate a sample of a multivariate Gaussian random vector $\mathbf{b}$ with mean $\mathbf{0}$ and covariance matrix $\mathbf{M}$, where $\mathbf{M}$ is the Gramian with respect to any fixed basis $\Phi_{h}=\{\phi_{j,h}\}_{j=1}^{N_{h}}$ of $V_{h}$, i.e., $M_{ij}:=(\phi_{i,h},\phi_{j,h})_{L_{2}(\mathcal{D})}$. This follows from the identical distribution of the GRFs $\mathcal{Z}^{0}_{h}$ and $\widetilde{\mathcal{Z}}^{0}_{h}$ colored by $\Pi_{h}$ and $\widetilde{\Pi}_{h}$, respectively, as well as from the chain of equalities

[TABLE]

which we obtain from Lemma 6.10 with $T_{h}:=\operatorname{Id}_{V_{h}}$. Since

[TABLE]

the random vector $\mathbf{Z}_{k}^{\beta}$, given by

[TABLE]

is then the vector of coefficients when expressing the $V_{h}$-valued sample of $\mathcal{Z}_{h,k}^{\beta}$ (or of $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$) with respect to the basis $\Phi_{h}$. Here, $\mathbf{L}\in\mathbb{R}^{N_{h}\times N_{h}}$ represents the action of the Galerkin operator $L_{h}$ in (6.1), i.e., $L_{ij}:=(L_{h}\phi_{j,h},\phi_{i,h})_{L_{2}(\mathcal{D})}$, and, for $\beta_{\star}\in(0,1)$, $\mathbf{Q}^{\beta_{\star}}_{k}\in\mathbb{R}^{N_{h}\times N_{h}}$ is the matrix analog of the operator $Q_{h,k}^{\beta_{\star}}$ from (6.15), i.e.,

[TABLE]

6.2. Error analysis

The errors $\mathcal{Z}^{\beta}-\mathcal{Z}_{h,k}^{\beta}$ and $\mathcal{Z}^{\beta}-\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ of the approximations in (6.16)–(6.17) compared to the true Whittle–Matérn field $\mathcal{Z}^{\beta}$ from (4.1) are GRFs colored (see Definition 3.2) by

[TABLE]

respectively. In order to perform the error analysis for $\mathcal{Z}_{h,k}^{\beta}$ and $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$, we split these operators as follows

[TABLE]

where $\widetilde{E}_{N_{h}}^{\beta}:=L^{-\beta}-L_{N_{h}}^{-\beta}$ is a dimension truncation error (recall the finite-rank operator $L_{N_{h}}^{-\beta}$ from (5.2)) which can be estimated with the results from Section 5 on spectral Galerkin approximations. Furthermore, we shall refer to

[TABLE]

as the Galerkin errors and as the quadrature errors, respectively.

In the following we provide error estimates for both approximations, $\mathcal{Z}_{h,k}^{\beta}$ and $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ in (6.16)–(6.17), with respect to the norm on $L_{q}(\Omega;H^{\sigma}(\mathcal{D}))$ as well as for its covariance functions $\varrho_{h,k}^{\beta}$, $\widetilde{\varrho}_{h,k}^{\beta}$ in the mixed Sobolev norm, cf. (3.14). By exploiting Theorem 6.6 the bounds for $\mathcal{Z}_{h,k}^{\beta}$ and $\varrho_{h,k}^{\beta}$ in Proposition 6.12 below will be sharp if a conforming finite element method with piecewise linear basis functions is used. However, to derive optimal rates for the case of finite elements of higher polynomial degree, a different approach will be necessary, cf. Remark 6.9. To this end, we perform an error analysis for $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ and $\widetilde{\varrho}_{h,k}^{\beta}$ based on spectral expansions, see Proposition 6.13. Since these arguments work only if the differential operator $L$ in (2.2) is at least $H^{2}(\mathcal{D})$-regular, both approaches and results are needed for a complete discussion of smooth vs. $H^{1+\alpha}(\mathcal{D})$-regular problems in Subsection 6.3. Finally, in Proposition 6.14, we use the approximation $\mathcal{Z}_{h,k}^{\beta}$ from (6.16) to formulate a convergence result with respect to the Hölder norm (1.8) in $L_{q}(\Omega)$-sense and with respect to the $L_{\infty}(\mathcal{D}\times\mathcal{D})$-norm for its covariance function $\varrho^{\beta}_{h,k}$.

We note that, at the cost of other assumptions (e.g., $\alpha>\nicefrac{{1}}{{2}}$) on the parameters involved, it is possible to circumvent the additional condition $\beta>1$ (instead of $\beta>\nicefrac{{3}}{{4}}$) needed in the following proposition for the $L_{q}(\Omega;H^{\sigma}(\mathcal{D}))$-estimate if $d=3$.

Proposition 6.12**.**

Suppose Assumptions 2.1.I–II, 2.3.I, 6.1.II–III, and let Assumption 6.1.I be satisfied with parameters $\theta_{0}\in(0,1)$ and $\theta_{1}\geq 1+\alpha$, where $0<\alpha\leq 1$ is as in (6.7). Assume furthermore that $\Pi_{h}$ is $H^{1}(\mathcal{D})$-stable, see (6.11), and that $d\in\{1,2,3\}$, $\beta>0$ and $0\leq\sigma\leq 1$ are such that $2\beta-\sigma>\nicefrac{{d}}{{2}}$. Let $\mathcal{Z}^{\beta}$ be the Whittle–Matérn field in (4.1) and, for $h,k>0$, let $\mathcal{Z}^{\beta}_{h,k}$ be the sinc-Galerkin approximation in (6.16), with covariance functions $\varrho^{\beta}$ and $\varrho_{h,k}^{\beta}$, respectively. Then, for every $q,\varepsilon>0$ and sufficiently small $h>0$,

[TABLE]

where, if $d=3$, for (6.24) to hold, we also suppose that $\beta>1$ and $\alpha\geq\nicefrac{{1}}{{2}}-\sigma$.

Proof.

We start with splitting the errors with respect to the norms on $\dot{H}^{\sigma}_{L}$, cf. (2.4),

[TABLE]

and on $\dot{H}^{\sigma,\sigma}_{L}$, see (3.15), respectively,

[TABLE]

which by (2.6) of Lemma 2.4 bound the errors (6.24)–(6.25) in the Sobolev norms. Here $\mathcal{Z}^{\beta}_{h}$ denotes a GRF colored by $L_{h}^{-\beta}\Pi_{h}$, with covariance function $\varrho^{\beta}_{h}$. Furthermore, we note the following: For $m\geq 0$, we have

[TABLE]

where the observation of Remark 6.2 was used in the last step. Thus, by the spectral asymptotics from Lemma 2.2 and by Assumption 6.1.III we have for $m\geq 0$, $m\neq\nicefrac{{d}}{{4}}$,

[TABLE]

For terms ($\text{A}_{\mathcal{Z}}$) and ($\text{B}_{\mathcal{Z}}$), we obtain with the definitions of the Galerkin and quadrature errors ${E}_{V_{h}}^{\beta},{E}_{Q}^{\beta}$ from (6.22)–(6.23) by (3.19) of Proposition 3.7 that

[TABLE]

For bounding term ($\text{A}_{\mathcal{Z}}$), we let $\gamma\in(0,\beta)$ and rewrite $E_{V_{h}}^{\beta}$ from (6.22) as follows,

[TABLE]

We first bound ($\text{A}_{\mathcal{Z}}$) for $d\in\{1,2\}$. To this end, let $\varepsilon_{0}>0$ be chosen sufficiently small such that $2\beta-\sigma-\nicefrac{{d}}{{2}}>4\varepsilon_{0}$ and choose $\gamma:=\nicefrac{{d}}{{4}}+\varepsilon_{0}$ in (6.28). We obtain thus \text{(\text{A}{\mathcal{Z}})}\lesssim_{q}\text{(\text{A}^{\prime}{\mathcal{Z}})}+\text{(\text{A}^{\prime\prime}_{\mathcal{Z}})}, where

[TABLE]

For ($\text{A}^{\prime}_{\mathcal{Z}}$), we find by (6.13) of Theorem 6.6 and by (6.27), applied for the parameters $\beta^{\prime}:=\beta-\nicefrac{{d}}{{4}}-\varepsilon_{0}$, $\sigma^{\prime}:=\sigma$, $\delta^{\prime}:=0$, and $m=\nicefrac{{d}}{{4}}+\varepsilon_{0}$, respectively,

[TABLE]

for any $\varepsilon^{\prime}>0$ and sufficiently small $h>0$.

After rewriting term ($\text{A}^{\prime\prime}_{\mathcal{Z}}$) we again apply (6.13) of Theorem 6.6, this time for the parameters $\beta^{\prime\prime}:=\nicefrac{{d}}{{4}}+\varepsilon_{0}>0$, $\sigma^{\prime\prime}:=0$, and $\delta^{\prime\prime}:=\min\{2\beta-\sigma-d-4\varepsilon_{0},1+\alpha\}$. Note that, due to the choice of $\varepsilon_{0}>0$ and since $d\in\{1,2\}$, we have $\delta^{\prime\prime}>-1$ and

[TABLE]

We thus find that, for any $\varepsilon^{\prime\prime}>0$ and sufficiently small $h>0$,

[TABLE]

The Hilbert–Schmidt norm $\|L^{-(\nicefrac{{d}}{{4}}+\varepsilon_{0})}\|_{\mathcal{L}_{2}^{0;0}}$ converges for any $\varepsilon_{0}>0$ due to the spectral asymptotics (2.3) of Lemma 2.2. In addition, since $1+\alpha>\nicefrac{{d}}{{2}}$ for $d\in\{1,2\}$, we find that $1+\alpha+\delta^{\prime\prime}>\min\{2\beta-\sigma-\nicefrac{{d}}{{2}}-4\varepsilon_{0},1+\alpha\}$, and we conclude that

[TABLE]

for sufficiently small $h>0$ and any $\varepsilon>0$ (by adjusting $\varepsilon_{0},\varepsilon^{\prime},\varepsilon^{\prime\prime}>0$).

If $d=3$, let $\varepsilon_{0}>0$ be such that $2\varepsilon_{0}<\min\{2\beta-\sigma-\nicefrac{{3}}{{2}},\beta-1\}$, and choose $\gamma:=\nicefrac{{3}}{{4}}-\nicefrac{{\sigma}}{{2}}+\varepsilon_{0}\in(0,\beta)$ in (6.28). We thus need to bound the terms

[TABLE]

This can be achieved similarly as for $d\in\{1,2\}$ by picking the parameters

[TABLE]

(recall that $\beta>1$ if $d=3$ and, thus, $\delta^{\prime\prime}>-1$). These choices result, for sufficiently small $h>0$, in the estimates

[TABLE]

for all $\varepsilon^{\prime},\varepsilon^{\prime\prime}>0$, where we also have used (6.12) and (6.27) for ($\text{A}^{\prime}_{\mathcal{Z}}$). Finally, since $\alpha\geq\nicefrac{{1}}{{2}}-\sigma$ if $d=3$, we again have $1+\alpha+\delta^{\prime\prime}\geq\min\{2\beta-\sigma-\nicefrac{{d}}{{2}}-4\varepsilon_{0},1+\alpha\}$. Thus, (6.29) also holds for $d=3$.

To estimate ($\text{B}_{\mathcal{Z}}$), we recall the convergence result of the sinc quadrature from [5, Lem. 3.4, Rem. 3.1 & Thm. 3.5]. For sufficiently small $k>0$, we have

[TABLE]

Next, by equivalence of the norms $\|\,\cdot\,\|_{\sigma}$, $\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}$ for $\sigma\in\{0,1\}$, see Lemma 2.4, and by the inverse inequality (6.4) from Assumption 6.1.II, we find, for $\sigma\in\{0,1\}$,

[TABLE]

where we have applied (6.27) with $m=n_{\beta}\in\mathbb{N}_{0}$, $m\neq\nicefrac{{d}}{{4}}$ for $d\in\{1,2,3\}$ in the last step. If $\sigma\in(0,1)$, a respective bound for ($\text{B}_{\mathcal{Z}}$) follows by interpolation.

We proceed with the derivation of (6.25) by estimating ($\text{A}_{\varrho}$) and ($\text{B}_{\varrho}$) in (6.26). By (3.20) of Proposition 3.7 we obtain

[TABLE]

To bound ($\text{A}_{\varrho}$), we let $\varepsilon_{0}>0$ be such that $2\beta-\sigma-\nicefrac{{d}}{{2}}>2\varepsilon_{0}$ and write

[TABLE]

and find therefore that \text{(\text{A}{\varrho})}\leq\text{(\text{A}^{\prime}{\varrho})}+\text{(\text{A}^{\prime\prime}_{\varrho})}, where

[TABLE]

For term ($\text{A}^{\prime}_{\varrho}$), we apply (6.13) of Theorem 6.6, for $\beta^{\prime}:=2\beta-\nicefrac{{\sigma}}{{2}}-\nicefrac{{d}}{{4}}-\varepsilon_{0}$, $\sigma^{\prime}:=\sigma$, and $\delta^{\prime}:=0$. We thus obtain that, for any $\varepsilon^{\prime}>0$ and sufficiently small $h>0$,

[TABLE]

Here, the arising Hilbert–Schmidt norm is bounded by a constant, since

[TABLE]

and boundedness follows from (6.10) and (6.27).

For term ($\text{A}^{\prime\prime}_{\varrho}$), we choose the parameters in (6.13) of Theorem 6.6, as follows: $\beta^{\prime\prime}:=\nicefrac{{\sigma}}{{2}}+\nicefrac{{d}}{{4}}+\varepsilon_{0}$, $\sigma^{\prime\prime}:=\sigma$, and $\delta^{\prime\prime}:=\min\{4\beta-2\sigma-d-4\varepsilon_{0},1+\alpha\}>0$. This gives, for any $\varepsilon^{\prime\prime}>0$ and sufficiently small $h>0$,

[TABLE]

since \bigl{\|}L^{-\nicefrac{{d}}{{4}}-\varepsilon_{0}}\bigr{\|}_{\mathcal{L}_{2}^{0;0}} is bounded due to the spectral asymptotics (2.3) of Lemma 2.2. We conclude that

[TABLE]

for every $\varepsilon>0$ and sufficiently small $h>0$.

Finally, we use the estimate

[TABLE]

as well as the inverse inequality (6.4) to conclude for term ($\text{B}_{\varrho}$) for $\sigma\in\{0,1\}$ that

[TABLE]

Combining the above estimate with (6.30) and stability of the operators

[TABLE]

which is uniform in $h$ and $k$ for sufficiently small $h,k>0$, shows that

[TABLE]

Interpolation for $\sigma\in(0,1)$ completes the proof of (6.25). ∎

Due to the similarity in the derivation with the proof of [2, Thm. 2.10], we have moved the proof of the following proposition to Appendix A.

Proposition 6.13**.**

Suppose Assumptions 2.1.I–II, 2.3.I, and 6.1.II–III. Let Assumption 6.1.IV be satisfied with parameters $r,s_{0},t>0$ such that $\nicefrac{{r}}{{2}}\geq t-1$ and $s_{0}\geq t$. Let $d\in\mathbb{N}$, $\beta>0$ and $0\leq\sigma\leq 1$ be such that $2\beta-\sigma>\nicefrac{{d}}{{2}}$. For $\tau\geq 0$, set

[TABLE]

Furthermore, define, for $0\leq\sigma\leq 1$,

[TABLE]

Let $\mathcal{Z}^{\beta}$ be the Whittle–Matérn field in (4.1) and, for $h,k>0$, let $\widetilde{\mathcal{Z}}^{\beta}_{h,k}$ denote the sinc-Galerkin approximation in (6.17), with covariance functions $\varrho^{\beta}$ and $\widetilde{\varrho}^{\beta}_{h,k}$, respectively. Then, for all $q>0$,

[TABLE]

hold for sufficiently small $h,k>0$, where

[TABLE]

and $\mathcal{P}:=\{C_{0},C_{\lambda},\sigma,\beta,A,\kappa,\mathcal{D}\}$.

Proposition 6.14**.**

Suppose Assumptions 2.1.I–II, 6.1.II–III, and let Assumption 6.1.I be satisfied with parameters $\theta_{0}\in(0,1)$ and $\theta_{1}\geq 1+\alpha$, where $0<\alpha\leq 1$ is as in (6.7). Assume furthermore that $\Pi_{h}$ is $H^{1}(\mathcal{D})$-stable, see (6.11), and that $d=1$, $\beta>0$ and $0<\gamma\leq\nicefrac{{1}}{{2}}$ are such that $2\beta>\gamma+\nicefrac{{1}}{{2}}$. Then, the Whittle–Matérn field $\mathcal{Z}^{\beta}$ in (4.1) and the sinc-Galerkin approximation $\mathcal{Z}_{h,k}^{\beta}$ in (6.16) can be taken as continuous random fields. Moreover, for every $\delta\in(0,\gamma)$, all $\varepsilon,q>0$ and sufficiently small $h>0$, we have

[TABLE]

Here, $\varrho^{\beta},\varrho^{\beta}_{h,k}$ denote the covariance functions of $\mathcal{Z}^{\beta}$ and $\mathcal{Z}^{\beta}_{h,k}$, respectively.

Proof.

Clearly, Q^{\beta_{\star}}_{h,k}L_{h}^{-n_{\beta}}\Pi_{h}\in\mathcal{L}\bigl{(}L_{2}(\mathcal{D});H^{\gamma+\nicefrac{{1}}{{2}}}(\mathcal{D})\bigr{)}, since $Q^{\beta_{\star}}_{h,k}L_{h}^{-n_{\beta}}\Pi_{h}$ is a finite-rank operator and $V_{h}\subset H^{1}_{0}(\mathcal{D})\subset H^{\gamma+\nicefrac{{1}}{{2}}}(\mathcal{D})$ by assumption. Thus, by Corollary 3.5 $\mathcal{Z}^{\beta}_{h,k}$ can be taken as a continuous GRF; and the same is true for the Whittle–Matérn field $\mathcal{Z}^{\beta}$ by Corollary 4.2. Then, $\mathcal{Z}^{\beta}-\mathcal{Z}^{\beta}_{h,k}$ is a continuous random field, colored by $E^{\beta}_{h,k}=E^{\beta}_{V_{h}}+E_{Q}^{\beta}$, see (6.22)–(6.23). Furthermore, by (3.9) and by Lemma 2.4, since $d=1$ and $\nicefrac{{1}}{{2}}<\gamma+\nicefrac{{1}}{{2}}\leq 1$, we have, for $\delta\in(0,\gamma)$ and $q\in(0,\infty)$,

[TABLE]

By (6.13) of Theorem 6.6 we then find, for any $\varepsilon>0$ and sufficiently small $h>0$,

[TABLE]

For $E_{Q}^{\beta}\in\mathcal{L}(V_{h})$ we use the inverse inequality (6.4) as well as the quadrature error estimate from [5, Lem. 3.4, Rem. 3.1 & Thm. 3.5] and obtain

[TABLE]

for sufficiently small $h>0$, which completes the proof of (6.37).

For the $L_{\infty}(\mathcal{D}\times\mathcal{D})$-estimate (6.38) of the covariance function, fix $\varepsilon\in(0,2)$. First, we recall the Sobolev embedding $H^{\nicefrac{{\varepsilon}}{{4}}+\nicefrac{{1}}{{2}}}(\mathcal{D})\hookrightarrow C^{\nicefrac{{\varepsilon}}{{4}}}(\overline{\mathcal{D}})$ as well as the equivalence of the spaces $H^{\nicefrac{{\varepsilon}}{{4}}+\nicefrac{{1}}{{2}}}(\mathcal{D})\cong_{(A,\kappa,\mathcal{D})}\dot{H}^{\nicefrac{{\varepsilon}}{{4}}+\nicefrac{{1}}{{2}}}_{L}$, see Lemma 2.4. We then conclude with (3.12) of Proposition 3.6(ii) that, for $\sigma:=\nicefrac{{1}}{{2}}+\nicefrac{{\varepsilon}}{{4}}\in(\nicefrac{{1}}{{2}},1)$,

[TABLE]

By (6.13) of Theorem 6.6 we have

[TABLE]

Furthermore, we find, similarly as in (6.31), that

[TABLE]

where we have used the inverse inequality (6.4) in the last step. Finally, since \bigl{\|}\bigl{(}E^{\beta}_{Q}\bigr{)}^{*}\bigr{\|}_{\mathcal{L}\left(\dot{H}^{-\sigma}_{L};\dot{H}^{0}_{L}\right)}=\bigl{\|}E^{\beta}_{Q}\bigr{\|}_{\mathcal{L}\left(\dot{H}^{0}_{L};\dot{H}^{\sigma}_{L}\right)}, the proof is completed by (6.39) combined with the uniform stability (6.32) of $L_{h}^{-\beta}$ and $Q_{h,k}^{\beta_{\star}}$. ∎

6.3. Application to finite element approximations

We now discuss different scenarios of

1. (i)

regularity of the second-order differential $L$ in (2.2), 2. (ii)

finite element (FE) discretizations satisfying Assumptions 6.1.I–IV for specific values of $0<\theta_{0}<\theta_{1}$ and of $r,s_{0},t>0$.

We then obtain explicit rates of convergence for the FE Galerkin approximations $\mathcal{Z}_{h,k}^{\beta}$, $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ in (6.16)–(6.17) from Propositions 6.12, 6.13 and 6.14.

Assumption 6.15** (FE discretization).**

Throughout this subsection, we suppose the following setting:

1. $\bullet$

the (minimal) Assumptions 2.1.I–II on the coefficients $A$, $\kappa$ of the operator $L$; 2. $\bullet$

Assumptions 2.3.I, i.e., $\mathcal{D}\subset\mathbb{R}^{d}$ is a bounded Lipschitz domain; 3. $\bullet$

$(\mathcal{T}_{h})_{h>0}$ is a quasi-uniform family of triangulations on $\overline{\mathcal{D}}$, indexed by the mesh width $h>0$; 4. $\bullet$

the basis functions of the finite-dimensional space $V_{h}\subset H^{1}_{0}(\mathcal{D})$ are continuous on $\overline{\mathcal{D}}$ and piecewise polynomial with respect to $\mathcal{T}_{h}$ of degree at most $p\in\mathbb{N}$.

All further assumptions on the operator $L$, on the domain $\mathcal{D}$, and on the FE spaces are explicitly specified for each case. Note that quasi-uniformity of $(\mathcal{T}_{h})_{h>0}$ already guarantees that Assumptions 6.1.II and 6.1.III are satisfied (6.1.III is obvious, for the inverse inequality 6.1.II see, e.g., [15, Cor. 1.141]).

In Subsection 6.3.1 we briefly comment on the situation of smooth coefficients and apply Proposition 6.13 to derive optimal convergence rates when $p\geq 1$. Afterwards, in Subsection 6.3.2 we focus on less regular problems and $p=1$ by using the results from Propositions 6.12 and 6.14.

6.3.1. The smooth case

The remaining crucial ingredient in order to derive explicit rates of convergence from Proposition 6.13 is to prove validity of Assumption 6.1.IV for the finite element spaces $(V_{h})_{h>0}$. For the case of a second-order elliptic differential operator $L$ with smooth coefficients, these results are well-known and we summarize them below.

Assumption 6.16** (smooth case).**

The domain $\mathcal{D}$ has a smooth $C^{\infty}$-boundary $\partial\mathcal{D}$, and the coefficients of $L$ in (2.2) are smooth, i.e., $A\in C^{\infty}(\overline{\mathcal{D}})^{d\times d}$ and $\kappa\in C^{\infty}(\overline{\mathcal{D}})$. Furthermore, the Rayleigh–Ritz projection $R_{h}\colon H^{1}_{0}(\mathcal{D})\to V_{h}$ in (6.2) satisfies the a-priori estimates

[TABLE]

Lemma 6.17**.**

Suppose Assumptions 6.15 and 6.16. Then, Assumption 6.1.IV is satisfied for $r=2p$ and $s_{0}=t=p+1$.

Proof.

See, e.g., [38, Thm. 6.1 & Thm. 6.2]. ∎

Theorem 6.18**.**

Suppose Assumptions 6.15 and 6.16. Let $d\in\mathbb{N}$, $\beta>0$, and $0\leq\sigma\leq 1$ be such that $2\beta-\sigma>\nicefrac{{d}}{{2}}$, let $\mathcal{Z}^{\beta}$ be the Whittle–Matérn field in (4.1) and, for $h,k>0$, let $\widetilde{\mathcal{Z}}^{\beta}_{h,k}$ be the sinc-Galerkin approximation in (6.17), and let $\varrho^{\beta}$, $\widetilde{\varrho}^{\beta}_{h,k}$ denote their covariance functions. Then we have, for sufficiently small $h>0$, sufficiently small $k=k(h)>0$, and all $q>0$,

[TABLE]

where $C_{\beta,h}^{\varrho},C_{\beta,h}^{\mathcal{Z}}$ and $\mathcal{P}$ are as in Proposition 6.13.

Proof.

By Lemma 6.17 we have $r=2p$, $s_{0}=t=p+1$ and, thus, for $\gamma\in\{0,1\}$, $\rho_{\gamma}(\tau)=\min\left\{p+1-\gamma,\,2\beta+\tau-\gamma-\nicefrac{{d}}{{2}}\right\}$ in (6.33). Finally,

[TABLE]

in (6.34), for any $0\leq\sigma\leq 1$, and the assertion holds by Proposition 6.13. ∎

Remark 6.19*.*

The convergence rates with respect to the $L_{2}(\mathcal{D})$-norms ($\sigma=0$)

[TABLE]

of the sinc-Galerkin FE approximation $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ and its covariance function $\widetilde{\varrho}_{h,k}^{\beta}$ reflect the higher regularity of the Whittle–Matérn field $\mathcal{Z}^{\beta}$ in (4.1) for large $\beta>0$ in (4.2). In particular, when the integer part does not vanish, $n_{\beta}\in\mathbb{N}$, a polynomial degree $p>1$ is meaningful, since thus higher order convergence rates can be achieved, cf. the numerical experiments in Section 7.

6.3.2. Less regularity

We now discuss convergence of FE discretizations when the operator $L$ in (2.2) has a coefficient $A$ which is not necessarily Lipschitz continuous or the domain $\mathcal{D}$ is not convex, i.e., the general case that $L$ is only $H^{1+\alpha}(\mathcal{D})$-regular. In the following definition we specify what we mean by this.

Definition 6.20**.**

Suppose Assumptions 2.1.I–II, 2.3.I, let $0<\alpha\leq 1$ and $L$ be the second-order differential operator in (2.2). We say that the elliptic problem associated with $L$ is $H^{1+\alpha}(\mathcal{D})$-regular if the restriction of $L\colon H^{1}_{0}(\mathcal{D})\to H^{1}_{0}(\mathcal{D})^{*}$ to $H^{1}_{0}(\mathcal{D})\cap H^{1+\alpha}(\mathcal{D})$ is a continuous map to $\dot{H}^{-1+\alpha}_{L}=(\dot{H}^{1-\alpha}_{L})^{*}$, see (2.4), and if additionally the data-to-solution map $L^{-1}\colon f\mapsto L^{-1}f$ is a bounded linear operator from $\dot{H}^{-1+\alpha}_{L}$ to $H^{1}_{0}(\mathcal{D})\cap H^{1+\alpha}(\mathcal{D})$.

We quote the following extension of the equivalence in (2.6) of Lemma 2.4 from [5, Prop. 4.1] to values $1\leq\sigma\leq 1+\alpha$ which holds if the elliptic problem associated with $L$ be $H^{1+\alpha}(\mathcal{D})$-regular.

Lemma 6.21**.**

Let the elliptic problem associated with $L$ be $H^{1+\alpha}(\mathcal{D})$-regular, see Definition 6.20. Then the equivalence in (6.7) holds for this parameter $0<\alpha\leq 1$.

Lemma 6.22**.**

Suppose Assumptions 6.15, 2.3.III (i.e., $\mathcal{D}$ is a Lipschitz polytope) and let $p=1$. Then, Assumption 6.1.I is satisfied for $\theta_{0}=\nicefrac{{1}}{{2}}$ and $\theta_{1}=2$.

Proof.

The operator $\mathcal{I}_{h}\colon H^{\theta}(\mathcal{D})\to V_{h}$ in Assumption 6.1.I can be taken as the Scott–Zhang interpolant, see [15, Lem. 1.130]. ∎

Theorem 6.23**.**

In addition to Assumptions 6.15, 2.3.III, suppose that the elliptic problem associated with $L$ is $H^{1+\alpha}(\mathcal{D})$-regular for some $0<\alpha\leq 1$ (see Definition 6.20) and let $p=1$. Assume further that $d\in\{1,2,3\}$, $\beta>0$ and $0\leq\sigma\leq 1$ are such that $2\beta-\sigma>\nicefrac{{d}}{{2}}$. Let $\mathcal{Z}^{\beta}$ be the Whittle–Matérn field in (4.1) and, for $h,k>0$, let $\mathcal{Z}^{\beta}_{h,k}$ be the sinc-Galerkin approximation in (6.16), with covariance functions $\varrho^{\beta}$ and $\varrho_{h,k}^{\beta}$. Then, for every $q,\varepsilon>0$ and sufficiently small $h>0$, $k=k(h)>0$,

[TABLE]

where, if $d=3$, for (6.24) to hold, we also suppose that $\beta>1$ and $\alpha\geq\nicefrac{{1}}{{2}}-\sigma$.

In addition, if $d=1$ and $0<\gamma\leq\nicefrac{{1}}{{2}}$ is such that $2\beta>\gamma+\nicefrac{{1}}{{2}}$, then

[TABLE]

for sufficiently small $h>0$, $k=k(h)>0$, every $\delta\in(0,\gamma)$ and $\varepsilon,q>0$.

Proof.

By Lemma 6.21 the equivalence in (6.7) holds. Furthermore, by Lemma 6.22 Assumption 6.1.I is satisfied for $\theta_{0}=\nicefrac{{1}}{{2}}<1$ and $\theta_{1}=2\geq 1+\alpha$. Finally, since we assume that the family of triangulations $(\mathcal{T}_{h})_{h>0}$ of $\overline{\mathcal{D}}\subset\mathbb{R}^{d}$ is quasi-uniform, the $L_{2}(\mathcal{D})$-orthogonal projection $\Pi_{h}$ is $H^{1}(\mathcal{D})$-stable, see [10] for $d\in\{1,2\}$ and [6] for arbitrary $d\in\mathbb{N}$. Thus, Propositions 6.12 and 6.14 are applicable and yield the assertions of this theorem. ∎

7. Numerical experiments

In the following numerical experiments we consider the original Whittle–Matérn field from (1.1) in Subsection 1.1, i.e., $L:=-\Delta+\kappa^{2}$, on the unit interval $\mathcal{D}=(0,1)$, augmented with homogeneous Dirichlet boundary conditions. We choose $\kappa:=0.5$ and apply a finite element discretization with continuous, piecewise polynomial basis functions of degree at most $p\in\{1,2\}$ to compute the sinc-Galerkin approximation $\mathcal{Z}^{\beta}_{h,k}$ (or $\widetilde{\mathcal{Z}}^{\beta}_{h,k}$) in (6.16)/(6.17). More precisely, we investigate

1. (i)

the empirical convergence to the Whittle–Matérn field $\mathcal{Z}^{\beta}$, see (4.1), with respect to the norms on $L_{2}(\Omega;L_{2}(\mathcal{D}))$, $L_{1}(\Omega;L_{\infty}(\mathcal{D}))$, and $L_{2}(\Omega;H^{1}_{0}(\mathcal{D}))$ for $\beta\in\{0.5,0.8,1.1,1.4,1.7\}$; 2. (ii)

the empirical convergence of the covariance function with respect to the norms on $L_{2}(\mathcal{D}\times\mathcal{D})$ and $L_{\infty}(\mathcal{D}\times\mathcal{D})$ for $\beta\in\{0.5,0.6,0.7,0.8,0.9,1\}$.

To this end, we generate an equidistant initial mesh on $\overline{\mathcal{D}}=[0,1]$ with $N_{0}:=9$ nodes (resp. $N_{0}:=17$ for the $L_{\infty}$-studies), of mesh size $h_{0}:=2^{-3}$ (resp. $h_{0}:=2^{-4}$). This initial mesh is $4$ times uniformly refined, so that on level $\ell\in\{0,\ldots,4\}$ the mesh is of width $h_{\ell}=h_{0}2^{-\ell}$. For $p\in\{1,2\}$, we use the MATLAB-based package ppfem [1] to assemble the matrices $\mathbf{M}$ and $\mathbf{L}$ in (6.20) and (6.21) with respect to the Babuška–Shen nodal basis $\{\phi_{j,h}\}_{j=1}^{N_{h}}$. On level $\ell$, the step size $k=k_{\ell}>0$ of the sinc quadrature is calibrated with the finite element mesh width via $k_{\ell}=-1/(\beta\ln h_{\ell})$.

The reference solutions for the field and the covariance function are generated based on an overkill Karhunen–Loève expansion of $\mathcal{Z}^{\beta}$ with $N_{\mathrm{KL}}=1000$ terms,

[TABLE]

where $\lambda_{j}=j^{2}\pi^{2}+\kappa^{2}$ and $e_{j}(x)=\sqrt{2}\sin(j\pi x)$ are the eigenvalues and eigenfunctions of $L=-\Delta+\kappa^{2}$ on $\mathcal{D}=(0,1)$. Here, for each of $100$ Monte Carlo runs, the same realization of the set of random variables $\{\xi_{1},\ldots,\xi_{N_{\mathrm{KL}}}\}$ is used to generate $\mathcal{Z}_{\mathrm{ref}}^{\beta}$ and the load vector $\mathbf{b}\sim\mathcal{N}(\mathbf{0},\mathbf{M})$ via

[TABLE]

For $d=1$, the operator $L$ does not have multiple eigenvalues and we can assemble the matrix $\mathbf{R}$, for each $h\in\{h_{0},\ldots,h_{4}\}$, by computing the discrete eigenfunctions $\{e_{j,h}\}_{j=1}^{N_{h}}$ and by adjusting their sign so that $e_{j,h}$ indeed approximates $e_{j}$ for each $j\in\{1,\ldots,N_{h}\}$. Note that we only have to assemble this matrix $\mathbf{R}$ to have comparable samples of the sinc-Galerkin approximation and the reference solution needed for the strong error studies. For the simulation practice, one could compute the Cholesky factor of the Gramian $\mathbf{M}$ or approximate the matrix square root $\sqrt{\mathbf{M}}$, e.g., as proposed in [23], in order to sample from $\mathbf{b}$. Since furthermore the dimension of the finite element spaces, even at the highest level $\ell=4$, is relatively small, we can assemble the covariance matrices of the sinc-Galerkin approximation directly, without Monte Carlo sampling, as

[TABLE]

cf. (6.20)–(6.21).

Note that the operator $L:=-\Delta+0.25$ has constant (and, thus, smooth) coefficients. Therefore, Theorem 6.18 provides (essentially) optimal convergence rates for the error of $\widetilde{\mathcal{Z}}_{h,k}^{\beta}$ in $L_{2}(\Omega;L_{2}(\mathcal{D}))$, $L_{2}(\Omega;H^{1}_{0}(\mathcal{D}))$ and of $\widetilde{\varrho}_{h,k}^{\beta}$ in $L_{2}(\mathcal{D}\times\mathcal{D})$. Furthermore, the convergence results of Theorem 6.23 on the $L_{1}(\Omega;L_{\infty}(\mathcal{D}))$-error are (essentially) sharp if $\beta\in(\nicefrac{{1}}{{4}},1)$ (resp. if $\beta\in(\nicefrac{{1}}{{4}},\nicefrac{{5}}{{8}})$ for the $L_{\infty}$-error of the covariance). For this smooth case, we have $\alpha>p+1$ in (6.7). For this reason, we expect the convergence rates listed in Table 1. The expected rates corresponding to the values of $\beta>\nicefrac{{1}}{{4}}$ used in our experiments are shown in parentheses in Table 2.

For every of the $100$ Monte Carlo samples, we approximate the integrals needed for computing the $L_{2}(\mathcal{D})$ and $H^{1}_{0}(\mathcal{D})$-errors by using MATLAB’s built-in function integral with tolerance 1e-6. For the $L_{\infty}$-studies we consider the largest error with respect to an equidistant mesh on on $\overline{\mathcal{D}}=[0,1]$ with $N_{\operatorname{ok}}=1001$ nodes, i.e.,

[TABLE]

where $x_{j}:=(j-1)10^{-3}$. Furthermore, to compute the $L_{2}(\mathcal{D}\times\mathcal{D})$-error, we approximate the distance of the covariances by a function which is piecewise constant on a regular lattice with $N_{\operatorname{ok}}^{2}$ nodes. Finally, the empirical convergence rates, also shown in Table 2, are obtained via a least-squares affine fit with respect to the data set $\{(\ln h_{\ell},\ln\operatorname{err}_{\ell}):2\leq\ell\leq 4\}$. Here, $\operatorname{err}_{\ell}$ denotes the error on level $\ell$ with respect to the norm used in the study and for the respective value of $\beta$ and $p$.

The resulting observed errors are displayed in Figure 1 for the fields and in Figure 2 for the covariances. Overall, the empirical results validate our theoretical outcomes fairly well, with a slight deviation for the $L_{\infty}$-studies which may be caused by a larger pre-asymptotic range.

8. Conclusion and discussion

We have identified necessary and sufficient conditions for square-integrability, Sobolev regularity, and Hölder continuity (in $L_{q}(\Omega)$-sense) for GRFs in terms of their color, as well as square-integrability, mixed Sobolev regularity, and continuity of their covariance functions, see Propositions 3.4, 3.6 and 3.7. Subsequently, we have applied these findings to generalized Whittle–Matérn fields, see $\mathcal{Z}^{\beta}$ in (4.1), where these conditions become assumptions on the smoothness parameter $\beta>0$, corresponding to the fractional exponent of the color $L^{-\beta}$, see Lemmata 4.1–4.3.

While these regularity results readily implied convergence of spectral Galerkin approximations, see Corollaries 5.1–5.3, significantly more work was needed to derive convergence for general Galerkin (such as finite element) approximations, for the following reason: It was unknown, how the deterministic fractional Galerkin error $L^{-\beta}g-L_{h}^{-\beta}g$ behaves in the Sobolev space $H^{\sigma}(\mathcal{D})$, for $0\leq\sigma\leq 1$, all possible exponents $\beta>0$, and sources $g\in H^{\delta}(\mathcal{D})$ of possibly negative regularity $\delta<0$. We have identified this behavior in Theorem 6.6 for the general situation that the second-order elliptic differential operator $L$ is $H^{1+\alpha}(\mathcal{D})$-regular for some $0<\alpha\leq 1$. This result could be exploited to show convergence of the sinc-Galerkin approximations and their covariances to the Whittle–Matérn field $\mathcal{Z}^{\beta}$ and to its covariance function $\varrho^{\beta}$, respectively, see Theorems 6.18 and 6.23.

The fact that the Rayleigh–Ritz projection and, thus, the deterministic Galerkin error $L^{-1}g-L_{h}^{-1}g$ converges at the rate $\min\{1+\alpha-\sigma,\,2\alpha\}$ in $H^{\sigma}(\mathcal{D})$, $0\leq\sigma\leq 1$, if $L$ is $H^{1+\alpha}(\mathcal{D})$-regular, cf. Lemma 6.3, and at the rate $p+1-\sigma$ if the problem is “smooth” and a conforming finite element discretization with piecewise polynomial basis functions of degree at most $p\in\mathbb{N}$ is used, combined with the low regularity of white noise in $\dot{H}^{-\nicefrac{{d}}{{2}}-\varepsilon}_{L}$, show that the Sobolev convergence rates of Theorems 6.18 and 6.23 are (essentially, up to $\varepsilon>0$) optimal. In addition, we believe that our results on Hölder convergence of the field and on $L_{\infty}$-convergence of the covariance function for $d=1$ in Theorem 6.23 are optimal

(i) if the problem is only $H^{1+\alpha}(\mathcal{D})$-regular for $\alpha\in(0,1)$ maximal, or

(ii) if the problem is smooth and $\beta\in(\nicefrac{{1}}{{4}},1)$ (resp. $\beta\in(\nicefrac{{1}}{{4}},\nicefrac{{5}}{{8}})$ for the covariance).

However, the deterministic $p$-FEM $L_{\infty}$-rate for $d=1$ is known to be $p+1$ if the problem is smooth, see [14]. Thus, our results will not be sharp in this case, see also our numerical experiments in Section 7.

Since the approach on deriving optimal $L_{\infty}$-rates involves non-Hilbertian regularity of the solution in $W^{p+1,\infty}(\mathcal{D})$, such a discussion was beyond the scope of this article and we leave this problem as well as the $C^{\delta}(\overline{\mathcal{D}})\,$/$\,L_{\infty}(\mathcal{D}\times\mathcal{D})$ error analysis of sinc-Galerkin approximations in dimension $d\in\{2,3\}$ as topics for future research.

Appendix A Proof of Proposition 6.13

The following lemma will be the main tool for the derivation of Proposition 6.13.

Lemma A.1**.**

Suppose Assumptions 2.1.I–II and 6.1.III. Let Assumption 6.1.IV be fulfilled with parameters $r,s_{0},t>0$ such that $\nicefrac{{r}}{{2}}\geq t-1$ and $s_{0}\geq t$. Let $d\in\mathbb{N}$, $\beta>0$, and $\rho_{0}(\,\cdot)$, $\rho_{1}(\,\cdot\,)$ be as in (6.33), i.e.,

[TABLE]

and define the exception set

[TABLE]

Then, for $\sigma\in\{0,1\}$, the Galerkin error $\widetilde{E}_{V_{h}}^{\beta}$ in (6.22) satisfies

[TABLE]

for sufficiently small $h>0$. Here, $\{(\lambda_{j},e_{j})\}_{j\in\mathbb{N}}$ are the $L_{2}(\mathcal{D})$-orthonormal, ordered eigenpairs of $L$ in (2.2) and we set $C_{\tau,h}:=1$ if $\tau\not\in\mathcal{E}_{\tau}$ and $C_{\tau,h}:=\ln(1/h)$ if $\tau\in\mathcal{E}_{\tau}$.

Proof.

Fix $\tau\geq 0$. The definitions of $\widetilde{E}_{V_{h}}^{\beta}$ in (6.22) and of $\widetilde{\Pi}_{h}$ in (6.18) yield

[TABLE]

By the mean value theorem, $\lambda_{j}^{-\beta}-\lambda_{j,h}^{-\beta}=\widetilde{\lambda}_{j}^{-\beta-1}(\lambda_{j,h}-\lambda_{j})$ for some $\widetilde{\lambda}_{j}\in(\lambda_{j},\lambda_{j,h})$. Thus, we can use (6.5) from Assumption 6.1.IV and the spectral behavior (2.3) from Lemma 2.2 combined with Assumption 6.1.III to bound the first sum in (A.2),

[TABLE]

where we also have used that $r\geq 2(t-1)$ by assumption. For the second sum in (A.2) we distinguish the cases $\sigma=0$ and $\sigma=1$. If $\sigma=0$, we can apply (6.6) of Assumption 6.1.IV and obtain

[TABLE]

since $s_{0}\geq t$. For $\sigma=1$, we first note that (6.5)–(6.6) of Assumption 6.1.IV imply the following estimate with respect to the norm on $\dot{H}^{1}_{L}$,

[TABLE]

Here, we have used the identity $(e_{j},e_{j,h})_{1}=\lambda_{j}(e_{j},e_{j,h})_{0}$. Thus, if $\sigma=1$, we can bound the second sum in (A.2) as follows,

[TABLE]

since $s_{0}\geq t$ and $\nicefrac{{r}}{{2}}\geq t-1$ by assumption. Combining (A.2), (A.3), (A.4) and (A.5) completes the proof. ∎

Proof of Proposition 6.13.

In order to derive (6.35), we start with splitting the error in the norm $\|\,\cdot\,\|_{\sigma}$ on $\dot{H}^{\sigma}_{L}$, cf. (2.4), which by (2.6) of Lemma 2.4 implies an upper bound for the Sobolev norm:

[TABLE]

Here, $\mathcal{Z}^{\beta}_{N_{h}}$ is the spectral Galerkin approximation from (5.1) and $\widetilde{\mathcal{Z}}^{\beta}_{h}$ denotes a GRF colored by $L_{h}^{-\beta}\widetilde{\Pi}_{h}$. We readily obtain a bound for ($\text{A}_{\mathcal{Z}}$) from (5.3) of Corollary 5.1, combined with Assumption 6.1.III. This gives

[TABLE]

Note that it suffices to estimate the terms ($\text{B}_{\mathcal{Z}}$) and ($\text{C}_{\mathcal{Z}}$) for $\sigma\in\{0,1\}$. The respective bounds for $\sigma\in(0,1)$ then follow by interpolation. By definition of the Galerkin and the quadrature error, $\widetilde{E}_{V_{h}}^{\beta}$, $\widetilde{E}_{Q}^{\beta}$, in (6.22)–(6.23) and by Proposition 3.7,

[TABLE]

Since we have to consider these terms only for $\sigma\in\{0,1\}$, the first term can be bounded by (A.1) of Lemma A.1 (with $\tau:=0$),

[TABLE]

where $C_{\beta,h}^{\mathcal{Z}}>0$ is defined as in the statement of Proposition 6.13. To estimate ($\text{C}_{\mathcal{Z}}$), we first apply the convergence result of the sinc quadrature from [5, Lem. 3.4, Rem. 3.1, Thm. 3.5]. Thus, for sufficiently small $k>0$ and all $1\leq j\leq N_{h}$,

[TABLE]

Again by equivalence of the norms $\|\,\cdot\,\|_{\sigma}$, $\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}$ for $\sigma\in\{0,1\}$, see Lemma 2.4, and by the inverse inequality (6.4) from Assumption 6.1.II, we then find

[TABLE]

where we have used the spectral behavior (2.3) from Lemma 2.2 and Assumptions 6.1.III-IV in the last step. This completes the proof of (6.35).

We now proceed with the derivation of (6.36). To this end, we consider the error with respect to the norm $\|\,\cdot\,\|_{\sigma,\sigma}$, see (3.15), since the embedding in (2.6) implies that $\dot{H}^{\sigma,\sigma}_{L}\hookrightarrow H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})$. We again partition the error in three terms,

[TABLE]

where $\widetilde{\varrho}^{\beta}_{h}$ denotes the covariance function of the above-introduced GRF $\widetilde{\mathcal{Z}}^{\beta}_{h}$ colored by $L_{h}^{-\beta}\widetilde{\Pi}_{h}$. A bound for the truncation error is given by (5.4) in Proposition 5.1,

[TABLE]

where we also used Assumption 6.1.III. We bound the remaining terms ($\text{B}_{\varrho}$) and ($\text{C}_{\varrho}$) for $\sigma\in\{0,1\}$. Since \bigl{[}\dot{H}^{0,0}_{L},\dot{H}^{1,1}_{L}\bigr{]}_{\sigma}=\dot{H}^{\sigma,\sigma}_{L}, see [41, Thm. 16.1], we may again interpolate these results for $\sigma\in(0,1)$. To this end, we first exploit (3.20) from Proposition 3.7 and (6.31) to derive for ($\text{B}_{\varrho}$) that

[TABLE]

By Lemma A.1 (for $\tau:=2\beta-\sigma>\nicefrac{{d}}{{2}}>0$ in (A.1)) we have, for $\sigma\in\{0,1\}$ and for $C_{\beta,h}^{\varrho}>0$ as in the statement of Proposition 6.13,

[TABLE]

Next, we use the identity \bigl{(}L_{h}^{-\beta}\widetilde{\Pi}_{h}\bigr{)}^{*}e_{j}=\sum_{\ell=1}^{N_{h}}\lambda_{\ell,h}^{-\beta}(e_{j},e_{\ell,h})_{L_{2}(\mathcal{D})}\,e_{\ell}, the orthogonality $(e_{k,h},e_{\ell,h})_{\sigma}=\delta_{k\ell}\lambda_{k,h}^{\sigma}$ (here, $\delta_{k\ell}$ denotes the Kronecker delta), which holds for $\sigma\in\{0,1\}$, and the relation $\lambda_{j}\leq\lambda_{j,h}$ from Assumption 6.1.IV. With these steps, we obtain again by (A.1) of Lemma A.1 (with $\tau:=2\beta-\sigma$) a bound for ($\text{B}^{\prime\prime}_{\varrho}$),

[TABLE]

In conclusion, \bigl{\|}\varrho_{N_{h}}^{\beta}-\widetilde{\varrho}_{h}^{\beta}\bigr{\|}_{\sigma,\sigma}\leq\text{(\text{B}^{\prime}{\varrho})}+\text{(\text{B}^{\prime\prime}{\varrho})}\lesssim_{(C_{0},C_{\lambda},\sigma,\beta,A,\kappa,\mathcal{D})}\bigl{(}C_{\beta,h}^{\varrho}\bigr{)}^{2}h^{\rho_{\sigma}(2\beta-\sigma)} for $\sigma\in\{0,1\}$. For ($\text{C}_{\varrho}$), we derive with the equivalence of the norms $\|\,\cdot\,\|_{\sigma}$, $\|\,\cdot\,\|_{H^{\sigma}(\mathcal{D})}$, the inverse inequality (6.4) from Assumption 6.1.II, and the convergence result for the sinc quadrature [5, Lem. 3.4, Rem. 3.1, Thm. 3.5] the following, if $\sigma\in\{0,1\}$,

[TABLE]

Since $L_{h}^{-n_{\beta}}\widetilde{\Pi}_{h}\widetilde{\Pi}_{h}^{*}(L_{h}^{-n_{\beta}})^{*}e_{\ell,h}=\lambda_{\ell,h}^{-2n_{\beta}}e_{\ell,h}$ for all $\ell\in\{1,\ldots,N_{h}\}$, this shows that

[TABLE]

Next, again by the inverse inequality (6.4) we find

[TABLE]

Here, we have used the uniform stability of $L^{-\beta_{\star}}$, $Q_{h,k}^{\beta_{\star}}$ with respect to $h$ and $k$, see (6.32), as well as (2.3) from Lemma 2.2 and Assumption 6.1.I. Combining the bounds for ($\text{A}_{\varrho}$), ($\text{B}_{\varrho}$), ($\text{C}_{\varrho}$) completes the proof. ∎