Weak convergence of Galerkin approximations for fractional elliptic stochastic PDEs with spatial white noise

The numerical approximation of the solution to a stochastic partial differential equation with additive spatial white noise on a bounded domain is considered. The differential operator is assumed to be a fractional power of an integer order elliptic differential operator. The solution is approximated by means of a finite element discretization in space and a quadrature approximation of an integral representation of the fractional inverse from the Dunford-Taylor calculus. For the resulting approximation, a concise analysis of the weak error is performed. Specifically, for the class of twice continuously Fr\'echet differentiable functionals with second derivatives of polynomial growth, an explicit rate of weak convergence is derived, and it is shown that the component of the convergence rate stemming from the stochasticity is doubled compared to the corresponding strong rate. Numerical experiments for different functionals validate the theoretical results.


Introduction
The representation of Gaussian random fields as solutions to stochastic partial differential equations (SPDEs) has become a popular approach in spatial statistics in recent years. It was observed already in [21] and [22] that a Gaussian random field u on R d with a covariance function of Matérn type [13] solves an SPDE of the form (κ 2 − ∆) β u = W. Here, W is Gaussian white noise, κ > 0 is a parameter determining the practical correlation range of the field, and β > d/4 controls the smoothness parameter ν of the Gaussian Matérn field via the equality ν = 2β −d/2.
Later, this relation was the incentive to consider the SPDE for Gaussian random field approximations of Matérn fields on bounded domains D R d . On the boundary ∂D, the operator κ 2 − ∆ is augmented with, e.g., homogeneous Dirichlet or Neumann boundary conditions. In [12] it was shown that by restricting the value of β to 2β ∈ N and by solving the stochastic problem (1.1) by means of a finite element method, the computational costs of many operations, which are needed for statistical inference, such as sampling and likelihood evaluations can be significantly reduced. This decrease in computing time is one of the main reasons for the popularity of the SPDE approach in spatial statistics. In addition, it facilitates various extensions of the Matérn model which are difficult to formulate using a covariance-based approach, see, for instance, [2,5,10,12,20]. However, the constraint 2β ∈ N imposed by [12] restricts the value of the smoothness parameter ν, which is the most important parameter when the model is used for prediction [17]. In [4] we showed that this restriction can be avoided by combining a finite element discretization in space with a quadrature approximation based on an integral representation of the inverse fractional power operator from the Dunford-Taylor calculus. We furthermore derived an explicit rate of convergence for the strong mean-square error of the proposed approximation for a class of fractional elliptic stochastic equations including (1.1).
In practice, it is often not only necessary to sample from the solution u to (1.1), but also to estimate the expected value E[ϕ(u)] of a certain real-valued quantity of interest ϕ(u). The aim of this work is to provide a concise analysis of the weak error |E[ϕ(u)] − E[ϕ(u Q h,k )]| for the approximation u Q h,k proposed in [4]. This analysis includes the derivation of an explicit weak convergence rate for twice continuously Fréchet differentiable real-valued functions ϕ, whose second derivatives are of polynomial growth. Functions of this form occur in many applications, e.g., when integral means of the solution with respect to a certain subdomain of D are of interest, or when a transformation of the model is used as a component in a hierarchical model. An example of the latter situation is to consider logit or probit transformed Gaussian random fields for binary regression models [16, §4.3.3].
We prove that, compared to the convergence rate of the strong error formulated in [4], the component of the weak convergence rate stemming from the stochasticity of the problem is doubled. To this end, two time-dependent stochastic processes are introduced, which at time t = 1 have the same probability distribution as the exact solution u and the approximation u Q h,k , respectively. The weak error is then bounded by introducing an associated Kolmogorov backward equation on the interval [0, 1] and applying Itô calculus.
The structure of this article is as follows: In §2 we formulate the equation of interest in a Hilbert space setting similarly to [4] and state our main result on weak convergence of the approximation in Theorem 2.4. A detailed proof of Theorem 2.4 is given in §3. For validating the theoretical result in practice, we describe the outcomes of several numerical experiments in §4. Finally, §5 concludes the article with a discussion.

