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

We introduce the novel multivariate decomposition finite element method (MDFEM) for solving elliptic PDEs with lognormal diffusion coefficients, that is, when the diffusion coefficient a has the form a=(Z) where Z is a Gaussian random field defined as Z(y) = ∑_j ≥ 1 y_j ϕ_j with y_j ∼N(0,1) and a sequence of functions {ϕ_j}_j ≥ 1. We estimate the expected value of some linear functional of the solution as an infinite-dimensional integral over the parameter space. The proposed algorithm combines the multivariate decomposition method (MDM), to compute infinite-dimensional integrals, with the finite element method (FEM), to solve different instances of the PDE. This allows us to apply higher-order multivariate quadrature methods for integration over the Euclidean space with respect to the Gaussian distribution, and, hence, considerably improves upon existing results which only use first order cubature rules. We develop multivariate quadrature methods for integration over the finite-dimensional Euclidean space with respect to the Gaussian distribution. By linear transformations of interlaced polynomial lattice rules from a 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. These cubature rules are then used in the MDFEM algorithm. Under appropriate conditions, applying higher-order cubature rules we achieve higher-order convergence rates in term of error vs cost, i.e., the computational cost to achieve an accuracy of O(ϵ) is O(ϵ^-(1+d/τ)p^*/(1-p^*) -d/τ) where d is the physical dimension, τ is the convergence rate of the finite element method and p^* represents the sparsity of {ϕ_j}_j ≥ 1.

Authors

12/13/2021

Goal-oriented adaptive finite element method for semilinear elliptic PDEs

We formulate and analyze a goal-oriented adaptive finite element method ...
08/28/2021

Scaled lattice rules for integration on ℝ^d achieving higher-order convergence with error analysis in terms of orthogonal projections onto periodic spaces

We show that to approximate the integral ∫_ℝ^d f(x) dx one can simply us...
12/18/2020

A hybrid MGA-MSGD ANN training approach for approximate solution of linear elliptic PDEs

We introduce a hybrid "Modified Genetic Algorithm-Multilevel Stochastic ...
09/30/2019

H^1-norm error estimate for a nonstandard finite element approximation of second-order linear elliptic PDEs in non-divergence form

This paper establishes the optimal H^1-norm error estimate for a nonstan...
11/02/2020

The effect of quadrature rules on finite element solutions of Maxwell variational problems. Consistency estimates on meshes with straight and curved elements

We study the effects of numerical quadrature rules on error convergence ...
12/01/2019

A Multigrid Method for a Nitsche-based Extended Finite Element Method

We present a tailored multigrid method for linear problems stemming from...
03/25/2021

Near-optimal approximation methods for elliptic PDEs with lognormal coefficients

This paper studies numerical methods for the approximation of elliptic P...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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

 −∇⋅(a(x,y)∇u(x,y)) =f(x), for x in D, (1) u(x,y) =0, for x on ∂D.

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

 a(x,y):=exp(Z(x,y)), (2)

with

 Z(x,y):=∑j≥1yjϕj(x),yj∈Ω=R,yj∼N(0,1), (3)

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.,

 E[G(u)]:=I(G(u)):=∫RNG(u(⋅,y))dμ(y), (4)

where

 dμ(y):=∏j≥1ρ(yj)dyj,ρ(y):=1√2πe−y22.

The weak (or variational) formulation of problem (1) is to find for a given the solution such that

 B(y;u,v):=∫Da(x,y)∇u(x,y)⋅∇v(x)dx=∫Df(x)v(x)dx,∀v∈V. (5)

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

 κ:=∥∥ ∥∥∑j≥1|ϕj|bj∥∥ ∥∥L∞(D)<∞, (6)

and

 {bj}j≥1∈ℓp∗(N) for some p∗∈(0,∞). (7)

Let us define for some space and

The following result is implied from [1, Theorem 2.1 and Remark 2.2].

Proposition 1.

Under conditions (6) and (7) it holds that for all .

For -a.s.

we define two random variables

 amin(y):=missingessinfx∈Da(x,y),amax(y):=supx∈D|a(x,y)|.

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.

 amin(y)≥exp(−∥Z(⋅,y)∥L∞(D))=1∥a(⋅,y)∥L∞(D)>0, (8)

and

 amax(y)≤exp(∥Z(⋅,y)∥L∞(D))=∥a(⋅,y)∥L∞(D)<∞. (9)

 (10)

For any and belonging to using Hölder’s inequality we have

 |B(y;u,v)|=∣∣∣∫Da(x,y)∇u(x,y)⋅∇v(x)dx∣∣∣≤amax(y)∥u(⋅,y)∥V∥v∥V.

Taking in (5) yields

 B(y;u,u)=∫Da(x,y)∇u(x,y)⋅∇u(x,y)dx≥amin(y)∥u(⋅,y)∥2V.

