Uncertainty quantification for random domains using periodic random variables

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, Kuo, and Sloan (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.


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 (Ω, F, P) be a probability space, where Ω is the set of possible outcomes, the σ-algebra F ⊂ 2 Ω is the set of all events, and P is a probability measure.We are interested in uncertainty quantification for the Poisson problem −∆u(x, ω) = f (x), x ∈ D(ω), u(x, ω) = 0, x ∈ ∂D(ω), for P-almost every ω ∈ Ω, (1.1) where the domain D(ω) ⊂ R d , d ∈ {2, 3}, is assumed to be uncertain.In the domain mapping framework, one starts with a reference domain D ref ⊂ R d .For example, in a manufacturing application, D 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 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 ref needs only to be Lipschitz, while the random perturbations V (ω) have to be C 2smooth.This is in contrast to the perturbation approach mentioned above, where both the reference domain and the random perturbations need to be 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 V (ω) 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 ρ(z) = . The density ρ 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, for more discussion.
The paper is organized as follows.The problem is formulated in Section 2, after which Section 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 Section 5, and some conclusions are drawn in Section 6.The appendix A contains several technical results required for the parametric regularity analysis.We use multi-index notation throughout this paper.The set of all finitely supported multi-indices is denoted by j=1 of real numbers and any ν ∈ 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 y.
For any matrix M , let σ(M ) denote the set of all singular values of M and let ∥M ∥ 2 := max σ(M ) denote the matrix spectral norm.For a function v on D ref we define where we apply the vector 2-norm (Euclidean norm) or matrix 2-norm (spectral norm) depending on whether v(x) is a vector or a matrix.Similarly, we define where ∇v(x) is the gradient vector if v(x) is a scalar, and v ′ (x) is the Jacobian matrix if v(x) is a vector.We will also need the standard Sobolev norm for a scalar function v : We also define the norm where D ⊂ 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 ≥ m ≥ 0, except for S(0, 0) := 1.

Parameterization of domain uncertainty
with stochastic fluctuations ψ i : Denoting the Jacobian matrix of ψ i by ψ ′ i , the Jacobian matrix J : We assume that the family of admissible domains {D(y)} y∈U is parameterized by and define the hold-all domain by setting The convention of explicitly writing down the factor 1/ √ 6 in (2.1) is not a crucial part of the analysis, but it ensures that the mean and covariance of the vector field V match those of an affine and uniform parameterization for representing domain uncertainty using the same sequence of stochastic fluctuations (ψ i ) ∞ i=1 , 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: (A2) For some C > 0, there holds where σ(J(x, y)) denotes the set of all singular values of the Jacobian matrix.
Remark.Under assumption (A3), there holds This follows from the continuity of the determinant and det J(x, 0) = 1.

The variational formulation on the reference domain
The variational formulation of the model problem (1.1) can be stated as follows: for y ∈ U , find u(•, y) ∈ H 1 0 (D(y)) such that where f ∈ C ∞ (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 y ∈ U , find for all v ∈ H 1 0 (D 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 y → u(•, y) of (2.7) with respect to the parametric variable, we need to derive bounds on the partial derivatives ∂ ν y u(•, y) in the Sobolev norm ∥ • ∥ H 1 0 (D ref ) .The analysis presented in this section is based on [15] and proceeds as follows: • We derive bounds for ∂ ν y B −1 (x, y) and ∂ ν y det J(x, y) for all ν ∈ F in Lemmata 3.1 and 3.2.These results are then used to derive a bound on the partial derivatives of the coefficient ∂ ν y A(x, y) in Lemma 3.3.
• Since the right-hand side of (2.7) depends on the parametric variable y ∈ U , we develop the derivative bounds on ∂ ν y f ref (x, y) in Lemma 3.5, with the aid of the derivative bound developed for ∂ ν y f (x, y) in Lemma 3.4.
• An a priori bound is developed for the the solution u(x, 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 ∂ ν y u(x, y), ν ∈ F , making use of the aforementioned results.

Parametric regularity of the transport coefficient
with the sequence Proof.We will first derive a regularity bound for follows by applying implicit differentiation to the identity From the definition of J it is easy to see that for any multi-index m ∈ F and any x ∈ D ref and y ∈ U we have Taking the matrix 2-norm, we obtain and hence with (2.2) we obtain Now the Leibniz product rule yields for ν ∈ F , and after taking the matrix 2-norm on both sides and then the L ∞ -norm over x we obtain The bounds in (3.3) indicate that only the derivatives with |supp(ν)| ≤ 2 will survive.Indeed, for ν = k e j with k ≥ 1, only m = ℓ e j for ℓ = 0, . . ., k remain, while for ν = k e j + k ′ e j ′ with k, k ′ ≥ 1 and j ̸ = j ′ , only m = k e j and m = k ′ e j ′ remain.On the other hand, for ν with |supp(ν)| ≥ 3, in each term at least one of the factors must vanish.We conclude that where we used b j ≤ ξ b for all j ∈ N.
Our assumption (A3) immediately yields for all y ∈ U that which proves the lemma for the case ν = 0. Now we consider ν ∈ F \ {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 m = 0 term, post-multiplying both sides by B −1 (x, y), taking the matrix 2-norm followed by the L ∞ -norm over x, and applying (3.4), we obtain The assertion follows by applying Lemmata A.1 and A.
, and β j = b j .Lemma 3.2.Under assumptions (A1)-(A4), there holds for all y ∈ U and ν ∈ F that Proof.Let y ∈ U .The proof is carried out by induction over all the minors (subdeterminants) of the matrix J(x, y) =: denote the q × q submatrix specified by the indices ℓ := {ℓ 1 , . . ., ℓ q } and ℓ We will prove by induction on q that Similarly to (3.2), we have for the derivatives of matrix elements 2), we obtain for the (1 × 1)minors Suppose (3.5) holds for all submatrices of size up to (q − 1) × (q − 1), and now we consider the case for a q×q submatrix.For arbitrary t ′ ≤ 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.
Proof.Let ν ∈ F .By the Leibniz product rule, we have from (2.6) that Taking the matrix 2-norm followed by the L ∞ -norm over x, and then applying Lem-mata 3.1 and 3.2, we obtain where A w := |w|! a |w| w! b w and B µ := |µ|! b µ , and we overestimated some multiplying factors to get The last equality is due to Lemma A.3.Using a |w| ≤ a |m| and w! ≤ m!, we have where the last sum over w can be rewritten as thus completing the proof.

Parametric regularity of the transported source term
Lemma 3.4.Let f ∈ C ∞ (D) be analytic.Then there exist constants Furthermore, under the assumptions (A1)-(A4), there holds for all y ∈ U and ν ∈ F that the derivatives of f (x, y) = f (V (x, y)) are bounded by is a consequence of the Cauchy integral formula for analytic functions of several variables (cf., e.g., [20 so we may assume in the following that ν ∈ F \{0}.We will make use of the multivariate Faà di Bruno formula [29, Theorem 3.1 and Remark 3.3] where we define κ ν,λ (x, y) for general ν ∈ F and λ ∈ Z d in a recursive manner as follows: κ ν,0 (x, y) := δ ν,0 , κ ν,λ (x, y) := 0 if |ν| < |λ| or λ ̸ ≥ 0 (i.e., λ contains negative entries), and otherwise From the definition of V in (2.1) it is easy to see that Thus we obtain from (3.9) the recursion By Lemma A.4 with c = 2π, we have for all ν ∈ F and λ ∈ N d 0 that This yields the desired result since . The m ̸ = 0 condition for the above sum can be dropped since the m = 0 term is actually zero when ν ̸ = 0.This allows us to then extend the formula to cover also the case ν = 0. Lemma 3.5.Under the assumptions of Lemma 3.4, there holds for all y ∈ U and ν ∈ 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 ∞ -norm over x and using Lemmata 3.2 and 3.4 yields where now where the last sum over w can be rewritten as This completes the proof.

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 y ∈ U that where C P > 0 is the Poincaré constant associated with the domain D ref .
Proof.We take v(x) = u(x, y) as test function in (2.7) to obtain , where we used the Cauchy-Schwarz inequality and the Poincaré inequality , with the constant C P > 0 depending only on the domain D 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 σ max = 1 + ξ b , see line below (2.2).The proof is completed by combining the upper and lower bounds, and applying the bound for ∥f ref (•, y)∥ L ∞ (D ref ) by taking ν = 0 in Lemma 3.5.
Theorem 3.7.Under the assumptions of Lemma 3.4, there holds for all y ∈ U and ν ∈ F \ {0} that where Proof.We prove the theorem by induction on the order of ν.The case ν = 0 holds by Lemma 3.6.Let y ∈ U and ν ∈ F \ {0}, and suppose that the result holds for all multi-indices of order less than |ν|.
We differentiate both sides of (2.7) and apply the Leibniz product rule to obtain for Taking v(x) = ∂ ν y u(x, y) and separating out the m = ν 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 ≤ (1 + √ 3) k (which follows easily from (3.1)), In particular, we used |ν − m| ≥ 1 for the exponent in ( 4π(1+ξ b ) σ 2 min ) |ν−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 ν ̸ = 0 Redefining β j to incorporate the factor 2 + √ 2 completes the proof.

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.

Dimension truncation error
We are interested in finding a rate for , where u denotes the solution to (2.7) as before and the dimensionally-truncated PDE solution is defined by u s (•, y) := u(•, (y 1 , . . ., y s , 0, 0, . ..)) for y ∈ U, s ∈ N.
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 dimensionallytruncated solution to (2.7) satisfies the error bounds where the implied coefficients are independent of the dimension s.
Proof.We start by defining where the fluctuations (∥ψ i ∥ W 1,∞ ) i≥1 ∈ ℓ p for the same p as in (A5).Let us consider an auxiliary PDE problem where the domain realizations D(y) := V (D ref , y) are defined instead via an affine parametric random perturbation field with (∥ ψ i ∥ W 1,∞ ) i≥1 ∈ ℓ p for the same p as in (A5).The random field 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 b ∈ ℓ p , where the p is the same as above.
Let G ∈ H −1 (D ref ) be arbitrary and define Then developing the Taylor series of F about the point (y ≤s , 0) := (y 1 , . . ., y s , 0, 0, . ..) yields Let k ≥ 2 be an undetermined integer which will be specified later.By defining θ(y) := (sin(2πy i )) i≥1 , we observe that where u per (•, y) is the smooth periodic extension of the solution to (2.7) for all y ∈ . We obtain by a simple change of variable that Since the last equality holds pointwise for all y ∈ [−1, 1] N , we can integrate the equation on both sides over U ⊂ [−1, 1] N .We note that the integral 1 0 sin(2πy i ) ν i dy i takes the value between 0 and 1, and gives the value 0 when ν i = 1.Thus the summands corresponding to ν ∈ F with ν i = 1 for at least one i > s vanish in the first term.Hence , 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 and since the sequence is b is non-increasing we defined τ k := where we used the well-known Stechkin's lemma (cf., e.g., [25,Lemma 3.3]): for 0 < p ≤ q < ∞ and any sequence (a i ) i≥1 of real numbers ordered such that Therefore The term T 2 in (4.2) can be estimated similarly to the approach taken in [9, Theorem 4.1].The trivial inequality |ν|! ν! ≥ 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 = ⌈1/(1 − p)⌉.
Finally, the Poincaré inequality implies that and the L 1 error can be obtained by using the Cauchy-Schwarz inequality .
This completes the proof.

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 1 0 (D 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 ref by recursive uniform bisection of simplices.
We define the dimensionally-truncated finite element solution for y ∈ U as the solution for all v h ∈ V h , where we set A s (x, y) := A(x, (y 1 , . . ., y s , 0, 0, . ..)) as well as f ref,s (x, y) := f ref (x, (y 1 , . . ., y s , 0, 0, . ..)) for all y ∈ U .Since the reference domain D ref was assumed in (A7) to be a bounded, convex polyhedron, it follows that elements in W 1,∞ (D ref ) have Lipschitz continuous representatives.Hence we may apply Lemma 3.5 and [10, Theorem 3.2.1.2]to obtain where and the constant C > 0 only depends on d and the diameter of D ref .
We have by Céa's lemma that for a constant C ′ > 0 independent of s, h, and y, which can be seen by letting v h ∈ V h be arbitrary and deriving where the penultimate step follows from Galerkin orthogonality.The well-known approximation property (cf., e.g., [8,26]) 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 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 dimensionallytruncated finite element solution that where the generic constant is independent of s, h, and y.
Testing against v = u s (•, y) − u s,h (•, y) and letting v h ∈ V h be arbitrary, it follows from Galerkin orthogonality of the finite element solution that and, in consequence, Then there holds for all y ∈ U that By (4.6) and (4.4), we obtain inf Plugging this inequality and (4.7) into (4.8)yields the desired result.

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 ∈ N.An n-point QMC rule is an equal weight cubature rule of the form where the cubature nodes y (1) , . . ., y (n) ∈ [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 z ∈ {1, . . ., 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] 1 where the dual lattice Λ ⊥ := {ℓ ∈ Z s : ℓ • z ≡ 0 (mod n)} depends on the lattice generating vector z.For integrands belonging to the weighted Korobov space endowed with the norm where α > 1 is the smoothness parameter, γ := (γ u ) u⊆{1:s} are positive weights, and supp(ℓ) := {1 ≤ j ≤ s : ℓ j ≥ 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 (y) = u s,h (x, y) for x ∈ D 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 α := ⌊1/p⌋ + 1 with p from Assumption (A5), and choose the weights γ ∅ := 1 and where 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 α and the weights γ u as arbitrary and will specify their values later in the proof.We will assume that α ≥ 2 is an integer, and we will express the integration error in term of the mixed derivatives of F (y) = u s,h (x, y) of order up to α in each coordinate.For u ⊆ {1 : s} we denote by ν = (α u ; 0) the multi-index whose component ν j is equal to α if j ∈ u and is 0 otherwise.First we observe by differentiating the Fourier series of F that We are ready to derive the integration error for u s,h .Using (4.9) we have error =

.14)
Denoting u = supp(ℓ), we obtain using the definition of r α,γ (ℓ) 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 , and C 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.
By choosing λ = p and α = ⌊1/p⌋ + 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 O(n −1/p ), with the implied constant independent of the dimension s when the weights are chosen according to (4.12).

Overall error estimate
Let u be the solution to (2.7), u s the corresponding dimensionally-truncated solution, and u s,h the dimensionally-truncated FE solution.It is easy to see that We combine the error bounds derived in Subsections 4.1-4.3.The error decomposition (4.17) allows us to deduce the following overall error estimate.Theorem 4.4.Let (y k ) n k=1 be the lattice cubature nodes in U s generated by a CBC algorithm using the construction detailed in Subsection 4.3.Suppose that for each lattice point, we solve (2.7) using one common finite element discretization in the domain D 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.

Numerical experiments
Let us consider the domain parameterization sin(2πy j )ψ j (x), x ∈ [0, 1] and y ∈ U s ( with c > 0 and ψ j (x) := j −θ cos(jπx), θ > 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(y) = V (D(0), y) with the vector-valued expansion where .
We note that cos 2 (jπx 1 ) + j 2 π 2 x 2 2 sin 2 (jπx 1 ) = jπ. (5. 3) 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πx 1 )+j 2 π 2 x 2 2 sin 2 (jπx 1 ) = 1+sin 2 (jπx 1 )(j 2 π 2 x 2 2 − 1) is maximized by x 1 = 1 2j and x 2 = 1.The elements of the sequence (b j ) ∞ j=1 needed for the construction of the QMC weights now simplify to Moreover, we can take The expected rate of QMC convergence is O(n −θ+1 ).relative error vs. n with dimension s We consider the Poisson problem (1.1) equipped with f (x) := x 2 and the random domain described above for θ ∈ {2.1, 2.5, 3.0} and c = 3/2.The input random field is truncated to s = 100 terms.We construct a uniform regular triangulation for the reference domain D ref = D(0) with mesh size h = 2 −5 , which is then mapped to each realization of the random domain using V .The PDE is solved using a first order finite element method.

Source problem
For each realization of y ∈ U s , the finite element solution is computed in D(y) and then mapped back onto the reference domain D 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 σ min = ρ = 1 and c = 10 −6 for numerical stability.The results are displayed in the left-hand side image of Figure 1.In this case, the observed rates σ θ are for θ = 2.1: σ 2.1 = 0.95, θ = 2.5: σ 2.5 = 1.15, θ = 3: σ 3 = 1.72.The theoretically predicted convergence rates are clearly observed in the numerical results: using a higher value for the decay rate θ 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 Figure 1.The observed rates σ θ are for θ = 2.1: σ 2.1 = 0.89, θ = 2.5: σ 2.5 = 1.09, θ = 3: σ 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.

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 (5.7) As the QoI, we consider the capacity cap(D(y)) := where u is the solution to (5.6).We also define its conjugate as cap(D(y)) := where v is the solution to (5.7).It is well-known that (cf., e.g., [1,) cap(D(y)) cap(D(y)) = 1 for all y ∈ U s , ( provided that x → a(x, y) specifies a Jordan arc such that a(x, y) > 0 for all x ∈ [0, 1] and y ∈ 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.where cap(D(y)) and cap(D(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(y).One realization of D(y) with respective solutions u(x, y) and v(x, y) is shown in Figure 2.
The numerical experiments were carried out as follows: (i) We set s = 10, c = 3/2, and let θ ∈ {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 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 Figure 3a.
(ii) We set c = 3/2, θ = 2.1, and let s ∈ {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 Figure 3b.
The QMC point sets for the above experiments were computed using the fast CBC algorithm using the same specifications as in Subsection 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 σ θ are for θ = 2.1: σ 2.1 = 0.7, θ = 2.5: σ 2.5 = 1.3, θ = 3: σ 3 = 2.4 (see Figure 3).The rate does not appear to depend on the dimension s.The cone of convergence for θ = 2.5 shown in Figure 3c, is bounded by the rates σ 2.5A = 2.3 (above) and σ 2.5B = 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 Section 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.

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.
Let m ∈ F .It is interesting to note that if |supp(m)| ≤ 2, for example, m u = ke i +ℓe j with k, ℓ ∈ 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 (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 can be represented using the recurrence a see, e.g., entry A080599 in the On-line Encyclopedia of Integer Sequences (OEIS ) and references therein.We also have We prove the claim (A.2) by induction with respect to the order of the multi-indices m ∈ F while utilizing the recursive characterization of sequence (a ′ k ) ∞ k=0 .We immediately see that (A.2) holds with equality when |m| ≤ 1. Fix m ∈ F such that |m| > 1 and suppose that the claim holds for all multi-indices with order less than |m|.Then In the third equality we introduced a change of variables µ ′ = w + µ.In the fourth equality we swapped the order of the sums.Finally we relabel µ ′ to m.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 ν ∈ F .The case ν = 0 is resolved immediately by observing that Υ 0,0 = 1 and Υ 0,λ = 0 if λ ̸ = 0 as desired.Fix ν ∈ F and λ ∈ N d 0 , and suppose that the claim holds for all pairs of multi-indices, with first multi-index of order ≤ |ν| and second multi-index < λ.If λ = 0, then the claim is trivially true, so we may assume in the following that λ ̸ = 0.By letting j ≥ 1 be arbitrary, we find that Υ ν+e j ,λ where we used the change of variable m j = m j + 1 in conjunction with the property S(ν + 1, 0) = 0 for all ν ∈ N 0 .Recognizing the above expression as β m S(ν j + 1, m j ) i̸ =j S(ν i , m i ) concludes the proof.We remark that if the inequality for the recursion on Υ ν,λ is replaced by equality then all steps of the proof are equalities and so is the final formula.
Lemma A.5.Let C > 0 be a constant and suppose that β = (β j ) j≥1 is a sequence of non-negative real numbers.Suppose that Υ 0 ≤ c 0 and where we used Lemma A.6 in the first equality.This completes the proof.It remains to obtain an estimate on the (τ k ) k≥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 ≥ 1, then the sequence admits the following characterization (see the OEIS entry for this sequence and references therein): from which the claim readily follows.
Let U := [0, 1] N denote a set of parameters.Let us fix a bounded domain D ref ⊂ R d with Lipschitz boundary and d ∈ {2, 3} as the reference domain.

∀i≤s b ν ,
where C ′ = C k! and C ′′ = C (k + 1)!.Recalling the definition of F and taking the supremum over all G ∈ H −1 (D ref ) with ∥G∥ H −1 (D ref ) ≤ 1 yields for the left-hand side that sup
7) are not computed on the reference domain D ref ; for each realization of y ∈ U s , the finite element solution is computed in D(y).Numerically, we can use the reciprocal relation (5.8) to define an a posteriori error estimate err(y) := |1 − cap(D(y)) cap(D(y))|, (5.9)