Weak approximations
The subject of our investigations is the fractional order equation considered in [4], for β ∈ (0, 1), where W denotes Gaussian white noise defined on a complete probability space (Ω, A, P) with values in a separable Hilbert space H. Here and below, (in-)equalities involving random terms are meant to hold P-almost surely, if not specified otherwise. Furthermore, we use the notation X d = Y to indicate that two random variables X and Y have the same probability distribution.
Similarly to [4], we make the following assumptions: L : D(L) ⊂ H → H is a densely defined, self-adjoint, positive definite operator and has a compact inverse L −1 : H → H. In this case, −L generates an analytic strongly continuous semigroup (S(t)) t≥0 on H. The H-orthonormal eigenvectors of L are denoted by {e j } j∈N and the corresponding eigenvalues by {λ j } j∈N . These values are listed in nondecreasing order and we assume that there exist constants α, c λ , C λ > 0 such that The action of the fractional power operator L β in (2.1) is well-defined oṅ where L h denotes the Galerkin discretization of the operator L with respect to V h , i.e., We then consider the following numerical approximation of the solution u to (2.1) proposed in [4,Eq. (2.18)]. It is based on the following two components: (a) The operator Q β h,k is the quadrature approximation for L −β h of [6]: The quadrature nodes {y = k : ∈ Z, −K − ≤ ≤ K + } are equidistant with distance k > 0 and we set K − := π 2 4βk 2 and K + := j=1 is any basis of the finite element space V h . The vector ξ = (ξ 1 , . . . , ξ N h ) T is multivariate Gaussian distributed with mean zero and covariance matrix M −1 , where M denotes the mass matrix with respect to the basis Φ, i.e., The main outcome of [4] is strong convergence of the approximation u Q h,k in (2.3) to the solution u of (2.1) at an explicit rate. Subsequently, this work focusses on weak approximations based on u Q h,k , i.e., we investigate the error for continuous functions ϕ : H → R.
where z ∼ N (0, I), G is the Cholesky factor of M = GG T , and the vector g has entries g j := (g, φ j,h ) H .

