Numerical solution of fractional elliptic stochastic PDEs with spatial white noise

TL;DR

This paper develops an efficient numerical method for solving fractional elliptic stochastic PDEs with spatial white noise, using integral representations, quadrature, finite element discretization, and noise approximation, with proven convergence rates.

Contribution

It introduces a novel numerical approach combining integral representation, quadrature, and finite element methods for fractional elliptic SPDEs driven by white noise, with explicit error analysis.

Findings

Numerical experiments confirm theoretical convergence rates.

Method effectively approximates solutions in 1D, 2D, and 3D.

Error analysis provides explicit strong mean-square convergence rates.

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 non-local. 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 power operator $L^{-\beta}$ is approximated by a weighted sum of non-fractional resolvents $(I + t_j^2 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 the theoretical results.

Tables

Table 1. Table 1. Theoretical strong and weak type convergence rates

	calibration	rate of convergence
strong error (2.20)	$k\le -\frac{{\pi }^2}{2d\alpha \beta \ln \left(h\right)}$	$\min \{d\left(\alpha \beta -1/2\right),r,s\}$
weak type error (3.9)	$k\le -\frac{{\pi }^2}{2K_{\alpha ,\beta }\left(h\right)}$

Table 2. Table 2. Numbers of finite element basis functions on the considered meshes for d = 1 , 2 , 3 𝑑 1 2 3 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

Table 3. Table 3. Observed (resp. theoretical) rates of convergence for the strong errors shown in Figure  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)

Equations

$$C({\mathbf{x}}_1,{\mathbf{x}}_2)=\frac{2^{1-\nu }{\sigma }^2}{\Gamma (\nu )}(\kappa \Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert {)}^{\nu }{K}_{\nu }(\kappa \Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert ),{\mathbf{x}}_1,{\mathbf{x}}_2\in \mathcal{D},$$

$$({\kappa }^2-\Delta {)}^{\beta }u(\mathbf{x})=\mathcal{W}(\mathbf{x}),\mathbf{x}\in \mathcal{D}.$$

$$0<{\lambda }_1\le {\lambda }_2\le \dots \le {\lambda }_j\le {\lambda }_{j+1}\le \dots ,\underset{j\to \infty }{\lim }{\lambda }_j=\infty .$$

$$c_{\lambda }{j}^{\alpha }\le {\lambda }_j\le C_{\lambda }{j}^{\alpha }\forall j\in \mathbb{N}.$$

$$L^{\beta }\varphi :=\sum _{j\in \mathbb{N}}{\lambda }_j^{\beta }(\varphi ,e_j{)}_H{e}_j.$$

$$(\varphi ,\psi {)}_{2\beta }:=(L^{\beta }\varphi ,L^{\beta }\psi {)}_H,\Vert \varphi {\Vert }_{2\beta }^2:=\Vert L^{\beta }\varphi {\Vert }_H^2=\sum _{j\in \mathbb{N}}{\lambda }_j^{2\beta }(\varphi ,e_j{)}_H^2.$$

$$\displaystyle\|g\|_{-s}=\sup\limits_{\phi\in\dot{H}^{s}\setminus\{0\}}\frac{\langle g,\phi\rangle}{\|\phi\|_{s}}=\biggl{(}\sum_{j\in\mathbb{N}}\lambda_{j}^{-s}\langle g,e_{j}\rangle^{2}\biggr{)}^{1/2},$$

$${\dot{H}}^s\hookrightarrow {\dot{H}}^r\hookrightarrow {\dot{H}}^0:=H\cong {\dot{H}}^{-0}\hookrightarrow {\dot{H}}^{-r}\hookrightarrow {\dot{H}}^{-s},0\le r\le s.$$

$$\Vert L^{\beta }\varphi {\Vert }_{s-2\beta }^2$$

$$⟨L^{\beta }\varphi ,\psi ⟩=(L^{\beta }\varphi ,\psi {)}_H=\sum _{j\in \mathbb{N}}{\lambda }_j^{\beta }(\varphi ,e_j{)}_H(\psi ,e_j{)}_H\le \Vert \varphi {\Vert }_s\Vert \psi {\Vert }_{2\beta -s}.$$