The Lax–Milgram lemma implies that there is a unique solution of the weak form (5), moreover,

 ∥u(⋅,y)∥V≤∥f∥V∗amin(y). (11)

Taking the norm of both sides of (11) and applying (10) we have

 ∥u∥Lp,ρ(RN,V)≤∥∥∥1amin∥∥∥Lp,ρ(RN)∥f∥V∗<∞,

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

 ∥v∥Cm(¯¯¯¯D):=∑0≤|τ|≤msupx∈¯¯¯¯D|∂τxv(x)|,

for with . For any real such that we set and define the following norm

 ∥v∥Cr(¯¯¯¯D):=∥v∥C[r](¯¯¯¯D)+max|τ|=[r]supx,x′∈¯¯¯¯D:x≠x′∣∣∂τxv(x)−∂τx′v(x′)∣∣∥x−x′∥{r},

where is the Euclidean norm. Both the cardinality of a set and the

norm of a vector are denoted by

but 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:

 −∇⋅(a(x,yu)∇u(x,yu))=f(x) for x in D,u(x,yu)=0 for x∈∂D,

where . The weak formulation of the -projected problem involves finding for a given such that the following equation holds

 B(yu;u,v):=∫Da(x,yu)∇u(x,yu)⋅∇v(x)dx=∫Df(x)v(x)dx,∀v∈V.

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

 B(yu;uh,v):=∫Da(x,yu)∇uh(x,yu)⋅∇v(x)dx=∫Df(x)v(x)dx,∀v∈Vh.

The MDFEM algorithm relies on the multivariate decomposition, see, e.g., , of the solution

 u(⋅,y)=∑|u|<∞uu(⋅,yu),

where the sum is over all finite subsets , and we take the anchored decomposition with anchor at to obtain the components

 uu(⋅,yu):=∑v⊆u(−1)|u|−|v|u(⋅,yv), (12)

with .

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

Lemma 1.

For any the function depends only on the variables listed in and satisfies

 uu(⋅,yu)=0 when any component of yu equals 0. (13)

Moreover, if the solution has derivatives of arbitrarily high order with respect to the variable , see Proposition 4 below, then

 (∂ωuyuuu)(⋅,yu)=(∂ωuyuu(⋅,⋅u))(⋅,yu) for any ωu∈N|u| and any yu∈Ru, (14)

and

 (∂ωwywuu)(⋅,0u)=0 for any ωw∈N|w| such that w⊂u. (15)

Since the functional is linear we write

 G(u(⋅,y))=∑|u|<∞G(uu(⋅,yu))=∑|u|<∞∑v⊆u(−1)|u|−|v|G(u(⋅,yv)). (16)

We remark that also because is a linear operator the function has the same properties (13), (14) and (15) as .

Under the conditions (6) and (7) the decomposition (16) is well-defined. This is deduced using similar arguments as in [27, Remark 5]. Therefore, we can interchange integral and sum to obtain

 I(G(u)) =∫RN∑|u|<∞G(uu(⋅,yu))dμ(y)=∑|u|<∞∫RNG(uu(⋅,yu))dμ(y) =∑|u|<∞∫R|u|G(uu(⋅,yu))dμu(yu)=:∑|u|<∞Iu(G(uu)), (17)

where .

The MDFEM algorithm to approximate (2) takes the form

 Q(G(u))=∑u∈UQu,nu(G(uhuu)):=∑u∈Unu−1∑i=0w(i)uG(uhuu(⋅,y(i)u)) (18)

where are cubature rules with being their cubature nodes and respective weights, and

 G(uhuu(⋅,yu)):=∑v⊆u(−1)|u|−|v|G(uhu(⋅,yv)). (19)

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

 cost(Q):=∑u∈Unu×cost of evaluating G(uhuu). (20)

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

 |I(G(u))−Q(G(u))| ≤∣∣ ∣∣∑u∉UIu(G(uu))∣∣ ∣∣ +∣∣ ∣∣∑u∈U(Iu(G(uu))−Iu(G(uhuu)))+(Iu−Qu,nu)(G(uhuu))∣∣ ∣∣. (21)

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

 Is(F):=∫RsF(y)ρ(y)dy,

where with an abuse of notation we omit the subscript referring to the dimension of and write

 ρ(y):=s∏j=1ρ(yj)dyj,ρ(y)=1√2πe−y22.

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

 ITs(F):=∫[−T,T]sF(y)ρ(y)dy=(2T)s∫[0,1]sF(T(y))ρ(T(y))dy,

for some and the mapping is defined as

 T(y)=(2Ty1−T,…,2Tys−T).