Weak convergence.
For bounding the error in (2.5), we start by introducing some more notation and assumptions. Let E := {e j,h } N h j=1 ⊂ V h be the Horthonormal eigenvectors of the discrete operator L h with corresponding eigenvalues {λ j,h } N h j=1 listed in nondecreasing order. In addition, the strongly continuous semigroup on V h generated by −L h is denoted by (S h (t)) t≥0 .
Here and below, using the Riesz representation theorem, we identify the first two Fréchet derivatives Dϕ and D 2 ϕ of ϕ with functions taking values in H and in L(H), respectively. Furthermore, we say that the second derivative has polynomial growth of degree p ∈ N, if there exists a constant K > 0 such that All the properties of the finite element discretization, of the operator L, and of the function ϕ, which are of importance for our analysis of the weak error (2.5), are summarized in the assumption below.
The following example shows that Assumptions 2.2(i)-(iv) are satisfied, e.g., for the motivating problem (1.1) related to approximations of Matérn fields, if β > d/4, when using continuous piecewise linear finite element bases. Example 2.3. For κ ≥ 0 and a bounded, convex, polygonal domain D ⊂ R d , consider the stochastic model problem (1.1), i.e., the fractional order equation (2.1) for g = 0 and L = κ 2 − ∆ on H = L 2 (D). Furthermore, we assume that the differential operator L is augmented with homogeneous Dirichlet boundary conditions on ∂D. In this case, the eigenvalues {λ j } j∈N of L satisfy (2.2) for α = 2/d (see [8,Ch. VI.4] for D = (0, 1) d , the result for more general domains as above follows from the min-max principle). Consequently, the first inequality of Assumption 2.2(iii) In addition, if (V h ) h∈(0,1) ⊂Ḣ 1 = H 1 0 (D) are finite element spaces with continuous piecewise linear basis functions defined with respect to a quasi-uniform family of triangulations, Assumption 2.2(i) holds and Assumptions 2.2(ii), (iv) are satisfied for r = s = q = 2, see [18,Thm. 6 was provided, showing convergence to zero with the rate min{d(2αβ − 1), r, s}, see [4,Cor. 3.4]. In particular, the term d(2αβ − 1) stemming from the stochasticity is doubled compared to the strong rate in (2.7).
In the following, we generalize this result to weak errors of the form (2.5) for functions ϕ : H → R, which are twice continuously Fréchet differentiable and have a second derivative of polynomial growth. The bound of the weak error in Theorem 2.4 is our main result.
, and set θ = 0 otherwise. Then, for g ∈Ḣ θ and for sufficiently small h ∈ (0, h 0 ) and k ∈ (0, k 0 ), the weak error in (2.5) admits the bound Remark 2.5. In the derivation of the strong convergence rate (2.7), we balanced the error terms caused by the quadrature and by the finite element method by choosing the quadrature step size k sufficiently small with respect to the finite element mesh width h, namely e −π 2 /(2k) ∝ h dαβ , see [4, Table 1]. For calibrating the terms in the weak error estimate (2.8), we distinguish the cases αβ < 1, αβ = 1, and αβ > 1. If αβ < 1, then dαβ > d(2αβ − 1) and we let k > 0 be such that e −π 2 /(2k) ∝ h dαβ . With this choice, the error estimate (2.8) simplifies to For αβ > 1 (αβ = 1) we achieve the same bound if k and h are calibrated such that Note that the calibration for αβ < 1 coincides with the one for the strong error and that the term d(2αβ − 1) in the derived weak convergence rate min{d(2αβ − 1), r, s} is doubled compared to the first term of the strong convergence rate (2.7).
Remark 2.6. We emphasize that (under the same assumptions) both the strong and weak convergence rates remain the same when approximating the solution u to where σ > 0 is a constant parameter which scales the variance of u. This can be seen from the equality σ −1 L β = L β σ for L σ := σ −1/β L, combined with the fact that the eigenvalues of the operator L σ satisfy the growth assumption (2.2) with the same exponent α > 0 as the eigenvalues of L.
However, the constants c λ , C λ > 0 in (2.2) and the constants in the error estimates change. For instance, if ϕ(u) := u p * H for p * ∈ N, then the constant C > 0 in (2.8) will depend linearly on σ p * .
Note that one has to consider a problem of the form when approximating a Matérn field with variance σ 2 * .
Here and in what follows, Γ( · ) denotes the Gamma function.
Remark 2.7. We also comment on how the error bound in (2.8) is used. If there exists a function E : N → R ≥0 such that lim n→∞ E(n) = 0 as well as a constant C > 0, independent of h and n, such that is an immediate consequence of the arguments in our proof that a bound of the weak error for the approximation u R h,n : An example of such a family {R β h,n } n∈N are the approximations of L −β h proposed in [3], which are based on rational approximations of the function x −β of different degrees n ∈ N.

The derivation of Theorem 2.4
The main idea in our derivation of the weak error estimate (2.8) is to introduce two time-dependent stochastic processes with the property that their (random) values at time t = 1 have the same distribution as the solution u to (2.1) and its approximation u Q h,k in (2.3), respectively. We then use an associated Kolmogorov backward equation and Itô calculus to estimate the difference between these values.
3.1. The extension to time-dependent processes. Recall the eigenvalue-eigenvector pairs {(λ j , e j )} j∈N of L as well as the parameter α > 0 determining the growth of the eigenvalues via (2.2). In what follows, we assume that g ∈ H and 2αβ > 1 so that the solution u to (2.1) satisfies u ∈ L 2 (Ω; H). With the aim of introducing the time-dependent processes mentioned above, we start by defining where {B j } j∈N is a sequence of independent real-valued Brownian motions adapted to a filtration F := (F t , t ≥ 0). Owing to this construction, (W β (t), t ≥ 0) is an F-adapted H-valued Wiener process with covariance operator L −2β , which is of trace-class if 2αβ > 1. Since the random variables {B j (1)} j∈N are independent and identically N (0, 1)-distributed, the spatial white noise W satisfies The stochastic process Y := (Y (t), t ∈ [0, 1]) defined as the (strong) solution to the stochastic partial differential equation therefore takes the following random value in H at time t = 1, Its Gaussian distribution implies the existence of all moments, as shown in the following lemma. 1], and Y be the strong solution of (3.1). Then the p-th moment of Y (t) exists and, for p ≥ 2, it admits the following bound: Proof. For p = 2, the bound in (3.3) follows from the Itô isometry [15, Thm. 8.7(i)]: In order to define a another stochastic process Y : We now consider the following stochastic partial differential equation Note that the reproducing kernel Hilbert space of W β isḢ 2β . The finite rank of the operator Q β h,k P β h : H → V h implies that it is a Hilbert-Schmidt operator fromḢ 2β to H. For this reason, existence and uniqueness of a (strong) solution Y to (3.5) is evident. Furthermore, the solution process Y satisfies define the deterministic matrix R and the random vector B 1 by i.e., B 1 is the vector of the first N h Brownian motions at time t = 1. Due to are equal in L 2 (Ω; H). In particular, their first and second moments coincide.
Since W E h and W Φ h are Gaussian random variables, their distributions are uniquely characterized by their first two moments and we conclude that (3.8) Since ϕ : H → R is twice continuously Fréchet differentiable, we can furthermore express the first two derivatives of w with respect to x in terms of ϕ and Y by Let Y be the solution to (3.5). The application of Itô's lemma [7] to the stochastic process (w(t, Y (t)), t ∈ [0, 1]) yields where, for T ∈ L(H), the H-adjoint operator is denoted by T * . To simplify the second term in (3.11), we define the operator Π h : H → V h by Note that in contrast to the H-orthogonal projection Π h , the operator Π h is neither self-adjoint ( Π * h = Π h ) nor a projection ( Π 2 h = Π h ). We then use the following relation between Π h and P β h from (3.4), and express (3.11) as an integral equation for t = 1. Taking the expectation on both sides of this equation yields . As a final step in this subsection, we relate the quantity of interest E[ϕ(u)] with the expected value of w(1, Y (1)) and similarly for the approximation E[ϕ(u Q h,k )] and w(1, Y (1)). For this purpose, we extend the equalities in (3.8)-(3.10) to the case that x = ξ is a an H-valued random variable in the following lemma.
Owing to Lemma 3.2 and the tower property for conditional expectations, the stochastic process (w(t, Y (t)), t ∈ [0, 1]) has no drift, i.e., (3.14) Furthermore, it follows with (3.2) and (3.6) that Summing up the observations in (3.13)-(3.16), we find that the difference between the quantity of interest E[ϕ(u)] and the expected value of the approximation ϕ(u Q h,k ) can be expressed by This equality implies that the weak error (2.5) admits the following upper bound where we set Q β h,k := Q β h,k Π h and L −β h := L −β h Π h . The following subsections are structured as follows: In §3.3 we bound the deterministic error (L −β − L −β h Π h )g H caused by the finite element discretization. This result is essential for estimating the first error term (I) in (3.17). Secondly, we investigate the terms (II) and (III) stemming from applying the quadrature operator Q β h,k instead of the discrete fractional inverse L −β h in §3.4. Finally, in §3.5 we estimate the trace in (IV) and combine all our results to prove Theorem 2.4.

