Numerical solution of fractional elliptic stochastic PDEs with spatial white noise

Abstract

The numerical approximation of solutions to stochastic partial differential equations with additive spatial white noise on bounded domains in |$\mathbb{R}^d$| is considered. The differential operator is given by the fractional power |$L^\beta $|⁠, |$\beta \in (0,1)$| of an integer-order elliptic differential operator |$L$| and is therefore nonlocal. Its inverse |$L^{-\beta }$| is represented by a Bochner integral from the Dunford–Taylor functional calculus. By applying a quadrature formula to this integral representation the inverse fractional-order operator |$L^{-\beta }$| is approximated by a weighted sum of nonfractional resolvents |$( I + \exp(2 y_\ell) \, L )^{-1}$| at certain quadrature nodes |$t_j> 0$|⁠. The resolvents are then discretized in space by a standard finite element method. This approach is combined with an approximation of the white noise, which is based only on the mass matrix of the finite element discretization. In this way an efficient numerical algorithm for computing samples of the approximate solution is obtained. For the resulting approximation the strong mean-square error is analyzed and an explicit rate of convergence is derived. Numerical experiments for |$L=\kappa ^2-\Delta $|⁠, |$\kappa> 0$| with homogeneous Dirichlet boundary conditions on the unit cube |$(0,1)^d$| in |$d=1,2,3$| spatial dimensions for varying |$\beta \in (0,1)$| attest to the theoretical results.

stochastic partial differential equations, Gaussian white noise, fractional operators, finite element methods, Matérn covariances, spatial statistics

1. Introduction

A real-valued Gaussian random field |$u$| defined on a spatial domain |$\mathcal{D} \subset \mathbb{R}^d$| is called a Gaussian Matérn field if its covariance function |$C: \mathcal{D}\times \mathcal{D} \to \mathbb{R}$| is given by

$$\begin{align} C(\mathbf{x}_1, \mathbf{x}_2) = \frac{2^{1-\nu}\sigma^2}{\varGamma(\nu)} (\kappa \| \mathbf{x}_1 - \mathbf{x}_2\| )^{\nu} K_{\nu}(\kappa\|\mathbf{x}_1 - \mathbf{x}_2\|), \quad \mathbf{x}_1, \mathbf{x}_2 \in \mathcal{D}, \end{align}$$

(1.1)

where |$\|\cdot \|$| is the Euclidean norm on |$\mathbb{R}^d$|⁠, and |$\varGamma $|⁠, |$K_{\nu }$| denote the gamma function and the modified Bessel function of the second kind, respectively. Via the positive parameters |$\sigma $|⁠, |$\nu $| and |$\kappa $| the most important characteristics of the random field |$u$| can be controlled: its variance, smoothness and correlation range. Due to this flexibility, Gaussian Matérn fields are often used for modeling in spatial statistics (see, e.g., Stein, 1999). However, a major drawback of this traditional covariance-based representation of Matérn fields is its high computational effort. For instance, sampling from |$u$| at |$n$| locations |$\mathbf{x}_1,\ldots ,\mathbf{x}_n \in \mathcal{D}$| requires a matrix factorization of an |$n \times n$| covariance matrix and, thus, in general, |$\mathcal{O}(n^3)$| arithmetic operations.

There are several approaches attempting to cope with this computational problem (see, e.g., Furrer et al., 2006; Banerjee et al., 2008; Sun et al., 2012; Nychka et al., 2015). One of them is based on the insight that a Gaussian Matérn field on |$\mathcal{D}:= \mathbb{R}^d$| with parameters |$\sigma ,\nu $| and |$\kappa>0$| can be interpreted as the statistically stationary solution |$u$| to the stochastic partial differential equation (SPDE):

$$\begin{equation} ( \kappa^2 - \varDelta )^\beta \, u(\mathbf{x}) = \mathcal{W}(\mathbf{x}), \qquad \mathbf{x}\in\mathcal{D}. \end{equation}$$

(1.2)

Here |$\varDelta $| denotes the Laplacian, |$4\beta = 2\nu + d$|⁠, and |$\mathcal{W}$| is Gaussian white noise on |$\mathbb{R}^d$|⁠. The marginal variance of |$u$| is then given by |$\sigma ^2 = \varGamma (2\beta - d/2)\varGamma (2\beta )^{-1}(4\pi )^{-d/2}\kappa ^{d - 4\beta }$|⁠.

This relation between Matérn fields and SPDEs had already been noticed by Whittle (1954, 1963). Later on, based on a finite element discretization of (1.2), where the differential operator |$\kappa ^2-\varDelta $| is augmented with Neumann boundary conditions, Markov random field approximations of Matérn fields on bounded domains |$\mathcal{D}\subsetneq \mathbb{R}^d$| were introduced by Lindgren et al. (2011). Owing to the computational savings of this approach compared to the covariance-based approximations, these SPDE-based models have become very popular in spatial statistics (see, e.g., Bhatt et al., 2015; Lai et al., 2015), and they are still the subject of current research, mainly for the following reason: the SPDE formulation (1.2) facilitates various generalizations of approximations of Gaussian Matérn fields involving (a) other differential operators (Bolin & Lindgren, 2011; Lindgren et al., 2011; Fuglstad et al., 2015), (b) more general domains, such as the sphere (Bolin & Lindgren, 2011; Lindgren et al., 2011) and (c) non-Gaussian driving noise (Bolin, 2014; Wallin & Bolin, 2015). However, a considerable drawback of the finite element approximation proposed by Lindgren et al. (2011) is that it is computable only if |$2\beta \in \mathbb{N}$|⁠. This limits flexibility and, in particular, it implies that the method cannot be applied to the important special case of exponential covariance (⁠|$\nu = 1/2$|⁠) in |$\mathbb{R}^2$|⁠, which corresponds to |$\beta = 3/4$|⁠.

With the objective of extending the approach of Lindgren et al. (2011) to all admissible values |$\nu> 0$| (i.e., |$\beta> d/4)$| and the generalizations mentioned above we consider (1.2) in the more general framework of fractional-order elliptic equations driven by white noise. We propose an explicit numerical scheme for generating samples of an approximation to the Gaussian solution process, which allows for any fractional power |$\beta>d/4$|⁠. This method is based on (i) a standard finite element discretization in space, (ii) a quadrature approximation of the inverse fractional elliptic differential operator proposed by Bonito & Pasciak (2015) for deterministic equations and (iii) an approximation of the noise term on the right-hand side, whose covariance matrix after discretization is equal to the finite element mass matrix.

Due to (iii) no explicit knowledge of the eigenfunctions of the differential operator is needed, in contrast to approximations based on truncated Karhunen–Loève expansions of the noise term (e.g., Zhang et al., 2016).

While Zhang et al. (2016) remarked that any orthonormal basis can be used in the Karhunen–Loève expansion of the noise term before projecting it onto the finite element space, their error analysis strongly benefits from the fact that the eigenfunctions of the differential operator are used. In fact, if a general orthonormal basis was used, it would be difficult to obtain explicit rates of convergence with respect to the truncation level after truncating the expansion at a finite level. Only if the constructed basis has certain smoothness with respect to the differential operator is it possible to obtain explicit rates of convergence (see, e.g., Kovács et al., 2011, addressing this kind of problem). Constructing such an orthonormal basis requires a lot of computational effort, in particular, for complicated domains and |$d=2,3$|⁠. In contrast the present approach is based on an expansion with respect to the standard (nonorthonormal) finite element basis for the discretized problem, which is readily available even for complex geometries.

This work follows a series of investigations of numerical methods for deterministic fractional-order equations (Gavrilyuk, 1996; Gavrilyuk et al., 2004, 2005; Roop, 2006; Caffarelli & Silvestre, 2007; Baeumer et al., 2015; Bonito & Pasciak, 2015; Jin et al., 2015; Nochetto et al., 2015; Bonito et al., 2017) and for nonfractional elliptic SPDEs with random forcing (e.g., Du & Zhang, 2003; Babuska et al., 2004; Gyöngy & Martínez, 2006; Cao et al., 2007; Zhang et al., 2016). An essentially similar idea to our approach, which combines (iii) with Lanczos and Krylov subspace methods, was pursued by Simpson (2008) for the differential operator |$\kappa ^2-\varDelta $| in (1.2). However, neither weak nor strong convergence have been proven and the empirical results show generally poor performance of the proposed scheme. A weak error estimate for |$\kappa ^2 - \varDelta $| and the nonfractional case |$\beta =1$| in |$d=2$| spatial dimensions was derived by Simpson et al. (2012), showing quadratic weak convergence for that particular case. Since the first two moments uniquely determine the distribution of the Gaussian solution process, an alternative approach is to consider the problem of approximating the covariance function of the solution instead of solving for the solution process itself (Dölz et al., 2017). Note that this method cannot be generalized to non-Gaussian models.

The structure of this article is as follows: in Section 2 we present the fractional-order equation with Gaussian white noise in a Hilbert space setting, together with the necessary assumptions, and comment on existence and uniqueness of a solution to this problem. Furthermore, we introduce the numerical approximation of the solution process and state our main result in Theorem 2.10: strong mean-square convergence of this approximation at an explicit rate. This theorem is proven in Section 3 by partitioning the strong mean-square error into several terms and estimating these terms one by one. In addition a weak-type convergence result is obtained in Corollary 3.4. In Section 4 the SPDE (1.2) is considered for numerical experiments on the unit cube |$\mathcal{D} = (0,1)^d$| with continuous, piecewise linear finite element basis functions in |$d=1,2,3$| spatial dimensions. The outcomes of the paper are summarized in Section 5.

2. Model problem and main result

In the following let |$H$| denote a separable Hilbert space and |$L: \mathcal{D}(L) \subset H \to H$| be a densely defined, self-adjoint, positive definite linear operator with a compact inverse. In this case there exists an orthonormal basis |$\{e_j\}_{j\in \mathbb{N}}$| of |$H$| consisting of eigenvectors of |$L$|⁠. The eigenvalue–eigenvector pairs |$\{(\lambda _j, e_j)\}_{j\in \mathbb{N}}$| can be arranged such that the positive eigenvalues |$\{\lambda _j\}_{j\in \mathbb{N}}$| are in nondecreasing order:

$$\begin{align*} 0 < \lambda_1 \leqslant \lambda_2 \leqslant \cdots \leqslant \lambda_j \leqslant \lambda_{j+1}\leqslant \cdots, \qquad \lim_{j\to\infty} \lambda_j = \infty. \end{align*}$$

We assume that the growth of the eigenvalues is given by |$\lambda _j \propto j^\alpha $| for an exponent |$\alpha> 0$|⁠, i.e., there exist constants |$c_\lambda , C_\lambda> 0$| such that

$$\begin{align} c_\lambda \, j^{\, \alpha} \leqslant \lambda_j \leqslant C_\lambda \, j^{\, \alpha} \qquad \forall \,j\in\mathbb{N}. \end{align}$$

