Uncertainty quantification using periodic random variables

TL;DR

This paper introduces a novel periodic random variable model for uncertainty quantification in PDEs, achieving higher convergence rates with lattice cubature rules, and compares it to traditional affine models.

Contribution

The paper proposes a periodic random field model with the same mean and covariance as the affine model, improving convergence in uncertainty quantification.

Findings

Higher order cubature convergence rate of O(n^{-1/p})

Periodic model performs comparably or better than affine model

Numerical examples validate the theoretical convergence improvements

Abstract

Many studies in uncertainty quantification have been carried out under the assumption of an input random field in which a countable number of independent random variables are each uniformly distributed on an interval, with these random variables entering linearly in the input random field (the so-called affine model). In this paper we consider an alternative model of the random field, in which the random variables have the same uniform distribution on an interval, but the random variables enter the input field as periodic functions. The field is constructed in such a way as to have the same mean and covariance function as the affine random field. Higher moments differ from the affine case, but in general the periodic model seems no less desirable. The new model of the random field is used to compute expected values of a quantity of interest arising from an elliptic PDE with random coefficients. The periodicity is shown to yield a higher order cubature convergence rate of $\mathcal{O}(n^{-1/p})$ independently of the dimension when used in conjunction with rank-1 lattice cubature rules constructed using suitably chosen smoothness-driven product and order dependent weights, where $n$ is the number of lattice points and $p$ is the summability exponent of the fluctuations in the series expansion of the random coefficient. We present numerical examples that assess the performance of our method.

Equations

$$A(\mathbf{x},\omega )=\overline{a}(\mathbf{x})+\sum _{j\ge 1}{Y}_j(\omega ){\psi }_j(\mathbf{x}),\mathbf{x}\in D,\omega \in \Omega ,$$

$$\sum _{j\ge 1}\Vert {\psi }_j{\Vert }_{L_{\infty }}<\infty$$

$$\mathbb{E}[A(\mathbf{x},\cdot )]=\overline{a}(\mathbf{x})+\sum _{j\ge 1}\mathbb{E}[Y_j]{\psi }_j(\mathbf{x})=\overline{a}(\mathbf{x}),$$

$${\int }_{-1/2}^{1/2}\Theta (y)\mathrm{d}y=0\text{and}{\int }_{-1/2}^{1/2}{\Theta }^2(y)\mathrm{d}y=\frac{1}{12},$$

$$A(\mathbf{x},\omega )=\overline{a}(\mathbf{x})+\sum _{j\ge 1}\Theta (Y_j(\omega )){\psi }_j(\mathbf{x}),\mathbf{x}\in D,\omega \in \Omega .$$

$$\mathrm{cov}(A)(\mathbf{x},{\mathbf{x}}^{\prime }):=$$

$$=$$

$$\mathrm{var}(A)(\mathbf{x})=\frac{1}{12}\sum _{j\ge 1}{\psi }_j^2(\mathbf{x}),$$

$$\Theta (y)=\frac{1}{\sqrt{6}}\sin (2\pi y),y\in [-\frac{1}{2},\frac{1}{2}];$$

$$A(\mathbf{x},\omega )=\overline{a}(\mathbf{x})+\frac{1}{\sqrt{6}}\sum _{j\ge 1}\sin (2\pi Y_j(\omega )){\psi }_j(\mathbf{x}),\mathbf{x}\in D,\omega \in \Omega .$$

$$\mathbb{E}[(A(\mathbf{x},\cdot )-\overline{a}(\mathbf{x}){)}^4]=\frac{1}{96}\sum _{j\ge 1}{\psi }_j^4(\mathbf{x})+\frac{1}{24}\sum _{j\ge 1}\sum _{k\ge j+1}{\psi }_j^2(\mathbf{x}){\psi }_k^2(\mathbf{x})$$

$$\displaystyle\begin{array}[]{rll}-\nabla\cdot(A({\boldsymbol{x}},\omega)\nabla u({\boldsymbol{x}},\omega))&\!\!\!\!=f({\boldsymbol{x}}),&{\boldsymbol{x}}\in D,\\ u({\boldsymbol{x}},\omega)&\!\!\!\!=0,&{\boldsymbol{x}}\in\partial D,\end{array}$$

$${\int }_{D}K(\mathbf{x},{\mathbf{x}}^{\prime }){\psi }_j({\mathbf{x}}^{\prime })\mathrm{d}{\mathbf{x}}^{\prime }={\lambda }_j{\psi }_j(\mathbf{x}),\mathbf{x}\in D,$$

$${\int }_D|{\psi }_j(\mathbf{x}){|}^2\mathrm{d}\mathbf{x}=12{\lambda }_j\text{for }j\ge 1.$$

$$\sum _{j\ge 1}\Vert {\psi }_j{\Vert }_{L_{\infty }}^p<\infty \text{for some }0<p<1,$$

$$\mathcal{O}(s^{-2/p+1}+h^2+n^{-1/p})$$

$$U:=[-\frac{1}{2},\frac{1}{2}{]}^{\mathbb{N}}\text{and}{U}_s:=[-\frac{1}{2},\frac{1}{2}{]}^s\text{for }s\in \mathbb{N}.$$

$${\int }_{U}F(\mathbf{y})\mathrm{d}\mathbf{y}:=\underset{t\to \infty }{\lim }{\int }_{U_t}F(y_1,\dots ,y_t,0,0,\dots )\mathrm{d}{y}_1\cdots \mathrm{d}{y}_t.$$

$$\displaystyle\int_{\overline{U}_{s}}F\big{(}{\boldsymbol{y}}_{\{s+1:\infty\}}\big{)}\,{\rm d}{\boldsymbol{y}}_{\{s+1:\infty\}}=\lim_{t\to\infty}\int_{U_{t}}F(y_{s+1},\ldots,y_{s+t},0,0,\ldots)\,{\rm d}y_{s+1}\cdots{\rm d}y_{s+t}$$

$$\mathcal{I}:=\{\mathbf{m}\in {\mathbb{N}}_0^{\infty }:|\mathrm{supp}(\mathbf{m})|<\infty \},$$

$$|\mathbf{m}|:=\sum _{j\ge 1}{m}_j$$

$${\mathbf{x}}^{\mathbf{\alpha }}:=\prod _{j\in \mathrm{supp}(\mathbf{\alpha })}{x}_j^{{\alpha }_j}.$$

$$\Vert w{\Vert }_{H_0^1}:=\Vert \nabla w{\Vert }_{L_2},w\in H_0^1(D).$$

$${\int }_{D}a(\mathbf{x},\mathbf{y})\nabla u(\mathbf{x},\mathbf{y})\cdot \nabla \varphi (\mathbf{x})\mathrm{d}\mathbf{x}=⟨f,\varphi ⟩\forall \varphi \in H_0^1(D),$$

$$a(\mathbf{x},\mathbf{y})=\overline{a}(\mathbf{x})+\frac{1}{\sqrt{6}}\sum _{j\ge 1}\sin (2\pi y_j){\psi }_j(\mathbf{x}),\mathbf{x}\in D,\mathbf{y}\in U,$$