The deterministic finite element error.
In this subsection we focus on the deterministic error (L −β − L −β h Π h )g H caused by the inhomogeneity g. More precisely, we derive an explicit rate of convergence depending on theḢ θ -regularity of g in Lemma 3.3 below. Subsequently, in Lemma 3.5, we apply this result to bound the first term of (3.17). for all g ∈Ḣ θ and sufficiently small h ∈ (0, h 0 ).
Remark 3.4. We note that by letting σ 1 = σ 2 := s, θ 1 := s − 2β + , and θ 2 := 0 in the proof of Lemma 3.3 the optimal convergence rate for the deterministic error, , can be derived. The error estimate (3.18) is formulated in such a way that the smoothness of g ∈Ḣ θ is minimal for convergence with the rate min{d(2αβ − 1), s}, which will dominate the overall weak error, stemming from the term (IV) in the partition (3.17), see Lemma 3.9.
Having bounded the error between L −β g and L −β h Π h g, an estimate of the first error term (I) in (3.17) is an immediate consequence of the fundamental theorem of calculus and the chain rule for Fréchet derivatives. This bound is formulated in the next lemma.
Proof. Since the mapping x → w(0, x) is Fréchet differentiable, we obtain by the fundamental theorem of calculus and the Cauchy-Schwarz inequality A bound for the first term is given by (3.18) in Lemma 3.3. For the second term, we use (3.9), Y (0) = L −β g, and the polynomial growth (2.6) of D 2 ϕ to estimate where the constants C, C > 0 depend only on β and the smallest eigenvalue of L.
In the following, we use this error estimate of the quadrature approximation Q β h,k for bounding the second term of (3.17) in Lemma 3.7 as well as the trace occurring in the third term of (3.17) in Lemma 3.8.
Lemma 3.7. Suppose that Assumption 2.2(v) is satisfied and that 2αβ > 1. Then there exists a constant C > 0, independent of h and k, such that for all g ∈ H and sufficiently small h ∈ (0, h 0 ) and k ∈ (0, k 0 ).
Proof. As in the proof of Lemma 3.5, we apply the fundamental theorem of calculus and the chain rule for Fréchet derivatives. By (3.9) and Lemma 3.6 we then find . Again, the proof is completed by (3.3) and the fact that tr(L −2β ) < ∞.
Here, the function f α,β is defined as in Theorem 2.4.
Proof. By the definition of Π h in (3.12) we have Therefore, the trace of interest simplifies to a finite sum, where the second equality follows from the self-adjointness of T ∈ L(H).
The application of the Cauchy-Schwarz inequality and of Lemma 3.6 to the first sum yield the following upper bound By Assumption 2.2(i) we thus have |S 1 | e − π 2 k h −d T L(H) . The second sum can be bounded by Finally, due to the approximation property of the discrete eigenvalues λ j,h in Assumption 2.2(ii) as well as the growth (2.2) of the exact eigenvalues λ j we obtain λ −β j,h ≤ λ −β j ≤ c −β λ j −αβ and, for αβ = 1, we find where we have used Lemma 3.6 and Assumption 2.2(i). If αβ = 1, we instead estimate |S 2 | e −π 2 /(2k) (1 + | ln(h)|) T L(H) . This completes the proof.
Proof. Similarly to (3.20) we use the self-adjointness of T and rewrite the trace as In order to estimate the terms S 1 and S 2 , we note that for j ∈ {1, . . . , N h } By the mean value theorem, the existence of λ j ∈ (λ j , λ j,h ) satisfying Owing to (3.19) the series S 1 simplifies to the finite sum Using (3.21) as well as Assumptions 2.2(i)-(iii), this sum can be bounded by since dα(q − 1) ≤ r and dαq/2 ≤ s by Assumption 2.2(iii). For the second term we find since L −β h e j = 0 for j > N h by (3.19). Therefore, the application of (3.21) yields and |S 2 | h min{d(2αβ−1),r,s} T L(H) follows from Assumptions 2.2(i), (iii).
Lemma 3.10. Suppose that Assumptions 2.2(i)-(iii) are satisfied. Let p ∈ N, t ∈ [0, 1], and Y be the strong solution of (3.5). Then the p-th moment of Y (t) exists and, for p ≥ 2, it admits the following bound: where the constant C > 0 is independent of h and k.
, we obtain by Lemma 3.6, for any p ≥ 2, that where, again, µ p := E[|Z| p ] denotes the p-th absolute moment of Z ∼ N (0, 1) and the constant C > 0 is independent of h, k, and p. Furthermore, using 0 < λ j ≤ λ j,h of Assumption 2.2(ii) and applying the Hölder inequality gives where tr(L −2β ) < ∞ by Assumption 2.2(iii). Thus, we obtain for the solution Y of (3.5) that for any t ∈ [0, 1] the bound holds. Finally, the assertion follows by the boundedness of Q β h,k which is uniform in h and k, the finiteness of tr(L −2β ), and Assumption 2.2(i).
Proof of Theorem 2.4. Owing to the partition (3.17) and the estimates of the error terms (I)-(IV) in Lemmata 3.5 and 3.7-3.9 we can bound the weak error as follows since w xx (t, x) ∈ L(H) is self-adjoint for every t ∈ [0, 1] and x ∈ H. The application of Lemma 3.2 and of the tower property for conditional expectations yield By the polynomial growth (2.6) of D 2 ϕ and the boundedness of the p-th moments of Y (t) and Y (t) in Lemmata 3.1 and 3.10, respectively, we obtain that since tr(L −2β ) < ∞. This completes the proof of the weak error estimate (2.8).
Remark 3.11. Note that, if the first and second Fréchet derivatives of ϕ are bounded, the estimates of the Lemmata 3.1 and 3.10 are not needed and the weak error estimate in (2.8) simplifies to The calibration of the discretization parameters k and h remains as described in Remark 2.5.

