Numerical methods for fractional SPDEs with applications to spatial statistics

The numerical approximation of solutions to stochastic partial differential equations (SPDEs) with additive spatial white noise on bounded domains is considered. The differential operator is given by the fractional power of an integer order elliptic differential operator and is therefore non-local. Its inverse operator 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 is approximated by a weighted sum of non-fractional resolvents at certain quadrature nodes. These resolvents are then discretized in space by a standard finite element method. By combining this approach with approximate realizations of the white noise, which are based only on the mass matrix of the finite element discretization, an efficient numerical algorithm for computing samples of the approximate solution is obtained. For the resulting approximation, explicit rates of strong and weak convergence are derived. A key property of the presented scheme is that it does not require the knowledge of the eigenfunctions of the differential operator, which is necessary, e.g., for approximations based on truncated spectral Karhunen-Loève expansions of the noise term. For this reason, the method is particularly interesting for real-world applications in spatial statistics, such as to employ solutions to fractional order SPDEs as approximations of Gaussian Matérn fields. This application is taken up in numerical experiments to illustrate and attest the theoretical results.