(2.1)

For |$\beta> 0$| and |$\phi \in \mathcal{D}( L^{\beta } ) := \{ \psi \in H: \sum _{j\in \mathbb{N}} \lambda _j^{2\beta } (\psi , e_j)_{H}^2 < \infty \}$| the action of the fractional power operator |$L^{\beta }: \mathcal{D}( L^{\beta }) \to H$| is defined by

$$\begin{align} L^{\beta} \phi:= \sum_{j\in\mathbb{N}} \lambda_j^{\beta} (\phi, e_j)_{H} \, e_j. \end{align}$$

(2.2)

The subspace |$\dot{H}^{2\beta }:= \mathcal{D}(L^{\beta }) \subset H$| is itself a Hilbert space with respect to the inner product and corresponding norm given by

$$\begin{align*} (\phi, \psi)_{2\beta}:= \bigl( L^{\beta} \phi, L^{\beta} \psi \bigr)_{H}, \qquad \|\phi\|_{2\beta}^2 := \bigl\| L^{\beta} \phi \bigr\|_{H}^2 = \sum_{j\in\mathbb{N}} \lambda_j^{2\beta} (\phi, e_j)_{H}^2. \end{align*}$$

Its dual space after identification via the inner product on |$H$|⁠, which is continuously extended as a duality pairing, is denoted by |$\dot{H}^{-2\beta }$|⁠. For |$s \geqslant 0$| the norm on the dual space |$\dot{H}^{-s}$| enjoys the useful representation

$$\begin{align} \|g\|_{-s} = \sup\limits_{\phi\in\dot{H}^{s}\setminus\{0\}} \frac{\langle g, \phi\rangle}{\|\phi\|_{s}} = \left( \sum_{j\in\mathbb{N}} \lambda_j^{-s} \langle g, e_j\rangle^2 \right)^{1/2}, \end{align}$$

(2.3)

where |$\langle \,\cdot \,, \,\cdot \,\rangle $| denotes the duality pairing between |$\dot{H}^{-s}$| and |$\dot{H}^{s}$| (Thomée, 2006, Proof of Lemma 5.1). We obtain the following scale of densely and continuously embedded Hilbert spaces:

$$\begin{align} \dot{H}^{s} \hookrightarrow \dot{H}^{r} \hookrightarrow \dot{H}^{0} := H \cong \dot{H}^{-0} \hookrightarrow \dot{H}^{-r} \hookrightarrow \dot{H}^{-s}, \qquad 0\leqslant r \leqslant s. \end{align}$$

(2.4)

It is an immediate consequence of these definitions that |$L^\beta $| is an isometric isomorphism from |$\dot{H}^{s}$| to |$\dot{H}^{s-2\beta }$| for |$s\geqslant 2\beta $|⁠, since for |$\phi \in \dot{H}^{s}$| we have

$$\begin{align} \bigl\| L^\beta \phi \bigr\|_{s-2\beta}^2 &= \left\| \sum_{j\in\mathbb{N}} \lambda_j^\beta (\phi, e_j)_{H} \, e_j \right\|_{s-2\beta}^2 = \sum_{j\in\mathbb{N}} \lambda_j^{s-2\beta} \lambda_j^{2\beta} (\phi, e_j)_{H}^2 = \|\phi\|_{s}^2. \end{align}$$

(2.5)

The following lemma states that |$L^\beta $| can be extended to a bounded linear operator between |$\dot{H}^{s}$| and |$\dot{H}^{s-2\beta }$| for all |$s\in \mathbb{R}$|⁠.

 

Lemma 2.1

For |$s\in \mathbb{R}$| there exists a unique continuous extension of |$L^\beta $| defined in (2.2) to an isometric isomorphism |$L^\beta : \dot{H}^{s} \to \dot{H}^{s-2\beta }$|⁠.

 

Proof.

For |$s \geqslant 2\beta $| the isometry property, and thus injectivity, has already been observed in (2.5). Surjectivity readily follows, since for any |$g \in \dot{H}^{s-2\beta }$| the vector |$\phi :=\sum _{j\in \mathbb{N}} \lambda _j^{-\beta } (g, e_j)_{H} \, e_j$| is an element of |$\dot{H}^{s}$| with |$\|\phi \|_{s}=\|g\|_{s-2\beta }$| and |$L^\beta \phi = g$|⁠.

Assume now that |$s < 2\beta $|⁠. We obtain for |$\phi \in \dot{H}^{2\beta }$| and |$\psi \in \dot{H}^{2\beta -s}$|⁠,

$$\begin{align*} \bigl\langle L^\beta \phi, \psi \bigr\rangle = \bigl( L^\beta \phi, \psi \bigr)_{H} = \sum_{j\in\mathbb{N}} \lambda_j^{\beta} (\phi, e_j)_{H} (\psi, e_j)_{H} \leqslant \|\phi|_{s} \|\psi\|_{2\beta - s}. \end{align*}$$

By density of |$\dot{H}^{2\beta }\hookrightarrow \dot{H}^{s}$| there exists a unique linear continuous extension |$L^\beta : \dot{H}^{s} \to \dot{H}^{s-2\beta }$|⁠. Moreover, it is an isometry, since the above estimate attains equality for |$\phi \in \dot{H}^{s}$| and |$\psi = L^{s-\beta } \phi \in \dot{H}^{2\beta -s}$|⁠. Surjectivity follows, similarly to the case |$s\geqslant 2\beta $| from the representation of the dual norm in (2.3) applied to |$\dot{H}^{s-2\beta }$|⁠.

 

Example 2.2

For |$\kappa \geqslant 0$| and a bounded, convex, polygonal domain |$\mathcal{D}\subset \mathbb{R}^d$| consider the eigenvalue value problem

$$\begin{align*} \left( \kappa^2 - \varDelta \right) e &= \lambda e && \textrm{in} \ \mathcal{D}, \\ e &= 0 && \textrm{on} \ \partial \mathcal{D}, \end{align*}$$

i.e., the operator |$L = \kappa ^2 - \varDelta $| with homogeneous Dirichlet boundary conditions on |$H=L_2(\mathcal{D})$|⁠. For this case we have |$\dot{H}^{2}=\mathcal{D}(L) = H^2(\mathcal{D}) \cap H^1_0(\mathcal{D})$|⁠, |$\dot{H}^{1} = H^1_0(\mathcal{D})$|⁠, as well as |$\dot{H}^{-1} = H^1_0(\mathcal{D})^* = H^{-1}(\mathcal{D})$|⁠, where|$\,\,^*$| denotes the dual after identification via the inner product on |$L_2(\mathcal{D})$|⁠. For |$s\in (0,1)$| one obtains the intermediate spaces

$$\begin{align*} \dot{H}^{s} &= H^s_0(\mathcal{D}), & H^s_0(\mathcal{D}) &:= \big[ L_2(\mathcal{D}), H^1_0(\mathcal{D}) \big]_{s,2}, \\ \dot{H}^{-s} &= H^{-s}(\mathcal{D}), & H^{-s}(\mathcal{D}) &:= \big[ H^{-1}(\mathcal{D}), L_2(\mathcal{D}) \big]_{1-s,2} = H^s_0(\mathcal{D})^*, \end{align*}$$

where |$[\,\cdot \,, \,\cdot \,]_{s,q}$| denotes the real |$K$|-interpolation method (see, e.g., Thomée, 2006, Chapter 19). Furthermore, one can show (Bonito & Pasciak, 2015, Proposition 4.1) that for |$1\leqslant s \leqslant 2$|⁠,

$$\begin{align*} \dot{H}^{s} = \left[ \dot{H}^{1}, \dot{H}^{2} \right]_{s-1,2} = \left[ H^1_0(\mathcal{D}), H^2(\mathcal{D}) \cap H^1_0(\mathcal{D}) \right]_{s-1,2} = H^s(\mathcal{D}) \cap H^1_0(\mathcal{D}). \end{align*}$$

If |$\mathcal{D}=(0,1)^d$| is the |$d$|-dimensional unit cube, then it is well known (e.g., Courant & Hilbert, 1953, Chapter VI.4) that the above eigenvalue problem has the following eigenvectors:

$$\begin{align} e_{\boldsymbol{\theta}}(\mathbf{x}) = e_{(\theta_1,\ldots,\theta_d)}(x_1,\ldots,x_d) = \prod_{i=1}^{d} \left( \sqrt{2} \sin(\pi\theta_i x_i) \right), \end{align}$$

(2.6)

where |$\boldsymbol{\theta }=(\theta _1,\ldots ,\theta _d)\in \mathbb{N}^d$| is a |$d$|-dimensional multi-index. The corresponding eigenvalues are given by

$$\begin{align} \lambda_{\boldsymbol{\theta}} = \kappa^2 + \pi^2 |\boldsymbol{\theta}|^2 = \kappa^2 + \pi^2 \sum_{i=1}^{d} \theta_i^2. \end{align}$$

(2.7)

These eigenvalues satisfy (2.1) for |$\alpha = 2/d$|⁠. Note that only the values of |$c_\lambda $| and |$C_\lambda $| in (2.1) change when considering any other bounded domain |$\mathcal{D}\subset \mathbb{R}^d$| with smooth or polygonal boundary. The value |$\alpha =2/d$| is the same by the min–max principle.

2.1 Fractional-order equation

Motivated by (1.2) we consider for |$g\in H$| the fractional-order equation

$$\begin{align} L^\beta u = g + \mathcal{W}, \end{align}$$

(2.8)

where |$\mathcal{W}$| denotes Gaussian white noise defined on a complete probability space |$(\varOmega , \mathcal{A}, \mathbb{P})$| with values in the separable Hilbert space |$H$|⁠. We assume that |$\beta \in (0,1)$| and refer to Remark 3.6 for a discussion of the generalization to |$\beta> 0$|⁠. Equation (2.8), as well as all the following equalities involving noise terms, is understood to hold |$\mathbb{P}$|-almost surely (⁠|$\mathbb{P}$|-a.s.).

Note that the white noise |$\mathcal{W}$| can formally be represented by the Karhunen–Loève expansion with respect to the orthonormal eigenbasis |$\{e_j\}_{j\in \mathbb{N}} \subset H$| of |$L$|⁠,

$$\begin{align} \mathcal{W} = \sum_{j\in\mathbb{N}} \xi_j \, e_j, \end{align}$$

(2.9)

where |$\{\xi _j\}_{j\in \mathbb{N}}$| is a sequence of independent real-valued standard normally distributed random variables. As we will show in Proposition 2.3, this formally defined series converges in |$L_2(\varOmega ; \dot{H}^{-s})$| for any |$s>\tfrac{1}{\alpha }$|⁠, which implies that realizations of |$\mathcal{W}$| are elements of |$\dot{H}^{-\frac{1}{\alpha }-\epsilon }$| for any |$\epsilon>0$||$\mathbb{P}$|-a.s. Related to this representation we introduce for |$N\in \mathbb{N}$| the truncated white noise