$$({\kappa }^2-\Delta )e$$

$$e$$

$${\dot{H}}^s$$

$${\dot{H}}^{-s}$$

$${\dot{H}}^s=[{\dot{H}}^1,{\dot{H}}^2{]}_{s-1,2}=[H_0^1(\mathcal{D}),H^2(\mathcal{D})\cap H_0^1(\mathcal{D}){]}_{s-1,2}=H^s(\mathcal{D})\cap H_0^1(\mathcal{D}).$$

$$e_{\mathbf{\theta }}(\mathbf{x})=e_{({\theta }_1,\dots ,{\theta }_d)}(x_1,\dots ,x_d)=\prod _{i=1}^d\left(\sqrt{2}\sin (\pi {\theta }_i{x}_i)\right),$$

$${\lambda }_{\mathbf{\theta }}={\kappa }^2+{\pi }^2|\mathbf{\theta }{|}^2={\kappa }^2+{\pi }^2\sum _{i=1}^d{\theta }_i^2.$$

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

$$\mathcal{W}=\sum _{j\in \mathbb{N}}{\xi }_j{e}_j,$$

$${\mathcal{W}}_N:=\sum _{j=1}^N{\xi }_j{e}_j.$$

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

$$\displaystyle\mathbb{E}\bigl{[}\|\mathcal{W}\|_{-\frac{1}{\alpha}-\epsilon}^{2}\bigr{]}\leq c_{\lambda}^{-\frac{1}{\alpha}-\epsilon}\left(1+\tfrac{1}{\epsilon\alpha}\right).$$

$$\displaystyle\mathbb{E}[\|\mathcal{W}_{N}\|_{-s}^{2}]\leq c_{\lambda}^{-s}\sum_{j=1}^{N}j^{-\alpha s}\leq c_{\lambda}^{-s}\Bigl{(}1+\int_{1}^{N}x^{-\alpha s}\,\mathrm{d}x\Bigr{)}=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}$$

$$-\Delta u(\mathbf{x})+f(u(\mathbf{x}))$$

$$u(\mathbf{x})$$

$$⟨L{v}_h,{\varphi }_h⟩=(L_h{v}_h,{\varphi }_h{)}_H=({\Pi }_{h}g,{\varphi }_h{)}_H=(g,{\varphi }_h{)}_H\forall {\varphi }_h\in V_h.$$

$$(L_h{e}_{j,h},{\varphi }_h{)}_H={\lambda }_{j,h}(e_{j,h},{\varphi }_h{)}_H\forall {\varphi }_h\in V_h,1\le j\le N_h.$$

$$0<{\lambda }_{1,h}\le {\lambda }_{2,h}\le \dots \le {\lambda }_{N_h,h}.$$

$${\lambda }_j\le {\lambda }_{j,h}$$

$$\Vert e_j-e_{j,h}{\Vert }_H^2$$

$$⟨Lv,v⟩\ge \gamma \Vert v{\Vert }_{H^m(\mathcal{D})}^2\forall v\in V,$$

$$r=2(p+1-m),s=\min \{p+1,2(p+1-m)\},\text{and}q=\frac{p+1}{m}.$$

Full text

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 the theoretical results.

Acknowledgement. This work was supported in part by the Swedish Research Council (grant Nos. 2016-04187, 2017-04274), and the Knut and Alice Wallenberg Foundation (KAW 20012.0067). The authors thank Stefan A. Funken for his contribution to the FEM implementations in Matlab.

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\colon\mathcal{D}\times\mathcal{D}\to\mathbb{R}$ is given by

[TABLE]

