Low-Dimensional Spatial Embedding Method for Shape Uncertainty Quantification in Acoustic Scattering

1 Introduction

Considerable effort has been devoted in recent years to develop robust and efficient computational strategies for the simulation of physical phenomena, that take into account shape uncertainty. Often, the problem is formulated as an elliptic partial differential equation (PDE) whose domain boundaries are uncertain. Such problems arise due to imperfections in manufacturing processes, e.g., in nano-optics where the production of nano particles is, often, inaccurate relatively to nano-scale electromagnetic wave lengths bejan2000shape . Other examples arise in the context of inverse problems, such as tomography where the visual representation of some hidden object is constructed by partial, and possibly noisy, measurements sethian1999level .

The common practice for quantifying uncertainty in computational models, is to employ the parametric method. The theoretical basis for the method was laid down by Wiener gPC-Wiener-HomoChaos . The method itself was initially developed by Ghanem and Spanos gPC-GS-StochasticFE , and later generalized by Xiu and Karniadakis gPC-WXK-FlowSimulations ; gPC-XK-SteadyStateDiffusion ; gPC-XK-Wiener-Askey

. In this approach the uncertain parameters are replaced by random quantities, and the problem is recast as a system with random input. The solution is estimated via global expansion of the random variables into a basis of uncorrelated functions. Thus, the stochastic problem is transformed into a deterministic system in higher dimension. The most popular expansions employed are the

Karhunen-Loève expansion, and the generalized Polynomial Chaos (gPC) expansion.

The parametric approach, often, demonstrates superior performance in terms of computational effort over other traditional methods, see gPC-Xiu-Book for a detailed review. However, when the physical domain of the problem is uncertain the quantification by the parametric approach becomes much more challenging. The main difficulty stems from the fact that the problem is not characterized by smooth coefficients whose dependence on the random parameters is known. Thus, an accurate discretization which captures desired features of the solution for any realization of the random shape is not readily available.

The stochastic collocation method and Monte Carlo sampling, which rely on samplings of the random parameters and the solution of each realization deterministically are well established. However, for random domain problems each realization is, essentially, characterized by a different geometry and requires a custom discretization scheme. Generally, when the variations of the random domain are large and undergo complicated changes as a function of the random parameters, these methods become extremely expensive to employ with prohibitive computational costs.

To overcome the difficulties associated with the quantification of a random shape or domain, various techniques have been proposed. Typically, these are classified as one of the following: perturbation, fictitious domain, level-set or random domain mapping. Perturbation techniques

Shape-Perturbation are straightforward and simple to apply, however, their applicability is limited to small shape deformations. The fictitious domain Shape-Fictitious and level-set methods Shape-LevelSet-Nouy20084663 ; Shape-LevelSet-XFEM are based on embedding the random domain in a larger, deterministic domain containing all possible realizations. These methods are capable of handling very irregular non-smooth geometries. However, the embedding introduces non-smoothness in the spatial region, that intersects with the random boundary. Thus, high-order convergence is only partially ensured in the entire computational domain.

The random domain mapping method Shape-TX2006b ; Shape-TX2006a is the most common tool used for solving PDEs on uncertain domains. The method is based on a realization-dependent coordinate transformation uniformly mapping all the realizations of the domain to a fixed, reference configuration. The variational formulation of the PDE on the random domain can then be posed on the reference domain, reducing the problem to a PDE on a fixed domain with stochastic coefficients. The transformed PDE whose domain is fixed is solved using standard techniques. However, the method is highly sensitive to the non-linear dependence of the problem on the random boundary. In case of complex evolution of the shape, the random coefficients are difficult to obtain and typically exhibit highly varying behavior. The common practice to overcome this difficulty is to impose a highly accurate discretization grid, often combined with dimensionality reduction techniques, e.g., sparse grids, to ensure reasonable computational effort. See Shape-DomainMapping-Regularity ; Shape-Harbrecht2016 ; Shape-Hptmr2015 for further details.

In this work an alternative method that attempts to mitigate the difficulties associated with the more standard techniques for PDEs on uncertain domains is proposed. The method is of boundary integral type BIE-CK-InverseScattering-Book and, thus, can handle large shape deformations. The analysis is based on two observations. First, that as a function of the random boundary of the domain the solution is piece-wise smooth in the spatial domain. Second, that in practice we seek to approximate the outcome of a predetermined set of linear output functionals operating on the random boundary. The key idea is of this work is to construct an integration grid whose associated weight function encompasses the irregularities and non-smoothness imposed by the random boundary. Thus, the outcome of the functionals can be evaluated accurately with relatively low number of integration gridpoints. This idea is similar to certain classic numerical techniques for estimating integrals of highly oscillatory functions, which rely on oscillatory weighted Gaussian integration formulae.

The proposed method constructs a discretization grid of the random surface for all possible realizations in two stages. In the first stage a spatial low-dimensional embedding of the family of random surfaces is constructed via the Coarea formula Federer . The embedding, essentially, captures any irregular behavior of the random surface and a discretization is applied only on a compact region in the spatial domain. In the second stage a parametric grid corresponding to the low-dimensional spatial grid is imposed. A sparse or hierarchical parametric grid can be applied for dealing with high dimensionality, while the spatial grid effectively ensures that the bulk variation of output functionals defined on the boundary is captured. In general, the method allows the handling of non-trivial geometries without the loss of accuracy in the region intersecting with the random interface.

Since this is a first case study and for the ease of presentation, the discussion has been limited to time-harmonic wave scattering by star-shaped obstacles. In the analysis and numerical study a $2 D$ scattering object and low-dimensional parametric space are assumed. More complicated examples in higher spatial and parametric dimensions are discussed. However, in-depth study of this topic is deferred to future work. For its simplicity, acoustic fluid-structure interaction has been chosen as the physical application. In that case, the solution represents small oscillations of pressure in a compressible ideal fluid. The method and ideas presented in this work can also be applied to electrodynamics and elastodynamics.