$$\begin{align} \mathcal{W}_N:= \sum_{j=1}^N \xi_j \, e_j. \end{align}$$

(2.10)

The following proposition specifies mean-square regularity of the white noise with respect to the |$\dot{H}$|-spaces in (2.4).

 

Proposition 2.3

For all |$s\geqslant 0$| there exists a constant |$C> 0$| depending only on |$\alpha $|⁠, |$c_\lambda $| in (2.1), and |$s$|⁠, such that the truncated white noise |$\mathcal{W}_N$| in (2.10) satisfies

$$\begin{align*} \mathbb{E}\bigl[ \|\mathcal{W}_N\|_{-s}^2 \bigr] \leqslant C \begin{cases} 1 + N^{1-\alpha s}, & s \neq \alpha^{-1}, \\ 1 + \ln(N), & s = \alpha^{-1}. \end{cases} \end{align*}$$

Furthermore, |$\mathcal{W} \in L_2\bigl (\varOmega ; \dot{H}^{-\frac{1}{\alpha } - \epsilon }\bigr )$| holds for all |$\epsilon>0$| with

$$\begin{align*} \mathbb{E}\left[ \|\mathcal{W}\|_{-\frac{1}{\alpha} - \epsilon}^2 \right] \leqslant c_\lambda^{-\frac{1}{\alpha} - \epsilon} \left( 1 + \tfrac{1}{\epsilon\alpha} \right). \end{align*}$$

 

Proof.

The orthonormality of the vectors |$\{e_j\}$|⁠, along with the distribution of the random variables |$\{\xi _j\}$| in the expansion of |$\mathcal{W}_N$| and the representation (2.3) of the dual norm, yields |$\mathbb{E}\bigl [ \|\mathcal{W}_N\|_{-s}^2 \bigr ] = \sum _{j=1}^N \lambda _j^{-s}$|⁠. Thus, we conclude with (2.1) that

$$\begin{align*} \mathbb{E}\Bigl[ \|\mathcal{W}_N\|_{-s}^2 \Bigr] \leqslant c_\lambda^{-s} \sum_{j=1}^N j^{-\alpha s} \leqslant c_\lambda^{-s} \left(\!1 + \int_1^N x^{-\alpha s} \, \mathrm{d}\, x \right) = c_\lambda^{-s} \begin{cases} \frac{N^{1-\alpha s} - \alpha s}{1-\alpha s}, & s\neq \alpha^{-1}, \\ 1 + \ln(N), & s = \alpha^{-1}. \end{cases} \end{align*}$$

Choosing |$s=\alpha ^{-1} + \epsilon $| and taking the limit |$N\to \infty $| shows the |$L_2\bigl (\varOmega ; \dot{H}^{-\frac{1}{\alpha }-\epsilon }\bigr )$| regularity of |$\mathcal{W}$|⁠.

 

Remark 2.4

Proposition 2.3 implies that |$\mathcal{W} \in \dot{H}^{-\frac{1}{\alpha }-\epsilon }$| holds also |$\mathbb{P}$|-a.s. for any |$\epsilon>0$|⁠. Therefore, existence and uniqueness of a solution |$u$| to (2.8) with regularity |$u\in \dot{H}^{2\beta - \frac{1}{\alpha } - \epsilon }$| (⁠|$\mathbb{P}$|-a.s.) follow from Lemma 2.1. In addition the above results show that |$u\in L_2\bigl (\varOmega ; \dot{H}^{2\beta - \frac{1}{\alpha } - \epsilon }\bigr )$|⁠. In particular |$u\in L_2(\varOmega ; H)$| holds if |$2\alpha \beta> 1$|⁠.

 

Remark 2.5

The above results are in accordance with Zhang et al. (Lemma 4.2 and Theorem 4.3 2016), where the following semilinear elliptic stochastic boundary value problem is considered on |$\mathcal{D}:=(0,1)^d$| for |$d\in \{1,2,3\}$|⁠, |$g\in L_2(\mathcal{D})$| and |$\,f:\mathbb{R}\to \mathbb{R}$|⁠:

$$\begin{align*} -\varDelta u(\mathbf{x}) + f(u(\mathbf{x})) &= g(\mathbf{x}) + \mathcal{W\,}(\mathbf{x}), && \mathbf{x}\in \mathcal{D}, \\ u(\mathbf{x}) &= 0, && \mathbf{x}\in\partial \mathcal{D}, \end{align*}$$

and mean-square regularity in the Sobolev space |$H^s(\mathcal{D})$| is proven under appropriate assumptions on the nonlinearity |$f$| for (i) the spatial white noise |$\mathcal{W}$| and |$s < -d/2$| and (ii) the solution |$u$| and |$s < 2 - d/2$|⁠. Note that in the linear case when |$f(u) = \kappa ^2 u$|⁠, |$\kappa \geqslant 0$|⁠, this corresponds to problem (2.8) with |$\beta =1$|⁠, |$L=\kappa ^2-\varDelta $|⁠, and the exponent |$\alpha $| of the eigenvalue growth in (2.1) is given by |$\alpha = 2/d$|⁠; see Example 2.2.

2.2 Finite element approximation

In the following we introduce a numerical method based on a finite element discretization for solving the fractional-order equation (2.8) approximately. For this purpose we consider the family |$(V_h)_{h\in (0,1)}$| of subspaces of |$\dot{H}^{1}$| with finite dimensions |$N_h:= \dim (V_h) < \infty $|⁠. The Galerkin discretization of the operator |$L: \dot{H}^{1} \to \dot{H}^{-1}$| is denoted by |$L_h: V_h \to V_h$|⁠, i.e., |$(L_h\psi _h,\phi _h)_{H} = \langle L\psi _h,\phi _h\rangle $| for all |$\psi _h,\phi _h\in V_h$|⁠. For |$g\in H$| the finite element approximation of |$v = L^{-1}g$| is then given by |$v_h = L_h^{-1} (\varPi _h g)$|⁠, where |$\varPi _h: H \to V_h$| denotes the |$H$|-orthogonal projection onto the finite element space |$V_h$|⁠, i.e.,

$$\begin{align*} \langle L v_h, \phi_h\rangle = (L_h v_h, \phi_h)_{H} = (\varPi_h g, \phi_h)_{H} = (g, \phi_h)_{H} \qquad \forall \,\phi_h \in V_h. \end{align*}$$

The eigenvalues |$\{\lambda _{j,h}\}_{j=1}^{N_h}$|⁠, as well as the corresponding |$H$|-orthonormal eigenvectors |$\mathcal{E}:= \{e_{j,h}\}_{j=1}^{N_h}$| of |$L_h$|⁠, satisfy the variational equalities

$$\begin{align} ( L_h e_{j,h}, \phi_h )_{H} = \lambda_{j,h} (e_{j,h}, \phi_h)_{H} \qquad \forall \,\phi_h \in V_h, \quad 1\leqslant j \leqslant N_h. \end{align}$$

(2.11)

These eigenvalues are again arranged in nondecreasing order:

$$\begin{align*} 0 < \lambda_{1,h} \leqslant \lambda_{2,h} \leqslant \cdots \leqslant \lambda_{N_h,h}. \end{align*}$$

Further assumptions on the approximation properties of finite element spaces are specified below.

 

Assumption 2.6

The family |$(V_h)_{h\in (0,1)}$| of finite-dimensional subspaces of |$\dot{H}^{1}$| satisfies the following:

• (i) there exists |$d\in \mathbb{N}$| such that |$N_h = \dim (V_h) \propto h^{-d}$| for all |$h> 0$|⁠;

• (ii) there exist constants |$C_1, C_2> 0$| and |$h_0\in (0,1)$|⁠, as well as exponents |$r,s> 0$| and |$q>1$| such that |$\{\lambda _{j,h}\}$| and |$\{e_{j,h}\}$| in (2.11) satisfy

  $$\begin{align} \qquad\qquad\;\;\lambda_j \leqslant \lambda_{j,h} &\leqslant \lambda_j + C_1 h^r \lambda_j^q, \end{align}$$

  (2.12)

  $$\begin{align} \|e_j - e_{j,h}\|_{H}^2 &\leqslant C_2 h^{2s} \lambda_j^q,\qquad\qquad \end{align}$$

  (2.13)

  for all |$h\in (0,h_0)$| and |$j\in \{1,\ldots ,N_h\}$|⁠.

The following example illustrates that Assumption 2.6 is, in general, satisfied for |$d$|-dimensional elliptic linear differential operators.

 

Example 2.7

Let |$\mathcal{D}\subset \mathbb{R}^d$| be a bounded, convex, polygonal domain domain and assume that |$L:\mathcal{D}(L)\subset L_2(\mathcal{D}) \to L_2(\mathcal{D})$| is a (strongly) elliptic linear differential operator of order |$2m\in \mathbb{N}$| with smooth coefficients, i.e., there exists a constant |$\gamma> 0$| such that

$$\begin{align*} \langle Lv, v\rangle \geqslant \gamma\, \|v\|_{H^m(\mathcal{D})}^2 \qquad \forall \,v\in V, \end{align*}$$

where |$V = H^m(\mathcal{D})$| or |$V=H_0^1(\mathcal{D}) \cap H^m(\mathcal{D})$|⁠, and |$\langle \,\cdot \,, \,\cdot \,\rangle $| denotes the duality pairing between |$V$| and its dual |$V^*$| after identification via the inner product on |$L_2(\mathcal{D})$|⁠. Assume that |$V_h \subset V$| is an admissible finite element space of polynomial degree |$p\in \mathbb{N}$| with respect to a regular mesh on |$\bar{\mathcal{D}}$|⁠. In this case Assumption 2.6 is satisfied (Strang & Fix, 2008, Theorems 6.1 and 6.2) for |$H=L_2(\mathcal{D})$|⁠, |$\dot{H}^{1} = V$| and the exponents

$$\begin{align*} r = 2(\,p+1-m), \quad s = \min\{\,p+1, 2(\,p+1-m)\} \quad \textrm{and} \quad q = \tfrac{p+1}{m}. \end{align*}$$

In particular, if the family of finite element spaces |$V_h$| is quasi-uniform and has continuous, piecewise linear basis functions, we have for |$L= \kappa ^2 - \varDelta $| with homogeneous Dirichlet boundary conditions in Example 2.2 that |$r=s=q=2$|⁠.

In order to approximate the white noise |$\mathcal{W}$| on the finite element space |$V_h$| we introduce the following |$V_h$|-valued random variables:

• (i) an expansion with respect to the discrete eigenbasis |$\mathcal{E}$|⁠:

  $$\begin{align} \mathcal{W\,}_{h}^{\mathcal{E}}:= \sum_{j=1}^{N_h} \xi_j \, e_{j,h}, \end{align}$$

  (2.14)

  where |$\boldsymbol{\xi }:= \left ( \xi _1,\ldots ,\xi _{N_h} \right )^{\textrm{T}}$| is the vector of the first |$N_h$| independent standard normally distributed random variables in the expansion (2.9) of |$\mathcal{W}$|⁠;