where $\|\cdot\|$ is the Euclidean norm on $\mathbb{R}^{d}$ and $\Gamma$, $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., Banerjee et al., 2008; Furrer et al., 2006; Nychka et al., 2015; Sun et al., 2012). One of them is based on the insight that a Gaussian Matérn field on $\mathcal{D}:=\mathbb{R}^{d}$ with parameters $\sigma,\nu,\kappa>0$ can be seen as the statistically stationary solution $u$ to the stochastic partial differential equation (SPDE)

[TABLE]

Here, $\Delta$ 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}=\Gamma(2\beta-d/2)\Gamma(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}-\Delta$ 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 subject of current research, mainly because of 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; Fuglstad et al., 2015; Lindgren et al., 2011),

(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 only computable 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, it is 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 (non-orthonormal) 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 (Baeumer et al., 2015; Bonito et al., 2017; Bonito & Pasciak, 2015; Caffarelli & Silvestre, 2007; Gavrilyuk, 1996; Gavrilyuk et al., 2004, 2005; Jin et al., 2015; Nochetto et al., 2015; Roop, 2006), and for non-fractional elliptic SPDEs with random forcing (e.g., Babuska et al., 2004; Cao et al., 2007; Du & Zhang, 2003; Gyöngy & Martínez, 2006; Zhang et al., 2016). An essentially similar idea to our approach, which combines (iii) with Lanczos and Krylov subspace methods, was persued by Simpson (2008) for the differential operator $\kappa^{2}-\Delta$ in (1.2). However, neither weak nor strong convergence have been proven and the empirical results show a generally poor performance of the proposed scheme. A weak error estimate for $\kappa^{2}-\Delta$ and the non-fractional case $\beta=1$ in $d=2$ spatial dimensions has been 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 §2 we present the fractional order equation with Gaussian white noise in a Hilbert space setting, 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 §3 by partitioning the strong mean-square error in several terms and estimating these terms one by one. In addition, a weak type convergence result is obtained in Corollary 3.4. In §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 §5.

2. Model problem and main result

In the following, let $H$ denote a separable Hilbert space and $L\colon\mathscr{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:

[TABLE]

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

[TABLE]

For $\beta>0$ and $\phi\in\mathscr{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}\colon\mathscr{D}(L^{\beta})\to H$ is defined by

[TABLE]

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

[TABLE]

Its dual space after identification via the inner product on $H$ is denoted by $\dot{H}^{-2\beta}$. For $s\geq 0$, the norm on the dual space $\dot{H}^{-s}$ enjoys the useful representation

[TABLE]

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

[TABLE]

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\geq 2\beta$, since for $\phi\in\dot{H}^{s}$ it holds

[TABLE]

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}\colon\dot{H}^{s}\to\dot{H}^{s-2\beta}$.

Proof.

For $s\geq 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}$

[TABLE]

By density of $\dot{H}^{2\beta}\hookrightarrow\dot{H}^{s}$, it follows that there exists a unique linear continuous extension $L^{\beta}\colon\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 $s\geq 2\beta$ from the representation of the dual norm in (2.3) applied to $\dot{H}^{s-2\beta}$. ∎

Example 2.2**.**

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

[TABLE]