This work employs the null-field approach NFM-Martin-Book ; NFM-Waterman-MatrixFormulation ; NFM-Waterman-T-Matrix ; NFM-Wriedt-Review , which in contrary to the better known boundary element method (BEM) BIE-sauter2010boundary and the Nyström method kress2013linear , does not involve singular integrals. Null-field methods are fast and much easier to implement compared to BEM and the Nyström method. Their applicability range is, however, more limited. The null-field reconstruction technique NFM-EM-Random-Shape-1D ; NFM-Optimal-Reconstruction is inherently stable, admits a-priori error evaluation, and facilitates the extraction of features of interest without prior estimation of the entire solution. The method enables us to perform analysis from a purely geometric point of view, which avoids the additional complications associated with integration of weakly singular kernels. Combining low-dimensional surface embedding with BEM and Nyström method can be foreseen in a future study.

The paper is organized as follows. The fundamentals of the null-field reconstruction method for the time-harmonic wave scattering problem is presented in Section 2. Section 3 reviews the procedure of optimal reconstruction from a numerical linear algebra point of view. Section 4 consists of the main theoretical results of this work and includes the formulation of the problem. In Section 5 the proposed method is applied to a class of randomly shaped polygonal cylinders, as a proof of concept that the suggested method can, indeed, handle complex non-smooth shapes. Summary of the results, conclusions, and suggestions for applying the method in more complicated scenarios are given in Section 6.

2 Null-Field Reconstruction for Time Harmonic Wave Scattering

In this section a brief review on the null-field reconstruction method for the time-harmonic wave scattering problem is given. The time-harmonic acoustic scattering problem is presented, followed by a review of the fundamental theory of null-field methods. The main idea of the null-field reconstruction technique for time-harmonic wave scattering is presented in the concluding subsection.

2.1 Acoustic Scattering by Impenetrable Obstacles

Let $B$ denote a bounded domain in $R^{2}$ representing an impenetrable obstacle with boundary $S$ . We denote by $¯ ¯¯ ¯ B = B \cup S$ the closure of $B$ . Let $R^{2} ∖ ¯ ¯¯ ¯ B$ be the unbounded exterior region occupied by a uniform medium. Let $r \in R^{2}$ denote a general spatial point,

r = (r cos θ, r sin θ),

For an incident time-harmonic field $u^{i n c} (r)$ ’illuminating’ the obstacle, the scattered field $u^{s c a} (r)$ satisfies the following exterior boundary value problem:

Δ u^{s c a} + κ^{2} u^{s c a} = 0 \forall r \in R^{2} ∖ ¯ ¯¯ ¯ B,

(1)

u^{s c a} = - u^{i n c} or \partial_{ν} u^{s c a} = - \partial_{ν} u^{i n c} on S,

(2)

(3)

where $i = \sqrt{- 1}$ . Equation (1) is known as the Helmholtz equation, where $Δ$ is the Laplacian and $κ$ is the wavenumber. Equation (2) specifies the boundary condition, depending on the physical problem: adopting the acoustic terminology, it is sound-soft for Dirichlet problems and sound-hard for Neumann problems. Here $ν$ is the unit outward normal to $S$ and $\partial_{ν} u$ is the normal derivative of $u$ . The last condition (3), known as the Sommerfeld radiation condition, ensures that the scattered field propagates from the obstacle to infinity. The solution of the exterior scattering problem is unique. A solution to the Helmholtz equation is called a wavefunction. A wavefunction satisfying the Sommerfeld condition (3) is called an outgoing wavefunction.

2.2 Null-Field Theory Fundamentals

Null-field methods for the acoustic scattering problem (1) are based on Green’s second theorem

\iint_{B} (w Δ v - v Δ w) d V = \int_{S} (w \frac{\partial v}{\partial ν} - v \frac{\partial w}{\partial ν}) d S,

(4)

which holds for any bounded domain $R$ with a Lipschitz piecewise smooth boundary $S$ , where $v$ and $w$ are scalar fields, and $\partial v / \partial ν$ and $\partial w / \partial ν$ denote corresponding normal derivatives.

Let $ψ$ be an outgoing wavefunction and let $u$ denote the total field, $u^{i n c} + u^{s c a}$ . Assuming $ψ$ is analytic in $R^{2} ∖ ¯ ¯¯ ¯ B$ , it can be shown by (4) that

Thus, for a sound-soft obstacle ( $u = 0$ on $S$ )

\int_{S} ψ \frac{\partial u}{\partial ν} d S = - \int_{S} (u^{i n c} \frac{\partial ψ}{\partial ν} - ψ \frac{\partial u^{i n c}}{\partial ν}) d S,

(5)

while for a sound-hard obstacle ( $\partial_{ν} u = 0$ on $S$ )

\int_{S} u \frac{\partial ψ}{\partial ν} d S = \int_{S} (u^{i n c} \frac{\partial ψ}{\partial ν} - ψ \frac{\partial u^{i n c}}{\partial ν}) d S .

(6)

Using (5) or (6), an infinite set of equations can be produced from which $u$ or $\partial_{ν} u$ on $S$ are approximated. In practice, one chooses a finite subset of equations of the form of (5) and (6) which are employed to optimally reconstruct the scattered field without an explicit estimation of $u$ or $\partial_{ν} u$ on $S$ . The core idea of reconstruction by functionals is presented in the next subsection, while the numerical procedure for its practical implementation is covered in Section 3.

2.3 Reconstruction of Surface Functionals

Typically we are interested in estimating features of interest which are expressed by the unknown surface density, $u$ or $\partial_{ν} u$ on $S$ . Often such features are the outcomes of functionals in an appropriate Hilbert space. Indeed, let $h$ denote the complex conjugate of $u$ or $\partial_{ν} u$ on $S$ . Then for a general non-smooth surface $S$ , the surface density $h$ belongs to the complex Hilbert space $L_{d S}^{2} (S)$ whose inner-product and norm are defined by

{⟨ g, h ⟩}_{L_{d S}^{2}} = \int_{S} h^{*} (r) g (r) d S, {∥ h ∥}_{L_{d S}^{2}} = \sqrt{{⟨ h, h ⟩}_{L_{d S}^{2}}},