• (ii) an expansion with respect to any basis |$\varPhi :=\{\phi _{j,h}\}_{j=1}^{N_h}$| of |$V_h$|⁠:

  $$\begin{align} \mathcal{W\,}_{h}^{\varPhi}:= \sum_{j=1}^{N_h} \widetilde{\xi_j} \, \phi_{j,h}, \end{align}$$

  (2.15)

  where the random vector |$\widetilde{\boldsymbol{\xi }}:= ( \widetilde{\xi _1},\ldots , \widetilde{\xi }_{N_h} )^{\textrm{T}}$| is given by

  $$\begin{align*} \widetilde{\boldsymbol{\xi}} = \mathbf{R}^{-1} \boldsymbol{\xi}, \qquad R_{ij}:= (e_{i,h}, \phi_{j,h})_{H}, \qquad 1\leqslant i,j \leqslant N_h. \end{align*}$$

  The vector |$\widetilde{\boldsymbol{\xi }}$| is therefore Gaussian distributed with mean zero and covariance matrix given by |$\mathbf{R}^{-1}(\mathbf{R}^{-1})^{\textrm{T}} = (\mathbf{R}^{\textrm{T}} \mathbf{R})^{-1} = \mathbf{M}^{-1}$|⁠, where |$\mathbf{M} = ( (\phi _{i,h}, \phi _{j,h})_{H} )_{i,j=1}^{N_h}$| is the mass matrix with respect to the basis |$\varPhi $| of the finite element space |$V_h$|⁠.

The following lemma shows that the above approximations of the white noise are equal in a mean-square sense.

 

Lemma 2.8

The noise approximations |$\mathcal{W\,}_h^{\mathcal{E}}$| and |$\mathcal{W\,}_h^{\varPhi} $| in (2.14)–(2.15) are equivalent in |$L_2(\varOmega ; H)$|⁠, i.e., |$\|\mathcal{W\,}_h^{\mathcal{E}} - \mathcal{W\,}_{h}^{\varPhi} \|_{L_2(\varOmega ; H)} = 0$|⁠.

 

Proof.

Inserting the definitions (2.14)–(2.15) of |$\mathcal{W\,}_h^{\mathcal{E}}$| and |$\mathcal{W\,}_h^{\varPhi} $| yields

$$\begin{align*} \mathbb{E}\Big[ \|\mathcal{W\,}_h^{\mathcal{E}} - \mathcal{W\,}_h^\varPhi\|_{H}^2 \Big] &= \sum_{i=1}^{N_h} \sum_{j=1}^{N_h} \Bigl( \mathbb{E}\big[ \xi_i \xi_j \big] ( e_{i,h}, e_{j,h} )_{H} - 2 \, \mathbb{E}\left[\widetilde{\xi_i} \xi_j \right] ( \phi_{i,h}, e_{j,h} )_{H} + \mathbb{E}\left[\widetilde{\xi_i} \widetilde{\xi_j} \right] ( \phi_{i,h}, \phi_{j,h} )_{H} \Bigr) . \end{align*}$$

By definition of the matrices |$\mathbf{R}$| and |$\mathbf{M}$| and due to the distribution of the random vectors |$\boldsymbol{\xi }$| and |$\widetilde{\boldsymbol{\xi }}$| we have for |$i,j\in \{1,\ldots ,N_h\}$|⁠,

$$\begin{align*} \qquad \mathbb{E}[\xi_i \xi_j] &= \delta_{ij}, & ( e_{i,h}, e_{j,h} )_{H} &= \delta_{ij}, \qquad \\ \qquad \mathbb{E}\left[\widetilde{\xi_i} \widetilde{\xi_j}\right]&= \left[ \mathbf{M}^{-1} \right]_{ij}, & ( \phi_{i,h}, \phi_{j,h} )_{H} &= M_{ji}, \qquad \\ \qquad \mathbb{E}\left[\widetilde{\xi_i} \xi_j\right] &= \left[ \mathbf{R}^{-1} \right]_{ij}, & ( \phi_{i,h}, e_{j,h} )_{H} &= R_{ji}, \qquad \end{align*}$$

where |$\delta _{ij}$| is the Kronecker delta. Thus, in terms of the trace |$\operatorname{tr}$| of matrices we have

$$\begin{align*} \mathbb{E}\bigl[ \|\mathcal{W\,}_h^{\mathcal{E}} - \mathcal{W\,}_h^\varPhi\|_{H}^2 \bigr] &= \operatorname{tr}(\mathbf{I}) - 2 \operatorname{tr}\big( \mathbf{R}^{-1} \mathbf{R} \big) + \operatorname{tr}\big( \mathbf{M}^{-1} \mathbf{M} \big) = 0, \end{align*}$$

which proves the equivalence of |$\mathcal{W\,}_h^{\mathcal{E}}$| and |$\mathcal{W\,}_h^{\varPhi} $| in |$L_2(\varOmega ;H)$|⁠.

Our numerical approach to cope with the fractional-order equation (2.8) will be based on the following representation of the inverse |$L^{-\beta }$| from the Dunford–Taylor calculus due to Balakrishnan (1960) (see also Yosida, 1995, §IX.11, equation (4)):

$$\begin{align*} L^{-\beta} &= \frac{\sin(\pi\beta)}{\pi} \int_0^{\infty} \lambda^{-\beta} (\lambda I + L)^{-1} \, \mathrm{d} \,\lambda = \frac{2\sin(\pi\beta)}{\pi} \int_{-\infty}^{\infty} e^{2\beta y} \bigl( I + e^{2y} L \bigr)^{-1} \, \mathrm{d}\, y. \end{align*}$$

We choose an equidistant grid |$\{y_\ell = \ell k: \ell \in \mathbb{Z}, -K^{-} \leqslant \ell \leqslant K^+\}$| with step size |$k>0$| for |$y$| and replace the differential operator |$L$| with its discrete version |$L_h$| to formulate the following quadrature for |$L_h^{-\beta }$| proposed by Bonito & Pasciak (2015):

$$\begin{align} Q^\beta_{h,k} := \frac{2 k \sin(\pi\beta)}{\pi} \sum_{\ell=-K^{-}}^{K^{+}} e^{2\beta y_\ell} \left( \mathrm{Id}_{V_h} + e^{2 y_\ell} L_h \right)^{-1}. \end{align}$$

(2.16)

Exponential convergence of order |$\mathcal{O}\bigl (e^{-\pi ^2/(2k)} \bigr )$| to |$L_h^{-\beta }$| with respect to the norm

$$\begin{align} \|T\|_{\mathcal{L}(V_h)}:= \sup\limits_{\phi_h \in V_h\setminus\{0\}} \frac{\|T\phi_h\|_{H}}{\|\phi_h\|_{H}} \end{align}$$

(2.17)

on the space |$\mathcal{L}(V_h):= \{ T: V_h \to V_h \textrm{ linear}\}$| has been proven (Bonito & Pasciak, 2015, Lemma 3.4, Remark 3.1, Theorem. 3.5) for the choice

$$\begin{align*} K^-:= \left\lceil \frac{\pi^2}{4\beta k^2} \right\rceil, \qquad K^+:= \left\lceil \frac{\pi^2}{4(1-\beta)k^2} \right\rceil. \end{align*}$$

Note that the quadrature (2.16) involves only nonfractional resolvents of the discrete operator |$L_h$|⁠. Thus, the corresponding numerical method is readily implementable.

With the notions of the |$H$|-orthogonal projection |$\varPi _h$| on |$V_h$|⁠, the noise approximation |$\mathcal{W\,}_h^{\varPhi} $| in (2.15) and the quadrature |$Q_{h,k}^\beta $| in (2.16) at hand we can now introduce the numerical approximation |$u_{h,k}^{Q}$| of the solution |$u$| to (2.8) as

$$\begin{align} u_{h,k}^{Q}:= Q_{h,k}^\beta \left( \varPi_h g + \mathcal{W\,}_h^{\varPhi} \right). \end{align}$$

(2.18)

 

Remark 2.9

We emphasize the construction of the noise term |$\mathcal{W\,}_h^{\varPhi} $| on the right-hand side of (2.18). For discretizing the white noise |$\mathcal{W}$| it is common to project the truncated Karhunen–Loève expansion |$\mathcal{W}_N$| in (2.10) onto the finite element space |$V_h$| and use the noise approximation |$\varPi _h \mathcal{W}_N \in L_2(\varOmega ; V_h)$| (see, e.g., Zhang et al., 2016). Instead we define |$\mathcal{W\,}_h^{\varPhi} $| in (2.15) in such a way that it is mean-square equivalent to the noise approximation |$\mathcal{W\,}_h^{\mathcal{E}}$| in (2.14), which can be interpreted as |$V_h$|-valued white noise expanded with respect to the |$H$|-orthonormal eigenvectors of |$L_h$|⁠. Note that the noise approximation |$\mathcal{W\,}_h^{\mathcal{E}}$| is needed only for the error analysis but not for the actual implementation of the numerical algorithm.

This approach has the following two major advantages in practice:

• (i) Samples of the truncated white noise |$\mathcal{W}_N$| are not needed and, thus, neither are the exact eigenvectors |$\{e_j \}$| of the operator |$L$|⁠.

• (ii) For the computation of the approximation |$u_{h,k}^{Q}$| in (2.18) in practice one has to sample from the load vector |$\mathbf{b}$| with entries

  $$\begin{align*} b_i:= (\varPi_h g + \mathcal{W\,}_h^\varPhi, \phi_{i,h})_{H}, \qquad 1\leqslant i \leqslant N_h, \end{align*}$$

  where |$\{\phi _{j,h}\}_{j=1}^{N_h}$| is a basis of the finite element space |$V_h$|⁠. We emphasize that this is computationally feasible if the basis |$\varPhi $| in the noise approximation |$\mathcal{W\,}_h^{\varPhi} $| is chosen as the same, since then |$\mathbf{b} \sim \mathcal{N}(\mathbf{g}, \mathbf{M})$|⁠, where |$g_i = (g, \phi _{i,h})_{H}$| and |$\mathbf{M}$| is the mass matrix with respect to the basis |$\varPhi $|⁠, which is usually sparse. Hence, samples of |$\mathbf{b}$| can be generated from |$\mathbf{b} = \mathbf{g} + \mathbf{G}\mathbf{z}$|⁠, where |$\mathbf{z} \sim \mathcal{N}(\mathbf{0},\mathbf{I})$| and |$\mathbf{G}$| is the Cholesky factor of |$\mathbf{M} = \mathbf{G}\mathbf{G}^{\textrm{T}} $|⁠.

The following estimate of the strong mean-square error between the exact solution |$u$| to the fractional-order equation (2.8) and the numerical approximation |$u^Q_{h,k}$| in (2.18) is our main result.

 

Theorem 2.10