i.e., the operator $L=\kappa^{2}-\Delta$ with homogeneous Dirichlet boundary conditions on $H=L_{2}(\mathcal{D})$. For this case, we have $\dot{H}^{2}=\mathscr{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

[TABLE]

where $[\cdot,\cdot]_{s,q}$ denotes the real $K$-interpolation method (see, e.g., Thomee, 2007, Ch. 19). Furthermore, one can show (Bonito & Pasciak, 2015, Prop. 4.1) that for $1\leq s\leq 2$ it holds that

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where $\mathcal{W}$ denotes Gaussian white noise defined on a complete probability space $(\Omega,\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 following equalities involving noise terms are 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$:

[TABLE]

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}(\Omega;\dot{H}^{-s})$ for any $s>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

[TABLE]

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\geq 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

[TABLE]

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

[TABLE]

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 yield \mathbb{E}\bigl{[}\|\mathcal{W}_{N}\|_{-s}^{2}\bigr{]}=\sum_{j=1}^{N}\lambda_{j}^{-s}. Thus, we conclude with (2.1)

[TABLE]

Choosing $s=\alpha^{-1}+\epsilon$ and taking the limit $N\to\infty$ shows the L_{2}\bigl{(}\Omega;\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{(}\Omega;\dot{H}^{2\beta-\frac{1}{\alpha}-\epsilon}\bigr{)}. In particular, $u\in L_{2}(\Omega;H)$ holds if $2\alpha\beta>1$.

Remark 2.5*.*

The above results are in accordance with Lemma 4.2 and Theorem 4.3 by Zhang et al. (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\colon\mathbb{R}\to\mathbb{R}$:

[TABLE]

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\equiv\kappa^{2}\geq 0$, this corresponds to Problem (2.8) with $\beta=1$, $L=\kappa^{2}-\Delta$, 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\colon\dot{H}^{1}\to\dot{H}^{-1}$ is denoted by $L_{h}\colon 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}(\Pi_{h}g)$, where $\Pi_{h}\colon H\to V_{h}$ denotes the $H$-orthogonal projection onto the finite element space $V_{h}$, i.e.,

[TABLE]

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

[TABLE]

These eigenvalues are again arranged in nondecreasing order:

[TABLE]

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:

1. (i)

there exists $d\in\mathbb{N}$ such that $N_{h}=\dim(V_{h})\propto h^{-d}$ for all $h>0$; 2. (ii)

there exist constants $C_{1},C_{2}>0$, $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

[TABLE]

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 and assume that $L\colon\mathscr{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}$, i.e., there exists a constant $\gamma>0$ such that

[TABLE]

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, Thms. 6.1 & 6.2) for $H=L_{2}(\mathcal{D})$, $\dot{H}^{1}=V$, and the exponents

[TABLE]

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}-\Delta$ 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:

1. (i)

an expansion with respect to the discrete eigenbasis $\mathcal{E}$:

[TABLE]

where $\boldsymbol{\xi}:=\left(\xi_{1},\ldots,\xi_{N_{h}}\right)^{T}$ is the vector of the first $N_{h}$ independent standard normally distributed random variables in expansion (2.9) of $\mathcal{W}$; 2. (ii)

an expansion with respect to any basis $\Phi:=\{\phi_{j,h}\}_{j=1}^{N_{h}}$ of $V_{h}$:

[TABLE]

where the random vector $\widetilde{\boldsymbol{\xi}}:=(\widetilde{\xi_{1}},\ldots,\widetilde{\xi}_{N_{h}})^{T}$ is given by

[TABLE]

It is therefore Gaussian distributed with zero mean and covariance matrix $\mathbf{R}^{-1}(\mathbf{R}^{-1})^{T}=(\mathbf{R}^{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 $\Phi$ of the finite element space $V_{h}$.

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

Lemma 2.8**.**

The noise approximations $\mathcal{W}_{h}^{\mathcal{E}}$ and $\mathcal{W}_{h}^{\Phi}$ in (2.14)–(2.15) are equivalent in $L_{2}(\Omega;H)$, i.e., $\|\mathcal{W}_{h}^{\mathcal{E}}-\mathcal{W}_{h}^{\Phi}\|_{L_{2}(\Omega;H)}=0$.

Proof.

Inserting the definitions (2.14)–(2.15) of $\mathcal{W}_{h}^{\mathcal{E}}$ and $\mathcal{W}_{h}^{\Phi}$ yields

[TABLE]

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}\}$:

[TABLE]

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

[TABLE]

which proves the equivalence of $\mathcal{W}_{h}^{\mathcal{E}}$ and $\mathcal{W}_{h}^{\Phi}$ in $L_{2}(\Omega;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, Eq. (4)):

[TABLE]