The final integral is approximated using an -point QMC rule of the form

 Qs,n(F):=(2T)snn−1∑i=0F(T(y(i)))ρ(T(y(i))), (22)

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

 ⟨F,G⟩Hα,0,ρ(R):=α−1∑τ=1F(τ)(0)G(τ)(0)+∫RF(α)(y)G(α)(y)ρ(y)dy. (23)

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

 Kα,0,ρ(x,y)=α−1∑τ=1xττ!yττ!+1{xy>0}(xy)∫+∞0(|x|−t)α−1+(α−1)!(|y|−t)α−1+(α−1)!1ρ(t)dt, (24)

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

 Kα,0,ρ(y,y) =α−1∑r=1(y2r(r!)2)+∫y0(y−t)2(α−1)((α−1)!)21ρ(t)dt≤α−1∑r=1(y2r(r!)2)+1ρ(y)∫y0(y−t)2(α−1)((α−1)!)2dt ≤α−1∑r=1(|y|2r(r!)2)+1ρ(y)|y|2α−1((α−1)!)2(2α−1).

The same bound holds for any . Thus, using we have

 ∫R(Kα,0,ρ(y,y))1/2ρ(y)dy≤∫R(α−1∑r=1(|y|rr!)+1√ρ(y)|y|α−1/2(α−1)!(2α−1)1/2)ρ(y)dy ≤α−1∑r=11r!2r/2Γ(r/2+1/2)√π+1(α−1)!(2α−1)1/22α+1/4Γ(α/2+1/4)π1/4=:M, (25)

where is the Gamma function.

The multivariate space is defined as the tensor product of the one dimensional spaces with the kernel given by

 Kα,0,ρ,s(x,y):=s∏j=1Kα,0,ρ(xj,yj).

The corresponding inner product is

 ⟨F,G⟩Hα,0,ρ,s(Rs):=∑τ∈{1:α}sv={j:τj=α}∫R|v|F(τ)(yv,0−v)G(τ)(yv,0−v)ρ(yv)dyv,

where , and is a vector of variables such that for and otherwise. In many references, see, e.g., [4, 20], the inner product is usually written as a double sum as follows

 ⟨F,G⟩Hα,0,ρ,s(Rs) (26)

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

 ⟨F,G⟩H′α([0,1]):=α−1∑τ=0(∫10F(τ)(y)dy)(∫10G(τ)(y)dy)+∫10F(α)(y)G(α)(y)dy,

with norm .

The multivariate space is the tensor product of the one dimensional space with inner product given by

 ⟨F,G⟩H′α,s([0,1]s) :=∑τ∈{0:α}sv={j:τj=α}∫[0,1]|v|(∫[0,1]s−|v|F(τ)(y)dy−v)(∫[0,1]s−|v|G(τ)(y)dy−v)dyv,

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.

Proposition 2.

For any and , the function belongs to and

 ∥(Fρ)∘T∥H′α,s([0,1]s)≤Cs1,αT(α−1/2)s∥F∥Hα,0,ρ,s(Rs),

where

 C1,α:=α!23α(α(α/2+1)1√2πΓ(2α)I0(1/4))1/2, (27)

and is the modified Bessel function of the first kind with .

3.2 Higher-order quasi-Monte Carlo for integration over Rs

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.

Theorem 1.

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

 |Is(F)−Qs,n(F)| ≤C2,α,λ,s∥F∥Hα,0,ρ,s(Rs)(ln(n))s(α/2+1/4)nλ, (28)

where is a constant independent of and defined below by (35).

Proof.

The error splits into two terms

 |Is(F)−Qs,n(F)| ≤∣∣∣Is(F)−∫[−T,T]sF(y)ρ(y)dy∣∣∣+∣∣∣(2T)s∫[0,1]sF(T(y))ρ(T(y))dy−Qs,n(F)∣∣∣. (29)

The first term is the domain truncation error. Using the reproducing property of and the Cauchy–Schwarz inequality we have

 ∣∣∣Is(F)−∫[−T,T]sF(y)ρ(y)dy∣∣∣ ≤∫Rs∖[−T,T]s|F(y)|ρ(y)dy =∫Rs∖[−T,T]s|⟨F(⋅),Kα,0,ρ,s(y,⋅)⟩|ρ(y)dy ≤∥F∥Hα,0,ρ,s(Rs)∫Rs∖[−T,T]s∥Kα,0,ρ,s(y,⋅)∥Hα,0,ρ,s(Rs)ρ(y)dy =∥F∥Hα,0,ρ,s(Rs)∫Rs∖[−T,T]s(Kα,0,ρ,s(y,y))1/2ρ(y)dy. (30)

The integral of the right hand side is bounded as

 ∫Rs∖[−T,T]s(Kα,0,ρ,s(y,y))1/2ρ(y)dy ≤2s∑j=1∫R⋯∫+∞T⋯∫Rs∏j=1(Kα,0,ρ,s(yj,yj))1/2ρ(yj)dy. (31)

We will estimate each of the above integrals. Using the same argument as in (3.1) we have for any

 ∫+∞T(Kα,0,ρ,s(y,y))1/2ρ(y)