In this paper we are concerned with the application of higher-order quasi-Monte Carlo (QMC) rules and multivariate decomposition methods (MDM) to elliptic PDEs with random diffusion coefficients. We focus on the lognormal diffusion coefficient, the logarithm of which is a Gaussian random field. The goal is to compute the expected value of some functional of the solution. This method was motivated by the need for new techniques for elliptic PDEs with smooth lognormal diffusion coefficients where the classical approaches fail.
A theoretical analysis of (higher-order) QMC rules and the MDM applied to this type of model problem but with uniform diffusion was recently introduced in . By using suitable higher-order QMC rules and FE approximations the MDM takes less computational cost to achieve a comparable approximation than the standard QMCFEM, see, e.g., . For lognormal diffusions we cope with a more challenging problem where the expectation or corresponding integration is taken with respect to the Gaussian distribution over an unbounded domain of the Euclidean space, where existing higher-order QMC algorithms are not directly applied. To solve this problem we propose to use a truncation method recently developed in , see also . Exploiting the fast decay of the Gaussian distribution toward infinity the Euclidean domain is truncated, the resulting integral is transformed to the unit cube using a linear transformation, and finally, suitable higher-order QMC rules are applied. The proposed algorithm allows us to achieve arbitrarily higher-order convergence for sufficiently smooth integrands.
Let be a bounded polygon domain in , with typically or , with boundary . We consider the following elliptic Dirichlet problem
|for in ,||(1)|
Here, the gradient operator is taken with respect to and with some .
We consider the case when is a sequence of parameters distributed on according to the product Gaussian measure , and the diffusion coefficient takes the form
where is a suitable system of real-valued, bounded, and measurable functions.
Let be a linear and bounded functional of the solution . We are interested in computing the expected value of
with respect to the probability distribution, i.e.,
The weak (or variational) formulation of problem (1) is to find for a given the solution such that
Under some assumptions on the system we will show the existence, the uniqueness and an a priori estimation of the solution of the weak formulation. The following assumptions are standard, see, e.g., [1, 2, 16, 19]: assume that there exists a positive sequence with for all such that the positive constant , given below, satisfies
Let us define for some space and
The following result is implied from [1, Theorem 2.1 and Remark 2.2].
we define two random variables
Under the assumption of Proposition 1 the bilinear form is continuous and coercive for -a.s. with constants and , respectively. Indeed, using Proposition 1 and by the definition of the parameter we have for -a.s.
Applying Proposition 1 leads to
For any and belonging to using Hölder’s inequality we have
Taking in (5) yields
The Lax–Milgram lemma implies that there is a unique solution of the weak form (5), moreover,
for any and .
To compute (4) we cope with three computational challenges: first, the infinite number of variables of the integrand; second, the integration domain is unbounded whereas most existing quasi-Monte Carlo methods integrate only over unit cubes; and, third, the integrand involves solutions of PDEs.
Let us discuss how to solve these problems in more detail. First, to approximate an infinite-dimensional integral we will exploit the multivariate decomposition method, which is developed based on the earlier changing dimension method, see, e.g., [24, 30, 21]. The goal is to decompose the infinite-dimensional problem into multiple finite-dimensional ones. By assessing the relative contribution of specific sets of variables, for a given desired error, the MDM will decide which ones to include in a so-called active set to approximate the infinite MDM sum. We argue, see Proposition 8, that provided certain conditions on are satisfied, sets in the active set have relatively low cardinalities, such that only relatively low-dimensional problems remain which can be solved at small cost. This is the reason why the MDM algorithm improves the complexity of computation.
Next, to compute integrals over the Euclidean space, we exploit the fast decay of the Gaussian distribution to truncate the unbounded domain into boxes. We then use a linear transformation to map the truncated integral into the unit cube. Finally, we apply existing higher-order quasi-Monte Carlo rules, see, e.g., [7, 27]
. This is in contrast to existing methods which use the inverse of the Gaussian cumulative distribution function to map the integral into the unit cube and then apply particular QMC rules, see, e.g.,[23, 28]. Using such nonlinear mappings might make the transformed integrand singular in the sense that their mixed derivatives might not exist or are unbounded. As a remedy, special function spaces and QMC rules were developed, however, only achieved first order convergence rates which are independent of the dimension of the integrand. The proposed QMC method in this paper by using the linear mapping avoids damaging the smoothness of the integrand. As a result, it allows to achieve arbitrarily higher-order convergence rates for sufficiently smooth integrands. Here, the constant might depend exponentially of the dimension, see Theorem 1, however, since the active set of the MDM algorithm consists of functions of few variables this exponential growth will be controlled, see Proposition 8.
Last, for each variable sampled by the QMC method the original stochastic PDE becomes a deterministic one. To solve such problem we use the finite element method (FEM). The combination of the multivariate decomposition method with the finite element method is called the multivariate decomposition finite element method or MDFEM in short.
The structure of this paper is as follows. Section 2 presents main steps of the MDFEM algorithm. Section 3 introduces higher-order QMC rules for multivariate integrals over the Euclidean space with respect to the Gaussian distributions. A novel anchored Gaussian Sobolev function space is introduced and QMC rules are developed for that specific space. Section 4 recalls the definition of interlaced polynomial lattice rules, studies and analyzes their error. Section 5 discusses the parametric regularity of the solution of the PDE. Section 6 previews the finite element method. Under some conditions on the system a novel result on the spatial Hölder smoothness of the solution of the PDE is derived. Section 7 presents the main contribution of this paper where the cost model, the construction and the complexity of the MDFEM algorithm are presented. A comparison of the MDFEM and the standard QMCFEM is also given which shows the out-performance of the MDFEM.
We will use some notations. We write and . For any we denote by the set of indices . For any let denote the Hölder space of all functions that are times continuously differentiable on with norm
for with . For any real such that we set and define the following norm
where is the Euclidean norm. Both the cardinality of a set and the
norm of a vector are denoted bybut it should be clear from the context whichever is meant.
2 Application of the MDM to PDEs: MDFEM
In this section we introduce the main steps of MDFEM algorithm. Let us first recall some definitions and notations from . For any , with , we denote by the vector such that for and otherwise, and by the -projected solution of (1) with , that is, the solution of the problem:
where . The weak formulation of the -projected problem involves finding for a given such that the following equation holds
An exact analytical solution of this problem is typically not available so a numerical approximation will be computed using the FEM. Let us define a finite dimensional subspace , where the is to be specified below, and such that for . The finite element approximation of the weak formulation of the -projected problem denoted by involves finding for a given such that the following equation holds
The MDFEM algorithm relies on the multivariate decomposition, see, e.g., , of the solution
where the sum is over all finite subsets , and we take the anchored decomposition with anchor at to obtain the components
For any multi-index let us denote with and by the value of such partial derivative of the function at . For any the function satisfies the following properties. The proof of this result is provided in the appendix
For any the function depends only on the variables listed in and satisfies
|when any component of equals .||(13)|
Moreover, if the solution has derivatives of arbitrarily high order with respect to the variable , see Proposition 4 below, then
Since the functional is linear we write
The MDFEM algorithm to approximate (2) takes the form
where are cubature rules with being their cubature nodes and respective weights, and
We remind that is the FE approximation with mesh width of the solution and emphasize here that to approximate for all we use the same .
The computational cost of the MDFEM algorithm is given as
There are three sources of error in the approximation (18): the error comes from truncating the infinite sum, the QMC error and the FEM error. They are gathered into two terms as follows
A sufficient condition to achieve an approximation error of is that both of these terms are less than . This forces us to construct the active set such that the first term is bounded by . For all the FE space and the cubature rule are chosen such that the second term is bounded by with optimal computational cost (20).
3 Higher-order quasi-Monte Carlo rules for finite dimensional integration with respect to the Gaussian distribution
In this section we consider quasi-Monte Carlo rules for approximating integrals over with respect to the Gaussian distribution. Particularly, we are interested in computing -dimensional integrals of the form
where with an abuse of notation we omit the subscript referring to the dimension of and write
We first truncate the Euclidean domain into a multidimensional box, then use a linear mapping to transform the truncated integral into one over the unit cube, which is finally approximated using suitable cubature rules. More precisely, the truncated and transformed integrals respectively have the following forms
for some and the mapping is defined as
The final integral is approximated using an -point QMC rule of the form
where are well chosen cubature points.
In the present for simplifying the notation we will drop the index in the derivative and write for any , the value of such derivative at is denoted by . The function is general and depends on variables, however, in further applications will be of the form for some such that . As discussed in the previous section the function have the same properties (13), (14) and (15) as . Currently only the properties (13) and (15) are needed, the property (14) will be used in further parts. As a result, we now only consider the integrand such that if any component of is equal to , and for any such that being a proper subset of . Here, is the smoothness parameter.
3.1 Reproducing kernel Hilbert space
We begin with introducing the one dimensional function space. The space consists of integrable functions over with respect to the Gaussian distribution having absolutely continuous derivatives up to order and square integrable derivative of order over with respect to the Gaussian distribution. The inner product is defined as
This space is called the anchored Gaussian Sobolev space with anchor at . The associated norm is given by which is actually a norm due to the fact that .
There exists yet another Gaussian Sobolev space whose norm as well involves derivatives, see, e.g., [18, 17, 7]. However, instead of anchoring the values of the function and its derivatives up to order at as in (23) that space takes their averages over against the Gaussian distribution. We refer to such space as the unanchored Gaussian Sobolev space. Since functions in there can be represented using Hermite polynomials the unanchored Gaussian Sobolev space is also called the Hermite function space. For more details on this space we refer to [18, 17, 7].
In this paper it suffices to consider the anchored Gaussian Sobolev space . This is also a reproducing kernel Hilbert space, the reader will find that this property is particularly useful for the error analysis of the MDFEM algorithm. The kernel is given by
where is the indicator function on the set , and , and the empty sum is defined as . Such formulas for with was given in [28, Section 3.3]. A slightly different kernel for another anchored Sobolev space over the Euclidean space of higher order smoothness, although, without the specific Gaussian distribution was given in [29, Section 11.5.1]. For completeness we provide the full derivation of this kernel in the appendix.
According to the setting of the MDM, see , we need to show that the square root of is measurable. Indeed, we have for any
The same bound holds for any . Thus, using we have
where is the Gamma function.
The multivariate space is defined as the tensor product of the one dimensional spaces with the kernel given by
The corresponding inner product is
where is a vector of variables such that for and otherwise. The corresponding norm is given by .
We will also need a function space over the unit cube. Let us define the unanchored Sobolev space over the unit cube . For the one dimensional case its inner product is given by
with norm .
The multivariate space is the tensor product of the one dimensional space with inner product given by
and norm .
The following result states the relation between norms in and in which is then needed for the QMC error analysis. The proof of this result is presented in the appendix.
For any and , the function belongs to and
and is the modified Bessel function of the first kind with .
3.2 Higher-order quasi-Monte Carlo for integration over
This subsection investigates the cubature error of approximating integrals over with respect to the Gaussian distribution using the truncation strategy and higher-order quasi-Monte Carlo rules. The main result will be given in Theorem 1. A similar analysis, however, for the unanchored Gaussian Sobolev funtion space was given in [7, Theorem 2, Corollary 1] where infinite-dimensional interlaced digital sequences with integer rates of convergence were used. Nevertheless, because the MDM algorithm needs cubature rules with possibly non-integer convergence rates, we use interlaced polynomial lattice rules. We will provide their error analysis in the next subsection. But first we state our result.
For any with such that and any , let be an interlaced polynomial lattice rule of order with points achieving the convergence rate of order as in Theorem 2 below. The cubature rule defined as in (22) with and cubature point set has error bounded by
where is a constant independent of and defined below by (35).
The error splits into two terms
The first term is the domain truncation error. Using the reproducing property of and the Cauchy–Schwarz inequality we have
The integral of the right hand side is bounded as
We will estimate each of the above integrals. Using the same argument as in (3.1) we have for any