We choose an equidistant grid $\{y_{\ell}=\ell k:\ell\in\mathbb{Z},-K^{-}\leq\ell\leq 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 method proposed by Bonito & Pasciak (2015):

[TABLE]

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

[TABLE]

on the space $\mathcal{L}(V_{h}):=\{T\colon V_{h}\to V_{h}\text{ linear}\}$ has been proven (Bonito & Pasciak, 2015, Lem. 3.4, Rem. 3.1, Thm. 3.5) for the choice

[TABLE]

Note that the quadrature (2.16) involves only non-fractional resolvents of the discrete operator $L_{h}$. Thus, the corresponding numerical method is readily implementable.

With the notions of the $H$-orthogonal projection $\Pi_{h}$ on $V_{h}$, the noise approximation $\mathcal{W}_{h}^{\Phi}$ 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

[TABLE]

Remark 2.9*.*

We emphasize the construction of the noise term $\mathcal{W}_{h}^{\Phi}$ 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 $\Pi_{h}\mathcal{W}_{N}\in L_{2}(\Omega;V_{h})$ (e.g., Zhang et al., 2016). Instead, we define $\mathcal{W}_{h}^{\Phi}$ 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:

1. (i)

Samples of the truncated white noise $\mathcal{W}_{N}$ are not needed and, thus, neither are the exact eigenvectors of the operator $L$. 2. (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

[TABLE]

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 $\Phi$ in the noise approximation $\mathcal{W}_{h}^{\Phi}$ 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 $\Phi$, 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 a Cholesky factor of $\mathbf{M}=\mathbf{G}\mathbf{G}^{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

[TABLE]

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}(\Omega;H)$ error between the solution $u$ of (2.8) and the approximation $u_{h,k}^{Q}$ in (2.18) is bounded by

[TABLE]

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 follows:

[TABLE]

and partition the strong error in (2.20) accordingly. Here, $u_{N_{h}}$, $u_{h}^{\mathcal{E}}$, and $u_{h}^{\Phi}$ 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 following equations:

[TABLE]

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}^{\Phi}$ vanishes identically in $L_{2}(\Omega;H)$:

[TABLE]

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) based on the truncated right-hand side $g_{N}+\mathcal{W}_{N}$.

Lemma 3.1**.**

Assume that the eigenvalues of the operator $L$ satisfy (2.1) and that $2\alpha\beta>1$. Let $u\in L_{2}(\Omega;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}(\Omega;\dot{H}^{2\beta})$ to the truncated equation (3.1) and it satisfies

[TABLE]

where the constant $C>0$ depends only on $\alpha,\beta$, and the constants in (2.1). In particular, under Assumption 2.6 it holds that

[TABLE]

Proof.

Existence and uniqueness of a solution $u_{N}$ in $L_{2}(\Omega;\dot{H}^{2\beta})$ 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}(\Omega;H)$ as well as the isomorphism property of $L^{\beta}\colon\dot{H}^{2\beta}\to H$ in Lemma 2.1. If $2\alpha\beta>1$, we obtain for any $N\in\mathbb{N}$:

[TABLE]

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 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}(\Omega;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}}$ has been 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}(\Omega;V_{h})$ satisfying (3.2), i.e.,

[TABLE]

and let $u_{N_{h}}\in L_{2}(\Omega;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

[TABLE]

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 $\Pi_{h}g$ can be expanded in terms of the orthonormal eigenvectors $\{e_{j,h}\}$ of $L_{h}$ by

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

for $h\in(0,h_{0})$, since since $d\alpha q\leq 2s$ by (2.19). For the second term we find

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where we used the relations $N_{h}\propto h^{-d}$ from Assumption 2.6(i) and $\alpha d(q-1)\leq 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 $\Pi_{h}g+\mathcal{W}_{h}^{\Phi}$. 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}^{\Phi}$.

Lemma 3.3**.**

Let $Q^{\beta}_{h,k}\colon V_{h}\to V_{h}$ be the operator in (2.16) approximating $L_{h}^{-\beta}$. Then it holds for sufficiently small $k\in(0,k_{0})$ and any $\phi_{h}\in V_{h}$ that

[TABLE]

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}^{\Phi}$ in (3.3) and $u_{h,k}^{Q}$ in (2.18), it holds

[TABLE]

Proof.

