MDFEM: Multivariate decomposition finite element method for elliptic PDEs with lognormal diffusion coefficients using higher-order QMC and FEM

¹ Faculty of Computer Science and Engineering, Ho Chi Minh City University of Technology, VNU-HCM, Ho Chi Minh City, Vietnam
² Department of Computer Science, KU Leuven, Celestijnenlaan 200A Box 2402, Leuven B-3001, Belgium

^* Corresponding author: dirk.nuyens@kuleuven.be

Received: 17 December 2019
Accepted: 24 June 2021

Abstract

We introduce the multivariate decomposition finite element method (MDFEM) for elliptic PDEs with lognormal diffusion coefficients, that is, when the diffusion coefficient has the form a = exp(Z) where Z is a Gaussian random field defined by an infinite series expansion Z(y) = ∑_j≥1y_jϕ_j with y_j ~ N(0,1) and a given sequence of functions {ϕ_j}j≥1. We use the MDFEM to approximate the expected value of a linear functional of the solution of the PDE which is an infinite-dimensional integral over the parameter space. The proposed algorithm uses the multivariate decomposition method (MDM) to compute the infinite-dimensional integral by a decomposition into finite-dimensional integrals, which we resolve using quasi-Monte Carlo (QMC) methods, and for which we use the finite element method (FEM) to solve different instances of the PDE. We develop higher-order quasi-Monte Carlo rules for integration over the finite-dimensional Euclidean space with respect to the Gaussian distribution by use of a truncation strategy. By linear transformations of interlaced polynomial lattice rules from the unit cube to a multivariate box of the Euclidean space we achieve higher-order convergence rates for functions belonging to a class of anchored Gaussian Sobolev spaces while taking into account the truncation error. These cubature rules are then used in the MDFEM algorithm. Under appropriate conditions, the MDFEM achieves higher-order convergence rates in terms of error versus cost, i.e., to achieve an accuracy of O(ε) the computational cost is O(ε^{-1/λ-d′/λ})=O(ε^{-(p*+d′/τ)/(1-p*})) where ε^−1/λ and ε^−d′/λ) are respectively the cost of the quasi-Monte Carlo cubature and the finite element approximations, with d′= d (1+δ′) for some δ′ ≥ 0 and d the physical dimension, and 0<p*≤(2+d′/τ)^-1 is a parameter representing the sparsity of {ϕ_j}j≥1.