where $h^{*}$ denotes the complex conjugate of $h$ and $d S$ is the induced volume form on the surface. Recall that by Riesz representation theorem any bounded linear functional $L_{d S}^{2} \to C$ operating on surface densities, is of the form ${⟨ f, h ⟩}_{L_{d S}^{2}}$ . We call such functionals surface functionals.

In this work we focus on the estimation of the scattering coefficients of the expansion of $u^{s c a}$ to cylinder harmonics in the $2 D$ case. These coefficients, denoted by $b_{m}$ , satisfy

u^{s c a} (r) = \infty \sum m = - \infty b_{m} H_{m} (κ r) e^{i m θ},

where $H_{m} (z)$ denotes the $m$ th-order Hankel function of the first kind. See BIE-CK-InverseScattering-Book for further details. The scattering coefficients are very useful features of the surface density, since they can easily express other important quantities such as the far-field pattern and the radar cross section NFM-Thesis .

Consider the sound-soft case (5). Using the Hilbert space notation, it follows that the scattering coefficients satisfy

\frac{i}{4} {⟨ J_{m} (κ r) e^{- i m θ}, \partial_{ν} u^{*} (r) ⟩}_{L_{% d S}^{2}}^{- i m θ} = b_{m}, \forall m \in Z,

(7)

where $J_{m} (z)$ denotes the $m$ th-order Bessel function of the first kind, and $\partial_{ν} u^{*} (r)$ denotes the complex conjugate of the surface density $\partial_{ν} u (r)$ . In practice only ${b_{m}}_{| m | \leq μ}$ satisfying

μ = ⌈ 3 κ r^{max} ⌉, r^{max} = inf {| r | ∣ ∣ r \in R^{2} ∖ ¯ ¯¯ ¯ B},

(8)

are required for an accurate description of $u^{s c a}$ , see NFM-Optimal-Reconstruction for more details. Similar expressions can be derived for the sound-hard case.

The core idea of the reconstruction procedure is to approximate the outcome of target functionals (7) and without producing an explicit approximation of the surface density $\partial_{ν} u$ on $S$ . This, generally, allows us to handle complex geometries as well as irregular or singular surface densities much more accurately . Explicitly, we approximate the each outcome (7) by a linear combination of the following form

\frac{i}{4} {⟨ J_{m} (κ r) e^{- i m θ}, ¯ ¯¯¯¯¯¯ ¯ \partial_{ν} u ⟩}_{L_{d% S}^{2}}^{- i m θ} \approx {⟨ \sum c_{m, l} ψ_{ℓ}, ¯ ¯¯¯¯¯¯ ¯ \partial_{ν} u ⟩}_{L_{d S}^{2}},

where ${ψ_{ℓ}}$ are predetermined sets of functionals, whose outputs are either known or can be calculated directly. We call such functionals the information functionals.

As shown in (5) the outcome of the information functionals are readily available if $ψ_{ℓ}$ are outgoing wavefunctions whose singularities are located in $¯ ¯¯ ¯ B$ . Hence, given the information

{⟨ ψ_{ℓ}, ¯ ¯¯¯¯¯¯ ¯ \partial_{ν} u ⟩}_{L_{d S}^{2}} = a_{ℓ},

the outcome of the target functional can be approximate by

\frac{i}{4} {⟨ J_{m} (κ r) e^{- i m θ}, ¯ ¯¯¯¯¯¯ ¯ \partial_{ν} u ⟩}_{L_{d% S}^{2}}^{- i m θ} \approx \sum c_{m, l} a_{ℓ}, \forall m \in Z .

Obtaining the coefficients ${c_{m, ℓ}}$ while ensuring measurable error bounds of the estimations of the scattering coefficients can be achieved by reconstruction kernel approximation which is the main topic of Section (3).

3 Optimal Reconstruction with A-priori Error Estimate

In this section we review the procedure for the recovery of target functionals by information functionals in general Hilbert space. Reconstruction problems often involve regularization parameters which govern the stability and accuracy of the procedure. The result of the optimization of the reconstruction with respect to the regularization parameters is referred to as optimal reconstruction. Optimal reconstruction can be traced back to the notion of optimal recovery OR-GW-OptimalApproximation ; OR-MR-OptimalEstimation . A more modern analysis from an inverse problem point of view can be found in Moment-Louis-FeatureReconstruction ; Moment-Louis-UnifiedApproach ; Moment-LM-MollifierMethod .

We begin with the definition of the reconstruction problem and the notion of reconstruction kernel. This is followed by a brief description of the numerical procedure including error analysis. The final part elaborates on proper numerical integration rules, that are needed for the error estimates. The method and error analysis presented here, as well as further technical details have been initially introduced in NFM-Optimal-Reconstruction .

3.1 The Reconstruction Problem and Reconstruction Kernels

Let $H$ be a complex Hilbert space, whose inner-product is denoted by ${⟨, ⟩}_{H}$ . The reconstruction problem is to approximate a finite set of target functionals

{⟨ g_{m}, h ⟩}_{H} = b_{m}, g_{m} \in H, m = 1, 2, \dots, M,

(9)

by a given finite set of information functionals,

{⟨ f_{ℓ}, h ⟩}_{H} = a_{ℓ}, f_{ℓ} \in H, l = 1, 2, \dots, L,

(10)

where the element $h \in H$ is unknown.

Definition 1

Let ${∥ ∥}_{H}$ denote the norm induced by ${⟨, ⟩}_{H}$ in $H$ , and let $C_{m}^{L}$ denote a closed convex subset of $C^{L}$ . A linear combination $\sum_{ℓ = 1}^{L} {ˆ c}_{m, ℓ} f_{ℓ}$ whose coefficients $({ˆ c}_{m, 1}, \dots, {ˆ c}_{m, L}) \subset C_{m}^{L}$ satisfy the minimality condition