The first assertion is proven in (Bonito & Pasciak, 2015, Lem. 3.4, Thm. 3.5). This bound, $\|\Pi_{h}g\|_{H}\leq\|g\|_{H}$, and $\mathbb{E}[\|\mathcal{W}_{h}^{\Phi}\|_{H}^{2}]=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 §§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 follows:

[TABLE]

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

[TABLE]

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

In §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:

[TABLE]

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}(\Omega;H)$ attains the following values for $u$ in (2.8), $u_{N_{h}}$ in (3.1), and $u_{h}^{\Phi}$ in (3.3):

[TABLE]

Again, we partition the error,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Remark 3.6*.*

In the non-fractional case $\beta=1$, the discretized problem (3.3) is also non-fractional. Therefore, realizations of its solution $u_{h}^{\Phi}=L_{h}^{-1}(\Pi_{h}g+\mathcal{W}_{h}^{\Phi})$ 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 $L^{\beta}=\widetilde{L}^{\widetilde{\beta}}$. Yet, the finite element theory for the operator $\widetilde{L}$ may not be trivial.

4. An application and numerical examples

In the following numerical experiment we take up the SPDE from (1.2) in §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:

[TABLE]

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 with $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 in §3.4, 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

1. (i)

of the strong error for $d=1,2,3$ and $\beta=\tfrac{2d+n}{8}$, $n\in\{1,\ldots,7-2d\}$; 2. (ii)

of the 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 corresponding 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.

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

[TABLE]

are generated and evaluated on a uniform overkill mesh of $\bar{\mathcal{D}}$ with $N_{\mathrm{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}_{\mathrm{ok}}$ corresponds to a truncated Karhunen–Loève expansion with $N_{\mathrm{ok}}^{d}$ terms. From this, samples of the overkill solution $u_{\mathrm{ok}}$ are obtained via

[TABLE]

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}_{\mathrm{ok}}$ and use the load vector $\widetilde{\mathbf{b}}$ with entries $\widetilde{b}_{i}=(\mathcal{W}_{\mathrm{ok}},\phi_{i,h})_{L_{2}(\mathcal{D})}$ instead of $\mathbf{b}\sim\mathcal{N}(\mathbf{0},\mathbf{M})$ from Remark 2.9, as $((\mathcal{W},\phi_{i,h})_{L_{2}(\mathcal{D})})_{i=1}^{N_{h}}\sim\mathcal{N}(\mathbf{0},\mathbf{M})$ and we treat $\mathcal{W}_{\mathrm{ok}}$ as the true white noise. The resulting approximation is denoted by $\widetilde{u}_{h,k}$. We choose $N_{\mathrm{ok}}=2^{18}+1$ for $d=1$, $N_{\mathrm{ok}}=2^{12}+1$ for $d=2$, and $N_{\mathrm{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}_{\mathrm{ok}}$ to generate samples $u_{\mathrm{ok}}^{(1)},\ldots,u_{\mathrm{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

[TABLE]

The results are shown in the first three panels of Figure 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.

For (ii) we study the weak type error $|\mathbb{E}[\|u\|_{L_{2}(\mathcal{D})}^{2}]-\mathbb{E}[\|u_{h,k}^{Q}\|_{L_{2}(\mathcal{D})}^{2}]|$ addressed in Corollary 3.4. Note that for this study no sample-wise comparison is needed and, thus, the load vector is sampled from $\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}[\|u\|_{L_{2}(\mathcal{D})}^{2}]=\sum_{\boldsymbol{\theta}}\lambda_{\boldsymbol{\theta}}^{-2\beta}$. The variance of $u_{h,k}^{Q}$ is approximated via Monte Carlo integration by

[TABLE]

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 Figure 1. The theoretical rate of convergence predicted by Corollary 3.4 is 0.5 for all three cases.

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}$ is 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 in 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 §§3.1–3.3. By means of similar techniques, we have proven a weak type error estimate in Corollary 3.4.

Finally, in §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 regarded 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 Figure 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.