Kristin Kirchner

David Bolin, Kristin Kirchner, Mihály Kovács

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^β$, $β\in(0,1)$, of an integer order elliptic differential operator $L$ and is therefore non-local. Its inverse $L^{-β}$ 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^{-β}$ 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=κ^2-Δ$, $κ> 0$, with homogeneous Dirichlet boundary conditions on the unit cube $(0,1)^d$ in $d=1,2,3$ spatial dimensions for varying $β\in(0,1)$ attest the theoretical results.

David Bolin, Kristin Kirchner, Mihály Kovács

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échet 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.

Kristin Kirchner, Fabio Nobile, Christoph Schwab et al.

We present a theoretical and numerical analysis of Monte Carlo methods for the estimation of statistical moments of random variables $X:Ω\rightarrow E$ taking values in a Banach space $E$. For practical computation, we consider finite-dimensional approximation subspaces ${(E_\ell)_{\ell\in\mathbb{N}}\subset E}$ of increasing dimension. We develop a refined error analysis that explicitly accounts for a dependence of the Rademacher type constants on the dimension of $E_\ell$, leading to novel complexity results for single- and multilevel Monte Carlo methods to estimate the mean and injective moments of arbitrary order, which are, in certain cases, sharper than those derived in [Kirchner, Schwab, J. Funct. Anal, 2024]. Moreover, we show that, in favorable cases, the resulting error-vs.-work bounds are independent of the Rademacher type of $E$. We then focus on $L^p(S)$-valued random variables for a $σ$-finite measure space satisfying certain approximation properties, and prove that for a random variable $X\in L^q(Ω;L^p(S))\cap L^p(S;L^q(Ω))$, with $q\in (1,\infty)$ and $p\in [1,\infty)$, the $L^q$-convergence rate of a Monte Carlo estimator is determined exclusively by the integrability parameter $\min\{q,2\}$, with no dependence on the Rademacher type $\min\{p,2\}$ of $L^p(S)$. We further investigate the impact of measuring the (multilevel) Monte Carlo error in the $L^q(Ω;L^p(S))$-norm while $X$ possesses additional regularity, $X\in L^{\tilde{q}}(Ω;L^p(S))\cap L^p(S;L^{\tilde{q}}(Ω))$ with $\tilde{q}\in [q,\infty)$. This analysis reveals an interplay between the sampling error and the strong approximation error, and leads to optimized error-vs.-work bounds for both single- and multilevel Monte Carlo methods. Numerical experiments confirm the sharpness of the analyses presented.

Kristin Kirchner

Numerical methods for stochastic partial differential equations typically estimate moments of the solution from sampled paths. Instead, we shall directly target the deterministic equations satisfied by the first and second moments, as well as the covariance. In the first part, we focus on stochastic ordinary differential equations. For the canonical examples with additive noise (Ornstein-Uhlenbeck process) or multiplicative noise (geometric Brownian motion) we derive these deterministic equations in variational form and discuss their well-posedness in detail. Notably, the second moment equation in the multiplicative case is naturally posed on projective-injective tensor product spaces as trial-test spaces. We construct Petrov-Galerkin discretizations based on tensor product piecewise polynomials and analyze their stability and convergence in these natural norms. In the second part, we proceed with parabolic stochastic partial differential equations with affine multiplicative noise. We prove well-posedness of the deterministic variational problem for the second moment, improving an earlier result. We then propose conforming space-time Petrov-Galerkin discretizations, which we show to be stable and quasi-optimal. In both parts, the outcomes are illustrated by numerical examples.

Sonja G. Cox, Kristin Kirchner

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