{∥ ∥ ∥ ∥ g_{m} - L \sum ℓ = 1 {ˆ c}_{m, ℓ} f_{ℓ} ∥ ∥ ∥ ∥}_{H} \leq {∥ ∥ ∥ ∥ g_{m} - L \sum ℓ = 1 c_{m, ℓ} f_{ℓ} ∥ ∥ ∥ ∥}_{H} \forall (c_{m, 1}, \dots, c_{m, L}) \in C_{m}^{L},

(11)

is called an optimal reconstruction kernel of the target functional $g_{m}$ by the information functionals ${f_{ℓ}}$ over $C_{m}^{L}$ .

Remark 1

Clearly, (11) is a projection on a convex set. The key point which is addressed later, is how to determine the convex set $C_{m}^{L}$ . Note that almost no prior knowledge on the element $h$ is assumed.

We will show in the next subsection, that obtaining (11) vastly exceeds our needs. In practice, it is sufficient to obtain an approximation satisfying

{∥ ∥ ∥ ∥ g_{m} - L \sum ℓ = 1 {ˆ c}_{m, ℓ} f_{ℓ} ∥ ∥ ∥ ∥}_{H} \leq ϵ {∥ g_{m} ∥}_{H} ({ˆ c}_{m, 1}, \dots, {ˆ c}_{m, L}) \in C_{m}^{L},

(12)

with respect to some predetermined threshold, $ϵ > 0$ . In that case the linear combination $\sum_{ℓ = 1}^{L} {ˆ c}_{m, ℓ} f_{ℓ}$ is simply called a reconstruction kernel (i.e., not optimal).

3.2 The Discrete Reconstruction Procedure with Error Analysis

For the evaluation of the reconstruction kernel, we assume a finite dimensional discretization satisfying the following definition.

Definition 2

Let $H_{B}$ be a bounded subset of a Hilbert space $H$ whose inner-product is denoted by ${⟨, ⟩}_{H}$ . A mapping,

T_{H_{B}} : H_{B} \to C^{1 \times P}, P \in N,

is called an inner-product preserving discretization of $H_{B}$ of accuracy $ϵ_{dis} > 0$ if

∣ ∣ ∣ \to f {\to g}^{*} - {⟨ f, g ⟩}_{H} ∣ ∣ ∣ < ϵ_{dis}, \forall f, g \in H_{B},

(13)

where $\to f = T_{H_{B}} (f)$ and $\to g = T_{H_{B}} (g)$

. The vectors

\to f

and

\to g

are called the corresponding inner-product preserving discretizations of

f

and

g

H_{B}

Let $_{ℓ},_{m} \in C^{1 \times P}$ denote inner-product preserving discretizations of some $f_{ℓ}, g_{ℓ} \in H_{B}$ , respectively. Let $η \in H$ denote the orthogonal projection of $h$ on the subspace spanned by $H_{B}$ . By definition (2) we obtain

∣ ∣ ∣_{ℓ} {\to η}^{*} - {⟨ f_{ℓ}, h ⟩}_{H} ∣ ∣ ∣, ∣ ∣_{m} {\to η}^{*} - {⟨ g_{m}, h ⟩}_{H} ∣ ∣ \leq ∥ η ∥ ϵ_{dis},

for all $ℓ \in {1, 2, \dots, L}$ and $m \in {1, 2, \dots, M}$ . Hence, given the information (10) and an approximation of $_{m}$ ,

{ˆ g}_{m} = L \sum ℓ = 1 {ˆ c}_{m, ℓ}_{ℓ},

(14)

we can reconstruct the unknown target coefficients (9) via

b_{m} = {⟨ g_{m}, h ⟩}_{H} \approx {ˆ g}_{m} {\to η}^{*} = L \sum ℓ = 1 {ˆ c}_{m, ℓ} \to f_{ℓ} {\to η}^{*} \approx L \sum ℓ = 1 {ˆ c}_{m, ℓ} a_{ℓ} .

(15)

To evaluate the error of the reconstruction (15) we denote for each $ℓ = 1, 2, \dots, L$ and $m = 1, 2, \dots, M$ the discretization errors

ϵ_{ℓ} = {⟨ f_{ℓ}, h ⟩}_{H} -_{ℓ} {\to η}^{*}, δ_{m} = {⟨ g_{m}, h ⟩}_{H} -_{m} {\to η}^{*},

and obtain the following estimate

b_{m} - L \sum ℓ = 1 {ˆ c}_{m, ℓ} a_{ℓ} = δ_{m} + (_{m} - {ˆ g}_{m}) {\to η}^{*} + L \sum ℓ = 1 {ˆ c}_{m, ℓ} ϵ_{ℓ},

where our assumption ensure that $| δ_{m} |, | ϵ_{ℓ} | \leq ∥ η ∥ ϵ_{dis}$ . Note that $(_{m} - {ˆ g}_{m}) {\to η}^{*}$ is the projection error which can not be reduced if the set of information functionals, ${f_{ℓ}}$ , is predetermined.

To control the error we impose the following regularization constraint

| {ˆ c}_{m, ℓ} | \leq ϵ_{evl} / ϵ_{dis}, ϵ_{evl} \geq ϵ_{dis},

(16)

where $ϵ_{evl}$ is a chosen or given evaluation error bound. Thus, we obtain

∣ ∣ ∣ ∣ b_{m} - L \sum ℓ = 1 {ˆ c}_{m, ℓ} a_{ℓ} ∣ ∣ ∣ ∣ \leq ∥ η ∥ \cdot ϵ_{dis} + ∣ ∣ (_{m} - {ˆ g}_{k}) {\to η}^{*} ∣ ∣ + L \cdot ∥ η ∥ \cdot ϵ_{evl} .

(17)

The regularization constraint (16) explicitly defines the convex set $C_{m}^{L}$ in (11) as

C_{m}^{L} = {(c_{1}, c_{2}, \dots, c_{L}) ∣ ∣ ∣ | c_{ℓ} | \leq \frac{ϵ_{evl}}{ϵ_{dis}}} .