Let the family |$(V_h)_{h\in (0,1)}$| of finite element spaces satisfy Assumption 2.6 and assume that the growth of the eigenvalues of the operator |$L$| is given by (2.1) for an exponent |$\alpha $| with

$$\begin{align} \tfrac{1}{2\beta} < \alpha \leqslant \min\left\{ \tfrac{r}{(q-1)d}, \tfrac{2s}{qd} \right\}, \end{align}$$

(2.19)

where the values of |$d\in \mathbb{N}$|⁠, |$r,s>0$| and |$q>1$| are the same as in Assumption 2.6. Then, for sufficiently small |$h\in (0,h_0)$| and |$k\in (0,k_0)$|⁠, the strong |$L_2(\varOmega ; H)$| error between the solution |$u$| of (2.8) and the approximation |$u_{h,k}^{Q}$| in (2.18) is bounded by

$$\begin{align} \bigl\| u - u_{h,k}^{Q} \bigr\|_{L_2(\varOmega; H)} \leqslant C \left(h^{\min\{d(\alpha\beta-1/2),\, r,\, s\}} + e^{-\pi^2/(2k)} h^{-d/2} \right) (1 + \|g\|_{H}), \end{align}$$

(2.20)

where the constant |$C>0$| is independent of |$h$| and |$k$|⁠.

3. Partition of the error and error estimates

In order to prove Theorem 2.10 we express the difference between the exact solution |$u$| and the approximation |$u_{h,k}^{Q}$| as

$$\begin{align*} u - u_{h,k}^{Q} = \bigl( u - u_{N_h} \bigr) + \bigl( u_{N_h} - u_h^{\mathcal{E}} \bigr) + \bigl( u_h^{\mathcal{E}} - u_h^\varPhi \bigr) + \bigl( u_h^\varPhi - u_{h,k}^{Q} \bigr), \end{align*}$$

and partition the strong error in (2.20) accordingly. Here |$u_{N_h}$|⁠, |$u_h^{\mathcal{E}}$| and |$u_h^\varPhi $| are defined in terms of the truncated white noise in (2.10) and the white noise approximations in (2.14)–(2.15) as the |$\mathbb{P}$|-almost sure solutions to the equations

$$\begin{align} L^{\beta} u_{N} &= g_{N} + \mathcal{W}_{N}, \end{align}$$

(3.1)

for |$N\in \mathbb{N}$|⁠, where |$g_{N}:= \sum _{j=1}^{N} (g, e_j)_{H} \, e_j$|⁠,

$$\begin{align} L_h^{\beta} u_h^{\mathcal{E}} &= \varPi_h g + \mathcal{W\,}_h^{\mathcal{E}}, \end{align}$$

(3.2)

$$\begin{align} L_h^{\beta} u_h^\varPhi &= \varPi_h g + \mathcal{W\,}_h^\varPhi, \end{align}$$

(3.3)

and |$u_{N_h}$| refers to the truncation with |$N=N_h=\dim (V_h)$| terms in (3.1). Note that by Lemma 2.8 the difference between |$u_h^{\mathcal{E}}$| and |$u_h^\varPhi $| vanishes identically in |$L_2(\varOmega ;H)$|⁠:

$$\begin{align} \bigl\| u_h^{\mathcal{E}} - u_h^\varPhi \bigr\|_{L_2(\varOmega;H)} \leqslant \bigl\| L_h^{-\beta} \bigr\|_{\mathcal{L}(V_h)} \bigl\| \mathcal{W\,}_h^{\mathcal{E}} - \mathcal{W\,}_h^\varPhi \bigr\|_{L_2(\varOmega;H)} = 0. \end{align}$$

(3.4)

In the following we address the three remaining terms separately: the truncation error, the error of the finite element discretization and the error caused by the quadrature approximation |$Q^\beta _{h,k}$| in (2.16) of the discrete fractional inverse |$L_h^{-\beta }$|⁠.

3.1 Truncation error

In the lemma below we bound the strong mean-square error between the exact solution |$u$| to the fractional-order equation (2.8) and the approximation |$u_N$| in (3.1) which is based on the truncation of the right-hand side to |$g_N + \mathcal{W}_N$|⁠.

 

Lemma 3.1

Suppose (2.1) for the eigenvalues of the operator |$L$| and that |$2\alpha \beta> 1$|⁠. Let |$u \in L_2(\varOmega ; H)$| be the solution to (2.8). Then, for any |$g\in H$|⁠, |$N\in \mathbb{N}$|⁠, there exists a unique solution |$u_N\in L_2\bigl ( \varOmega ; \dot{H}^{2\beta } \bigr )$| to the truncated equation (3.1) and it satisfies

$$\begin{align*} \mathbb{E}\left[ \|u-u_N\|_{H}^2 \right] \leqslant C N^{-(2\alpha\beta-1)} \left(1 + \|g\|_{H}^2\right), \end{align*}$$

where the constant |$C>0$| depends only on |$\alpha , \beta $| and the constants in (2.1). In particular, under Assumption 2.6(i) we have

$$\begin{align} \|u-u_{N_h}\|_{L_2(\varOmega;H)} \lesssim h^{d(\alpha\beta-1/2)} (1 + \|g\|_{H}). \end{align}$$

(3.5)

 

Proof.

Existence and uniqueness of a solution |$u_N$| in |$L_2\bigl ( \varOmega ; \dot{H}^{2\beta } \bigr )$| follows from the fact that for any |$g\in H$| the truncated right-hand side |$g_N + \mathcal{W}_N$| is an element of |$L_2(\varOmega ; H)$| as well as the isomorphism property of |$L^\beta :\dot{H}^{2\beta }\to H$| in Lemma 2.1. If |$2\alpha \beta> 1$| we obtain for any |$N\in \mathbb{N}$|⁠,

$$\begin{align*} \mathbb{E}\left[ \|u-u_N\|_{H}^2 \right] &= \mathbb{E}\left[ \bigl\| L^{-\beta} ( g - g_N + \mathcal{W} - \mathcal{W}_N ) \bigr\|_{H}^2 \right] = \sum_{j>N} \lambda_j^{-2\beta} \left( (g,e_j)_{H}^2 + 1 \right) \\ &\leqslant c_\lambda^{-2\beta} \left(1 + \|g\|_{H}^2 \right) \sum_{j>N} j^{-2\alpha\beta} \leqslant \frac{N^{1-2\alpha\beta}}{c_\lambda^{2\beta} (2\alpha\beta - 1)} \left(1 + \|g\|_{H}^2\right). \end{align*}$$

Under Assumption 2.6(i) we have |$N_h \propto h^{-d}$| and (3.5) readily follows.

3.2 Finite element discretization error

The next ingredient in the derivation of (2.20) in Theorem 2.10 is an upper bound for the error caused by introducing the finite element discretization. More precisely the solution |$u_N$| of the truncated problem (3.1) corresponds to the best |$\widetilde{V}_N$|-valued approximation of |$u\in L_2(\varOmega ; H)$|⁠. Here the finite-dimensional subspace |$\widetilde{V}_N \subset \dot{H}^{1}$| is the linear span of the first |$N$| eigenvectors of the operator |$L$|⁠. Subsequently the approximation |$u_h^{\mathcal{E}}$| was defined in (3.2), which takes values in another finite-dimensional subspace, namely in the finite element space |$V_h \subset \dot{H}^{1}$|⁠. The purpose of the following lemma is to bound the error between these two approximations when |$N=\dim (V_h)$|⁠.

 

Lemma 3.2

Let Assumption 2.6 be satisfied. Assume that the eigenvalue growth of the operator |$L$| is given by (2.1) for an exponent |$\alpha $| with (2.19). Let |$u_h^{\mathcal{E}}$| be the unique element in |$L_2(\varOmega ; V_h)$| satisfying (3.2), i.e.,

$$\begin{align*} \bigl( L_h u_h^{\mathcal{E}}, \phi_h \bigr)_{H} = \bigl( \varPi_h g + \mathcal{W\,}_h^{\mathcal{E}}, \phi_h \bigr)_{H} \qquad \forall \,\phi_h \in V_h, \quad \mathbb{P}\text{-a.s.,} \end{align*}$$

and let |$u_{N_h} \in L_2(\varOmega ;H)$| denote the solution to the truncated equation (3.1) with |$N=N_h=\dim (V_h)$| terms. Then their difference can be bounded by

$$\begin{align} \bigl\| u_{N_h} - u_h^{\mathcal{E}} \bigr\|_{L_2(\varOmega;H)} \leqslant C h^{\min\{ d(\alpha\beta-1/2),\, r,\, s\}} (1 + \|g\|_{H}), \end{align}$$

(3.6)

for sufficiently small |$h\in (0,h_0)$|⁠, where the constant |$C>0$| depends only on |$\alpha , \beta $| and the constants in (2.1), (2.12) and (2.13).

 

Proof.

In order to show the estimate in (3.6) we first note that |$\varPi _h g$| can be expanded in terms of the orthonormal eigenvectors |$\{e_{j,h}\}$| of |$L_h$| by

$$\begin{align*} \varPi_h g = \sum_{j=1}^{N_h} (\varPi_h g, e_{j,h})_{H} \, e_{j,h} = \sum_{j=1}^{N_h} (g, e_{j,h})_{H} \, e_{j,h}. \end{align*}$$

Thus, we can partition the difference in (3.6) as

$$\begin{align*} \bigl\| u_{N_h} - u_h^{\mathcal{E}} \bigr\|_{L_2(\varOmega;H)} &= \Bigl\| L^{-\beta} \bigl( g_{N_h} + \mathcal{W}_{N_h} \bigr) - L_h^{-\beta} \bigl( \varPi_h g + \mathcal{W\,}_h^{\mathcal{E}} \bigr) \Bigr\|_{L_2(\varOmega;H)} \\ &= \left\| \sum_{j=1}^{N_h} \lambda_j^{-\beta} ((g, e_j)_{H} + \xi_j) e_j - \sum_{j=1}^{N_h} \lambda_{j,h}^{-\beta} ((g, e_{j,h})_{H} + \xi_j) e_{j,h} \right\|_{L_2(\varOmega;H)} \\ &\leqslant \left\| \sum_{j=1}^{N_h} \lambda_j^{-\beta} ((g, e_j)_{H} + \xi_j) (e_j - e_{j,h}) \right\|_{L_2(\varOmega;H)} + \left\| \sum_{j=1}^{N_h} \lambda_j^{-\beta} (g, e_j - e_{j,h})_{H} \, e_{j,h} \right\|_{H} \\ &\quad + \left\| \sum_{j=1}^{N_h} \left(\lambda_j^{-\beta} - \lambda_{j,h}^{-\beta} \right) ( (g, e_{j,h})_{H} + \xi_j) e_{j,h} \right\|_{L_2(\varOmega;H)} =: \text{(I)} + \text{(II)} + \text{(III)}. \end{align*}$$