$$\mathbb{E}[G(u)]={\int }_{U}G(u(\mathbf{y}))\mathrm{d}\mathbf{y}.$$

$$\Vert v{\Vert }_{W^{1,\infty }}:=\max \{\Vert v{\Vert }_{L_{\infty }},\Vert \nabla v{\Vert }_{L_{\infty }}\};$$

$$u^s(\cdot ,\mathbf{y}):=u(\cdot ,(y_1,\dots ,y_s,0,0,\dots ))\forall \mathbf{y}\in U,$$

$$\Vert u(\cdot ,\mathbf{y}){\Vert }_{H_0^1}\le \frac{\Vert f{\Vert }_{H^{-1}}}{a_{\min }}$$

$${\partial }^{\mathbf{m}}a(\mathbf{x},\mathbf{y})=\{\begin{array}{ll}a(\mathbf{x},\mathbf{y})& \text{if }\mathbf{m}=0,\\ \frac{1}{\sqrt{6}}(2\pi {)}^k\sin \left(2\pi y_j+k\frac{\pi }{2}\right){\psi }_j(\mathbf{x})& \text{if }\mathbf{m}=k{\mathbf{e}}_j,k\ge 1,\\ 0& \text{otherwise,}\end{array}$$

Full text

1 Introduction

This paper is concerned with the development and use of specially designed random fields on a physical domain $D\subseteq\mathbb{R}^{d}$, where $d=1$, $2$, or $3$. For simplicity we assume that the boundary $\partial D$ is Lipschitz.

Many studies in uncertainty quantification are modeled by PDEs over the domain $D$, in which one or more of the coefficients is a random field over $D$. In particular, many recent papers (including [3, 5, 6, 9, 16, 17, 23]) have used an “affine” model of the random field, taking the form

[TABLE]

where $(\Omega,A,\mathbb{P})$ is a probability space, $Y_{1},Y_{2},\ldots$ are independently and identically distributed random variables uniformly distributed on $[-\tfrac{1}{2},\tfrac{1}{2}],$ and $(\psi_{j})_{j\geq 1}$ are real-valued $L_{\infty}$ functions on $D$ satisfying

[TABLE]

and otherwise for the moment not specified. With these definitions the sum in (1.1) converges uniformly on $D$ for all values of the $Y_{j}$, and the random field is pointwise well defined. The expected value at ${\boldsymbol{x}}\in D$ is

[TABLE]

since $\mathbb{E}[Y_{j}]=\int_{-1/2}^{1/2}y\,{\mathrm{d}}y=0$ for $j\geq 1.$

The fact that the random variable in (1.1) occurs linearly seems to be a result of history, rather than something imposed by modeling assumptions. Suppose instead that we replace $Y_{j}$ in (1.1) by $\Theta(Y_{j})$, where $\Theta\!:[-\tfrac{1}{2},\tfrac{1}{2}]\to\mathbb{R}$ is a continuous function with the properties

[TABLE]

both of which are satisfied by the special affine choice $\Theta(y)=y$. Thus we replace (1.1) by

[TABLE]

It has exactly the same mean as (1.1) because $\mathbb{E}[\Theta(Y_{j})]=0$ and also has the same covariance

[TABLE]

because $\mathbb{E}[\Theta(Y_{j})\,\Theta(Y_{j^{\prime}})]$ vanishes for $j\neq j^{\prime}$ by the independence of $Y_{j}$ and $Y_{j^{\prime}}$ (and hence of $\Theta(Y_{j})$ and $\Theta(Y_{j^{\prime}})$), and $\mathbb{E}[\Theta^{2}(Y_{j})]=\frac{1}{12}=\mathbb{E}[Y_{j}^{2}]$ from (1.3). In particular, the variance, obtained by setting ${\boldsymbol{x}}={\boldsymbol{x}}^{\prime}$, is

[TABLE]

independently of the choice of the function $\Theta$.

In this paper we explore a periodic choice of $\Theta$ satisfying (1.3), namely,

[TABLE]

thus our model of the random field, instead of (1.1), becomes

[TABLE]

An ensemble of randomly generated realizations of a pair of affine and periodic random fields is illustrated in Figure 1  with the choices $\overline{a}(x):=2$ and $\psi_{j}(x):=j^{-3/2}\sin((j-\tfrac{1}{2})\pi x)$ for $x\in[0,1]$ and $j\in\mathbb{N}$. While the individual realizations of the affine and periodic fields are obviously different, the statistical moments of the fields coincide up to second order, but all higher moments will be different; see Figure 2. We note that the fourth central moment (or nonstandardized kurtosis) of the random field is

[TABLE]

for the periodic model, whereas the factor $1/96$ is replaced by $1/80$ for the affine model. We are not aware of any modeling reason to prefer the affine model over the periodic model.

To give further insight into the proposed new model, we note that we are in effect replacing each of the countably infinite number of uniformly distributed random variables $Y_{j}$ by new random variables $\sin(2\pi Y_{j})$ (when renormalized to lie between $-1$ and $+1$). It is easily seen that the probability density function for each of these new random variables is $1/(\pi\sqrt{1-y^{2}})$ for $-1\leq y\leq 1$, which in the context of orthogonal polynomials is the weight function associated with Chebyshev polynomials of the first kind.

This observation means that there is a close connection between our new model and generalized polynomial chaos (GPC) (see [28]) in which the dependence of the solution on the stochastic variables is expressed as a linear combination of multivariate basis functions of orthogonal polynomials with respect to a variety of weight functions. In the GPC setting we have merely changed from the uniform probability distribution to the one associated with Chebyshev polynomials of the first kind: a popular choice (see [1, 2, 12, 22, 25]) because of the attractive simplicity of the Chebyshev polynomials. Nevertheless, it should be emphasized that the proposed approximation scheme is not the same as GPC, since the input field (and hence also the solution) is here made to be periodic, and the natural approximation is by trigonometric rather than algebraic polynomials.

We use the random field (1.6) to model the uncertain diffusion coefficient of the following PDE problem: find $u\!:D\times\Omega\to\mathbb{R}$ that satisfies

[TABLE]

for almost all events $\omega\in\Omega$. Then we approximate $\mathbb{E}[G(u)]$, where $G$ is a bounded, linear functional of the solution to (1.9). The motivation for the choice (1.6) of the random field is that the random field $u$ is now a 1-periodic function of the random variable $\omega$, and periodic integrands are known to be especially advantageous in the context of so-called lattice cubature rules [26]. By using the periodic model of random fields instead of the affine model, it is possible to carry out lattice rule calculations of expected values in high dimensions with higher order convergence rates, instead of being restricted, as, for example, in [16, 17, 23], to a convergence rate of at best $\mathcal{O}(n^{-1})$.

In this paper we leave open the choice of the fluctuations $(\psi_{j})_{j\geq 1}$ but note that if the covariance function $K({\boldsymbol{x}},{\boldsymbol{x}}^{\prime}):={\rm cov}(A)({\boldsymbol{x}},{\boldsymbol{x}}^{\prime})$ of the field $A({\boldsymbol{x}},\omega)$ is specified, then the appropriate choice is to take the $\psi_{j}$ to be suitably normalized eigenfunctions of the integral operator with kernel $K$:

[TABLE]

where $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$ are the eigenvalues of the integral operator, and the eigenfunctions are orthogonal with respect to the $L_{2}$-inner product $\langle\cdot,\cdot\rangle_{L_{2}}$ and normalized by

[TABLE]

In this case (1.5) becomes Mercer’s theorem for the covariance function, and (1.4) is a version of the Karhunen–Loève expansion; see [11, 18, 24]. The only nonstandard point in the proof of the Karhunen–Loève theorem is the occurrence of the function $\Theta$, but the only properties that are needed are those in (1.3).

The main result in this paper is as follows. We show that if in addition to (1.2), the fluctuation operators satisfy

[TABLE]

as well as certain regularity assumptions to be made precise later, the overall error of the discretized PDE problem (1.9) with the diffusion coefficient truncated to the first $s$ terms is given by

[TABLE]

using a first order finite element solver with mesh size $h$ and $n$ rank-1 lattice cubature points in $[-\tfrac{1}{2},\tfrac{1}{2}]^{s}$ generated using the component-by-component (CBC) algorithm. The error term (1.10) consists of the dimension truncation error, first order finite element discretization error, and the cubature error, respectively, and the implied coefficient is independent of the truncated dimension $s$ as well as $h$ and $n$. In particular, we note that we are able to obtain a higher order lattice cubature convergence rate $\mathcal{O}(n^{-1/p})$ beating the $\mathcal{O}(n^{-\min\{1/p-1/2,1-\delta\}})$, $\delta>0$, rate for randomly shifted lattice rules in, e.g., [16, 23]. The same rate $\mathcal{O}(n^{-1/p})$ has been obtained for the affine model in, e.g., [6] but with interlaced polynomial lattice rules, which are more complicated than rank-1 lattice rules in their construction. Moreover, the dimension truncation error rate $\mathcal{O}(s^{-2/p+1})$ in (1.10) matches the recent result for affine-parametric operator equations [8]. We also discuss the case $p=1$. Higher order convergence for the finite element error can potentially be obtained by using higher order elements.

This paper is structured as follows. We present the notations and discuss the preliminaries in Subsection 1.1. The periodic parametric mathematical model is introduced in Section 2. We assess the regularity of this model with respect to the parametric variable in Subsection 2.1 and consider the dimension truncation error and finite element discretization errors in Subsection 2.2. The quasi-Monte Carlo (QMC) method as it applies to the periodic framework is discussed in Section 3, and we show in Subsection 3.1 that the use of rank-1 lattice rules in the periodic setting yields a higher order convergence rate for our model provided that appropriate smoothness-driven product and order dependent (SPOD) weights are used in the lattice rule construction. To this end, Subsection 3.2 contains a description of the fast CBC algorithm for rank-1 lattice cubature rules using SPOD weights. The overall error estimate for the discretized PDE problem is presented in Section 4. We present numerical experiments in Section 5 that assess the cubature convergence rate. We end this paper with some conclusions on our results.

1.1 Notations and preliminaries

We follow the convention $\mathbb{N}=\{1,2,3,\ldots\}$ and use $\mathbb{N}_{0}$ to denote the set of natural numbers including zero. Moreover, we use the shorthand notation $\{m:n\}=\{m,m+1,\ldots,n\}$ for integers such that $m\leq n$ and set

[TABLE]

We define the integral over the set $U$ by

[TABLE]

For fixed $s\in\mathbb{N}$, we introduce the set $\overline{U}_{s}:=\{(y_{j})_{j\geq s+1}:y_{j}\in[-\tfrac{1}{2},\tfrac{1}{2}],~{}j\geq s+1\}$ and write

[TABLE]

to mean integration over the variables $(y_{j})_{j\geq s+1}$.

Let the set of all multi-indices with finite support be denoted by

[TABLE]

where we define the support of a multi-index by ${\rm supp}(\boldsymbol{m}):=\{j\in\mathbb{N}:m_{j}\neq 0\}$, and $|{\rm supp}(\boldsymbol{m})|$ is the cardinality of the support. Here and throughout this manuscript, we refer to the $j{\rm th}$ component of a multi-index $\boldsymbol{m}$ as $m_{j}$. Moreover, we define

[TABLE]

for multi-indices $\boldsymbol{m}\in\mathcal{I}$. Let $\boldsymbol{x}=(x_{j})_{j\geq 1}$ be a sequence and $\boldsymbol{\alpha}\in\mathcal{I}$. We denote

[TABLE]

In addition, we use the notation $\boldsymbol{\alpha}\leq\boldsymbol{\beta}$ to signify that $\alpha_{j}\leq\beta_{j}$ for all $j\geq 1$.

We assume in the sequel that $D\subseteq\mathbb{R}^{d}$, $d\in\{1,2,3\}$, is a bounded domain with a Lipschitz regular boundary. This assumption also justifies us taking the Sobolev norm of the space $H_{0}^{1}(D)$ to be

[TABLE]

The duality pairing between $H^{-1}(D)$ of $H_{0}^{1}(D)$ is denoted by $\langle\cdot,\cdot\rangle$.

We establish the following notations and assumptions regarding the finite element approximation of $w\in H_{0}^{1}(D)$. Let us assume that $D$ is a convex and bounded polyhedron with plane faces. We denote by $\{V_{h}\}_{h}$ a family of finite element subspaces $V_{h}\subset H_{0}^{1}(D)$, parametrized by the mesh size $h>0$, which are spanned by continuous, piecewise linear finite element basis functions such that each $V_{h}$ is obtained from an initial, regular triangulation of $D$ by recursive, uniform bisection of simplices. We use the notation $w_{h}\in V_{h}$ to denote the finite element approximation of $w$ in the finite element space $V_{h}$.

2 Parametric weak formulation

The parametric weak formulation of (1.9) is, for ${\boldsymbol{y}}\in U$, to find $u(\cdot,{{\boldsymbol{y}}})\in H_{0}^{1}(D)$ such that

[TABLE]

where $f\in H^{-1}(D)$, and the diffusion coefficient is assumed to have the form

[TABLE]

consistently with (1.6). Furthermore, let $G\!:H_{0}^{1}(D)\to\mathbb{R}$ be a bounded, linear mapping. As the quantity of interest, we consider the expectation of $\boldsymbol{y}\mapsto G(u(\cdot,\boldsymbol{y}))$ taken over the parametric space:

[TABLE]

We state the following assumptions which are the same as the assumptions in [16]:

• (A1)

$\overline{a}\in L_{\infty}(D)$ and $\sum_{j\geq 1}\|\psi_{j}\|_{L_{\infty}}<\infty$;

• (A2)

there exist positive constants $a_{\rm max}$ and $a_{\rm min}$ such that $0<a_{\rm min}\leq a({\boldsymbol{x}},{{\boldsymbol{y}}})\leq a_{\rm max}<\infty$ for all ${\boldsymbol{x}}\in D$ and ${\boldsymbol{y}}\in U$;

• (A3)

$\sum_{j\geq 1}\|\psi_{j}\|_{L_{\infty}}^{p}<\infty$ for some $0<p\leq 1$;

• (A4)

$\overline{a}\in W^{1,\infty}(D)$ and $\sum_{j\geq 1}\|\psi_{j}\|_{W^{1,\infty}}<\infty$, where

[TABLE]

• (A5)

$\|\psi_{1}\|_{L_{\infty}}\geq\|\psi_{2}\|_{L_{\infty}}\geq\cdots$;

• (A6)

the physical domain $D\subseteq\mathbb{R}^{d}$, $d\in\{1,2,3\}$, is a convex and bounded polyhedron with plane faces.

We refer to these assumptions as they are needed.

For convenience, we introduce the following notation to mean the dimensionally truncated exact solution to (2.1):

[TABLE]

and we define $u_{h}^{s}(\cdot,{\boldsymbol{y}}):=u_{h}(\cdot,(y_{1},\ldots,y_{s},0,0,\ldots))\in V_{h}$ for all ${\boldsymbol{y}}\in U$ to mean the dimensionally truncated finite element solution to (2.1).

2.1 Parametric regularity of the solution

We proceed to derive a regularity estimate for the problem (2.1) with respect to the parametric variable ${\boldsymbol{y}}$. The approach we take here follows the argument of [16], where a uniform affine model of the uncertain diffusion coefficient was considered.

We begin by remarking that a straightforward application of the Lax–Milgram lemma ensures that (2.1) is uniquely solvable over the whole parametric domain and that the solution can be bounded a priori.

Lemma 2.1**.**

Under the assumptions (A1) and (A2), the weak formulation (2.1) has a unique solution $u(\cdot,\boldsymbol{y})\in H_{0}^{1}(D)$ for any $\boldsymbol{y}\in U$ such that

[TABLE]

for any source term $f\in H^{-1}(D)$.

Let $\boldsymbol{m}\in\mathcal{I}$ be a multi-index. It is easy to see that the mixed partial derivatives of (2.2) with respect to ${\boldsymbol{y}}$ are

[TABLE]

where $\mathbf{e}_{j}\in\mathcal{I}$ denotes the multi-index whose $j{\rm th}$ component is $1$ and all other components are [math]. This is due to the dependence of $a$ on each $y_{j}$ being in separate additive terms: if we differentiate once or more with respect to $y_{j}$, then we obtain an expression depending only on $y_{j}$ and $\psi_{j}$, and if we then differentiate with respect to a different component of the ${\boldsymbol{y}}$ variable, we get [math].

Let $\boldsymbol{\nu}\in\mathcal{I}$ be a multi-index with $\boldsymbol{\nu}\neq\mathbf{0}$. We differentiate (2.1) on both sides to get

[TABLE]

which, after an application of the Leibniz product rule, yields

[TABLE]

Plugging in (2.4) and separating out the case $\boldsymbol{m}=\mathbf{0}$, we obtain

[TABLE]

for all $\phi\in H_{0}^{1}(D)$. In particular, we can choose to test this formula against $\phi=(\partial^{\boldsymbol{\nu}}u)(\cdot,{\boldsymbol{y}})$. By applying the ellipticity assumption $a({\boldsymbol{x}},{\boldsymbol{y}})\geq a_{\rm min}$ on the left-hand side and $|\psi_{j}({\boldsymbol{x}})|\leq\|\psi_{j}\|_{L_{\infty}}$ as well as the Cauchy–Schwarz inequality on the right-hand side, we obtain

[TABLE]

Eliminating the common factor on both sides and using (1.12) yields for ${\boldsymbol{\nu}}\neq\mathbf{0}$

[TABLE]

where we set

[TABLE]

This differs from the definition of $b_{j}$ in [16] by the factor $1/\sqrt{6}$.

Our goal is to use the recurrence (2.5) to derive an explicit upper bound on the term $\|\partial^{\boldsymbol{\nu}}u(\cdot,{\boldsymbol{y}})\|_{H_{0}^{1}}$ for all $\boldsymbol{\nu}\in\mathcal{I}$. It turns out that Stirling numbers of the second kind (or Stirling partition numbers) play a large role in the forthcoming analysis; they are defined by

[TABLE]

for $n\geq k\geq 0$ except for $S(0,0):=1$. The bound on $\|\partial^{\boldsymbol{\nu}}u(\cdot,{\boldsymbol{y}})\|_{H_{0}^{1}}$ (see Theorem 2.3 below) follows from the following result, stated in a general form in case it is useful in other contexts.

Lemma 2.2**.**

Let $B,c>0$, and let $(\mathbb{A}_{\boldsymbol{\nu}})_{\boldsymbol{\nu}\in\mathcal{I}}$ and $(\Upsilon_{j})_{j\in\mathbb{N}}$ be sequences of non-negative real numbers that satisfy the recurrence

[TABLE]

Then

[TABLE]

Moreover, if equalities hold in the formulae (2.7), then there is equality in (2.8).

Proof.

We prove this result by carrying out an induction argument on $|\boldsymbol{\nu}|$ based on the recurrence (2.7). The base step $\boldsymbol{\nu}=\mathbf{0}$ is resolved immediately. For arbitrary $\boldsymbol{\nu}\in\mathcal{I}\setminus\{\boldsymbol{0}\}$, suppose that the claim holds for all multi-indices of order $<|\boldsymbol{\nu}|$. In particular, if $\nu_{j}\geq k\geq 1$ for some $j\geq 1$, then the induction hypothesis gives

[TABLE]

Applying the recursion (2.7) in conjunction with the inequality above yields

[TABLE]

For given $\boldsymbol{m}$, ${\boldsymbol{\nu}}$, $\boldsymbol{\Upsilon}$ and an index $j$, we define $\boldsymbol{m^{\prime}}=(m_{1},\ldots,m_{j-1},m_{j+1},\ldots)$, $\boldsymbol{\nu^{\prime}}=(\nu_{1},\ldots,\nu_{j-1},\nu_{j+1},\ldots)$, and $\boldsymbol{\Upsilon^{\prime}}=(\Upsilon_{1},\ldots,\Upsilon_{j-1},\Upsilon_{j+1},\ldots)$, respectively. Then we may write the $j{\rm th}$ term in the outer sum from (2.9) as

[TABLE]

where we swapped the order of the sums over $k$ and $m_{j}$. Furthermore, it holds that

[TABLE]

which can be verified either by direct calculation based on the definition of $S(n,k)$ or as a consequence of [21, equation (9.25)]. Thus (2.1) becomes

[TABLE]

and together with (2.9) this yields

[TABLE]

Since $S(k,0)=0$ for all $k\geq 1$, a straightforward computation shows that

[TABLE]

which simplifies the upper bound into

[TABLE]

completing the proof. ∎

The desired result can be obtained as an immediate corollary to Lemma 2.2 using Lemma 2.2 and (2.5).

Theorem 2.3**.**

Under the assumptions (A1) and (A2), for any ${\boldsymbol{y}}\in U$, let $u(\cdot,{\boldsymbol{y}})\in H_{0}^{1}(D)$ be the solution of the problem (2.1) with the source term $f\in H^{-1}(D)$, and let $\boldsymbol{b}=(b_{j})_{j\geq 1}$ be the sequence defined by (2.6). Then for any multi-index $\boldsymbol{\nu}\in\mathcal{I}$ we have

[TABLE]

The result also holds for the dimension-truncated finite element solution $u_{h}^{s}(\cdot,{\boldsymbol{y}})\in V_{h}$ for all $s\in\mathbb{N}$, ${\boldsymbol{y}}\in U$.

2.2 Dimension truncation and finite element discretization errors

In practice, it is generally only possible to solve the problem (2.1) approximately using, e.g., the finite element method and with the series (2.2) truncated to finitely many terms. In this section, we discuss the approximation errors caused by the finite element discretization and dimension truncation.

In the affine setting, the fundamental dimension truncation error bound has already been discussed in [16] leading to an error bound of the order $\mathcal{O}(s^{-2/p+2})$. While this analysis can also be applied to the periodic setting with only minuscule changes to the argument, Gantner [8] has recently proved an improved bound of order $\mathcal{O}(s^{-2/p+1})$ in the context of affine-parametric operator equations. In the following, we prove an analogous result for the problem (2.1)–(2.3). While the proof technique we use is the same as in [8], we present the proof for completeness in order to highlight that the result holds also in the periodic framework. The following proof also differs from [8] insofar as we do not need to put a restriction on the size of the sum, e.g., $\sum_{j\geq 1}b_{j}<\sqrt{6}$.

Lemma 2.4** (cf. [8, Theorem 1]).**

Under the assumptions (A1)–(A3) and (A5), for any ${\boldsymbol{y}}\in U$, let $u(\cdot,{\boldsymbol{y}})\in H_{0}^{1}(D)$ denote the solution to the problem (2.1) with the source term $f\in H^{-1}(D)$, and let $G\in H^{-1}(D)$. If $0<p<1$, then for any $s\in\mathbb{N}$ there exists a constant $C>0$ such that

[TABLE]

If $p=1$, then

[TABLE]

In both cases, $C>0$ denotes a generic constant that does not depend on $s$, $f$, or $G$.

Proof.

We define the operators $A({\boldsymbol{y}})\!:H_{0}^{1}(D)\to H^{-1}(D)$ for ${\boldsymbol{y}}\in U$, and $A_{j}\!:H_{0}^{1}(D)\to H^{-1}(D)$ for $j\in\mathbb{N}$, by setting $\langle A({\boldsymbol{y}})w,\phi\rangle:=\langle a(\cdot,{\boldsymbol{y}})\nabla w,\nabla\phi\rangle_{L_{2}}$ for all $\phi\in H_{0}^{1}(D)$, and $\langle A_{j}w,\phi\rangle:=\langle\frac{1}{\sqrt{6}}\psi_{j}\nabla w,\nabla\phi\rangle_{L_{2}}$ for all $\phi\in H_{0}^{1}(D)$, respectively. Moreover, we define $A^{s}({\boldsymbol{y}}):=A((y_{1},\ldots,y_{s},0,0,\ldots))$ and denote $u({\boldsymbol{y}}):=u(\cdot,{\boldsymbol{y}})$ and $u^{s}({\boldsymbol{y}}):=u^{s}(\cdot,{\boldsymbol{y}})$ for all $s\in\mathbb{N}$, ${\boldsymbol{y}}\in U$. These definitions lead to the identity

[TABLE]

Let $w\in H_{0}^{1}(D)$. Lemma 2.1 and (1.12) together with

[TABLE]

imply that both operators $A({\boldsymbol{y}})$ and $A^{s}({\boldsymbol{y}})$ are boundedly invertible linear maps for all ${\boldsymbol{y}}\in U$. Furthermore, we obtain

[TABLE]

where the sequence $(b_{j})_{j\geq 1}$ is defined as in (2.6). In consequence, this yields

[TABLE]

In what follows, we omit the argument ${\boldsymbol{y}}$ and denote the operator norm by $\|\cdot\|=\|\cdot\|_{\mathscr{L}(H_{0}^{1}(D))}$ for brevity.

Since the sequence $(b_{j})_{j\geq 1}$ is summable, there exists $s^{\prime}\in\mathbb{N}$ such that for all $s\geq s^{\prime}$ the upper bound in (2.12) is at most $1/2$. Let us assume that $s\geq s^{\prime}$. For future reference, we note that this implies for all $s\geq s^{\prime}$

[TABLE]

It follows from (2.12) and our assumption $s\geq s^{\prime}$ that the Neumann series

[TABLE]

is well defined. Moreover, we have the representation

[TABLE]

For each $k\in\mathbb{N}$, we note that the integrand in (2.14) can be expanded as

[TABLE]

where the product symbol is assumed to respect the order of the noncommutative operators. Using the independence of the components of ${\boldsymbol{y}}\in U$ and (1.11), the integral over $U$ in (2.14) can be written as a product of integrals

[TABLE]

where $I_{1}\geq 0$ because it can be written as a product of univariate integrals of the form $\int_{-1/2}^{1/2}\sin(2\pi y_{j})^{m}\,{\mathrm{d}}y_{j}$ for $m\in{\mathbb{N}}$, which take values between [math] and $1$ (importantly, this expression is zero when $m=1$), while we can estimate $I_{2}$ by

[TABLE]

Thus

[TABLE]

where we have used the multinomial theorem together with $\boldsymbol{\nu}!:=\prod_{i\geq 1}\nu_{i}!$ for $\boldsymbol{\nu}\in\mathcal{I}$, Lemma 2.1, and the bound (2.11). The key observation is that this term vanishes whenever any component of $\boldsymbol{\nu}$ is equal to $1$, and consequently the term vanishes when $k=1$.

We may now estimate (2.14) by splitting the sum into the $k\geq k^{\prime}$ terms and the $k<k^{\prime}$ terms for a value of $k^{\prime}$ to be specified later. We obtain

[TABLE]

Consider first the case $0<p<1$. The $k\geq k^{\prime}$ terms can be bounded using the geometric series as

[TABLE]

where we used the inequality $\sum_{j\geq s+1}b_{j}\leq(\sum_{j\geq 1}b_{j}^{p})^{1/p}s^{-1/p+1}$ (see [16, Theorem 5.1]), and the ensuing constant $C_{1}:=2(\sum_{j\geq 1}b_{j}^{p})^{k^{\prime}/p}$ is independent of $s$, $f$, and $G$. On the other hand, for each $2\leq k<k^{\prime}$ we use the estimate

[TABLE]

where we used both inequalities in (2.13), the inequalities ${\rm e}^{x}\leq 1+({\rm e}-1)x$ for all $x\in[0,1]$ and $\sum_{j\geq s+1}b_{j}^{2}\leq\frac{1}{2/p-1}(\sum_{j\geq 1}b_{j}^{p})^{1/p}s^{-2/p+1}$, and the resulting constant $C_{2}:=\frac{2({\rm e}-1)}{2/p-1}(\sum_{j\geq 1}b_{j}^{p})^{1/p}$ is independent of $s$, $f$, and $G$. Hence we conclude that

[TABLE]

We therefore choose $k^{\prime}:=\lceil(2-p)/(1-p)\rceil$ to balance the two terms. This proves the assertion for $s\geq s^{\prime}$ after a trivial adjustment of the constant factors. The result can be extended to all $s\in\mathbb{N}$ by noticing that

[TABLE]

holds for all $1\leq s<s^{\prime}$, and the claim follows by a trivial adjustment of all of the constants involved.

For $p=1$ we amend the above argument slightly to obtain

[TABLE]

where $C>0$ is a constant independent of $s$, $f$, and $G$. ∎

Regarding the finite element approximation error, it is clear that an analogous result to the one presented in [16] holds.

Lemma 2.5** (cf. [16, Theorem 5.1]).**

Under assumptions (A1), (A2), (A4), and (A6), for any ${\boldsymbol{y}}\in U$, let $u(\cdot,{\boldsymbol{y}})\in H_{0}^{1}(D)$ denote the solution to (2.1) with the source term $f\in H^{-1+t}(D)$ such that $0\leq t\leq 1$, and let $G\in H^{-1+t^{\prime}}(D)$ with $0\leq t^{\prime}\leq 1$. Then the finite element approximations satisfy the following asymptotic convergence estimate as $h\to 0$:

[TABLE]

where $0\leq t+t^{\prime}\leq 2$ and the constant $C>0$ is independent of $h$ and $\boldsymbol{y}$.

Remark. We note that the limiting case $t=t^{\prime}=1$ in Lemma 2.5 corresponds to taking $f\in L_{2}(D)$ and $G\in L_{2}(D)$, where the dual of $L_{2}(D)$ is identified with itself, resulting in a convergence rate of $\mathcal{O}(h^{2})$.

3 QMC in the periodic setting

QMC methods are a class of numerical methods designed to approximate multivariate integrals such as

[TABLE]

for a continuous integrand $F$ by using an equal weight cubature formula of the form

[TABLE]

where $\boldsymbol{y}_{0},\ldots,\boldsymbol{y}_{n-1}\in[0,1]^{s}$ are prescribed cubature nodes.

We consider rank-1 lattice rules, where the QMC nodes $\Lambda:=\{\boldsymbol{y}_{0},\ldots,\boldsymbol{y}_{n-1}\}$ are taken to be of the form

[TABLE]

where $\{\boldsymbol{x}\}$ denotes taking the componentwise fractional part of $\boldsymbol{x}\in\mathbb{R}^{s}$ and $\boldsymbol{z}\in\mathbb{N}^{s}$ is called the generating vector of a lattice rule. It is well known that the lattice rule error for functions with absolutely convergent Fourier series is precisely [27]

[TABLE]

where $\hat{F}(\boldsymbol{h}):=\int_{[0,1]^{s}}F({\boldsymbol{y}}){\rm e}^{-2\pi{\rm i}{\boldsymbol{y}}\cdot\boldsymbol{h}}\,{\rm d}{\boldsymbol{y}}$ for $\boldsymbol{h}\in\mathbb{Z}^{s}$ and we denote the dual lattice by $\Lambda^{\perp}=\Lambda^{\perp}(\boldsymbol{z})=\{\boldsymbol{h}\in\mathbb{Z}^{s}:\boldsymbol{h}\cdot\boldsymbol{z}\equiv 0~{}({\rm mod}~{}n)\}$, which is defined with respect to the generating vector $\boldsymbol{z}$ of the rank-1 lattice rule.

Let $F\in C([0,1)^{s})$ be a 1-periodic function with respect to each of its variables, and set

[TABLE]

where ${\rm supp}(\boldsymbol{h}):=\{j\in\{1:s\}:h_{j}\neq 0\}$ and ${\boldsymbol{\gamma}}=(\gamma_{\mathfrak{u}})_{\mathfrak{u}\subseteq\{1:s\}}$ denotes a collection of nonnegative weights. Using the error formula (3.1), we obtain

[TABLE]

where the factor depending only on the QMC nodes is defined by

[TABLE]

and the norm is given by

[TABLE]

Since the inequality (3.2) is sharp, we see that $P_{\alpha}({\boldsymbol{\gamma}},\boldsymbol{z})$ is the worst-case error in the space with $\|F\|_{\alpha}\leq 1$. The quantity $P_{\alpha}({\boldsymbol{\gamma}},\boldsymbol{z})$ is well known in classical lattice rule literature (at least for the unweighted case $\gamma_{\mathfrak{u}}\equiv 1$, see [26]) and coincides with the squared error term in the Hilbert space setting considered in the paper [7], leading us to conclude the following.

Lemma 3.1**.**

Let $s\in\mathbb{N}$ and prime $n$, and let ${\boldsymbol{\gamma}}=(\gamma_{\mathfrak{u}})_{\mathfrak{u}\subseteq\{1:s\}}$ be a collection of nonnegative weights. Let $F\in C([0,1)^{s})$ be a 1-periodic function with respect to each of its variables such that $\|F\|_{\alpha}<\infty$. Then a generating vector $\boldsymbol{z}\in\mathbb{N}^{s}$ can be constructed by the CBC algorithm such that

[TABLE]

for $\lambda\in(1/\alpha,1]$. Here, $\zeta(x):=\sum_{k\geq 1}k^{-x}$ denotes the Riemann zeta function for $x>1$.

Proof.

It can be readily verified that $F$ has an absolutely convergent Fourier series given that $\|F\|_{\alpha}<\infty$. The claim then follows from the previous discussion in conjunction with [7, Theorem 5]. ∎

The result can be extended to nonprime $n$ by replacing $n-1$ with Euler’s totient function $\varphi_{\rm tot}(n):=|\{m\in\mathbb{\{}1:n-1\}:{\rm gcd}(m,n)=1\}|$. In particular, $1/\varphi_{\rm tot}(n)\leq 2/n$ if $n$ is a prime power.

When $\alpha\geq 2$ is an integer, it can be shown that

[TABLE]

provided that $F$ has mixed partial derivatives of order $\alpha$. Furthermore, when $\alpha$ is even, we can write

[TABLE]

where

[TABLE]

and $B_{\alpha}$ denotes the Bernoulli polynomial of degree $\alpha$.

3.1 Higher order convergence in the PDE context

In this section, we let the assumptions (A1)–(A3) be in effect. We are interested in the expectation of the functional $F(\boldsymbol{y}):=G(u(\cdot,\boldsymbol{y}-\tfrac{\mathbf{1}}{\mathbf{2}}))$, where $G$ denotes a bounded, linear functional $G\!:H_{0}^{1}(D)\to\mathbb{R}$, $u(\cdot,\boldsymbol{y}-\tfrac{\mathbf{1}}{\mathbf{2}})\in H_{0}^{1}(D)$ is the solution to the weak formulation (2.1), and we let $\boldsymbol{y}\in[0,1]^{\mathbb{N}}$.

For an integer $\alpha\geq 2$, we estimate the norm as follows:

[TABLE]

We thus obtain, using (3.3) and Theorem 2.3,

[TABLE]

since $S(\alpha,0)=0$ for $\alpha\neq 0$. We now choose the weights to be

[TABLE]

which ensures that $\|F\|_{\alpha}$ is bounded. These weights have a very specific form: they are SPOD weights, first seen in [6]. We then observe that the bound for the error term in Lemma 3.1 becomes

[TABLE]

where

[TABLE]

for $\lambda\in(1/\alpha,1]$ and $n$ a prime power.

Finally, we need to choose $\lambda$ in such a way that $C(s,\alpha,\lambda)$ is bounded independently of $s$. First applying the inequality (cf. [13, Theorem 19])

[TABLE]

to the inner sum of $C(s,\alpha,\lambda)$ and denoting $S_{\max}(\alpha):=\max_{k\in\{1:\alpha\}}S(\alpha,k)$ yields

[TABLE]

where we have set $\beta_{j}:=\max\{1,S_{\max}(\alpha)(2\zeta(\alpha\lambda))^{1/\lambda}\}b_{j}$. We recast the double sum as a sum over multi-indices ${\boldsymbol{\nu}}$:

[TABLE]

Let us define the sequence $d_{j}=\beta_{\lceil j/\alpha\rceil}$, $j\geq 1$. In concrete terms, this means that

[TABLE]

We relate this definition to $C(s,\alpha,\lambda)$ by observing that

[TABLE]

The final inequality holds because $(\sum_{j\geq 1}d_{j}^{\lambda})^{\ell}$ includes all the products of the form $\prod_{j\in\mathfrak{v}}d_{j}^{\lambda}$ with $|\mathfrak{v}|=\ell$, but since the order in which the terms in the product appear does not matter, we can divide this by $\ell!$.

We now choose $\lambda=p$ and verify that the last expression is finite with this choice of $\lambda$. Our assumption that $(\|\psi_{j}\|_{L_{\infty}})_{j\geq 1}\in\ell^{p}$ for some $p\in(0,1]$ implies that $\sum_{j\geq 1}b_{j}^{p}<\infty$. For the inner sum, we now have

[TABLE]

provided that $\alpha p>1$. If additionally $p<1$, then the ratio test implies convergence of the outer sum since

[TABLE]

If $p=1$, then the sum is geometric and converges if and only if $T<1$. An equivalent condition is that

[TABLE]

Since the condition $1/\alpha<p\leq 1$ needs to be in effect, we conclude that by choosing $\alpha:=\lfloor 1/p\rfloor+1$ we obtain $\mathcal{O}(n^{-1/p})$ convergence with an implied constant independent of $s$. If $p=1$, then we assume additionally that (3.6) holds.

3.2 CBC construction with SPOD weights

We describe the CBC construction of lattice rules for SPOD weights of the general form

[TABLE]

which is specified by a smoothness degree $\sigma\in{\mathbb{N}}$, a sequence $(\Gamma_{\ell})_{\ell\geq 0}$, plus a sequence $(\gamma_{j,m})_{j\geq 1}$ for every $m=1,\ldots,\sigma$. Note that for ${\mathrm{\mathfrak{u}}}=\varnothing$, we use the convention that the empty product is one, and we interpret the sum over $\boldsymbol{m}_{\varnothing}$ as a sum with a single term ${\boldsymbol{0}}$, so that $\gamma_{\varnothing}=\Gamma_{0}$ (which in turn is typically set to $1$).

The choice of weights (3.5) corresponds to the specific case $\sigma=\alpha$, $\Gamma_{\ell}=\ell!$, and $\gamma_{j,m_{j}}=b_{j}^{m_{j}}S(\sigma,m_{j})$. We consider a generic search criterion $P({\boldsymbol{z}})$ which takes the same form as (3.4) but with a generic function $\omega\!:[0,1]\to\mathbb{R}$. Substituting in the weights, we can write

[TABLE]

Next we find a recursive definition for $p_{s,\ell}(k)$. By considering whether or not $\nu_{s}$ is zero, we can write

[TABLE]

Thus we have

[TABLE]

Note that the first term in (3.2) is exactly the value of $P(z_{1},\ldots,z_{s-1})$ in the first $s-1$ dimensions, but this is irrelevant for the construction.

Let ${\mathbb{Z}}_{n}:=\{0,1,\ldots,n-1\}$ denote the set of the integers modulo $n$, and let $\mathbb{U}_{n}:=\{u\in{\mathbb{Z}}_{n}:\gcd(u,n)=1\}$ denote the multiplicative group of integers modulo $n$ with $|\mathbb{U}_{n}|=\varphi_{\rm tot}(n)$. We define the matrix

[TABLE]

and the vectors

[TABLE]

where the entries $p_{s,\ell}(k)$ are defined recursively by (3.2) together with $p_{s,0}(k):=1$ for all $k$.

At step $s$, we see from (3.2) that the CBC algorithm should pick the value of $z_{s}\in\mathbb{U}_{n}$ which corresponds to the smallest entry in the matrix-vector product

[TABLE]

Then it is clear from (3.2) that the vectors ${\boldsymbol{p}}_{s,\ell}$ for the next iteration can be obtained recursively via

[TABLE]

where ${\boldsymbol{\Omega}}_{n}(z_{s})$ denotes the row of ${\boldsymbol{\Omega}}_{n}$ corresponding to the chosen $z_{s}$, and the operator $.*$ denotes the elementwise vector multiplication. Since the vectors ${\boldsymbol{p}}_{s-1,\ell}$ are no longer needed in the next iteration, we can simply overwrite ${\boldsymbol{p}}_{s-1,\ell}$ with ${\boldsymbol{p}}_{s,\ell}$. Hence, starting with the vectors ${\boldsymbol{p}}_{0,\ell}:={\boldsymbol{1}}_{n}$ requires $\mathcal{O}(\sigma s\,n)$ storage overall.

The fast implementation is based on ordering the indices $z\in\mathbb{U}_{n}$ and $k\in{\mathbb{Z}}_{n}$ in (3.9) and (3.10) to allow fast matrix-vector multiplication using FFT; see [4, 19, 20] for details. The overall CBC construction cost is $\mathcal{O}(s\,n\log n+\sigma^{2}s^{2}\,n)$ operations.

4 Combined error analysis

The overall error of the PDE problem (2.1) is a combination of the dimension truncation error, finite element discretization error, and QMC cubature error as

[TABLE]

where $({\boldsymbol{y}}_{i})_{i=0}^{n-1}$ are QMC nodes in $[0,1]^{s}$, $u$ denotes the solution to (2.1), $u^{s}$ and $u_{h}^{s}$ denote the dimension-truncated solution and the corresponding finite element solution, and $G\!:H_{0}^{1}(D)\to\mathbb{R}$ is a bounded, linear functional.

We can combine the results of the previous sections to produce the following overall error bound.

Theorem 4.1**.**

For any ${\boldsymbol{y}}\in U$, let $u(\cdot,{\boldsymbol{y}})\in H_{0}^{1}(D)$ denote the solution to (2.1) with the source term $f\in H^{-1+t}(D)$ for some $0\leq t\leq 1$, and let $G\in H^{-1+t^{\prime}}(D)$ for some $0\leq t^{\prime}\leq 1$. Let $({\boldsymbol{y}}_{k})_{k=0}^{n-1}$ be the lattice cubature nodes in $[0,1]^{s}$ generated by the CBC construction detailed in Subsection 3.2 for any prime power $n$, and for each lattice point we solve the approximate elliptic problem (2.1) using one common finite element discretization in the domain $D$. If $p=1$, then we assume in addition that (3.6) holds. Under the assumptions (A1)–(A6), we have the combined error estimate

[TABLE]

where $0\leq t+t^{\prime}\leq 2$, $h$ denotes the mesh size of the piecewise linear finite element mesh, $C>0$ is a constant independent of $s$, $h$, $f$, and $G$, and

[TABLE]

5 Numerical experiments

We solve (2.1) in the two-dimensional physical domain $D=(0,1)^{2}$ with the source term $f({\boldsymbol{x}})=x_{2}$ and the periodic diffusion coefficient (2.2), denoted below by $a_{\rm per}({\boldsymbol{x}},{\boldsymbol{y}})$, where $\overline{a}({\boldsymbol{x}})=2$ and

[TABLE]

For the numerical experiments, we truncate the parametric dimension to $s=100$ and use a first order finite element solver to compute solutions to (2.1) numerically by using a regular finite element mesh of the square domain $D$ with the one-dimensional mesh width $h=2^{-7}$. We use lattice rules generated by the fast CBC algorithm detailed in Subsection 3.2 with

[TABLE]

nodes and choose $\sigma=\alpha=\beta\in\{2,4\}$. Moreover, all computations have been carried out using three different values for the scaling parameter $c\in\{1,0.5,0.1\}$ to allow us to vary the difficulty of the resulting integration problem. The reference solution was computed using a rank-1 lattice rule with $n=128\,021$ nodes.

In addition, we compare the convergence rates obtained in the periodic setting to the rates obtained using interlaced polynomial lattice rules generated for the problem (2.1), equipped instead with the affine diffusion coefficient

[TABLE]

which has the same mean field and covariance as the periodic field $a_{\rm per}$ when the fluctuations are chosen as in (5.1). To generate interlaced polynomial lattice rules tailored for the affine diffusion coefficient, we used the QMC4PDE toolbox [14, 15] with the interlacing factors chosen to be equal to $\beta$ and $n=2^{k}$, $k\in\{4,\ldots,16\}$. In this case, the reference solution was computed using a corresponding interlaced polynomial lattice rule with $n=2^{17}$ nodes.

The quantity of interest in the first numerical experiment is the expectation $\mathbb{E}[G(u)]$ of the linear functional

[TABLE]

where the value of this integral can be calculated exactly when the integrand is a finite element solution. The results obtained using rank-1 lattice rules for the periodic model (2.2) are displayed on the left-hand sides of Figures 3 and 4 for the decay rates $\beta=2$ and $\beta=4$, respectively. The corresponding results obtained using interlaced polynomial lattices for the affine model (5.2) are displayed on the right-hand sides of Figures 3 and 4. The expected rates of convergence are $\mathcal{O}(n^{-2})$ and $\mathcal{O}(n^{-4})$, respectively. We observe that the solutions computed using the periodic diffusion coefficient $a_{\rm per}$ appear to converge at a rate at least as good as the expected rate. When the scaling parameter is set to $c=1$, it is notable that the rank-1 lattice rules used in conjunction with the periodic model appear to outperform the solution computed using interlaced polynomial lattice rules within the affine framework. It is apparent from Figure 3 that the observed rate of convergence is actually slightly better than the expected rate. This may be attributed to the fact that the dependence of the solution $u$ to the parametric variable ${\boldsymbol{y}}$ is analytic. For the solutions obtained using the affine model in Figure 3, the interlacing factor $2$ actually acts as a bottleneck, capping the convergence rate at $\mathcal{O}(n^{-2})$. Similar numerical behavior for interlaced polynomial lattice rules has been previously reported, e.g., in [10].

In our second experiment, we consider the problem of approximating

[TABLE]

The parameters and weights used for the construction of the rank-1 lattice rules with the periodic diffusion coefficient (2.2) as well as the interlaced polynomial lattice rules generated for the affine diffusion coefficient (5.2) are exactly the same as in the first numerical experiment. The results are displayed in Figures 5 and 6. We find that the general trend of the results matches that of the first numerical experiment, with the observed rates being at least as good as the expected rates with the scaling parameters $c\in\{0.5,0.1\}$, while the results obtained for the periodic model with $c=1$ and $\beta=4$ appear to remain in the preasymptotic regime.

Remark. Since the higher order moments of the input random fields (2.2) and (5.2) are in general different, so are the corresponding solutions to the respective integration problems. Making a direct numerical comparison of the values obtained in either setting is therefore not sensible.

Conclusions

From a modeling point of view, there does not seem to be a reason to prefer an affine expansion of a random field over a periodic expansion. Yet in the context of uncertainty quantification for PDEs with uncertain coefficients, we have seen that the model chosen for the random coefficient can make all the difference between obtaining essentially linear convergence with the affine model on the one hand, and on the other hand higher order convergence with the periodic model using rank-1 lattice cubature rules for the task of approximating the response statistics of the system. Higher order convergence can also be obtained with the affine model using, for example, interlaced polynomial lattice rules, but the overwhelming simplicity of constructing rank-1 lattice cubature rules makes the periodic framework a very enticing model for solving PDE problems equipped with uncertain coefficients. We have also presented numerical experiments that assess the QMC error derived in this work, in which the results are at least as good as those for a comparable affine model with interlaced polynomial lattice rules.

Acknowledgements

We gratefully acknowledge the financial support from the Australian Research Council (DP180101356). We are also grateful to Fabio Nobile and Yoshihito Kazashi for collaboration on a related joint project that led to the inception of this paper. Frances Kuo thanks the participants from the Oberwolfach Workshop 1911 on Uncertainty Quantification for stimulating discussions about this work.