The error estimate (17) implies that it is sufficient to obtain an approximation (14) satisfying $∥ ∥ (_{m} - {ˆ g}_{k}) ∥ ∥ \leq ϵ_{evl}$ . Indeed, in that case (17) reduces to

∣ ∣ ∣ ∣ b_{m} - L \sum ℓ = 1 {ˆ c}_{m, ℓ} a_{ℓ} ∣ ∣ ∣ ∣ \leq ∥ η ∥ \cdot ϵ_{dis} + (L + 1) \cdot ∥ η ∥ \cdot ϵ_{% evl} .

Often, the summation of evaluation errors $\sum_{ℓ = 1}^{L} {ˆ c}_{m, ℓ} ϵ_{ℓ}$ is not cumulative. Thus, the overall error is typically $O (ϵ_{evl})$ assuming $∥ η ∥ = O (1)$ . This facilitates an efficient approximation technique which is based performing successive singular value decompositions on subsets of information functionals. The technique was presented in NFM-Optimal-Reconstruction and demonstrated high stability and good convergence properties. Further details including different variants of the technique can be found in NFM-Thesis .

3.3 Inner-Product Preserving Discretization and Numerical Integration

Obtaining inner-product preserving discretizations of surface functionals is a fundamental issue. Let us focus on the case, as in this work, where ${f_{ℓ}}$ and ${g_{m}}$ are smooth functions in some $H = L^{2}$ space with an inner-product,

{⟨ f, g ⟩}_{H} = \int_{D} f g^{*} ω d s,

(18)

where $D \subset R^{d}$ is compact and Jordan measurable, $g^{*}$ is the complex conjugate of $g$ and $ω$ is a proper weight function.

To numerically compute the integrals 18, we observe that it is sufficient to employ an integration rule which is accurate on the finite dimensional subspace of smooth functions spanned by ${f_{ℓ}} \cup {g_{m}}$ . Thus, we assume the availability of a standard rule of the following form

{⟨ f, g ⟩}_{H} \approx N \sum i = 1 f (σ^{(i)}) g^{*} (σ^{(i)}) ω^{(i)} \forall f, g \in span ({f_{ℓ}} \cup {g_{m}}),

with integration nodes $σ^{(1)}, \dots, σ^{(N)}$ contained in $D$ and real positive weights $ω^{(1)}, \dots, ω^{(N)}$ . Discretizing an element $f \in H$ as a weighted gridfunction

\to f = (f (σ^{(1)}) \sqrt{ω^{(1)}}, \dots, f (σ^{(N)}) \sqrt{ω^{(N)}}),

(19)

essentially, satisfies the inner-product preserving assumption (13) if $N \in N$ is sufficiently large. The number of elements $N$ required for an effective inner-product preserving discretization depends on the convergence rate of the numerical integration formula and, typically, under some smoothness assumption of the integrands. Indeed, if the weight function $ω$ encompasses all the singularities while $f (s)$ and $g (s)$ are analytic, a Gaussian numerical integration rule with respect to $ω$ ensures exponential convergence. Note that the weights of Gaussian rules are always positive and uniformly bounded. See davis2007methods for more details.

4 Surface Embedding of $2 d$ Random Star-Shaped Obstacles

In this section the main theoretical contribution of this paper is presented. The first two subsections cover the setting of the problem, where Subsection (4.1) defines the random shape properties, and Subsection (4.2) covers relevant components of the generalized Polynomial Chaos (gPC) expansion theory. The chosen framework leads to a reconstruction problem in a Hilbert space. A concise discussion on the disadvantages of naive discretization of the reconstruction problem concludes Subsection (4.2).

In Subsection (4.3) we present an analytic approach for overcoming the difficulties associated with the naive discretization approach. Using the Coarea formula we construct a low-dimensional spatial embedding within the family of random surfaces, which facilitates a natural choice for setting a cubature rule in a compact region of $R^{2}$

. The chosen integration weight function is a strictly positive minimal variance quantity encompassing the irregularities of the family of random surfaces.

In Subsections (4.4) and (4.5) we focus on the case of a single random variable describing the randomness of the object. Using the implicit function theorem we obtain explicit formulas including full characterization of the singular behavior of the integration weight function. The usage of the single random variable formulation as a building block for the more general case of multiple random variables is considered and discussed in Section (6).

Subsections (4.6) and (4.7) are devoted to the demonstration of the preceding theoretical parts on a model problem of a randomly oriented elliptic cylinder. The random orientation problem is a very simple ’toy’ problem. However, it allows us us to demonstrate in an affable fashion the implementation of the theory.

4.1 The $2 d$ Random Shape Setting

For brevity, we focus on the sound-soft case and assume that $B$ represents a star-shaped obstacle in $R^{2}$ whose boundary, $S$

, depends smoothly on a real valued vector of mutually independent and continuous random variables

Z = (Z_{1}, \dots, Z_{P}), P \in N .

The boundary $S$ is, however, not assumed to be uniformly smooth in the spatial domain. We assume that each random variable $Z_{p}$

has finite even moments

E [{Z_{p}}^{2 n}] = \int_{I_{Z_{p}}} z_{p}^{2 n} \frac{d F_{Z_{p}}}{d z_{p}} d z_{p} < \infty, n \in {0, 1, \dots, N},

(20)

where $I_{Z_{p}}$ is the support of $Z_{p}$ and $\frac{d F_{Z_{p}} (z_{p})}{d z_{p}}$ is the probability density function of $Z_{p}$ . Property (20) effectively ensures the existence of surface functionals suitable for the reconstruction of the scattering coefficients.

Our assumption that the obstacle is star-shaped for any realization of the random vector $Z$ , ensures that its boundary possesses a polar representation,

S (Z) = {(ρ (θ; Z) \cdot cos θ, ρ (θ; Z) \cdot sin θ) | θ \in [0, 2 π]},

(21)

and the existence of two positive radial bounds, $r^{max}$ and $r^{min}$ , satisfying

0 < r^{min} = inf θ, z ρ (θ; z) < sup θ, z ρ (θ; z) = r^{max} < \infty .