Since the random variables |$\{\xi _j\}$| are independent and standard normally distributed we obtain for the first term by the Cauchy–Schwarz inequality for sums,

$$\begin{align*} \text{(I)}^2 &= \left\| \sum_{j=1}^{N_h} \lambda_j^{-\beta} (g, e_j)_{H} (e_j - e_{j,h}) \right\|_H^2 + \sum_{j=1}^{N_h} \lambda_j^{-2\beta} \|e_j-e_{j,h}\|_{H}^2 \\ &\leqslant \left( 1 + \sum_{j=1}^{N_h} (g, e_j)_{H}^2 \right) \sum_{j=1}^{N_h} \lambda_j^{-2\beta} \|e_j-e_{j,h}\|_{H}^2 \leqslant \left(1 + \|g\|_{H}^2\right) \sum_{j=1}^{N_h} \lambda_j^{-2\beta} \|e_j-e_{j,h}\|_{H}^2. \end{align*}$$

Assumption 2.6(i), (2.1) and (2.13) complete the estimation of the first term:

$$\begin{align*} \text{(I)}^2 &\lesssim h^{2s} \left(1 + \|g\|_{H}^2\right) \sum_{j=1}^{N_h} \lambda_j^{q-2\beta} \lesssim h^{2s} \left(1 + \|g\|_{H}^2\right) \sum_{j=1}^{N_h} j^{\alpha(q-2\beta)} \\ &\lesssim h^{2s} \left(1 + N_h^{1 + \alpha(q-2\beta)} \right) \left(1 + \|g\|_{H}^2\right) \lesssim \left( h^{2s} + h^{2d(\alpha\beta-1/2)} \right) \left(1 + \|g\|_{H}^2\right), \end{align*}$$

for |$h\in (0,h_0)$|⁠, since |$d \alpha q \leqslant 2s$| by (2.19). For the second term we find

$$\begin{align*} \text{(II)}^2 &= \sum_{j=1}^{N_h} \lambda_j^{-2\beta} (g, e_j - e_{j,h})_{H}^2 \leqslant \|g\|_{H}^2 \sum_{j=1}^{N_h} \lambda_j^{-2\beta} \|e_j-e_{j,h}\|_{H}^2, \end{align*}$$

and conclude as for (I) that |$\text{(II)}^2 \lesssim \bigl ( h^{2s} + h^{2d(\alpha \beta -1/2)} \bigr ) \|g\|_{H}^2$|⁠. For (III) we use again the independence and distribution of the random variables |$\{\xi _j\}$| and obtain

$$\begin{align*} \text{(III)}^2 &= \sum_{j=1}^{N_h} \left( \lambda_j^{-\beta} - \lambda_{j,h}^{-\beta} \right)^2 \left( (g, e_{j,h})_{H}^2 + 1 \right) \leqslant \left( 1 + \|g\|_{H}^2 \right) \sum_{j=1}^{N_h} \left(\lambda_j^{-\beta} - \lambda_{j,h}^{-\beta} \right)^2. \end{align*}$$

By the mean value theorem there exists |$\widetilde{\lambda }_j \in (\lambda _j, \lambda _{j,h})$| such that

$$\begin{align*} \lambda_j^{-\beta} - \lambda_{j,h}^{-\beta} = \beta \widetilde{\lambda}_j^{-\beta-1} (\lambda_{j,h} - \lambda_j ) \leqslant \beta \lambda_j^{-\beta-1} (\lambda_{j,h} - \lambda_j ) \qquad \forall \,j \in \{1,\ldots,N_h\}. \end{align*}$$

Therefore, we can bound the third term by (2.1) and (2.12) as

$$\begin{align*} \text{(III)}^2 &\lesssim \left(1 + \|g\|_{H}^2\right) \sum_{j=1}^{N_h} \lambda_j^{-2\beta-2} (\lambda_{j,h} - \lambda_j)^2 \lesssim h^{2r} \left(1 + \|g\|_{H}^2\right) \sum_{j=1}^{N_h} \lambda_j^{2(q-\beta-1)} \\[8pt] &\lesssim h^{2r} \left(1 + \|g\|_{H}^2\right) \sum_{j=1}^{N_h} j^{2\alpha(q-\beta-1)} \lesssim \left(h^{2r} + h^{2d(\alpha\beta - 1/2)} \right) \left(1 + \|g\|_{H}^2\right), \end{align*}$$

where we used the relations |$N_h \propto h^{-d}$| from Assumption 2.6(i) and |$\alpha d (q-1) \leqslant r$| from (2.19) in the last estimate.

3.3 Quadrature approximation error

As the final step for proving the strong convergence result of Theorem 2.10 we investigate the error caused by applying the quadrature |$Q^\beta _{h,k}$| in (2.16) instead of the discrete fractional inverse |$L_h^{-\beta }$| to the right-hand side |$\varPi _h g + \mathcal{W\,}_h^{\varPhi} $|⁠. The difference between these two operators in the norm (2.17) on the space |$\mathcal{L}(V_h)$| has been bounded by Bonito & Pasciak (2015). The following lemma is a consequence of that result and the distribution of the noise approximation |$\mathcal{W\,}_h^{\varPhi} $|⁠.

 

Lemma 3.3

Let |$Q^\beta _{h,k}: V_h \to V_h$| be the operator in (2.16) approximating |$L_h^{-\beta }$|⁠. Then, for sufficiently small |$k\in (0,k_0)$| and any |$\phi _h \in V_h$|⁠, the estimate

$$\begin{align} \bigl\| \bigl( Q^\beta_{h,k} - L_h^{-\beta} \bigr) \phi_h \bigr\|_{H} \leqslant C e^{-\pi^2/(2k)} \|\phi_h\|_{H} \end{align}$$

(3.7)

holds, where the constant |$C>0$| depends on |$\beta $| and grows linearly in the largest eigenvalue |$\lambda _{1,h}^{-1}$| of |$L_h^{-1}$|⁠. In particular, for |$u_h^\varPhi $| in (3.3) and |$u_{h,k}^Q$| in (2.18), we have

$$\begin{align} \bigl\| u_h^\varPhi - u_{h,k}^Q \bigr\|_{L_2(\varOmega; H)} \lesssim e^{-\pi^2/(2k)} \left( h^{-d/2} + \|g\|_{H} \right). \end{align}$$

(3.8)

 

Proof.

The first assertion is proven in Bonito & Pasciak (2015, Lemma 3.4, Theorem. 3.5). This bound, |$\|\varPi _h g\|_{H} \leqslant \|g\|_{H}$|⁠, and |$\mathbb{E}\bigl [ \|\mathcal{W\,}_h^{\varPhi} \|_{H}^2 \bigr ] = N_h \lesssim h^{-d}$| imply (3.8).

3.4 Proof of Theorem 2.10 and some remarks

After having bounded the truncation, the discretization and the quadrature errors in Sections 3.1–3.3 the strong convergence result of Theorem 2.10 is an immediate consequence.

 

Proof of Theorem 2.10.

As outlined at the beginning of the section we partition the mean-square error as

$$\begin{align*} \bigl\| u - u_{h,k}^Q \bigr\|_{L_2(\varOmega; H)} &\leqslant \bigl\| u - u_{N_h} \bigr\|_{L_2(\varOmega; H)} + \bigl\| u_{N_h} - u_h^{\mathcal{E}} \bigr\|_{L_2(\varOmega; H)} \\ &\quad + \bigl\| u_h^{\mathcal{E}} - u_h^\varPhi \bigr\|_{L_2(\varOmega; H)} + \bigl\| u_h^\varPhi - u_{h,k}^Q \bigr\|_{L_2(\varOmega; H)} \\ &=: \text{(I)} + \text{(II)} + \text{(III)} + \text{(IV)}. \end{align*}$$

By (3.5), (3.6) and (3.8) of Lemmas 3.1–3.3 we have

$$\begin{align*} \text{(I)} &\lesssim h^{d(\alpha\beta-1/2)} (1 + \|g\|_{H}), \\ \text{(II)} &\lesssim h^{\min\{d(\alpha\beta-1/2),\, r,\, s\}} (1 + \|g\|_{H}), \\ \text{(IV)} &\lesssim e^{-\pi^2/(2k)} \left( h^{-d/2} + \|g\|_{H} \right) \lesssim e^{-\pi^2/(2k)} h^{-d/2} (1 + \|g\|_{H}), \end{align*}$$

and by (3.4) the third term vanishes: |$\text{(III)} = 0$|⁠. Thus, the assertion is proven.

In Section 4 we will not only verify the above derived rate of strong convergence by means of numerical experiments but also investigate weak-type errors. The following result is proven similarly to Theorem 2.10.

 

Corollary 3.4

Let the assumptions of Theorem 2.10 be satisfied. Then there exists a constant |$C>0$| independent of |$h\in (0,h_0)$| and |$k\in (0,k_0)$| such that the following weak-type error estimate between the approximation |$u_{h,k}^Q$| in (2.18) and the solution |$u$| to (2.8) holds:

$$\begin{align} \Bigl| \|u\|_{L_2(\varOmega; H)}^2 - \bigl\| u_{h,k}^Q \bigr\|_{L_2(\varOmega;H)}^2 \Bigr| &\leqslant C \left( h^{\min\{ d(2\alpha\beta-1),\, r,\, s\}} + e^{- \pi^2/(2k)} \right) \|g\|_{H}^2 \\ &\quad + C \left( h^{\min\{ d(2\alpha\beta-1),\, r \}} + e^{- \pi^2 / k } h^{-d} + e^{- \pi^2 / (2k)} + e^{- \pi^2 / (2k)} f_{\alpha,\beta}(h) \right), \notag \end{align}$$

(3.9)

where |$\,f_{\alpha ,\beta }(h):= h^{d(\alpha \beta -1)}$|⁠, if |$\alpha \beta \neq 1$|⁠, and |$\,f_{\alpha ,\beta }(h):= |\ln (h)|$| if |$\alpha \beta = 1$|⁠.

 

Proof.

First we note that the norm on |$L_2(\varOmega ;H)$| attains the following values for |$u$| in (2.8), |$u_{N_h}$| in (3.1) and |$u_h^\varPhi $| in (3.3):

$$\begin{align*} \|u\|_{L_2(\varOmega;H)}^2 &= \sum_{j\in\mathbb{N}} \lambda_j^{-2\beta} \left( 1 + (g,e_j)_{H}^2 \right), \\ \bigl\| u_{N_h} \bigr\|_{L_2(\varOmega;H)}^2 &= \sum_{j=1}^{N_h} \lambda_j^{-2\beta} \left( 1 + (g,e_j)_{H}^2 \right), \\ \bigl\| u_h^\varPhi \bigr\|_{L_2(\varOmega;H)}^2 &= \sum_{j=1}^{N_h} \lambda_{j,h}^{-2\beta} \left( 1 + (g,e_{j,h})_{H}^2 \right) = \bigl\| L_h^{-\beta} \varPi_h g \bigr\|_{H}^2 + \bigl\| L_h^{-\beta} \mathcal{W\,}_h^{\mathcal{E}} \bigr\|_{L_2(\varOmega;H)}^2. \end{align*}$$

