## 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 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 denote a bounded domain in representing an impenetrable obstacle with boundary . We denote by the closure of . Let be the unbounded exterior region occupied by a uniform medium. Let denote a general spatial point,

For an incident time-harmonic field ’illuminating’ the obstacle, the scattered field satisfies the following exterior boundary value problem:

(1) |

(2) |

(3) |

where .
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
and is the normal derivative of .
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

(4) |

which holds for any bounded domain with a Lipschitz piecewise smooth boundary , where and are scalar fields, and and denote corresponding normal derivatives.

Let be an outgoing wavefunction and let denote the total field, . Assuming is analytic in , it can be shown by (4) that

Thus, for a sound-soft obstacle ( on )

(5) |

while for a sound-hard obstacle ( on )

(6) |

Using (5) or (6), an infinite set of equations can be produced from which or on 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 or on . 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, or on .
Often such features are the outcomes of functionals in an
appropriate *Hilbert space*.
Indeed, let denote the complex conjugate of
or on . Then for a general non-smooth surface
, the surface density belongs to the complex Hilbert space
whose inner-product and norm are defined by

where denotes the complex conjugate of and is the
induced volume form on the surface.
Recall that by *Riesz representation theorem* any bounded linear
functional operating on surface densities,
is of the form . We call such
functionals surface functionals.

In this work we focus on the estimation of the *scattering coefficients*
of the expansion of to cylinder harmonics in
the case.
These coefficients, denoted by , satisfy

where denotes the 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

(7) |

where denotes the th-order
*Bessel function of the first kind*,
and denotes the complex conjugate
of the surface density .
In practice only satisfying

(8) |

are required for an accurate description of , 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 on . 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

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 . Hence, given the information

the outcome of the target functional can be approximate by

Obtaining the coefficients 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 be a complex Hilbert space, whose inner-product is denoted by . The reconstruction problem is to approximate a finite set of target functionals

(9) |

by a given finite set of information functionals,

(10) |

where the element is unknown.

###### Definition 1

Let denote the norm induced by in , and let denote a closed convex subset of . A linear combination whose coefficients satisfy the minimality condition

(11) |

is called an *optimal reconstruction kernel* of the target functional
by the information functionals over .

###### 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 . Note that almost no prior knowledge on the element 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

(12) |

with respect to some predetermined threshold, .
In that case the linear combination
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 be a bounded subset of a Hilbert space whose inner-product is denoted by . A mapping,

is called an inner-product preserving discretization of of accuracy if

(13) |

where and

. The vectors

and are called the corresponding inner-product preserving discretizations of and on .Let denote inner-product preserving discretizations of some , respectively. Let denote the orthogonal projection of on the subspace spanned by . By definition (2) we obtain

for all and . Hence, given the information (10) and an approximation of ,

(14) |

we can reconstruct the unknown target coefficients (9) via

(15) |

To evaluate the error of the reconstruction (15) we denote for each and the discretization errors

and obtain the following estimate

where our assumption ensure that . Note that is the projection error which can not be reduced if the set of information functionals, , is predetermined.

To control the error we impose the following regularization constraint

(16) |

where is a chosen or given evaluation error bound. Thus, we obtain

(17) |

The regularization constraint (16) explicitly defines the convex set in (11) as

The error estimate (17) implies that it is sufficient to obtain an approximation (14) satisfying . Indeed, in that case (17) reduces to

Often, the summation of evaluation errors
is not cumulative. Thus, the overall error is typically
assuming .
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 and are smooth functions in some space with an inner-product,

(18) |

where is compact and Jordan measurable,
is the complex conjugate of
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 . Thus, we assume the availability of a standard rule of the following form

with *integration nodes*
contained in and real positive
weights .
Discretizing an element as a weighted
*gridfunction*

(19) |

essentially, satisfies the inner-product preserving assumption
(13) if is sufficiently large.
The number of elements 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
and 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 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

. 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 Random Shape Setting

For brevity, we focus on the sound-soft case and assume that represents a star-shaped obstacle in whose boundary,

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

The boundary is, however, not assumed to be uniformly smooth in the spatial domain. We assume that each random variable

has finite even moments

(20) |

where is the support of
and is the *probability density function* of .
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 , ensures that its boundary possesses a polar representation,

(21) |

and the existence of two positive radial bounds, and , satisfying

(22) |

Thus, as illustrated in Figure 1, is confined to the transition region,

(23) |

### 4.2 Random Shape and Generalized Polynomial Chaos Expansion

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

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

(24) |

where is a multi-index and is an orthogonal basis of the inner-product space induced by the probability density function of ,

whose support is . Hence, the expansion coefficients are readily available by the orthogonality via

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

(25) |

Using the polar form representation (21) whose associated induced volume form on the surface is , the general representation (25) can be explicitly written as

(26) |

where the normalized metric tensor in polar coordinates is given by

In principle, we need to devise a discretization scheme for (26) and apply the optimal reconstruction procedure of Section (3). However, 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 . 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, . 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

(27) |

where the weight function , is proportional to
the conditional expectation of given the information
. Thus, has minimal variance
while, essentially, encompassing the irregularities of the family of random surfaces, .
The term is a linear functional uniformly bounded in
operating on . A high-order numerical integration rule
with respect to 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 be an open *Jordan measurable* subset of
where and is a non-negative integer.
Let be a piecewise smooth
*Lipschitz* function, such that the level set,

is a piecewise smooth -dimensional manifold in . Then for any integrable function, , we have

(28) |

where is the Jacobian of , and denotes surface measure of .

###### Remark 2

The coarea formula expresses the integral of a function over
in terms of the level sets of
the function . The level sets, , are called
fibers of the domain . The formula is a kind of ”curvilinear” version of *Fubini’s theorem*.

Let us consider the function

By direct calculations we obtain

where

Now, for applying (28) on (26) with the chosen implicit function, , we can only consider spatial points whose associated level set,

contains at least one smooth -dimensional manifold in . Explicitly, these points satisfy where is non-empty. Thus, we obtain the following representation,

(29) | ||||

(30) |

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

(31) |

and the domain of integration of (30) is defined by the subset of irregular points of the set ,

(32) |

Note that or (for certain values of ) 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, , in (29) is no longer composed with the boundary and does not inherit its irregular properties. The subset of irregular points, (32), defines portions of the random surface, (21), which are independet on ; i.e., non-random, thus reduces to an integral of the following general form,

whose discretization is straitforward.

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

(33) |

This integral can become infinite since and are, essentially, singular. Accordingly, using the following weight function

(34) |

we have that

(35) |

where denotes the irregular component (30) and

(36) |

is a linear functional operating on the surface density . If is bounded then the choice (34) implies that (33) is effectively desingularized.

To show that is a bounded linear functional, let us assume for simplicity that (32) is an empty set. In that case, we we observe that (34) is, in fact, proportional to the conditional expectation of given ,

where is the probability density function of . 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 as a function of . Employing a similar argument we obtain that the linear functional (36) is, in fact,

Hence, the *Cauchy-Schwarz* inequality implies that

which shows that the linear functional, , is, indeed, bounded. A similar argument can be applied to show that the functional remains bounded when 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 and corresponding cubature weights with respect to the weight function , regardless of . In the second stage we identify the integration grid in which corresponds for the evaluation of .

### 4.4 Explicit Representation for a Single Random Variable

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

for which . We also assume for brevity, that (32) is an empty set. Thus, the target functional (29,30) reduces to

(37) |

We will show that in this case the functional (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 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 the fiber set (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,

(38) |

where is the reduction of to a zero dimensional subset of ,

(39) |

Comments

There are no comments yet.