(22)

Thus, as illustrated in Figure 1, $S (Z)$ is confined to the transition region,

R^{t r a} = {r \in R^{2} ∣ ∣ r^{min} < | r | < r^{max}} \subset R^{2} .

(23)

Figure 1: The Transition Region $R^{t r a}$ .

4.2 Random Shape and Generalized Polynomial Chaos Expansion

Given our assumptions we observe that the scattering coefficients (7) are finite dimensional random fields,

b_{m} = b_{m} (Z) .

A common method to approximate these fields is to obtain their generalized Polynomial Chaos (gPC) expansions,

b_{m} (Z) \approx b_{m}^{N} (Z) = \sum | n | \leq N b_{m, n} P_{n} (Z),

(24)

where $n = (n_{1}, \dots, n_{P})$ is a multi-index and ${P_{n}}$ is an orthogonal basis of the inner-product space induced by the probability density function of $Z$ ,

{⟨ ϕ (Z), ψ (Z) ⟩}_{L_{d F_{Z}}^{2}} = \int_{I_{Z}} ϕ (z) ψ (z) d F_{Z} (z), d F_{Z} (z) = P \prod k = 1 d F_{Z_{p}} (z_{p}),

whose support is $I_{Z} = I_{Z_{1}} \times \dots \times I_{Z_{P}}$ . Hence, the expansion coefficients are readily available by the orthogonality via

b_{m, n} = \frac{1}{γ_{n}} {⟨ b_{m}, P_{n} ⟩}_{L_{d F_{Z}}^{2}}, γ_{n} = {⟨ P_{n}, P_{n} ⟩}_{L_{d F_{Z}}^{2}} .

For a randomly shaped obstacle each coefficient in the gPC expansion (24) is a target functional of the following form

b_{m, n} = \int_{I_{Z}} \int_{S (z)} g_{m} (r) h (z, ρ) d S (z) P_{n} (z) d F_{Z} (z) .

(25)

Using the polar form representation (21) whose associated induced volume form on the surface is $d S (θ; z) = \sqrt{ρ^{2} + ρ_{θ}^{2}} d θ$ , the general representation (25) can be explicitly written as

b_{m, n} = \int_{I_{Z}} \int_{[0, 2 π]} g_{m} (ρ) h (z, θ) s (θ; z) ρ d θ P_{n} (z) d F_{Z} (z),

(26)

where the normalized metric tensor in polar coordinates is given by

s (θ; z) = \frac{1}{ρ (θ; z)} \cdot \frac{\partial S (θ; z)}{\partial θ} = \sqrt{1 + {(\frac{1}{ρ} \cdot \frac{\partial ρ}{\partial θ})}^{2}} .

In principle, we need to devise a discretization scheme for (26) and apply the optimal reconstruction procedure of Section (3). However, $ρ \in S (Z)$ inherits any irregularity of family of surfaces; e.g., lack of smoothness and oscillatory behaviour, which often necessitates specialized high-order discretization of the surface $S (Z)$ . Additionally, discretizing the random surface integral with a grid of numerical integration nodes has to be realized for every grid point in the parameters domain, $I_{Z}$ . Hence, in general, the practical implementation of an inner-product preserving discretization satisfying (13) is a difficult task. An analytic approach for overcoming this fundamental difficulty is presented in the next subsection.

4.3 Random Surface Embedding and the Coarea Formula

In this subsection we present an analytic approach for producing inner-product preserving discretizations of functionals of the form of (26). The key idea is to apply a change of variables transforming (26) to the following equivalent representation

b_{m, n} = \iint_{R^{t r a}} g_{m} (r) ϕ_{n} (h) ω (r) r d θ d r,

(27)

where the weight function $ω (r) > 0$ , is proportional to the conditional expectation of $s (θ; z)$ given the information $ρ (θ; z) = r \in R^{t r a}$ . Thus, $ω (r)$ has minimal variance while, essentially, encompassing the irregularities of the family of random surfaces, $S (Z)$ . The term $ϕ_{n} (h)$ is a linear functional uniformly bounded in $R^{t r a}$ operating on $h$ . A high-order numerical integration rule with respect to $ω (r)$ would serve as a discretization satisfying (13). The transformed representation (27) is obtained by the so-called Coarea Formula Federer which allows us to express the surface integral in terms of the integral of the level sets of another function.

Theorem 4.1

(The Coarea Formula)
Let $D$ be an open Jordan measurable subset of $R^{d + δ d}$ where $d \in N$ and $δ d$ is a non-negative integer. Let $Φ : R^{d + δ d} \to R^{d}$ be a piecewise smooth Lipschitz function, such that the level set,

Φ^{- 1} (r) = {x \in D | Φ (x) = r},

is a piecewise smooth $δ d$ -dimensional manifold in $D \subset R^{d + δ d}$ . Then for any integrable function, $g : D \to R$ , we have

\int_{D} g (x) d x = \int_{R^{d}} (\int_{Φ^{- 1} (r)} \frac{g (x)}{J_{Φ} (x)} d S_{r} (x)) d r,

(28)

where $J_{Φ} (x) = \sqrt{det ({[D_{Φ} (x)]}_{Φ}^{T} \cdot [D_{Φ} (x)])}$ is the Jacobian of $Φ$ , and $d S_{r}$ denotes surface measure of $Φ^{- 1} (r)$ .

Remark 2

The coarea formula expresses the integral of a function $g$ over $D$ in terms of the level sets of the function $Φ$ . The level sets, $Φ^{- 1} (y)$ , are called fibers of the domain $D$ . The formula is a kind of ”curvilinear” version of Fubini’s theorem.

Let us consider the function

Φ (z, θ) = ρ (θ; z) = (ρ (θ; z) \cdot cos θ, ρ (θ; z) \cdot sin θ), Φ : I_{Z} \times [0, 2 π] \to R^{t r a} .

By direct calculations we obtain