Again we partition the error,

$$\begin{align*} \Bigl| \|u\|_{L_2(\varOmega; H)}^2 - \bigl\| u_{h,k}^Q \bigr\|_{L_2(\varOmega;H)}^2 \Bigr| &\leqslant \Bigl| \|u\|_{L_2(\varOmega; H)}^2 - \bigl\| u_{N_h} \bigr\|_{L_2(\varOmega;H)}^2 \Bigr| + \Bigl| \bigl\| u_{N_h} \bigr\|_{L_2(\varOmega; H)}^2 - \bigl\| u_h^\varPhi \bigr\|_{L_2(\varOmega;H)}^2 \Bigr| \\ &\quad + \Bigl| \bigl\| u_h^\varPhi \bigr\|_{L_2(\varOmega; H)}^2 - \bigl\| u_{h,k}^Q \bigr\|_{L_2(\varOmega;H)}^2 \Bigr| \\ &=: \text{(I)} + \text{(II)} + \text{(III)}, \end{align*}$$

and bound every term separately. By applying similar steps as in the proofs of Lemmas 3.1–3.2 we obtain for the first two terms

$$\begin{align*} \text{(I)} &= \sum_{j>N_h} \lambda_j^{-2\beta} \left( 1 + (g,e_j)_{H}^2 \right) \lesssim h^{d(2\alpha\beta-1)} \left( 1 + \|g\|_{H}^2 \right), \end{align*}$$

as well as

$$\begin{align*} \text{(II)} &= \sum_{j=1}^{N_h} \left( \lambda_j^{-2\beta} - \lambda_{j,h}^{-2\beta} \right) \left( 1 + (g,e_{j,h})_{H}^2 \right) + \sum_{j=1}^{N_h} \lambda_j^{-2\beta} (g,e_{j} - e_{j,h})_{H} (g,e_j + e_{j,h})_{H} \\ &\leqslant 2\beta \left( 1 + \|g\|_{H}^2\right) \sum_{j=1}^{N_h} \lambda_j^{-2\beta-1} (\lambda_{j,h} - \lambda_j) + 2 \, \|g\|_{H}^2 \sum_{j=1}^{N_h} \lambda_j^{-2\beta} \|e_{j} - e_{j,h}\|_{H} \\ &\lesssim h^r \left( 1 + \|g\|_{H}^2 \right) \sum_{j=1}^{N_h} \lambda_j^{q-2\beta-1} + h^s \|g\|_{H}^2 \sum_{j=1}^{N_h} \lambda_j^{q/2 - 2\beta} \\ &\lesssim h^{\min\{ d(2\alpha\beta-1),\, r \}} \left( 1 + \|g\|_{H}^2 \right) + h^{\min\{ d(2\alpha\beta-1), \, s\}} \|g\|_{H}^2, \end{align*}$$

since |$d \alpha (q-1) \leqslant r$| and |$d\alpha q/2 \leqslant s$| as assumed in (2.19). For the third term we find

$$\begin{align*} \text{(III)} &\leqslant \Bigl| \bigl\| Q_{h,k}^\beta \varPi_h g \bigr\|_{H}^2 - \bigl\| L_h^{-\beta} \varPi_h g \bigr\|_{H}^2 \Bigr| + \Bigl| \bigl\| Q_{h,k}^{\beta} \mathcal{W\,}_h^{\mathcal{E}} \bigr\|_{L_2(\varOmega;H)}^2 - \bigl\| L_h^{-\beta} \mathcal{W\,}_h^{\mathcal{E}} \bigr\|_{L_2(\varOmega;H)}^2 \Bigr| \\ &=: \text{(IIIa)} + \text{(IIIb)}. \end{align*}$$

Since |$\|L_h^{-\beta }\|_{\mathcal{L}(V_h)} = \lambda _{1,h}^{-\beta }$| and |$\|Q_{h,k}^\beta \|_{\mathcal{L}(V_h)} \lesssim \lambda _{1,h}^{-\beta }$| for sufficiently small |$k\in (0,k_0)$| we conclude with (3.7) of Lemma 3.3 that

$$\begin{align*} \text{(IIIa)} &= \Bigl| \bigl( \bigl( Q_{h,k}^\beta - L_h^{-\beta} \bigr) \varPi_h g, \bigl(Q_{h,k}^\beta + L_h^{-\beta} \bigr) \varPi_h g \bigr)_{H} \Bigr| \lesssim e^{-\pi^2/(2k)} \|g\|_{H}^2, \\[8pt] \text{(IIIb)} &= \Bigl| \bigl\| \bigl( Q_{h,k}^\beta - L_h^{-\beta} \bigr) \mathcal{W\,}_h^{\mathcal{E}} \bigr\|_{L_2(\varOmega;H)}^2 + 2 \, \bigl( \bigl( Q_{h,k}^\beta - L_h^{-\beta} \bigr) \mathcal{W\,}_h^{\mathcal{E}}, L_h^{-\beta} \mathcal{W\,}_h^{\mathcal{E}} \bigr)_{L_2(\varOmega; H)} \Bigr| \\[8pt] &\leqslant \sum_{j=1}^{N_h} \bigl\| \bigl( Q_{h,k}^\beta - L_h^{-\beta} \bigr) e_{j,h} \bigr\|_{H}^2 + 2 \, \Bigl| \sum_{j=1}^{N_h} \lambda_{j,h}^{-\beta} \bigl( \bigl(Q_{h,k}^\beta - L_h^{\beta} \bigr) e_{j,h}, e_{j,h} \bigr)_{H} \Bigr| \\[8pt] &\lesssim e^{-\pi^2/k} N_h + e^{-\pi^2/(2k)} \sum_{j=1}^{N_h} \lambda_{j,h}^{-\beta} \lesssim e^{-\pi^2/k} h^{-d} + e^{-\pi^2/(2k)} \left( 1 + f_{\alpha,\beta}(h) \right), \end{align*}$$

where we have used Assumption 2.6 in the last estimate. This proves (3.9).

 

Remark 3.5

The error estimates in (2.20) and (3.9) imply that the distance of the quadrature nodes |$k$| has to be adjusted to the finite element mesh size |$h$|⁠. Table 1 shows the calibration between |$h$| and |$k$| for error studies of strong and weak type as well as the corresponding theoretical convergence rates. For the calibration we set

$$\begin{align*} K_{\alpha,\beta}(h) := \begin{cases} d\alpha\beta \ln(h), & \alpha\beta < 1, \\[2pt] d \ln(h) - \max\{0, \ln(|\ln(h)|)\}, & \alpha\beta = 1, \\[2pt] d(2\alpha\beta - 1) \ln(h), & \alpha\beta> 1. \end{cases} \end{align*}$$

Table 1

Theoretical strong and weak-type convergence rates

	Calibration	Rate of convergence
Strong error (2.20)	|$k\leqslant -\tfrac{\pi ^2}{2d\alpha \beta \ln (h)}$|	|$\min \{d(\alpha \beta -1/2),r,s\}$|
Weak-type error (3.9)	|$k\leqslant -\tfrac{\pi ^2}{2 K_{\alpha ,\beta }(h)} $|	$$\begin{cases} \min \{d(2\alpha \beta -1),r\} \quad\quad \textrm{if} \ g=0 \\ \min \{d(2\alpha \beta -1),r,s\} \quad \ \textrm{otherwise} \end{cases}$$

 

Remark 3.6

In the nonfractional case |$\beta = 1$| the discretized problem (3.3) is also nonfractional. Therefore realizations of its solution |$u_h^\varPhi = L_h^{-1}(\varPi _h g + \mathcal{W\,}_h^{\varPhi} )$| can be computed directly and no quadrature is needed. If |$\beta \in (n,n+1]$| for some |$n\in \mathbb{N}$| one may apply the above-described method and error estimates to |$\widetilde{L}:= L^{n+1}$| and |$\widetilde{\beta }:= \tfrac{\beta }{n+1} \in (0,1]$|⁠. Since, however, the finite element theory for the operator |$\widetilde{L}$| may not be trivial, an alternative is to use an ‘iterated finite element method’, i.e., to define the approximation by |$Q_{h,k}^{\beta -n} L_h^{-n} (\varPi _h g + \mathcal{W\,}_h^{\varPhi} )$|⁠.

4. An application and numerical examples

In the following numerical experiment we take up the SPDE from (1.2) in Section 1. More precisely, with the objective of generating computationally efficient approximations of Gaussian Matérn fields on the unit cube |$\mathcal{D}:= (0,1)^d$| in |$d=1,2,3$| spatial dimensions, we consider the following problem:

$$\begin{align} \big(\kappa^2 - \varDelta\big)^{\beta} u(\mathbf{x}) &= \mathcal{W}(\mathbf{x}), && \mathbf{x} \in \mathcal{D}, \end{align}$$

(4.1a)

$$\begin{align} u(\mathbf{x}) &= 0, && \mathbf{x} \in \partial \mathcal{D}, \end{align}$$

(4.1b)

and study the above-presented numerical method generating the approximation |$u_{h,k}^Q$| in (2.18). As already observed in Example 2.2 the exponent of the eigenvalue growth is given by |$\alpha = 2/d$| in this case.

Furthermore, using a finite element discretization on uniform meshes with continuous, piecewise linear basis functions, Assumption 2.6 is satisfied for this problem in all three dimensions for the values |$r=s=q=2$|⁠; see Example 2.7. The condition in (2.19) of Theorem 2.10 becomes |$\beta> d/4$|⁠. We emphasize that this assumption is meaningful also from the statistical point of view; on all of |$\mathbb{R}^d$|⁠, |$\beta\!>\! d/4$| corresponds to a positive smoothness parameter |$\nu\!>\! 0$| of the Matérn covariance function (1.1).

Thus, if the quadrature step size |$k$| and the finite element mesh width |$h$| are calibrated as indicated in Table 1 the theoretical rates of convergence for |$\beta \in (d/4,1)$| are |$2\beta - d/2$| for the strong error and |$\min \{4\beta -d,2\}$| for the weak-type error according to Theorem 2.10 and Corollary 3.4.

For problem (4.1) with |$\kappa =0.5$| we investigate the empirical convergence rates of the following:

• (i) strong error for |$d=1,2,3$| and |$\beta = \tfrac{2d + n}{8}$|⁠, |$n\in \{1,\ldots ,7-2d\}$|⁠;

• (ii) weak-type error for |$d=1,2,3$| and |$\beta = \tfrac{2d + 1}{8}$|⁠.

In each dimension we use a finite element method in space with continuous, piecewise linear basis functions on uniform meshes with mesh diameter |$h$|⁠, mesh nodes |$\mathbf{x}_1,\ldots , \mathbf{x}_{N_h}$| and mass matrix |$\mathbf{M}$|⁠. We choose |$k = -1/(\beta \ln h)$|⁠. The numbers of finite element basis functions and the corresponding numbers of quadrature nodes depending on |$\beta $| used in the strong error study are shown in Table 2.