An application and numerical experiments
In this section we validate the theoretical results of the previous sections within the scope of a simulation study based on the model for Matérn approximations in (1.1) on the domain D = (0, 1) d for d = 1, 2, κ = 0.5, and u = 0 on ∂D, i.e., L = κ 2 − ∆ with homogeneous Dirichlet boundary conditions. In this case, the operator L has the following eigenvalue-eigenvector pairs [8,Ch. VI.4]: where j = (j 1 , . . . , j d ) ∈ N d is a d-dimensional multi-index. As already mentioned in Example 2.3, these eigenvalues satisfy (2.2) for α = 2/d. Note that for every x ∈ D the solution u satisfies u(x) ∼ N (0, σ(x) 2 ). Following a Karhunen-Loève expansion of u with respect to the eigenfunctions {e j } j∈N d in (4.1), the variance σ(x) 2 can be expressed explicitly in terms of the eigenvalues and eigenfunctions in (4.1) by where ξ j j∈N d are independent N (0, 1)-distributed random variables. Considering continuous evaluation functions ϕ : L 2 (D) → R of the form and the value of E[f (u(x))] can be derived analytically from u(x) ∼ N (0, σ(x) 2 ). More precisely, we choose f (u) = |u| p , p = 2, 3, 4, as well as f (u) = Φ(20(u − 0.5)), where Φ(·) denotes the cumulative distribution function for the standard normal distribution. The motivation of the latter function is given by its correspondence to a probit transform which is often used to approximate step functions (see, e.g., [1]), in this case 1(u > 0.5). These four functions satisfy Assumption 2.2(v) and we obtain for the quantity of interest, if f (u) = |u| p , and if f (u) = Φ(c(u − a)) for a ∈ R and c > 0. We truncate the series in (4.2) in order to approximate the variance σ(x) 2 , Here, we choose N ok = 1 + 2 18 for d = 1 and N ok = 1 + 2 11 for d = 2 so that, in both cases, N d ok N h for all considered finite element spaces with N h basis functions. This estimate of σ(x) is used at N d ok equally spaced locations x ∈ D, and the reference solution E[ϕ(u)] is then approximated by applying the trapezoidal rule in order to evaluate the integrals in (4.3) and (4.4) numerically.
We consider (1.1) for β = 0.6, 0.7, 0.8, 0.9 and use a finite element discretization based on continuous piecewise linear basis functions with respect to uniform meshes onD = [0, 1] d . We use four different mesh sizes h in each dimension d = 1, 2, and calibrate the quadrature step size k with h for each value of β by k = −1/(β ln h). This results in the numbers of basis functions and quadrature nodes shown in Table 1. As already pointed out in Example 2.3, the growth exponent of the    Figure 1. For each function ϕ and for each value of β, we compute the empirical convergence rate r by a least-squares fit of a line c + r ln h to the data set {h , err }. The results are shown in Table 2 and can be seen to validate the theoretical rates given in Theorem 2.4 for d = 1. For d = 2, the observed rates deviate slightly from the theoretical rates for β = 0.9, which is caused by the fact that we had to use coarser finite element meshes for d = 2 than for d = 1 in order to be able to assemble the dense matrices Q β h,k ∈ R N h ×N h for performing the simulation study without Monte Carlo simulations.   Table 2.