D_{Φ} = {(\begin{matrix} \nabla_{z} ρ cos θ & \nabla_{z} ρ sin θ \partial_{θ} (ρ cos θ) & \partial_{θ} (ρ sin θ) \end{matrix})}^{T} \in R^{2 \times (P + 1)}, J_{Φ} (z, θ) = ρ \cdot | \nabla_{z} ρ |,

where

\nabla_{z} ρ = {(\partial_{z_{1}} ρ, \dots, \partial_{z_{P}} ρ)}_{z_{1}}^{T}, {| \nabla_{z} ρ |}_{z}^{2} = P \sum p = 1 {(\partial_{z_{p}} ρ)}_{z_{p}}^{2} .

Now, for applying (28) on (26) with the chosen implicit function, $Φ$ , we can only consider spatial points $r \in R^{t r a}$ whose associated level set,

Φ^{- 1} (r) = {(z, θ) | Φ (z, θ) = r},

contains at least one smooth $(P - 1)$ -dimensional manifold in $R^{P + 1}$ . Explicitly, these points satisfy $J_{Φ} (z, θ) \neq 0$ where $Φ^{- 1} (r)$ is non-empty. Thus, we obtain the following representation,

	$b_{m, n}$	$= \iint_{R^{t r a}} g_{m} (r) \int_{I_{Z_{r}}} h (z, θ) P_{n} (z) \frac{s (θ; z)}{\| \nabla_{z} ρ \|} P \prod p = 1 \frac{d F_{Z_{p}}}{d z_{p}} d S_{r} r d θ d r$		(29)
		$+ \int_{E^{t r a}} g_{m} (ρ) h (z, θ) s (θ; z) ρ d θ P_{n} (z) d F_{Z} (z),$		(30)

where the domain of integration of the inner integral in (29) is given by

I_{Z_{r}} = {z \in I_{Z} | ρ (θ; z) = r \in R^{t r a}, | \nabla_{z} ρ | \neq 0} \subset I_{Z},

(31)

and the domain of integration of (30) is defined by the subset of irregular points of the set $Φ^{- 1} (r)$ ,

E^{t r a} = {(z, θ) \in I_{Z} \times [0, 2 π] | ρ (θ; z) = r \in R^{t r a}, | \nabla_{z} ρ | = 0} \subseteq Φ^{- 1} (r) .

(32)

Note that $E^{t r a}$ or $I_{Z_{r}}$ (for certain values of $r$ ) can be empty sets, and in that case the associated integral is taken to be zero.

The advantage of the representation (29,30) is that the wave function, $g_{m}$ , in (29) is no longer composed with the boundary and does not inherit its irregular properties. The subset of irregular points, $E^{t r a}$ (32), defines portions of the random surface, $S (Z)$ (21), which are independet on $Z$ ; i.e., non-random, thus reduces to an integral of the following general form,

\int_{E^{t r a}} [g_{m} (ρ (θ)) h (z, θ) s (θ) ρ (θ) d θ] P_{n} (z) d F_{Z} (z),

whose discretization is straitforward.

A major challenge is to efficiently evaluate the inner integral in (29),

\int_{I_{Z_{r}}} h (z, θ) P_{n} (z) \frac{s (θ; z)}{| \nabla_{z} ρ |} K \prod k = 1 \frac{d F_{Z_{p}}}{d z_{p}} d S_{r} .

(33)

This integral can become infinite since ${| \nabla_{z} ρ |}_{z}^{- 1}$ and $\prod_{k = 1}^{K} \frac{d F_{Z_{p}}}{d z_{p}}$ are, essentially, singular. Accordingly, using the following weight function

ω (r) = \int_{I_{Z_{r}}} \frac{s (θ; z)}{| \nabla_{z} ρ |} K \prod k = 1 \frac{d F_{Z_{p}}}{d z_{p}} d S_{r},

(34)

we have that

b_{m, n} = \iint_{R^{t r a}} g_{m} (r) ϕ_{n} (h) ω (r) r d θ d r + b_{m, n}^{% irregular},

(35)

where $b_{m, n}^{irregular}$ denotes the irregular component (30) and

ϕ_{n} (h) = ω^{- 1} (r) \cdot \int_{I_{Z_{r}}} h (z, θ) P_{n} (z) \frac{s (θ; z)}{| \nabla_{z} ρ |} K \prod k = 1 \frac{d F_{Z_{p}}}{d z_{p}} d S_{r},

(36)

is a linear functional operating on the surface density $h$ . If $ϕ_{n} (h)$ is bounded then the choice (34) implies that (33) is effectively desingularized.

To show that $ϕ_{n} (h)$ is a bounded linear functional, let us assume for simplicity that $E^{t r a}$ (32) is an empty set. In that case, we we observe that (34) is, in fact, proportional to the conditional expectation of $s (θ; z)$ given $ρ (θ; z) = r$ ,

ω (r) = E [s (θ; z) | ρ (θ, z) = r] \cdot f_{ρ} (z),

where $f_{ρ} (z)$ is the probability density function of $ρ (θ; Z)$ . It is well known that conditional expectation is a minimum variance predictor as a function of the given information. Hence, the choice (34) implies minimization of oscillatory behaviour of $s (θ; z)$ as a function of $ρ (θ; z)$ . Employing a similar argument we obtain that the linear functional (36) is, in fact,

ϕ_{n} (h) = \frac{E [h (z, θ) P_{n} (z) s (θ; z) | ρ (θ; z) = r]}{E [s (θ; z) | ρ (θ; z) = r]} .

Hence, the Cauchy-Schwarz inequality implies that

| ϕ_{n} (h) | \leq (max z \in I_{Z_{r}} | P_{n} (z) |) \cdot E [| h (z, θ) | | ρ (θ; z) = r],

which shows that the linear functional, $ϕ_{n}$ , is, indeed, bounded. A similar argument can be applied to show that the functional $ϕ_{n}$ remains bounded when $E^{t r a}$ is not an empty set.