Table 2

Numbers of finite element basis functions on the considered meshes for |$d=1,2,3$| as well as the corresponding numbers of quadrature nodes as a function of |$\beta $| for the strong error study

		|$\beta $|
	|$N_h$|	|$3/8$|	|$4/8$|	|$5/8$|	|$6/8$|	|$7/8$|
|$d=1$|	127	37	61	99	176	408
	255	48	77	129	229	533
	511	60	99	163	291	675
	1023	73	121	200	357	832
|$d=2$|	961	---	---	43	75	171
	3969	---	---	62	109	253
	16129	---	---	86	152	352
	65025	---	---	113	203	469
|$d=3$|	729	---	---	---	---	55
	6859	---	---	---	---	105
	59319	---	---	---	---	172

For (i), measuring the strong mean-square error between the exact solution |$u$| to our model problem (4.1) and the approximation |$u_{h,k}^Q$| in (2.18), we proceed as follows: first, samples of an overkill approximation of the white noise

$$\begin{align*} \mathcal{W}_{\textrm{ok}}:= \sum_{\theta_1=1}^{N_{\textrm{ok}}} \ldots \sum_{\theta_d=1}^{N_{\textrm{ok}}} \xi_{(\theta_1,\ldots,\theta_d)} e_{(\theta_1,\ldots,\theta_d)} \end{align*}$$

are generated and evaluated on a uniform overkill mesh of |$\bar{\mathcal{D}}$| with |$N_{\textrm{ok}}^d$| nodes. Here |$\{\xi _{\boldsymbol{\theta }}\}$| are independent standard normally distributed random variables and |$\{e_{\boldsymbol{\theta }}\}$| are the eigenfunctions in (2.6). The approximation |$\mathcal{W}_{\textrm{ok}}$| corresponds to a truncated Karhunen–Loève expansion with |$N_{\textrm{ok}}^d$| terms. From this samples of the overkill solution |$u_{\textrm{ok}}$| are obtained via

$$\begin{align*} u_{\textrm{ok}}:= \big( \kappa^2 - \varDelta \big)^{-\beta} \mathcal{W}_{\textrm{ok}} = \sum_{\theta_1=1}^{N_{\textrm{ok}}} \ldots \sum_{\theta_d=1}^{N_{\textrm{ok}}} \lambda_{(\theta_1,\ldots,\theta_d)}^{-\beta} \xi_{(\theta_1,\ldots,\theta_d)} e_{(\theta_1,\ldots,\theta_d)}, \end{align*}$$

where |$\{\lambda _{\boldsymbol{\theta }}\}$| are the eigenvalues from (2.7). For the sake of generating comparable samples of the approximation |$u_{h,k}^Q$| we consider the same realizations of |$\mathcal{W}_{\textrm{ok}}$| and use the load vector |$\widetilde{\mathbf{b}}$| with entries |$\widetilde{b}_i = (\mathcal{W}_{\textrm{ok}}, \phi _{i,h})_{L_2(\mathcal{D})}$| instead of |$\mathbf{b} \sim \mathcal{N}(\mathbf{0},\mathbf{M})$| from Remark 2.9, as |$\bigl ( (\mathcal{W},\phi _{i,h})_{L_2(\mathcal{D})} \bigr )_{i=1}^{N_h} \sim \mathcal{N}(\mathbf{0},\mathbf{M}),$| and we treat |$\mathcal{W}_{\textrm{ok}}$| as the true white noise. The resulting approximation is denoted by |$\widetilde{u}_{h,k}$|⁠. We choose |$N_{\textrm{ok}} = 2^{18} + 1$| for |$d=1$|⁠, |$N_{\textrm{ok}} = 2^{12} + 1$| for |$d=2$| and |$N_{\textrm{ok}} = 5 \cdot 2^{6} + 1$| for |$d=3$|⁠.

For each value of |$\beta $| and in every spatial dimension we use |$50$| samples of |$\mathcal{W}_{\textrm{ok}}$| to generate samples |$u_{\textrm{ok}}^{(1)}, \ldots , u_{\textrm{ok}}^{(50)}$| of the overkill solution and of the numerical approximation |$\widetilde{u}_{h,k}^{(1)}, \ldots ,\widetilde{u}_{h,k}^{(50)}$| for every mesh size |$h$|⁠. The observed strong errors are then computed as the average |$L_2$|-errors

$$\begin{align*} \mathrm{err}:= \tfrac{1}{50} \sum_{i=1}^{50} \sqrt{ \bigl(\mathbf{v}^{(i)}\bigr)^{\textrm{T}} \mathbf{M} \, \mathbf{v}^{(i)}}, \qquad v^{(i)}_j:= u_{\textrm{ok}}^{(i)}(\mathbf{x}_j) - \widetilde{u}_{h,k}^{(i)}(\mathbf{x}_j). \end{align*}$$

The results are shown in the first three panels of Fig. 1. The data set |$\{(h_\ell , \mathrm{err}_\ell )\}_\ell $| is then used to compute the observed rate of convergence |$\mathrm{r}$| as the least-squares solution to the linear regression |$\ln \mathrm{err} = \mathrm{c} + \mathrm{r} \ln h$| for each combination of |$(d,\beta )$|⁠. As shown in Table 3 the resulting observed rates of convergence are in accordance with the theoretical values predicted by Theorem 2.10.

Fig. 1.

The panels in the first row show the observed strong error for different values of |$\beta $| and |$d=1,2$|⁠. In the third graph the observed strong error for |$\beta = 7/8$| and |$d=1,2,3$| is displayed. The corresponding observed strong convergence rates are shown in Table 3. The final graph shows the observed average weak-type errors for three experiments in |$d=1,2,3$| dimensions. For all of them the theoretical rate of convergence is |$0.5$| and the observed rates are shown in the legend. All errors are shown as functions of the mesh size |$h$| used in the computations in a log-log scale.

For (ii) we study the weak-type error |$\bigl | \mathbb{E}\bigl [ \|u\|_{L_2(\mathcal{D})}^2 \bigr ] - \mathbb{E}\bigl [ \|u_{h,k}^Q\|_{L_2(\mathcal{D})}^2 \bigr ] \bigr |$| addressed in Corollary 3.4. Note that for this study no sample-wise comparison is needed and thus the load vector is sampled via |$\mathbf{b} \sim \mathcal{N\,}(\mathbf{0},\mathbf{M})$| as discussed in Remark 2.9. The variance of the exact solution can be computed directly from the known eigenvalues of the differential operator as |$\mathbb{E}\bigl [ \|u\|_{L_2(\mathcal{D})}^2 \bigr ] = \sum _{\boldsymbol{\theta }} \lambda _{\boldsymbol{\theta }}^{-2\beta }$|⁠. The variance of |$u_{h,k}^Q$| is approximated via Monte Carlo integration by

$$\begin{align*} \mathbb{E}\left[ \|u_{h,k}^Q\|_{L_2(\mathcal{D})}^2 \right] \approx \tfrac{1}{N_{\mathrm{MC}}} \sum_{i=1}^{N_{\mathrm{MC}}} \bigl( \mathbf{u}^{(i)} \bigr)^{\textrm{T}} \mathbf{M} \, \mathbf{u}^{(i)}, \qquad u_j^{(i)}:= u_{h,k}^{(i)}(\mathbf{x}_j), \end{align*}$$

where the number of Monte Carlo samples is |$N_{\mathrm{MC}} = 10^3$| and |$u_{h,k}^{(i)}$| denotes a realization of the numerical approximation |$u_{h,k}^Q$|⁠. The observed weak-type errors and the observed rates of convergence are displayed in the fourth panel of Fig. 1. The theoretical rate of convergence predicted by Corollary 3.4 is 0.5 for all three cases.

Table 3

Observed (respectively theoretical) rates of convergence for the strong errors shown in Fig. 1

	|$\beta $|
	|$3/8$|	|$4/8$|	|$5/8$|	|$6/8$|	|$7/8$|
|$d=1$|	0.25 (0.25)	0.50 (0.5)	0.75 (0.75)	1.00 (1)	1.21 (1.25)
|$d=2$|	---	---	0.29 (0.25)	0.51 (0.5)	0.74 (0.75)
|$d=3$|	---	---	---	---	0.26 (0.25)

5. Conclusion

We have considered the fractional-order equation (2.8) with Gaussian white noise in a Hilbert space setting. We have shown that the fractional operator |$L^\beta $| extends to an isometric isomorphism between the |$\dot{H}$|-spaces in (2.4). From this result and the mean-square regularity of the white noise with respect to the |$\dot{H}$|-spaces we have deduced existence and uniqueness of a solution |$u$| to (2.8) with a certain regularity.

We have proposed the approximation |$u_{h,k}^Q$| in (2.18) based on two numerical ingredients: (i). finite-dimensional subspaces |$V_h$| of |$\dot{H}^{1}$| and (ii) the quadrature approximation (2.16) of the inverse fractional operator. The most advantageous and novel properties of the corresponding numerical scheme are (a) that only solutions to integer-order (i.e., local) elliptic equations have to be computed and (b) that it does not require the knowledge of the eigenfunctions of the differential operator |$L$|⁠.

Our main result, Theorem 2.10, shows strong mean-square convergence of the approximation |$u_{h,k}^Q$| to the exact solution |$u$|⁠. If the quadrature step size |$k$| and the finite element mesh width |$h$| are calibrated appropriately, see Table 1, the resulting strong convergence rate depends only on the fractional-order |$\beta $|⁠, the dimension |$d$|⁠, the eigenvalue growth |$\alpha $| of the operator |$L$| and the approximation properties of the finite element spaces. In order to prove this result we have partitioned the strong error into three terms: the truncation error, the error caused by the finite element discretization and the error of the quadrature approximation. We have derived bounds for each of these error terms separately in Sections 3.1–3.3. By means of similar techniques we have proven a weak-type error estimate in Corollary 3.4.

Finally, in Section 4 we have applied the proposed numerical method to an explicit problem with relevance for spatial statistics; the solution |$u$| to (4.1) can be interpreted as an approximation of a Gaussian Matérn field on the unit cube |$(0,1)^d$|⁠. The performed numerical experiments with continuous, piecewise linear finite element basis functions in |$d=1,2,3$| spatial dimensions verify the derived theoretical strong and weak-type convergence rates; see Fig. 1 and Table 3.

We hope that these results and insights will prove valuable for applications in spatial statistics, which often require sampling from (approximations of) Gaussian Matérn fields and their various extensions.

Acknowledgements

The authors thank Stefan A. Funken for his contribution to the implementations of the finite element method in MATLAB.

Funding

Swedish Research Council (2016-04187 and 2017-04274); Knut and Alice Wallenberg Foundation (KAW 20012.0067).