Conclusion
Gaussian random fields are of great importance as models in spatial statistics. A popular method for reducing the computational cost for operations, which are needed during statistical inference, is to represent the Gaussian field as a solution to an SPDE. In this work, we have investigated a recent extension of this approach to Gaussian random fields with general smoothness proposed in [4]. The method considers the fractional order equation (2.1) and is based on combining a finite element discretization in space with the quadrature approximation (2.4) of the inverse fractional power operator. This yields an approximate solution u Q h,k of the SPDE, which in [4] was shown to converge to the solution u of (2.1) in the strong mean-square sense with rate (2.7).
In many applications one is mostly interested in a certain quantity of the random field u which can be expressed by ϕ(u) for some real-valued function ϕ. For this reason, the focus of the present work has been the weak error |E[ϕ(u)]−E[ϕ(u Q h,k )]|. The main outcome of this article, Theorem 2.4, shows convergence of this type of error to zero at an explicit rate for twice continuously Fréchet differentiable functions ϕ, which have a second derivative of polynomial growth. Notably, the component of the convergence rate stemming from the stochasticity of the problem is doubled compared to the strong convergence rate (2.7) derived in [4]. For proving this result, we have performed a rigorous error analysis in §3, which is based on an extension of the equation (2.1) to a time-dependent problem as well as an associated Kolmogorov backward equation and Itô calculus.
In order to validate the theoretical findings, we have performed a simulation study for the stochastic model problem (1.1) on the domain D = (0, 1) d for d = 1, 2 in §4. This model is highly relevant for applications in spatial statistics, since it is often used to approximate Gaussian Matérn fields. We have considered four different functions ϕ and the fractional orders β = 0.6, 0.7, 0.8, 0.9. The observed empirical weak convergence rates can be seen to verify the theoretical results. One of the considered functions ϕ is based on a transformation of the random field by a Gaussian cumulative distribution function. Quantities of this form are particularly important for applications to porous materials, as they are used to model the pore volume fraction of the material, see, e.g., [1]. Thus, we see ample possibilities for applying the outcomes of this work to problems in spatial statistics and related disciplines.