Setting the integration grid for (35) can be done in two stages. First we obtain a cubature rule with nodes ${r^{(i, j)}} \subset R^{t r a}$ and corresponding cubature weights ${ω^{(i, j)}}$ with respect to the weight function $ω (r)$ , regardless of $z$ . In the second stage we identify the integration grid in $I_{Z_{r}}$ which corresponds ${r^{(i, j)}}$ for the evaluation of $ϕ_{n} (h)$ .

4.4 Explicit Representation for a Single Random Variable

Let us consider a simplified case of one-dimensional random vector,

Z = Z \in I_{Z} \subset R .

for which $d F_{Z} = d F_{Z}$ . We also assume for brevity, that $E^{t r a}$ (32) is an empty set. Thus, the target functional (29,30) reduces to

b_{m, n} = \int_{I_{Z}} \int_{[0, 2 π]} g_{m} (ρ) h (z, θ) s (θ; z) ρ d θ P_{n} (z) d F_{Z} (z) .

(37)

We will show that in this case the functional $ϕ_{n} (h)$ (36) can be explicitly represented. This approach can serve as a building block for the case of a general random vector, which is discussed in Section (6). The assumption $E^{t r a} = \emptyset$ does not imply a loss of generality, since the discretization of the irregular part (30) is carried out directly without applying the coarea formula.

Our assumptions imply that for any $r \in R^{t r a}$ the fiber set $I_{Z_{r}}$ (31) is either an empty set or composed of a finite set of discrete points; i.e., a zero-dimensional sub-surface. Thus, we obtain the following explicit representation of (29) for a single random variable,

\iint_{R^{t r a}} g_{m} (r) ⎡ ⎣ \sum z_{r} \in I_{Z_{r}} h (z_{r}, θ) P_{n} (z_{r}) \frac{s (θ; z_{r})}{{| \partial_{z} ρ |}_{z = z_{r}}} \frac{d F_{Z}}{d z} (z_{r}) ⎤ ⎦ r d θ d r,

(38)

where $I_{Z_{r}}$ is the reduction of $I_{Z_{r}}$ to a zero dimensional subset of $I_{Z}$ ,

I_{Z_{r}} (r) = {z \in I_{Z} | ρ (θ; z) = r} .

(39)

To efficiently evaluate (38), we employ the implicit function theorem Dini which ensures the following identities

1 = \partial_{r} ρ (θ; z_{r}) = {[\partial_{z} ρ (θ; z)]}_{z = z_{r}} \partial_{r} z_{r},

\partial_{θ} ρ (θ; z_{r}) = {[\partial_{z} ρ (θ; z)]}_{z = z_{r}} \partial_{θ} z_{r} + {[\partial_{θ} ρ (θ; z)]}_{z = z_{r}} .

Thus,

{[\partial_{z} ρ (θ; z)]}_{z = z_{r}} = \frac{1}{\partial_{r} z_{r}}, {[\partial_{θ} ρ (}_{θ}

Low-Dimensional Spatial Embedding Method for Shape Uncertainty Quantification in Acoustic Scattering

Authors

A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem

Acoustic scattering by impedance screens/cracks with fractal boundary: well-posedness analysis and boundary element approximation

A higher order perturbation approach for electromagnetic scattering problems on random domains

Uncertainty Quantification in Three Dimensional Natural Convection using Polynomial Chaos Expansion and Deep Neural Networks

Helmholtz scattering by random domains: first-order sparse boundary elements approximation

Direct Volume Rendering with Nonparametric Models of Uncertainty

Scaled boundary cubature scheme for numerical integration over polytopes and curved solids. Part I: Two-dimensional domains

1 Introduction

2 Null-Field Reconstruction for Time Harmonic Wave Scattering

2.1 Acoustic Scattering by Impenetrable Obstacles

2.2 Null-Field Theory Fundamentals

2.3 Reconstruction of Surface Functionals

3 Optimal Reconstruction with A-priori Error Estimate

3.1 The Reconstruction Problem and Reconstruction Kernels

Definition 1

Remark 1

3.2 The Discrete Reconstruction Procedure with Error Analysis

Definition 2

3.3 Inner-Product Preserving Discretization and Numerical Integration

4 Surface Embedding of $2 d$ Random Star-Shaped Obstacles

4.1 The $2 d$ Random Shape Setting

4.2 Random Shape and Generalized Polynomial Chaos Expansion

4.3 Random Surface Embedding and the Coarea Formula

Theorem 4.1

Remark 2

4.4 Explicit Representation for a Single Random Variable

Low-Dimensional Spatial Embedding Method for Shape Uncertainty Quantification in Acoustic Scattering

Authors

A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem

Acoustic scattering by impedance screens/cracks with fractal boundary: well-posedness analysis and boundary element approximation

A higher order perturbation approach for electromagnetic scattering problems on random domains

Uncertainty Quantification in Three Dimensional Natural Convection using Polynomial Chaos Expansion and Deep Neural Networks

Helmholtz scattering by random domains: first-order sparse boundary elements approximation

Direct Volume Rendering with Nonparametric Models of Uncertainty

Scaled boundary cubature scheme for numerical integration over polytopes and curved solids. Part I: Two-dimensional domains

This week in AI

1 Introduction

2 Null-Field Reconstruction for Time Harmonic Wave Scattering

2.1 Acoustic Scattering by Impenetrable Obstacles

2.2 Null-Field Theory Fundamentals

2.3 Reconstruction of Surface Functionals

3 Optimal Reconstruction with A-priori Error Estimate

3.1 The Reconstruction Problem and Reconstruction Kernels

Definition 1

Remark 1

3.2 The Discrete Reconstruction Procedure with Error Analysis

Definition 2

3.3 Inner-Product Preserving Discretization and Numerical Integration

4 Surface Embedding of 2d Random Star-Shaped Obstacles

4.1 The 2d Random Shape Setting

4.2 Random Shape and Generalized Polynomial Chaos Expansion

4.3 Random Surface Embedding and the Coarea Formula

Theorem 4.1

Remark 2

4.4 Explicit Representation for a Single Random Variable

4 Surface Embedding of $2 d$ Random Star-Shaped Obstacles

4.1 The $2 d$ Random Shape Setting