Abstract
We consider uncertainty quantification for the Poisson problem subject to domain uncertainty. For the stochastic parameterization of the random domain, we use the model recently introduced by Kaarnioja et al. (SIAM J. Numer. Anal., 2020) in which a countably infinite number of independent random variables enter the random field as periodic functions. We develop lattice quasi-Monte Carlo (QMC) cubature rules for computing the expected value of the solution to the Poisson problem subject to domain uncertainty. These QMC rules can be shown to exhibit higher order cubature convergence rates permitted by the periodic setting independently of the stochastic dimension of the problem. In addition, we present a complete error analysis for the problem by taking into account the approximation errors incurred by truncating the input random field to a finite number of terms and discretizing the spatial domain using finite elements. The paper concludes with numerical experiments demonstrating the theoretical error estimates.
Similar content being viewed by others
1 Introduction
Modeling domain uncertainty is pertinent in many engineering applications, where the shape of the object may not be perfectly known. For example, one might think of manufacturing imperfections in the shape of products fabricated by line production, or shapes which stem from inverse problems such as tomography.
If the domain perturbations are small, one can apply the perturbation technique as in, e.g., [2, 13, 17]. In that approach, which is based on Eulerian coordinates, the shape derivative is used to linearize the problem under consideration around a nominal reference geometry. By Hadamard’s theorem, the resulting linearized equations for the first order shape sensitivities are homogeneous equations posed on the nominal geometry, with inhomogeneous boundary data only. By using a first order shape Taylor expansion, one can derive tensor deterministic partial differential equations (PDEs) for the statistical moments of the problem.
The approach we consider in this work is the domain mapping approach. It transfers the shape uncertainty onto a fixed reference domain and hence reflects the Lagrangian setting. Especially, it allows us to deal with large deformations, see [27, 31] for example. Analytic dependency of the solution on the random domain mapping in the case of the Poisson equation was shown in [3, 15] and in the case of linear elasticity in [14]. Related shape holomorphy for acoustic and electromagnetic scattering problems has been shown in [19, 21] and for Stokes and Navier–Stokes equations in [4]. Mixed spatial and stochastic regularity of the solutions to the Poisson equation on random domains, required for multilevel cubature methods, has been recently proved in [16].
In order to introduce the domain mapping approach, let \((\varOmega ,{\mathcal {F}},{\mathbb {P}})\) be a probability space, where \(\varOmega \) is the set of possible outcomes, the \(\sigma \)-algebra \({\mathcal {F}}\subset 2^{\varOmega }\) is the set of all events, and \({\mathbb {P}}\) is a probability measure. We are interested in uncertainty quantification for the Poisson problem
where the domain \(D(\omega )\subset {\mathbb {R}}^d\), \(d\in \{2,3\}\), is assumed to be uncertain. In the domain mapping framework, one starts with a reference domain \({D_{\textrm{ref}}}\subset {\mathbb {R}}^d\). For example, in a manufacturing application, \({D_{\textrm{ref}}}\) might be the planned domain provided by computer-aided design with no imperfections. It is assumed that a random perturbation field
is prescribed. The Poisson problem (1.1) on the perturbed domain is then transformed onto the fixed reference domain \({D_{\textrm{ref}}}\) by a change of variable. This results in a PDE on the reference domain, equipped with a random coefficient and a random source term, which are correlated by the random deformation, and complicated by the occurrence of the Jacobian of the transformation. Note that the dependence on the random perturbation field is nonlinear. Nonetheless, in the domain mapping approach, the reference domain \({D_{\textrm{ref}}}\) needs only to be Lipschitz, while the random perturbations \({{\varvec{V}}}(\omega )\) have to be \({\mathcal {C}}^2\)-smooth. This is in contrast to the perturbation approach mentioned above, where both the reference domain and the random perturbations need to be \({\mathcal {C}}^2\)-smooth.
The approach adopted in this paper is first to truncate the initially infinite number of scalar random variables to a finite (but possibly large) number; then to approximate the transformed and truncated PDE on the reference domain by a finite element method; and finally to approximate the expected value of the solution (which is an integral, possibly high dimensional) by a carefully designed quasi-Monte Carlo (QMC) method—i.e., an equal weight cubature rule.
In [3, 15], the perturbation field was represented by an affine, vector-valued Karhunen–Loève expansion, under the assumption of uniformly independent and identically distributed random variables. In this work, we instead expand the perturbation field \({{\varvec{V}}}(\omega )\) using periodic random variables, following the recent work [23]. The advantage of the periodic setting is that it permits higher order convergence of the QMC approximation to the expectation of the solution to (1.1). The periodic representation is equivalent in law to an affine representation of the input random field, where the random variables entering the series are i.i.d. with the Chebyshev probability density \(\rho (z)=\frac{1}{\pi \sqrt{1-z^2}}\), \(z\in (-1,1)\). The density \(\rho \) is associated with Chebyshev polynomials of the first kind, which is a popular family of basis functions in the method of generalized polynomial chaos (gPC). The choice of density is a modeling assumption, and it might be argued that in many applications choosing the Chebyshev density over the uniform density is a matter of taste rather than conviction, see [22, Sections 1–2] for more discussion.
The paper is organized as follows. The problem is formulated in Sect. 2, after which Sect. 3 establishes the regularity of the solution with respect to the stochastic variables; a more difficult task than in most such applications because of the nonlinear nature of the transformation, but a prerequisite for the design of good QMC rules. Section 4 considers error analysis, estimating separately the errors from dimension truncation, finite element approximation and QMC cubature. In the latter case, the error analysis leads to the design of appropriate weight parameters for the function space, ones that allow rigorous error bounds for the QMC contribution to the error. Numerical experiments that test the theoretical analysis are presented in Sect. 5, and some conclusions are drawn in Sect. 6. The Appendix A contains several technical results required for the parametric regularity analysis.
2 Problem formulation
2.1 Notations
Let \(U:={[0,1]^{\mathbb {N}}}\) denote a set of parameters. Let us fix a bounded domain \({D_{\textrm{ref}}}\subset {\mathbb {R}}^d\) with Lipschitz boundary and \(d\in \{2,3\}\) as the reference domain.
We use multi-index notation throughout this paper. The set of all finitely supported multi-indices is denoted by \({\mathscr {F}}:=\{{\varvec{\nu }}\in {\mathbb {N}}_0^{{\mathbb {N}}}:|\textrm{supp}({\varvec{\nu }})|<\infty \}\), where \(\textrm{supp}({\varvec{\nu }}):=\{j\in {\mathbb {N}}:\nu _j\ne 0\}\) is the support of a multi-index \({\varvec{\nu }}:=(\nu _1,\nu _2,\ldots )\). For any sequence \({{\varvec{b}}}=(b_j)_{j=1}^\infty \) of real numbers and any \({\varvec{\nu }}\in {\mathscr {F}}\), we define
where the product over the empty set is defined as 1 by convention. We also define the multi-index notation
for higher order partial derivatives with respect to variable \({{\varvec{y}}}\).
For any matrix M, let \(\sigma (M)\) denote the set of all singular values of M and let \(\Vert M\Vert _2:={\max \sigma (M)}\) denote the matrix spectral norm. For a function v on \({D_{\textrm{ref}}}\) we define
where we apply the vector 2-norm (Euclidean norm) or matrix 2-norm (spectral norm) depending on whether \(v({{\varvec{x}}})\) is a vector or a matrix. Similarly, we define
where \(\nabla v({{\varvec{x}}})\) is the gradient vector if \(v({{\varvec{x}}})\) is a scalar, and \(v'({{\varvec{x}}})\) is the Jacobian matrix if \(v({{\varvec{x}}})\) is a vector. We will also need the standard Sobolev norm
for a scalar function \(v:{D_{\textrm{ref}}}\rightarrow {\mathbb {R}}\) in \(H_0^1({D_{\textrm{ref}}})\). We also define the norm
where \(D\subset {\mathbb {R}}^d\) is a nonempty domain.
In this paper we will make use of Stirling numbers of the second kind defined by
for integers \(n\ge m\ge 0\), except for \(S(0,0):=1\).
2.2 Parameterization of domain uncertainty
Let \({{\varvec{V}}}\!:{\overline{{D_{\textrm{ref}}}}}\times U\rightarrow {\mathbb {R}}^d\) be a vector field such that
with stochastic fluctuations \({\varvec{\psi }}_i: {D_{\textrm{ref}}}\rightarrow {\mathbb {R}}^d\). Denoting the Jacobian matrix of \({\varvec{\psi }}_i\) by \({\varvec{\psi }}_i'\), the Jacobian matrix \(J: {D_{\textrm{ref}}}\rightarrow {\mathbb {R}}^{d\times d}\) of vector field \({{\varvec{V}}}\) is
We assume that the family of admissible domains \(\{D({{\varvec{y}}})\}_{{{\varvec{y}}}\in U}\) is parameterized by
and define the hold-all domain by setting
The convention of explicitly writing down the factor \(1/\sqrt{6}\) in (2.1) is not a crucial part of the analysis, but it ensures that the mean and covariance of the vector field \({{\varvec{V}}}\) match those of an affine and uniform parameterization for representing domain uncertainty using the same sequence of stochastic fluctuations \(({\varvec{\psi }}_i)_{i=1}^\infty \), with the uniform random variables supported on \([-1/2,1/2]\). See also the discussion in [23].
In order to ensure that the aforementioned domain parameterization is well-posed and to enable our subsequent analysis of dimension truncation and finite element errors, we state the following assumptions for later use:
-
(A1)
For each \({{\varvec{y}}}\in U\), \({{\varvec{V}}}(\cdot ,{{\varvec{y}}})\!:{\overline{{D_{\textrm{ref}}}}}\rightarrow {\mathbb {R}}^d\) is an invertible, twice continuously differentiable vector field.
-
(A2)
For some \(C>0\), there holds
$$\begin{aligned}\Vert {{\varvec{V}}}(\cdot ,{{\varvec{y}}})\Vert _{{\mathcal {C}}^2({\overline{{D_{\textrm{ref}}}}})}\le C\quad \text {and}\quad \Vert {{\varvec{V}}}^{-1}(\cdot ,{{\varvec{y}}})\Vert _{{\mathcal {C}}^2(\overline{D({{\varvec{y}}})})}\le C\end{aligned}$$for all \({{\varvec{y}}}\in U\).
-
(A3)
There exist constants \(0< \sigma _{\min } \le 1 \le \sigma _{\max }<\infty \) such that
$$\begin{aligned} \sigma _{\min }\le \min \sigma (J({{\varvec{x}}},{{\varvec{y}}})) \le \max \sigma (J({{\varvec{x}}},{{\varvec{y}}})) \le \sigma _{\max }~\text {for all}~{{\varvec{x}}}\in {D_{\textrm{ref}}},~{{\varvec{y}}}\in U, \end{aligned}$$where \(\sigma (J({{\varvec{x}}},{{\varvec{y}}}))\) denotes the set of all singular values of the Jacobian matrix.
-
(A4)
There holds \(\Vert {\varvec{\psi }}_i\Vert _{W^{1,\infty }({D_{\textrm{ref}}})}<\infty \) for all \(i\in {\mathbb {N}}\), and
$$\begin{aligned} \sum _{i=1}^\infty \Vert {\varvec{\psi }}_i\Vert _{W^{1,\infty }({D_{\textrm{ref}}})}<\infty . \end{aligned}$$ -
(A5)
For some \(p\in (0,1)\), there holds
$$\begin{aligned} \sum _{i=1}^\infty \Vert {\varvec{\psi }}_i\Vert _{W^{1,\infty }(D_{\textrm{ref}})}^p<\infty . \end{aligned}$$ -
(A6)
\(\Vert {\varvec{\psi }}_1\Vert _{W^{1,\infty }({D_{\textrm{ref}}})}\ge \Vert {\varvec{\psi }}_2\Vert _{W^{1,\infty }({D_{\textrm{ref}}})}\ge \cdots \).
-
(A7)
The reference domain \(D_{\textrm{ref}}\subset {\mathbb {R}}^d\) is a convex, bounded polyhedron.
For later convenience we define the sequence \({{\varvec{b}}}=(b_i)_{i=1}^\infty \) and the constant \(\xi _{{{\varvec{b}}}} \ge 0\) by setting
It follows that we can take \(\sigma _{\max }:= 1 + \xi _{{{\varvec{b}}}}\).
Remark. Under assumption (A3), there holds
This follows from the continuity of the determinant and \(\det J({{\varvec{x}}},{\textbf{0}})=1\).
2.3 The variational formulation on the reference domain
The variational formulation of the model problem (1.1) can be stated as follows: for \({{\varvec{y}}}\in U\), find \(u(\cdot ,{{\varvec{y}}})\in H_0^1(D({{\varvec{y}}}))\) such that
where \(f\in {\mathcal {C}}^\infty ({\mathcal {D}})\) is assumed to be an analytic function.
We can transport the variational formulation (2.4) to the reference domain by a change of variable. Let us first define the matrix-valued function
the transport coefficient
and the transported source term
Then we can recast the problem (2.4) on the reference domain as follows: for \({{\varvec{y}}}\in U\), find \({\widehat{u}}(\cdot ,{{\varvec{y}}})\in H_0^1({D_{\textrm{ref}}})\) such that
for all \({\widehat{v}}\in H_0^1({D_{\textrm{ref}}})\).
The solutions to problems (2.4) and (2.7) are connected to one another by
In the sequel, we focus on analyzing the problem (2.7).
3 Parametric regularity of the solution
In order to develop higher order rank-1 lattice cubature rules for the purpose of integrating the solution \({{\varvec{y}}}\mapsto {\widehat{u}}(\cdot ,{{\varvec{y}}})\) of (2.7) with respect to the parametric variable, we need to derive bounds on the partial derivatives \(\partial _{{\varvec{y}}}^{{\varvec{\nu }}}{\widehat{u}}(\cdot ,{{\varvec{y}}})\) in the Sobolev norm \(\Vert \cdot \Vert _{H_0^1({D_{\textrm{ref}}})}\).
The analysis presented in this section is based on [15] and proceeds as follows:
-
We derive bounds for \(\partial _{{\varvec{y}}}^{\varvec{\nu }}B^{-1}({{\varvec{x}}},{{\varvec{y}}})\) and \(\partial _{{\varvec{y}}}^{\varvec{\nu }}\det J({{\varvec{x}}},{{\varvec{y}}})\) for all \({\varvec{\nu }}\in {\mathscr {F}}\) in Lemmata 3.1 and 3.2. These results are then used to derive a bound on the partial derivatives of the coefficient \(\partial _{{\varvec{y}}}^{\varvec{\nu }}A({{\varvec{x}}},{{\varvec{y}}})\) in Lemma 3.3.
-
Since the right-hand side of (2.7) depends on the parametric variable \({{\varvec{y}}}\in U\), we develop the derivative bounds on \(\partial _{{\varvec{y}}}^{{\varvec{\nu }}}f_{\textrm{ref}}({{\varvec{x}}},{{\varvec{y}}})\) in Lemma 3.5, with the aid of the derivative bound developed for \(\partial _{{\varvec{y}}}^{\varvec{\nu }}{\widehat{f}}({{\varvec{x}}},{{\varvec{y}}})\) in Lemma 3.4.
-
An a priori bound is developed for the the solution \({\widehat{u}}({{\varvec{x}}},{{\varvec{y}}})\) of (2.7) in Lemma 3.6.
-
Finally, the main result of this section is stated in Theorem 3.7 which contains the partial derivative bound for \(\partial _{{\varvec{y}}}^{\varvec{\nu }}{\widehat{u}}({{\varvec{x}}},{{\varvec{y}}})\), \({\varvec{\nu }}\in {\mathscr {F}}\), making use of the aforementioned results.
3.1 Parametric regularity of the transport coefficient
Lemma 3.1
Let \({{\varvec{y}}}\in U\) and \({\varvec{\nu }}\in {\mathscr {F}}\). Under assumptions (A1)– (A4), the matrix-valued function B defined by (2.5) satisfies
with the sequence
Proof
We will first derive a regularity bound for \(\Vert \partial _{{{\varvec{y}}}}^{{\varvec{\nu }}}B(\cdot ,{{\varvec{y}}})\Vert _{L^\infty (D_{\textrm{ref}})}\), \({\varvec{\nu }}\in {\mathscr {F}}\). The bound for \(\Vert \partial _{{{\varvec{y}}}}^{{\varvec{\nu }}}B(\cdot ,{{\varvec{y}}})^{-1}\Vert _{L^\infty (D_{\textrm{ref}})}\) follows by applying implicit differentiation to the identity \(B^{-1}B=I\).
From the definition of J it is easy to see that for any multi-index \({{\varvec{m}}}\in {\mathscr {F}}\) and any \({{\varvec{x}}}\in {D_{\textrm{ref}}}\) and \({{\varvec{y}}}\in U\) we have
Taking the matrix 2-norm, we obtain
and hence with (2.2) we obtain
Now the Leibniz product rule yields for \({\varvec{\nu }}\in {\mathscr {F}}\),
and after taking the matrix 2-norm on both sides and then the \(L^\infty \)-norm over \({{\varvec{x}}}\) we obtain
The bounds in (3.3) indicate that only the derivatives with \(|\textrm{supp}({\varvec{\nu }})|\le 2\) will survive. Indeed, for \({\varvec{\nu }}= k\,{{\varvec{e}}}_j\) with \(k\ge 1\), only \({{\varvec{m}}}= \ell \,{{\varvec{e}}}_j\) for \(\ell =0,\ldots , k\) remain, while for \({\varvec{\nu }}= k\,{{\varvec{e}}}_j + k'{{\varvec{e}}}_{j'}\) with \(k,k'\ge 1\) and \(j\ne j'\), only \({{\varvec{m}}}= k\,{{\varvec{e}}}_j\) and \({{\varvec{m}}}= k'{{\varvec{e}}}_{j'}\) remain. On the other hand, for \({\varvec{\nu }}\) with \(|\textrm{supp}({\varvec{\nu }})|\ge 3\), in each term at least one of the factors must vanish. We conclude that
where we used \(b_j\le \xi _{{{\varvec{b}}}} \) for all \(j\in {\mathbb {N}}\).
Our assumption (A3) immediately yields for all \({{\varvec{y}}}\in U\) that
which proves the lemma for the case \({\varvec{\nu }}= {{\varvec{0}}}\). Now we consider \({\varvec{\nu }}\in {\mathscr {F}}\setminus \{{{\varvec{0}}}\}\). To obtain the derivative bounds on \(B^{-1}\) we start with \(B^{-1}B = I\) and apply the Leibniz product rule to obtain
Separating out the \({{\varvec{m}}}= {{\varvec{0}}}\) term, post-multiplying both sides by \(B^{-1}({{\varvec{x}}},{{\varvec{y}}})\), taking the matrix 2-norm followed by the \(L^\infty \)-norm over \({{\varvec{x}}}\), and applying (3.4), we obtain
The assertion follows by applying Lemmata A.1 and A.2 together with \(\varUpsilon _{{\varvec{\nu }}}=\big \Vert \partial ^{{\varvec{\nu }}}_{{\varvec{y}}}B^{-1}(\cdot ,{{\varvec{y}}})\big \Vert _{L^\infty ({D_{\textrm{ref}}})}\), \(c_0 = 1/\sigma ^2_{\min }\), \(c = 4\pi (1+\xi _{{{\varvec{b}}}})/\sigma _{\min }^2\), \(q = 2\), and \(\beta _j=b_j\). \(\square \)
Lemma 3.2
Under assumptions (A1)– (A4), there holds for all \({{\varvec{y}}}\in U\) and \({\varvec{\nu }}\in {\mathscr {F}}\) that
Proof
Let \({{\varvec{y}}}\in U\). The proof is carried out by induction over all the minors (subdeterminants) of the matrix \({J}({{\varvec{x}}},{{\varvec{y}}})=:[J_{\ell ,\ell '}({{\varvec{x}}},{{\varvec{y}}})]_{\ell ,\ell '=1}^d\), \({{\varvec{x}}}\in {D_{\textrm{ref}}}\). For any \(q\le d\), let \(J^{{\varvec{\ell }},{\varvec{\ell }}'}({{\varvec{x}}},{{\varvec{y}}}):= [J_{\ell _t,\ell '_{t'}}({{\varvec{x}}},{{\varvec{y}}})]_{t,t'=1}^q\) denote the \(q\times q\) submatrix specified by the indices \({\varvec{\ell }}:= \{\ell _1,\ldots ,\ell _q\}\) and \({\varvec{\ell }}':=\{\ell '_1,\ldots ,\ell '_q\}\) with \(1\le \ell _1<\cdots <\ell _q\le d\) and \(1\le \ell '_1<\cdots <\ell '_q\le d\). We will prove by induction on q that
Similarly to (3.2), we have for the derivatives of matrix elements
for \(t,t'\in \{1,\ldots ,q\}\). Using \(\big | [{\varvec{\psi }}_i'({{\varvec{x}}})]_{\ell _t,\ell '_{t'}}\big | \le \Vert {\varvec{\psi }}_i'({{\varvec{x}}})\Vert _2\) and (2.2), we obtain for the \((1\times 1)\)-minors
Suppose (3.5) holds for all submatrices of size up to \((q-1)\times (q-1)\), and now we consider the case for a \(q\times q\) submatrix. For arbitrary \(t'\le q\), the Laplace cofactor expansion yields
The Leibniz product rule then gives
Together with (3.6) and the induction hypothesis (3.5) we obtain
To simplify the last term, we obtain in complete analogy with the derivation presented in [23, Lemma 2.2] the identity
Plugging this expression into (3.7) and collecting terms yields (3.5), which in turn completes the proof of the lemma. \(\square \)
Lemma 3.3
Under assumptions (A1)– (A4), there holds for all \({{\varvec{y}}}\in U\) and \({\varvec{\nu }}\in {\mathscr {F}}\) that
where the sequence \((a_k)_{k=0}^\infty \) is defined in (3.1).
Proof
Let \({\varvec{\nu }}\in {\mathscr {F}}\). By the Leibniz product rule, we have from (2.6) that
Taking the matrix 2-norm followed by the \(L^\infty \)-norm over \({{\varvec{x}}}\), and then applying Lemmata 3.1 and 3.2, we obtain
where \({\mathbb {A}}_{{\varvec{w}}}:= |{{\varvec{w}}}|!\,a_{|{{\varvec{w}}}|}\,{{\varvec{w}}}!\,{{\varvec{b}}}^{{{\varvec{w}}}}\) and \({{\mathbb {B}}}_{\varvec{\mu }}:= |{\varvec{\mu }}|!\, b^{{\varvec{\mu }}}\), and we overestimated some multiplying factors to get
The last equality is due to Lemma A.3. Using \(a_{|{{\varvec{w}}}|}\le a_{|{{\varvec{m}}}|}\) and \({{\varvec{w}}}!\le {{\varvec{m}}}!\), we have
where the last sum over \({{\varvec{w}}}\) can be rewritten as
thus completing the proof. \(\square \)
3.2 Parametric regularity of the transported source term
Lemma 3.4
Let \(f\in {\mathcal {C}}^\infty ({\mathcal {D}})\) be analytic. Then there exist constants \(C_f>0\) and \(\rho \ge 1\) such that \(\Vert \partial _{{\varvec{x}}}^{{\varvec{\nu }}}f\Vert _{L^\infty ({\mathcal {D}})}\le C_f\, {\varvec{\nu }}!\,\rho ^{|{\varvec{\nu }}|}\) for all \({\varvec{\nu }}\in {\mathbb {N}}_0^d\). Furthermore, under the assumptions (A1)– (A4), there holds for all \({{\varvec{y}}}\in U\) and \({\varvec{\nu }}\in {\mathscr {F}}\) that the derivatives of \({\widehat{f}}({{\varvec{x}}},{{\varvec{y}}})=f({{\varvec{V}}}({{\varvec{x}}},{{\varvec{y}}}))\) are bounded by
Proof
The bound \(\Vert \partial _{{\varvec{x}}}^{{\varvec{\nu }}}f\Vert _{L^\infty ({\mathcal {D}})}\le C_f\, {\varvec{\nu }}!\,\rho ^{|{\varvec{\nu }}|}\) is a consequence of the Cauchy integral formula for analytic functions of several variables (cf., e.g., [20, Theorem 2.2.1]).
Let \({{\varvec{y}}}\in U\). Trivially, we see that \(\Vert \widehat{f}(\cdot ,{{\varvec{y}}})\Vert _{L^\infty (D_{\textrm{ref}})}\le C_f\), so we may assume in the following that \({\varvec{\nu }}\in {\mathscr {F}}\setminus \{{{\varvec{0}}}\}\). We will make use of the multivariate Faà di Bruno formula [29, Theorem 3.1 and Remark 3.3]
where we define \(\kappa _{{\varvec{\nu }},{\varvec{\lambda }}}({{\varvec{x}}},{{\varvec{y}}})\) for general \({\varvec{\nu }}\in {\mathscr {F}}\) and \({\varvec{\lambda }}\in {\mathbb {Z}}^d\) in a recursive manner as follows: \(\kappa _{{\varvec{\nu }},{{\varvec{0}}}}({{\varvec{x}}},{{\varvec{y}}}):= \delta _{{\varvec{\nu }},{{\varvec{0}}}}\), \(\kappa _{{\varvec{\nu }},{\varvec{\lambda }}}({{\varvec{x}}},{{\varvec{y}}}):= 0\) if \(|{\varvec{\nu }}|<|{\varvec{\lambda }}|\) or \({\varvec{\lambda }}\not \ge {{\varvec{0}}}\) (i.e., \({\varvec{\lambda }}\) contains negative entries), and otherwise
From the definition of \({{\varvec{V}}}\) in (2.1) it is easy to see that
which yields
Thus we obtain from (3.9) the recursion
By Lemma A.4 with \(c=2\pi \), we have for all \({\varvec{\nu }}\in {\mathscr {F}}\) and \({\varvec{\lambda }}\in {\mathbb {N}}_0^d\) that
which together with (3.8) yields for \({\varvec{\nu }}\in {\mathscr {F}}\setminus \{{{\varvec{0}}}\}\) that
This yields the desired result since \(\sum _{{\varvec{\lambda }}\in \mathbb N_0^d,\,|{\varvec{\lambda }}|=|{{\varvec{m}}}|}1 = \left( {\begin{array}{c}|{{\varvec{m}}}|+d-1\\ d-1\end{array}}\right) \). The \({{\varvec{m}}}\ne {{\varvec{0}}}\) condition for the above sum can be dropped since the \({{\varvec{m}}}={{\varvec{0}}}\) term is actually zero when \({\varvec{\nu }}\ne {{\varvec{0}}}\). This allows us to then extend the formula to cover also the case \({\varvec{\nu }}={{\varvec{0}}}\). \(\square \)
Lemma 3.5
Under the assumptions of Lemma 3.4, there holds for all \({{\varvec{y}}}\in U\) and \({\varvec{\nu }}\in {\mathscr {F}}\) that
Proof
The proof follows essentially the same steps as the proof of Lemma 3.3. By the Leibniz product rule, we have
Taking the \(L^\infty \)-norm over \({{\varvec{x}}}\) and using Lemmata 3.2 and 3.4 yields
where now \({\mathbb {A}}_{{\varvec{w}}}:= \left( {\begin{array}{c}|{{\varvec{w}}}|+d-1\\ d-1\end{array}}\right) \,|{{\varvec{w}}}|!\,\rho ^{|{{\varvec{w}}}|}\,{{\varvec{b}}}^{{{\varvec{w}}}}\), \({{\mathbb {B}}}_{{\varvec{\mu }}}:= |{\varvec{\mu }}|!\,{{\varvec{b}}}^{{\varvec{\mu }}}\), and \(c_k:= C_f\, d!\,(2+\xi _{{{\varvec{b}}}})^d\, (2\pi )^k\). Using \(\rho ^{|{{\varvec{w}}}|}\le \rho ^{|{{\varvec{m}}}|}\), we have
where the last sum over \({{\varvec{w}}}\) can be rewritten as
This completes the proof. \(\square \)
3.3 Parametric regularity of the transported PDE
We have the following a priori bound.
Lemma 3.6
Under the assumptions of Lemma 3.4, there holds for all \({{\varvec{y}}}\in U\) that
where \(C_P>0\) is the Poincaré constant associated with the domain \({D_{\textrm{ref}}}\).
Proof
We take \({\widehat{v}}({{\varvec{x}}}) = {\widehat{u}}({{\varvec{x}}},{{\varvec{y}}})\) as test function in (2.7) to obtain
where we used the Cauchy–Schwarz inequality and the Poincaré inequality \(\Vert \cdot \Vert _{L^2({D_{\textrm{ref}}})}\le C_P\,\Vert \cdot \Vert _{H_0^1({D_{\textrm{ref}}})}\), with the constant \(C_P>0\) depending only on the domain \({D_{\textrm{ref}}}\).
From the definition of A in (2.6), the left-hand side of (3.11) can be written as
where we used assumptions (A3) and (2.3). We recall that \(\sigma _{\max }=1+\xi _{{{\varvec{b}}}}\), see the line below (2.2). The proof is completed by combining the upper and lower bounds, and applying the bound for \(\Vert f_{\textrm{ref}}(\cdot ,{{\varvec{y}}})\Vert _{L^\infty ({D_{\textrm{ref}}})}\) by taking \({\varvec{\nu }}= {{\varvec{0}}}\) in Lemma 3.5. \(\square \)
Theorem 3.7
Under the assumptions of Lemma 3.4, there holds for all \({{\varvec{y}}}\in U\) and \({\varvec{\nu }}\in {\mathscr {F}}\setminus \{{{\varvec{0}}}\}\) that
where \(\beta _j:= (2+\sqrt{2})\max (1+\sqrt{3},\rho )\,b_j\).
Proof
We prove the theorem by induction on the order of \({\varvec{\nu }}\). The case \({\varvec{\nu }}= {{\varvec{0}}}\) holds by Lemma 3.6. Let \({{\varvec{y}}}\in U\) and \({\varvec{\nu }}\in {\mathscr {F}}\setminus \{{\textbf{0}}\}\), and suppose that the result holds for all multi-indices of order less than \(|{\varvec{\nu }}|\).
We differentiate both sides of (2.7) and apply the Leibniz product rule to obtain for all \({\widehat{v}}\in H_0^1({D_{\textrm{ref}}})\) that
Taking \({\widehat{v}}({{\varvec{x}}}) = \partial ^{\varvec{\nu }}_{{\varvec{y}}}{\widehat{u}}({{\varvec{x}}},{{\varvec{y}}})\) and separating out the \({{\varvec{m}}}={\varvec{\nu }}\) term, we obtain
Using similar steps as in the proof of Lemma 3.6, we then arrive at
Combining this with Lemmata 3.3 and 3.5, we obtain
With some further estimations, this can be expressed in the form of the recursive bound in Lemma A.5 with, using \(a_k\le (1+\sqrt{3})^k\) (which follows easily from (3.1)),
In particular, we used \(|{\varvec{\nu }}-{{\varvec{m}}}|\ge 1\) for the exponent in \((\frac{4\pi (1+\xi _{{{\varvec{b}}}})}{\sigma _{\min }^2})^{|{\varvec{\nu }}-{{\varvec{m}}}|}\), which allowed us to enlarge the constant as defined in c.
We then conclude from Lemmata A.5, A.7 and A.8 that for \({\varvec{\nu }}\ne {{\varvec{0}}}\)
Redefining \(\beta _j\) to incorporate the factor \(2+\sqrt{2}\) completes the proof. \(\square \)
4 Error analysis
In practice, solving (2.7) is only possible after truncating the input random field (2.1) to finitely many terms and discretizing the spatial domain using, e.g., finite elements. In this section, we analyze the errors incurred by dimension truncation, finite element discretization, and QMC cubature for the expected value of the dimensionally-truncated finite element solution of (2.7). We also discuss the construction of QMC cubatures with higher order convergence rates independently of the dimension, and present a combined error estimate for the problem.
4.1 Dimension truncation error
We are interested in finding a rate for
where \({{\widehat{u}}}\) denotes the solution to (2.7) as before and the dimensionally-truncated PDE solution is defined by
In order to derive a dimension-truncation error bound for this problem, we employ the Taylor series approach used in [9, 12] together with a change of variable technique. (The Neumann series approach used in, e.g., [7] is not applicable here because the transported problem cannot be recast in affine parametric form.)
Theorem 4.1
Suppose that the assumptions (A1)– (A6) hold. Then the dimensionally-truncated solution to (2.7) satisfies the error bounds
and
where the implied coefficients are independent of the dimension s.
Proof
We start by defining
where the fluctuations \((\Vert \varvec{\psi }_i\Vert _{W^{1,\infty }})_{i\ge 1}\in \ell ^p\) for the same p as in (A5). Let us consider an auxiliary PDE problem
where the domain realizations \({{\widetilde{D}}}({{\varvec{y}}}):=\widetilde{{\varvec{V}}}(D_{\textrm{ref}},{{\varvec{y}}})\) are defined instead via an affine parametric random perturbation field
with \((\Vert \widetilde{\varvec{\psi }}_i\Vert _{W^{1,\infty }})_{i\ge 1}\in \ell ^p\) for the same p as in (A5). The random field \(\widetilde{{\varvec{V}}}\) satisfies the hypotheses of [15, Theorem 5] (note that the random variables in [15] are supported on \([-1,1]\) instead of \([-1/2,1/2]\)) and the transported PDE solution corresponding to (4.1) satisfies
for some \({{\widetilde{{{\varvec{b}}}}}}\in \ell ^p\), where the p is the same as above.
Let \(G\in H^{-1}(D_{\textrm{ref}})\) be arbitrary and define
Clearly,
Fix \({{\varvec{y}}}\in [-1,1]^{{\mathbb {N}}}\). Then developing the Taylor series of \({{\widetilde{F}}}\) about the point \(({{\varvec{y}}}_{\le s},{\textbf{0}}):=(y_1,\ldots ,y_s,0,0,\ldots )\) yields
Let \(k\ge 2\) be an undetermined integer which will be specified later. By defining \(\varvec{\theta }({{\varvec{y}}}):=(\sin (2\pi y_i))_{i\ge 1}\), we observe that
where \({\widehat{u}}_{\textrm{per}}(\cdot ,{{\varvec{y}}})\) is the smooth periodic extension of the solution to (2.7) for all \({{\varvec{y}}}\in [-1,1]^{{\mathbb {N}}}\). Let us define \(F({{\varvec{y}}}):=\langle G,{{\widehat{u}}}_{\textrm{per}}(\cdot ,{{\varvec{y}}})\rangle _{H^{-1}(D_{\textrm{ref}}),H_0^1(D_{\textrm{ref}})}\). We obtain by a simple change of variable that
Since the last equality holds pointwise for all \({{\varvec{y}}}\in [-1,1]^{{\mathbb {N}}}\), we can integrate the equation on both sides over \(U\subset [-1,1]^{{\mathbb {N}}}\). We note that the integral \(\int _{0}^1 \sin (2\pi y_i)^{\nu _i}\,\textrm{d}y_i\) takes the value between 0 and 1, and gives the value 0 when \(\nu _i = 1\). Thus the summands corresponding to \({\varvec{\nu }}\in {\mathscr {F}}\) with \(\nu _i=1\) for at least one \(i>s\) vanish in the first term. Hence
where \(C'=C\,k!\) and \(C''=C\,(k+1)!\). Recalling the definition of F and taking the supremum over all \(G\in H^{-1}(D_\textrm{ref})\) with \(\Vert G\Vert _{H^{-1}(D_{\textrm{ref}})}\le 1\) yields for the left-hand side that
which yields the overall bound
where the generic constant is independent of s. The goal will be to balance \(T_1\) and \(T_2\) in (4.2) by choosing the value of k appropriately.
The term \(T_1\) in (4.2) can be bounded similarly to [7, Theorem 1]:
where we abuse notation slightly by interpreting the value of the term \(\frac{1-{{\widetilde{b}}}_i^{k-1}}{1-{{\widetilde{b}}}_i}\) as \(k-1\) if \(\widetilde{b}_i=1\), and since the sequence is \({{\widetilde{{{\varvec{b}}}}}}\) is non-increasing we defined \(\tau _k:= \frac{1-{\widetilde{b}}_1^{k-1}}{1-{\widetilde{b}}_1}\) if \({{\widetilde{b}}}_1\ne 1\) and \(\tau _k:=k-1\) if \({{\widetilde{b}}}_1=1\). Moreover, it is an easy consequence of (A6) that
where we used the well-known Stechkin’s lemma (cf., e.g., [25, Lemma 3.3]): for \(0<p\le q<\infty \) and any sequence \((a_i)_{i\ge 1}\) of real numbers ordered such that \(|a_1|\ge |a_2|\ge \cdots \), there holds
Therefore
The term \(T_2\) in (4.2) can be estimated similarly to the approach taken in [9, Theorem 4.1]. The trivial inequality \(\frac{|{\varvec{\nu }}|!}{{\varvec{\nu }}!}\ge 1\) and the multinomial theorem imply together that
where the final inequality is a consequence of (4.3). The exponent of s in the bound on \(T_2\) is dominated by that of \(T_1\) by choosing \(k=\lceil 1/(1-p)\rceil \).
Finally, the Poincaré inequality implies that
and the \(L^1\) error can be obtained by using the Cauchy–Schwarz inequality
This completes the proof. \(\square \)
4.2 Finite element error
Let the assumptions (A1)–(A4) and (A7) hold. Let us denote by \(\{V_h\}_h\) a family of finite element subspaces of \(H_0^1(D_{\textrm{ref}})\) indexed by the mesh size h. We assume that the finite element spaces are spanned by continuous, piecewise linear finite element basis functions and the mesh corresponding to each \(V_h\) is obtained from an initial, regular triangulation of \(D_{\textrm{ref}}\) by recursive uniform bisection of simplices.
We define the dimensionally-truncated finite element solution for \({{\varvec{y}}}\in U\) as the solution \(u_{s,h}(\cdot ,{{\varvec{y}}})\in V_h\) to
for all \({{\widehat{v}}}_h\in V_h\), where we set \(A_s({{\varvec{x}}},{{\varvec{y}}}):=A({{\varvec{x}}},(y_1,\ldots ,y_s,0,0,\ldots ))\) as well as \(f_{\textrm{ref,s}}({{\varvec{x}}},{{\varvec{y}}}):=f_{\textrm{ref}}({{\varvec{x}}},(y_1,\ldots ,y_s,0,0,\ldots ))\) for all \({{\varvec{y}}}\in U\). Since the reference domain \(D_{\textrm{ref}}\) was assumed in (A7) to be a bounded, convex polyhedron, it follows that elements in \(W^{1,\infty }(D_{\textrm{ref}})\) have Lipschitz continuous representatives. Hence we may apply Lemma 3.5 and [10, Theorem 3.2.1.2] to obtain
where \(\Vert {{\widehat{v}}}\Vert _{H^2(D_{\textrm{ref}})\cap H_0^1(D_{\textrm{ref}})}:=(\Vert {{\widehat{v}}}\Vert _{L^2(D_{\textrm{ref}})}^2+\Vert \varDelta {{\widehat{v}}}\Vert _{L^2(D_{\textrm{ref}})}^2)^{1/2}\) and the constant \(C>0\) only depends on d and the diameter of \(D_{\textrm{ref}}\).
We have by Céa’s lemma that
for a constant \(C'>0\) independent of s, h, and \({{\varvec{y}}}\), which can be seen by letting \({{\widehat{v}}}_h\in V_h\) be arbitrary and deriving
where the penultimate step follows from Galerkin orthogonality.
The well-known approximation property (cf., e.g., [8, 26])
holds for all \({{\widehat{v}}}\in H^2(D_{\textrm{ref}})\cap H_0^1(D_{\textrm{ref}})\) with a constant \(C''>0\) independent of h. The inequalities (4.4), (4.5), and (4.6) yield altogether that
where the constant \(C'''>0\) is independent of s, h, and \({{\varvec{y}}}\).
An \(L^1\) error estimate can be obtained from a standard Aubin–Nitsche duality argument.
Theorem 4.2
Under assumptions (A1)– (A4) and (A7), there holds for the dimensionally-truncated finite element solution that
where the generic constant is independent of s, h, and \({{\varvec{y}}}\).
Proof
Since
it suffices to prove the \(L^2\) error estimate.
Let \({{\widehat{g}}}\in L^2(D_{\textrm{ref}})\). For \({{\varvec{y}}}\in U\), denote by \({{\widehat{u}}}_{{{\widehat{g}}},s}(\cdot ,{{\varvec{y}}})\in H_0^1(D_{\textrm{ref}})\) the solution to
Testing against \({{\widehat{v}}}={{\widehat{u}}}_s(\cdot ,{{\varvec{y}}})-{{\widehat{u}}}_{s,h}(\cdot ,{{\varvec{y}}})\) and letting \({{\widehat{v}}}_h\in V_h\) be arbitrary, it follows from Galerkin orthogonality of the finite element solution that
and, in consequence,
Then there holds for all \({{\varvec{y}}}\in U\) that
Plugging this inequality and (4.7) into (4.8) yields the desired result. \(\square \)
4.3 Quasi-Monte Carlo cubature error and a choice of weights
Let us consider QMC cubature over the s-dimensional unit cube \(U_s:=[0,1]^s\) for finite \(s\in {\mathbb {N}}\). An n-point QMC rule is an equal weight cubature rule of the form
where the cubature nodes \({{\varvec{y}}}^{(1)},\ldots ,{{\varvec{y}}}^{(n)}\in [0,1]^s\) are deterministically chosen. In the subsequent analysis, we consider the family of rank-1 lattice rules, where the QMC nodes in \([0,1]^s\) are defined by
where \({\varvec{z}}\in \{1,\ldots ,n-1\}^s\) is called the generating vector.
First we briefly recall the theory of lattice rules for 1-periodic integrands with absolutely convergent Fourier series: for
the lattice rule error is given by [30]
where the dual lattice \(\varLambda ^\perp :=\{{\varvec{\ell }}\in {\mathbb {Z}}^s: {\varvec{\ell }}\cdot {{\varvec{z}}}\equiv 0~(\textrm{mod}~n)\}\) depends on the lattice generating vector \({{\varvec{z}}}\). For integrands belonging to the weighted Korobov space endowed with the norm
where \(\alpha >1\) is the smoothness parameter, \({\varvec{\gamma }}:=(\gamma _{\mathrm {{\mathfrak {u}}}})_{\mathrm {{\mathfrak {u}}}\subseteq \{1:s\}}\) are positive weights, and \(\textrm{supp}({\varvec{\ell }}):= \{1\le j\le s: \ell _j\ge 0\}\), we then have the integration error bound
It is known [6, Theorem 5] that the component-by-component (CBC) algorithm can be used to find a generating vector satisfying the bound
where n was assumed to be prime in [6], but the result also generalizes to non-primes by replacing \(n-1\) by the Euler totient function.
In the following we will adapt this theory to the particular integrand \(F({{\varvec{y}}}) = {{\widehat{u}}}_{s,h}({{\varvec{x}}},{{\varvec{y}}})\) for \({{\varvec{x}}}\in {D_{\textrm{ref}}}\). We derive smoothness-driven product-and-order dependent (SPOD) weights and show that the error can be bounded independent of s; weights of this form have previously appeared in [5, 23].
Theorem 4.3
Suppose that assumptions (A1)– (A5) and (A7) hold. Let \(\alpha := \lfloor 1/p\rfloor +1\) with p from Assumption (A5), and choose the weights \(\gamma _{\varnothing }:=1\) and
where \(\widetilde{C}:=\frac{2\,d!\,(2+\xi _{{{\varvec{b}}}})^d(1+\xi _{{{\varvec{b}}}})^3}{\sigma _{\min }^{d+4}}\) and \(\beta _j:=(2+\sqrt{2})\max (1+\sqrt{3},\rho )\,b_j\). Then a lattice rule can be constructed by the CBC algorithm such that
where the implied constant is independent of s.
Proof
We will begin by treating \(\alpha \) and the weights \(\gamma _\mathrm {{\mathfrak {u}}}\) as arbitrary and will specify their values later in the proof. We will assume that \(\alpha \ge 2\) is an integer, and we will express the integration error in term of the mixed derivatives of \(F({{\varvec{y}}}) = \widehat{u}_{s,h}({{\varvec{x}}},{{\varvec{y}}})\) of order up to \(\alpha \) in each coordinate. For \(\mathrm {{\mathfrak {u}}}\subseteq \{1:s\}\) we denote by \({\varvec{\nu }}= ({\varvec{\alpha }}_\mathrm {{\mathfrak {u}}};{{\varvec{0}}})\) the multi-index whose component \(\nu _j\) is equal to \(\alpha \) if \(j\in \mathrm {{\mathfrak {u}}}\) and is 0 otherwise. First we observe by differentiating the Fourier series of F that
from which we conclude that the Fourier coefficients of \(\partial ^{({\varvec{\alpha }}_\mathrm {{\mathfrak {u}}};{{\varvec{0}}})}_{{\varvec{y}}}F\) are given by
We are ready to derive the integration error for \({{\widehat{u}}}_{s,h}\). Using (4.9) we have
Denoting \(\mathrm {{\mathfrak {u}}}= \textrm{supp}({\varvec{\ell }})\), we obtain using the definition of \(r_{\alpha ,{\varvec{\gamma }}}({\varvec{\ell }})\) in (4.10) with (4.13)
where for the inequalities we used the unimodularity of the exponential function, changed the order of integration, and applied both Cauchy–Schwarz and Poincaré inequalities.
Combining (4.11), (4.14) and (4.15) together with Theorem 3.7, it follows that the CBC algorithm can be used to find a generating vector satisfying the bound
where \(\widetilde{C}:=\frac{2\,d!\,(2+\xi _{{{\varvec{b}}}})^d(1+\xi _{{{\varvec{b}}}})^3}{\sigma _{\min }^{d+4}}\), \(\beta _j=(2+\sqrt{2})\max (1+\sqrt{3},\rho )b_j\), and \(C_{\textrm{init}}\) is defined as in (3.10). We choose the weights by (4.12) to ensure that the factor on the second line is bounded by 1.
Finally, we show that with this choice of weights there is a constant \(C_{\alpha ,{\varvec{\gamma }}}(\lambda )\) independent of s such that \(C_{\alpha ,{\varvec{\gamma }}}(\lambda ,s)\le C_{\alpha ,{\varvec{\gamma }}}(\lambda )<\infty \). To this end, we plug the weights (4.12) into the expression for \(C_{\alpha ,{\varvec{\gamma }}}(\lambda ,s)\) in (4.11), define \(S_{\max }(\alpha ):=\max _{k\in \{1:\alpha \}}S(\alpha ,k)\), and use Jensen’s inequality (see, e.g., [18, Theorem 19])
to find that
where we define \(\mu _j:=\max (1,{{\widetilde{C}}}^\alpha \, \alpha !\,S_{\max }(\alpha )\,(2\zeta (\alpha \lambda ))^{1/\lambda })\,\beta _j\). By defining the auxiliary sequence \(d_j:=\beta _{\lceil j/\alpha \rceil }\), \(j\ge 1\), that is
the last sum in (4.16) is bounded from above by
where the final inequality is due to the fact that \(\big (\sum _{j=1}^\infty d_j^\lambda \big )^\ell \) includes all products of the form \(d_{j_1} d_{j_2}\ldots d_{j_\ell }\), not only those (as in the line above) with all \(j_i, \ldots , j_\ell \) different; and since the terms with all \(j_1, \ldots , j_\ell \) different are repeated \(\ell !\) times, in the last step we can divide by \(\ell !\).
By choosing \(\lambda =p\) and \(\alpha =\lfloor 1/p\rfloor +1\), it is easily seen by Assumption (A5) and an application of the d’Alembert ratio test that the upper bound in the above expression is finite. This implies that the QMC error is of order \({\mathcal {O}}(n^{-1/p})\), with the implied constant independent of the dimension s when the weights are chosen according to (4.12). \(\square \)
4.4 Overall error estimate
Let \({{\widehat{u}}}\) be the solution to (2.7), \({{\widehat{u}}}_s\) the corresponding dimensionally-truncated solution, and \(\widehat{u}_{s,h}\) the dimensionally-truncated FE solution. It is easy to see that
We combine the error bounds derived in Sects. 4.1–4.3. The error decomposition (4.17) allows us to deduce the following overall error estimate.
Theorem 4.4
Let \(({{\varvec{y}}}_k)_{k=1}^n\) be the lattice cubature nodes in \(U_s\) generated by a CBC algorithm using the construction detailed in Sect. 4.3. Suppose that for each lattice point, we solve (2.7) using one common finite element discretization in the domain \(D_{\textrm{ref}}\). Under the assumptions (A1)–(A7), we have the combined error estimate
where h denotes the mesh size of the piecewise linear finite element mesh and \(C>0\) is a constant independent of s, h, and n.
5 Numerical experiments
Let us consider the domain parameterization
where
with \(c>0\) and \(\psi _j(x):=j^{-\theta }\cos (j\pi x)\), \(\theta >1\). Thus the reference domain in this case is the unit square, and the uncertain boundary is confined to the upper edge of the square. It is possible to write \(D({{\varvec{y}}})=V(D({\textbf{0}}),{{\varvec{y}}})\) with the vector-valued expansion
where
We note that
and
The identity (5.2) can be easily verified from the definition of the matrix spectral norm, and (5.3) follows by observing that \(\cos ^2(j \pi x_1) + j^2 \pi ^2 x_2^2\sin ^2(j \pi x_1) = 1 + \sin ^2(j \pi x_1) (j^2 \pi ^2 x_2^2 - 1)\) is maximized by \(x_1=\frac{1}{2j}\) and \(x_2=1\). The elements of the sequence \((b_j)_{j=1}^\infty \) needed for the construction of the QMC weights now simplify to
Moreover, we can take
The expected rate of QMC convergence is \({\mathcal {O}}(n^{-\theta +1})\).
5.1 Source problem
We consider the Poisson problem (1.1) equipped with \(f({{\varvec{x}}}):=x_2\) and the random domain described above for \(\theta \in \{2.1,2.5,3.0\}\) and \(c=\sqrt{3/2}\). The input random field is truncated to \(s=100\) terms. We construct a uniform regular triangulation for the reference domain \(D_{\textrm{ref}}=D({\textbf{0}})\) with mesh size \(h=2^{-5}\), which is then mapped to each realization of the random domain using \({\varvec{V}}\). The PDE is solved using a first order finite element method. For each realization of \({{\varvec{y}}}\in U_s\), the finite element solution is computed in \(D({{\varvec{y}}})\) and then mapped back onto the reference domain \(D_{\textrm{ref}}\). This is equivalent to but simpler than computing the solution on the reference domain, see [15]. We solve the PDE problem over QMC point sets with increasing number of cubature nodes n and compute the relative errors
where \(Q_{s,n}\) denotes the s-dimensional QMC operator over n cubature nodes. As the reference solution, we use the approximate QMC solution obtained using \(n=1\,024\,207\) cubature nodes. The QMC point sets for this experiment were constructed using the fast CBC algorithm with SPOD weights of the form (4.12), where we set \(\sigma _{\min }=\rho =1\) and \(c=10^{-6}\) for numerical stability. The results are displayed in the left-hand side image of Fig. 1. In this case, the observed rates \(\sigma _{\theta }\) are for \(\theta = 2.1\): \(\sigma _{2.1}=0.95\), \(\theta = 2.5\): \(\sigma _{2.5}= 1.15\), \(\theta = 3\): \(\sigma _{3} = 1.72\). The theoretically predicted convergence rates are clearly observed in the numerical results: using a higher value for the decay rate \(\theta \) results in a faster rate of convergence.
In addition, we also consider a nonlinear quantity of interest (QoI)
and compute the relative errors
Again, the QMC solution obtained using \(n=1\,024\,207\) cubature nodes is used as the reference. The results, computed using the same QMC rules as above, are displayed in the right-hand side image of Fig. 1. The observed rates \(\sigma _{\theta }\) are for \(\theta = 2.1\): \(\sigma _{2.1}=0.89\), \(\theta = 2.5\): \(\sigma _{2.5}= 1.09\), \(\theta = 3\): \(\sigma _{3} = 1.74\). Even though the use of a nonlinear QoI results in more irregular convergence behavior, the general trends are still visible, i.e., faster decay of the input random field yields a faster cubature convergence rate.
5.2 Capacity problem
In this section, we consider an example which is beyond the setting considered in the previous sections to demonstrate that our method also works well for a broader class of domain uncertainty quantification problems.
We will consider the Dirichlet–Neumann problem
and its conjugate problem
As the QoI, we consider the capacity
where u is the solution to (5.6). We also define its conjugate as
where v is the solution to (5.7). It is well-known that (cf., e.g., [1, Theorem 4-5])
provided that \(x\mapsto a(x,{{\varvec{y}}})\) specifies a Jordan arc such that \(a(x,{{\varvec{y}}})>0\) for all \(x\in [0,1]\) and \({{\varvec{y}}}\in U_s\).
The solutions to (5.6) and (5.7) are calculated numerically using hp-FEM. As before, the finite element solutions to (5.6) and (5.7) are not computed on the reference domain \(D_{\textrm{ref}}\); for each realization of \({{\varvec{y}}}\in U_s\), the finite element solution is computed in \(D({{\varvec{y}}})\). Numerically, we can use the reciprocal relation (5.8) to define an a posteriori error estimate
where \(\textrm{cap}(D({{\varvec{y}}}))\) and \(\overline{\textrm{cap}}(D({{\varvec{y}}}))\) are computed using hp-FEM, and we tune the finite element approximation to ensure that (5.9) is approximately of the same order across all realizations of \(D({{\varvec{y}}})\). One realization of \(D({{\varvec{y}}})\) with respective solutions \(u({{\varvec{x}}},{{\varvec{y}}})\) and \(v({{\varvec{x}}},{{\varvec{y}}})\) is shown in Fig. 2.
Sample realisation at \(s=10\): Contour plots of \(u({{\varvec{x}}},{{\varvec{y}}})\) and \(v({{\varvec{x}}},{{\varvec{y}}})\). Combination of the contour plots is the image of the conformal map of the domain
Error criterion (5.10) with increasing n. a Varying decay rate \(\theta \in \{2.1,2.5, 3.0\}\), and fixed dimension \(s=10\). b Varying dimensions \(s\in \{10,40,100\}\), and fixed decay rate \(\theta =2.1\). c Detail of the case \(s=10\), \(\theta =2.5\)
The numerical experiments were carried out as follows:
-
(i)
We set \(s=10\), \(c=\sqrt{3/2}\), and let \(\theta \in \{2.1,2.5,3.0\}\) in (5.1). We solve the problems (5.6) and (5.7) over QMC point sets with increasing number of cubature nodes n and consider the error criteria
$$\begin{aligned} \textrm{Capacity}=\frac{|{\mathbb {E}}[\textrm{cap}(D(\cdot ))]-Q_{s,n}(\textrm{cap}(D(\cdot )))|}{|{\mathbb {E}}[\textrm{cap}(D(\cdot ))]|}, \end{aligned}$$(5.10)where \(Q_{s,n}\) denotes the s-dimensional QMC operator over n cubature points. As the reference solution, we use the approximate QMC solution computed using \(n=1\,024\,207\). (The reference result is computed using a tailored lattice. In two cases the general lattice reference is available and the effect on the convergence curves is negligible.) The results are displayed in Fig. 3a.
-
(ii)
We set \(c=\sqrt{3/2}\), \(\theta =2.1\), and let \(s\in \{10,40,100\}\). We compute the relative errors (5.10) for the problems (5.6) and (5.7) over QMC point sets with increasing number of cubature nodes n. Again, as the reference solution, we use the approximate QMC solution computed using \(n=1\,024\,207\). The results are displayed in Fig. 3b.
The QMC point sets for the above experiments were computed using the fast CBC algorithm using the same specifications as in Sect. 5.1. In order to utilize the a posteriori estimate introduced above, we remark that the boundary conditions of the numerical experiments are not the homogeneous Dirichlet zero boundary conditions analyzed in the earlier theory.
The observed rates \(\sigma _{\theta }\) are for \(\theta = 2.1\): \(\sigma _{2.1}=0.7\), \(\theta = 2.5\): \(\sigma _{2.5}= 1.3\), \(\theta = 3\): \(\sigma _{3} = 2.4\) (see Fig. 3). The rate does not appear to depend on the dimension s. The cone of convergence for \(\theta = 2.5\) shown in Fig. 3c, is bounded by the rates \(\sigma _{2.5\text {A}}= 2.3\) (above) and \(\sigma _{2.5\text {B}}= 0.6\) (below). The general convergence patterns are clearly observed. The irregularities in the convergence behavior may be explained by the nonlinear QoI (compare with the experiment in Sect. 5.1) as well as the more complicated boundary conditions: the gradient term in the expression for the capacity appears to be sensitive to obtuse angles at the Dirichlet–Neumann interface for certain realizations of the random domain.
6 Conclusions
We analyzed the Poisson problem subject to domain uncertainty, employing a domain mapping approach in which the shape uncertainty is transferred onto a fixed reference domain. As a result, the coefficient of the transported PDE problem on the reference domain becomes highly nonlinear with respect to the stochastic variables. To this end, we developed a novel parametric regularity analysis for the transported problem subject to such diffusion coefficients. The parametric regularity analysis is an important ingredient in the design of tailored QMC rules and we demonstrated that these methods exhibit higher order cubature convergence rates when the random perturbation field is parameterized using the model introduced in [23] in which an infinite sequence of independent random variables enter the random field as periodic functions. The numerical results display that faster convergence in the stochastic fluctuations leads to faster QMC convergence in the computation of the PDE response, which is consistent with the theory.
Another novel feature of our work is the dimension truncation error analysis for the problem. Many studies in uncertainty quantification for PDEs with random coefficients exploit the parametric structure of the problem, using a Neumann series expansion to obtain dimension truncation error bounds. However, there have been several recent studies suggesting a Taylor series approach which is based on the parametric regularity of the problem [9, 11, 12]. The Taylor series approach can be used to obtain dimension truncation rates for non-affine parametric PDE problems, but a limitation of the aforementioned papers is that the parametric regularity bound needs to be of product-and-order dependent (POD) form in order for the Taylor-based approach to yield useful results. This issue is circumvented in the present paper by the use of a change of variables technique in order to transform our problem into a form in which the Taylor series approach can be used to produce improved dimension truncation error rates.
The results developed in this paper constitute important groundwork for other applications involving domain uncertainty such as domain shape recovery problems.
References
Ahlfors, L.V.: Conformal Invariants: Topics in Geometric Function Theory. Higher Mathematics Series. McGraw-Hill, New York (1973)
Babuška, I., Chatzipantelidis, P.: On solving elliptic stochastic partial differential equations. Comput. Methods Appl. Mech. Eng. 191, 4093–4122 (2002)
Castrillón-Candás, J.E., Nobile, F., Tempone, R.F.: Analytic regularity and collocation approximation for elliptic PDEs with random domain deformations. Comput. Math. Appl. 71(6), 1173–1197 (2016)
Cohen, A., Schwab, C., Zech, J.: Shape holomorphy of the stationary Navier-Stokes equations. SIAM J. Math. Anal. 50(2), 1720–1752 (2018)
Dick, J., Kuo, F.Y., Le Gia, Q.T., Nuyens, D., Schwab, C.: Higher order QMC Petrov-Galerkin discretization for affine parametric operator equations with random field inputs. SIAM J. Numer. Anal. 52(6), 2676–2702 (2014)
Dick, J., Sloan, I.H., Wang, X., Woźniakowski, H.: Good lattice rules in weighted Korobov spaces with general weights. Numer. Math. 103(1), 63–97 (2006)
Gantner, R.N.: Dimension truncation in QMC for affine-parametric operator equations. In: Owen, A.B., Glynn, P.W. (eds.) Monte Carlo and Quasi-Monte Carlo Methods 2016, pp. 249–264. Springer, Stanford, CA (2018)
Gilbarg, D., Trudinger, N.S.: Elliptic Partial Differential Equations of Second Order, 2nd edn. Springer, Berlin (2001)
Gilbert, A.D., Graham, I.G., Kuo, F.Y., Scheichl, R., Sloan, I.H.: Analysis of quasi-Monte Carlo methods for elliptic eigenvalue problems with stochastic coefficients. Numer. Math. 142(4), 863–915 (2019)
Grisvard, P.: Elliptic Problems in Nonsmooth Domains. Pitman Publishing Inc, New York (1985)
Guth, P.A., Kaarnioja, V.: Generalized dimension truncation error analysis for high-dimensional numerical integration: lognormal setting and beyond (2022). Preprint at arXiv:2209.06176
Guth, P.A., Kaarnioja, V., Kuo, F.Y., Schillings, C., Sloan, I.H.: Parabolic PDE-constrained optimal control under uncertainty with entropic risk measure using quasi-Monte Carlo integration (2022). Preprint at arXiv:2208.02767
Harbrecht, H.: On output functionals of boundary value problems on stochastic domains. Math. Methods Appl. Sci. 33(1), 91–102 (2010)
Harbrecht, H., Karnaev, V., Schmidlin, M.: Quantifying domain uncertainty in linear elasticity. Tech. Rep. 2023-06, Fachbereich Mathematik, Universität Basel, Switzerland (2023)
Harbrecht, H., Peters, M., Siebenmorgen, M.: Analysis of the domain mapping method for elliptic diffusion problems on random domains. Numer. Math. 134(4), 823–856 (2016)
Harbrecht, H., Schmidlin, M.: Multilevel quadrature for elliptic problems on random domains by the coupling of FEM and BEM. Stoch. Partial Differ. Equ. Anal. Comput. 10(4), 1619–1650 (2022)
Harbrecht, H., Schneider, R., Schwab, C.: Sparse second moment analysis for elliptic problems in stochastic domains. Numer. Math. 109(3), 385–414 (2008)
Hardy, G.H., Littlewood, J.E., Pólya, G.: Inequalities. Cambridge University Press, Cambridge (1934)
Hiptmair, R., Scarabosio, L., Schillings, C., Schwab, C.: Large deformation shape uncertainty quantification in acoustic scattering. Adv. Comput. Math. 44(5), 1475–1518 (2018)
Hörmander, L.: An Introduction to Complex Analysis in Several Variables, 3rd edn. North-Holland (1990)
Jerez-Hanckes, C., Schwab, C., Zech, J.: Electromagnetic wave scattering by random surfaces: Shape holomorphy. Math. Models Methods Appl. Sci. 27(12), 2229–2259 (2017)
Kaarnioja, V., Kuo, F.Y., Sloan, I.H.: Lattice-based kernel approximation and serendipitous weights for parametric PDEs in very high dimensions. To appear in: A. Hinrichs, P. Kritzer, F. Pillichshammer (eds.). Monte Carlo and Quasi-Monte Carlo Methods (2022). Springer Verlag
Kaarnioja, V., Kuo, F.Y., Sloan, I.H.: Uncertainty quantification using periodic random variables. SIAM J. Numer. Anal. 58(2), 1068–1091 (2020)
Katana: Published online (2010). https://doi.org/10.26190/669X-A286
Kressner, D., Tobler, C.: Low-rank tensor Krylov subspace methods for parametrized linear systems. SIAM J. Matrix Anal. Appl. 32(4), 1288–1316 (2011)
Kuo, F.Y., Schwab, C., Sloan, I.H.: Quasi-Monte Carlo finite element methods for a class of elliptic partial differential equations with random coefficients. SIAM J. Numer. Anal. 50(6), 3351–3374 (2012)
Mohan, P.S., Nair, P.B., Keane, A.J.: Stochastic projection schemes for deterministic linear elliptic partial differential equations on random domains. Int. J. Numer. Methods Eng. 85(7), 874–895 (2011)
Olver, F.W.J., Olde Daalhuis, A.B., Lozier, D.W., Schneider, B.I., Boisvert, R.F., Clark, C.W., Miller, B.R., Saunders, B.V., Cohl, H.S., McClain, M.A., eds.: NIST Digital Library of Mathematical Functions. Release 1.1.6 of 2022-06-30., http://dlmf.nist.gov/
Savits, T.H.: Some statistical applications of Faa di Bruno. J. Multivariate Anal. 97(10), 2131–2140 (2006)
Sloan, I.H., Kachoyan, P.J.: Lattice methods for multiple integration: Theory, error analysis and examples. SIAM J. Numer. Anal. 24(1), 116–128 (1987)
Xiu, D., Tartakovsky, D.M.: Numerical methods for differential equations in random domains. SIAM J. Sci. Comput. 28(3), 1167–1185 (2006)
Acknowledgements
F. Y. Kuo and I. H. Sloan acknowledge the support from the Australian Research Council (DP210100831). This research includes computations using the computational cluster Katana supported by Research Technology Services at UNSW Sydney [24].
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Technical results
Technical results
Lemma A.1
Let \({\varvec{\nu }}\in {\mathscr {F}}\), \(c_0>0\), \(c>0\), and \(q\in {\mathbb {Z}}^+\). Let \({\varvec{\beta }}=(\beta _j)_{j\ge 1}\) and \((\varUpsilon _{{\varvec{\nu }}})_{{\varvec{\nu }}\in \mathscr {F}}\) be sequences of non-negative real numbers and suppose that
for all \({\varvec{\nu }}\in {\mathscr {F}}\setminus \{{{\varvec{0}}}\}\). Then
where we define, with \({{\varvec{e}}}_\mathrm {{\mathfrak {u}}}:= \sum _{j\in \mathrm {{\mathfrak {u}}}} {{\varvec{e}}}_j\),
The result is sharp in the sense that all inequalities can be replaced by equalities. The result also extends to \(q=\infty \), i.e., we can remove the upper bounds on \(|\textrm{supp}({{\varvec{m}}})|\) and \(|\mathrm {{\mathfrak {u}}}|\).
Proof
We prove the claim by induction with respect to the order of the multi-indices \({\varvec{\nu }}\in {\mathscr {F}}\). The claim is resolved immediately for \({\varvec{\nu }}={{\varvec{0}}}\), so let us fix \({\varvec{\nu }}\in {\mathscr {F}} \setminus \{{{\varvec{0}}}\}\) and suppose that the assertion holds for all multi-indices with order less than \(|{\varvec{\nu }}|\). Then
where we swapped the order of the sums over \({{\varvec{m}}}\) and \({{\varvec{w}}}\). The inner sum over \({{\varvec{m}}}\) can be rewritten as
which can be further simplified using the identity (cf., e.g., [28, formula 26.8.23])
Hence
where we swapped the order of the sums over \({{\varvec{w}}}\) and \(\mathrm {{\mathfrak {u}}}\) and then carried out the change of variable \({{\varvec{w}}}'= {{\varvec{w}}}+{{\varvec{e}}}_\mathrm {{\mathfrak {u}}}\). Changing the order of the sums again leads to
Relabeling \({{\varvec{w}}}'\) to \({{\varvec{m}}}\) and applying the definition of \(D_q({{\varvec{m}}})\) yields the desired result.
We remark that if the inequalities for \(\varUpsilon _0\) and the recursion on \(\varUpsilon _{\varvec{\nu }}\) are replaced by equalities then all steps of the proof are equalities and so is the final formula. We also remark that the proof does not make any use of the upper bounds \(|\textrm{supp}({{\varvec{m}}})|\le q\) and \(|\mathrm {{\mathfrak {u}}}|\le q\), so the result extends trivially to the case \(q=\infty \). \(\square \)
Let \({{\varvec{m}}}\in {\mathscr {F}}\). It is interesting to note that if \(|\textrm{supp}({{\varvec{m}}})| \le 2\), for example, \({{\varvec{m}}}_{\mathrm {{\mathfrak {u}}}}= k{{\varvec{e}}}_i + \ell {{\varvec{e}}}_j\) with \(k,\ell \in {\mathbb {N}}\), then the numbers (A.1) have the following representation in closed form:
These are so-called Delannoy numbers. We prove an upper bound for \(D_2({{\varvec{m}}})\) in general.
Lemma A.2
The numbers defined by (A.1) with \(q=2\) can be bounded by
with the sequences
Proof
Let us first observe that the sequence \((a_k')_{k=0}^\infty \) can be represented using the recurrence \(a_0'=a_1'=1\) and \(a_k'=ka_{k-1}'+\frac{k^2-k}{2}a_{k-2}'\), \(k\ge 2\); see, e.g., entry A080599 in the On-line Encyclopedia of Integer Sequences (OEIS) and references therein. We also have \(a_0=a_1=1\) and \(a_k = a_{k-1} + a_{k-2}/2\), \(k\ge 2\).
We prove the claim (A.2) by induction with respect to the order of the multi-indices \({{\varvec{m}}}\in {\mathscr {F}}\) while utilizing the recursive characterization of sequence \((a_k')_{k=0}^\infty \). We immediately see that (A.2) holds with equality when \(|{{\varvec{m}}}|\le 1\). Fix \({{\varvec{m}}}\in {\mathscr {F}}\) such that \(|{{\varvec{m}}}|>1\) and suppose that the claim holds for all multi-indices with order less than \(|{{\varvec{m}}}|\). Then
as desired. \(\square \)
Lemma A.3
Let \({\varvec{\nu }}\in {\mathscr {F}}\) and let \(({\mathbb {A}}_{\varvec{\nu }})_{{\varvec{\nu }}\in {\mathscr {F}}}\) and \(({{\mathbb {B}}}_{\varvec{\nu }})_{{\varvec{\nu }}\in {\mathscr {F}}}\) be arbitrary sequences of real numbers. Then
Proof
Moving the sum over \({{\varvec{m}}}\) further inside, we can write
as required. In the second equality we used the identity (cf., e.g., [28, formula 26.8.23])
In the third equality we introduced a change of variables \({\varvec{\mu }}'= {{\varvec{w}}}+{\varvec{\mu }}\). In the fourth equality we swapped the order of the sums. Finally we relabel \({\varvec{\mu }}'\) to \({{\varvec{m}}}\). \(\square \)
Lemma A.4
Let \(c>0\), \({\varvec{\beta }}=(\beta _j)_{j\ge 1}\) and \((\varUpsilon _{{\varvec{\nu }},{\varvec{\lambda }}})_{{\varvec{\nu }}\in {\mathscr {F}},{\varvec{\lambda }}\in \mathbb Z^d}\) be sequences of non-negative real numbers satisfying
Then
The result is sharp in the sense that all inequalities can be replaced by equalities.
Proof
We prove the claim by induction with respect to the multi-indices \({\varvec{\nu }}\in {\mathscr {F}}\). The case \({\varvec{\nu }}={{\varvec{0}}}\) is resolved immediately by observing that
as desired.
Fix \({\varvec{\nu }}\in {\mathscr {F}}\) and \(\varvec{\lambda }\in {\mathbb {N}}_0^d\), and suppose that the claim holds for all pairs of multi-indices, with first multi-index of order \(\le |{\varvec{\nu }}|\) and second multi-index \(< {\varvec{\lambda }}\). If \({\varvec{\lambda }}={{\varvec{0}}}\), then the claim is trivially true, so we may assume in the following that \({\varvec{\lambda }}\ne {{\varvec{0}}}\). By letting \(j\ge 1\) be arbitrary, we find that
We introduce \({{\varvec{m}}}':=(m_1,\ldots ,m_{j-1},m_{j+1},\ldots )\), \({\varvec{\nu }}':=(\nu _1,\ldots ,\nu _{j-1},\nu _{j+1},\ldots )\), and \(\varvec{\beta '}:=(\beta _1,\ldots ,\beta _{j-1},\beta _{j+1},\ldots )\). Then
Making use of the fact that (see [28, formula 26.8.25])
we obtain
where we used the change of variable \({\overline{m}}_j = m_j+1\) in conjunction with the property \(S(\nu +1,0)=0\) for all \(\nu \in {\mathbb {N}}_0\). Recognizing the above expression as
concludes the proof. We remark that if the inequality for the recursion on \(\varUpsilon _{{\varvec{\nu }},\lambda }\) is replaced by equality then all steps of the proof are equalities and so is the final formula. \(\square \)
Lemma A.5
Let \(C>0\) be a constant and suppose that \({\varvec{\beta }}=(\beta _j)_{j\ge 1}\) is a sequence of non-negative real numbers. Suppose that \(\varUpsilon _{{{\varvec{0}}}}\le c_0\) and
Then
where
The result is sharp in the sense that all inequalities can be replaced by equalities.
Proof
The proof is carried out by induction with respect to the multi-indices \({\varvec{\nu }}\in {\mathscr {F}}\). The base step \({\varvec{\nu }}={{\varvec{0}}}\) is trivially true, so let us assume that \({\varvec{\nu }}\in \mathscr {F}{\setminus }\{{{\varvec{0}}}\}\) is such that the claim holds for all multi-indices with order \(<|{\varvec{\nu }}|\). Then the recurrence together with the induction hypothesis yields that
with
To simplify our presentation we abbreviate temporarily \({\mathbb {A}}_{{\varvec{w}}}:= {{\mathbb {P}}}_{{{\varvec{w}}}}\,{{{\varvec{w}}}!}\,{\varvec{\beta }}^{{{\varvec{w}}}}\) and \({{\mathbb {B}}}_{\varvec{\mu }}:= {(|{\varvec{\mu }}|+1)!}\,{{\varvec{\mu }}!}\,{\varvec{\beta }}^{{\varvec{\mu }}}\). Then we can write
where we used Lemma A.3 in the third equality.
Reinserting the definitions of \({\mathbb {A}}_{{\varvec{w}}}\) and \({{\mathbb {B}}}_{{{\varvec{m}}}-{{\varvec{w}}}}\), the inner sum over \({{\varvec{w}}}\) in (A.6) becomes
Substituting (A.6) and (A.7) into (A.5), we obtain
The result now follows from the definition of the sequence \({{\mathbb {P}}}_{{\varvec{m}}}\). \(\square \)
We need to derive an upper bound for the sequences \(({{\mathbb {P}}}_{{{\varvec{m}}}})_{{{\varvec{m}}}\in {\mathscr {F}}}\). In anticipation of this, let us define the auxiliary sequence
This sequence is listed as A003480 in The On-Line Encyclopedia of Integer Sequences (OEIS). We observe that this sequence satisfies the following property.
Lemma A.6
The sequence \((\tau _k)_{k\ge 0}\) defined by (A.8) satisfies
Proof
Let \({{\varvec{m}}}\in {\mathscr {F}}\setminus \{{\textbf{0}}\}\). Then
as required. \(\square \)
Remark. Lemma A.6 may be viewed as an extension of the renewal property of sequence A003480; see corresponding entry in the OEIS for further details.
We are now ready to state an upper bound for the sequence \(({{\mathbb {P}}}_{{\varvec{m}}})_{{{\varvec{m}}}\in {\mathscr {F}}}\).
Lemma A.7
The sequence \(({{\mathbb {P}}}_{{{\varvec{m}}}})_{{{\varvec{m}}}\in {\mathscr {F}}}\) defined by (A.4) satisfies
where the sequence \((\tau _k)_{k\ge 0}\) is defined by (A.8).
Proof
The proof is carried out by induction with respect to the multi-indices \({{\varvec{m}}}\in {\mathscr {F}}\setminus \{{\textbf{0}}\}\). The base step \({{\varvec{m}}}={\varvec{e}}_j\), \(j\ge 1\), is resolved by observing that \(\tau _1 = 2\) and
Let \({{\varvec{m}}}\in {\mathscr {F}}{\setminus }\{{\textbf{0}}\}\) with \(|{{\varvec{m}}}|\ge 2\) and suppose that the assertion holds for all multi-indices \({{\varvec{w}}}\in \mathscr {F}{\setminus }\{{\textbf{0}}\}\) with \(|{{\varvec{w}}}|<|{{\varvec{m}}}|\). Then by separating out the \({{\varvec{w}}}={{\varvec{0}}}\) term in (A.4) and applying the induction hypothesis to \({{\mathbb {P}}}_{{\varvec{w}}}\) for \({{\varvec{w}}}\ne {{\varvec{0}}}\), we obtain
where we used Lemma A.6 in the first equality. This completes the proof. \(\square \)
It remains to obtain an estimate on the \((\tau _k)_{k\ge 0}\) sequence.
Lemma A.8
The sequence defined by (A.8) satisfies
Proof
The claim is clearly true with equality when \(k=0\). If \(k\ge 1\), then the sequence admits the following characterization (see the OEIS entry for this sequence and references therein):
from which the claim readily follows. \(\square \)
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Hakula, H., Harbrecht, H., Kaarnioja, V. et al. Uncertainty quantification for random domains using periodic random variables. Numer. Math. 156, 273–317 (2024). https://doi.org/10.1007/s00211-023-01392-6
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00211-023-01392-6