# An efficient numerical method for acoustic wave scattering in random media

## Abstract

This paper is concerned with developing efficient numerical methods for acoustic wave scattering in random media which can be expressed as random perturbations of homogeneous media. We first analyze the random Helmholtz problem by deriving some wave-number-explicit solution estimates. We then establish a multi-modes representation of the solution as a power series of the perturbation parameter and analyze its finite modes approximations. Based on this multi-modes representation, we develop a Monte Carlo interior penalty discontinuous Galerkin (MCIP-DG) method for approximating the mode functions, which are governed by recursively defined nearly deterministic Helmholtz equations. Optimal order error estimates are derived for the method and an efficient algorithm, which is based on the LU direct solver, is also designed for efficiently implementing the proposed multi-modes MCIP-DG method. It is proved that the computational complexity of the whole algorithm is comparable to that of solving one deterministic Helmholtz problem using the LU director solver. Numerical experiments are provided to validate the theoretical results and to gauge the performance of the proposed numerical method and algorithm.

## 1Introduction

Partial differential equations with random coefficients arise naturally in the modeling of many physical phenomena. This is due to the fact that some level of uncertainty is usually involved if the knowledge of the physical behavior is not complete or when noise is present in the experimental measurements. In recent years, substantial progress has been made in the numerical approximation of such PDEs due to the significant development in computational resources. We refer to [1] and references therein for more details.

In this paper, we consider the propagation of the acoustic wave in a medium where the wave velocity is characterized by a random process. More precisely, we study the approximation of the solution to the following Helmholtz problem:

where is the wavenumber, and is a convex bounded polygonal domain with boundary . Let be a probability space with sample space , algebra and probability measure . For each fixed , the refractive index is a real-valued random variable defined over . We assume that the medium is a small random perturbation of a uniform background medium in the sense that

Here represents the magnitude of the random fluctuation, and is some random process satisfying

For notation brevity we only consider the case that is real-valued. However, we note that the results of this paper are also valid for complex-valued . On the boundary , a radiation boundary condition is imposed to absorb incoming waves [6]. Here denotes the unit outward normal to , and stands for the normal derivative of . The boundary value problem – arises in the modeling of the wave propagation in complex environments, such as composite materials, oil reservoir and geological basins [9]. In such instances, it is of practical interest to characterize the uncertainty of the wave energy transport when the medium contains some randomness. In particular, we are interested in the computation of some statistics of the wave field, e.g, the mean value of the solution .

To solve stochastic (or random) partial differential equations (SPDEs) numerically, the simplest and most natural approach is to use the Monte Carlo method, where a set of independent identically distributed (i.i.d.) solutions are obtained by sampling the PDE coefficients, and the mean of the solution is calculated via a statistical average over all the sampling in the probability space [3]. An alternative is the stochastic Galerkin method, where the SPDE is reduced to a high dimensional deterministic equation by expanding the random field in the equation using the Karhunen-Loève or Wiener Chaos expansions. We refer the reader to [1] for detailed discussions. However, it is known that a brute-force Monte Carlo or stochastic Galerkin method applied directly to the Helmholtz equation with random coefficients is computationally prohibitive even for a moderate wavenumber , since a large number of degrees of freedom is involved in the spatial discretization. It is apparent that in such cases, the Monte Carlo method requires solving a PDE with many sampled coefficients, while the high dimensional deterministic equation associated with the stochastic Galerkin method will be too expensive to be solved.

In this paper, we propose an efficient numerical method for solving the Helmholtz problem – when the medium is weakly random defined by . A multi-modes representation of the solution is derived, where each mode is governed by a Helmholtz equation with deterministic coefficients and a random source. We develop a Monte Carlo interior penalty discontinuous Galerkin (MCIP-DG) method for approximating the mode functions. In particular, we take the advantage that the coefficients of the Helmholtz equation for all the modes are identical, hence the associated discretized equations share the same constant coefficient matrix. Using this crucial fact, it is observed that an LU direct solver for the discretized equations leads to a tremendous saving in the computational costs, since the LU decomposition matrices can be used repeatedly, and the solutions for all modes and all samples can be obtained in an efficient way by performing simple forward and backward substitutions. Indeed, it turns out that the computational complexity of the proposed algorithm is comparable to that of solving one deterministic Helmholtz problem using the LU direct solver.

The rest of the paper is organized as follow. A wave-number-explicit estimate for the solution of the random Helmholtz equation is established in Section 2. In Section 3, we introduce the multi-modes expansion of the solution as a power series of and analyze the error estimation for its finite-modes approximation. The Monte Carlo interior penalty discontinuous Galerkin method is presented in Section Section 4, where the error estimates for the approximation of each mode function is also obtained. In Section 5, a numerical procedure for solving – is described and its computational complexity is analyzed in detail. In addition, we derive an optimal order error estimates for the proposed procedure. Several numerical experiments are provided in Section 6 to demonstrate the efficiency of the method and to validate the theoretical results.

## 2PDE analysis

### 2.1Preliminaries

Standard function and space notations are adopted in this paper. For example, denotes the complex-valued Sobolev space and . stands for the standard inner product on the complex-valued space for any subset of . and denote generic constants which are independent of and the mesh parameter . We also define spaces

Without loss of generality, we assume that the domain . Throughout this paper we also assume that is a star-shaped domain with respect to the origin in the sense that there exists a positive constant such that

Let be a probability space on which all the random variables of this paper are defined. denotes the expectation operator. The abbreviation *a.s.* stands for *almost surely*.

As it will be needed in the late sections of the paper, in this section we analyze the boundary value problem for the Helmholtz equation with the following slightly more general nonhomogeneous boundary condition:

### 2.2Wave-number-explicit solution estimates

In this subsection we shall derive stability estimates for the solution of problem , which is defined in Definition ?. Our focus is to obtain explicit dependence of the stability constants on the wave number , such wave-number-explicit stability estimates will play a vital role in our convergence analysis in the later sections. We note that wave-number-explicit stability estimates also play a pivotal role in the development of numerical methods, such as finite element and discontinuous Galerkin methods, for deterministic reduced wave equations (cf. [7]). As a byproduct of the stability estimates, the existence and uniqueness of solutions to problem , can be conveniently established.

Setting in yields

Taking the real and imaginary parts and using the definition of , we get

Applying the Cauchy-Schwarz inequality to produces

Thus, holds. Applying Cauchy-Schwarz to yields

To this one can apply with and obtain . The proof is complete.

follows immediately from applying the divergence theorem to

and the fact that . To show , we first recall the following differential identities [4]:

Then follows from adding the above two identities, integrating the sum over and applying the divergence theorem on the left-hand side of the resulting equation.

We are now ready to state and prove our wave-number-explicit estimate for solutions of problem , defined in Definition ?.

To avoid some technicalities, below we only give a proof for the case . For the general case, needs be replaced by its mollification at the beginning of the proof and followed by taking the limit after the integration by parts is done.

Setting in yields

Using and after taking the real part of and regrouping we get

It then follows from Schwarz inequality and the “star-shape” condition and the facts that for and a.s. that

At this point, we note that implies . Setting , and denoting , using Lemma ? we can bound right-hand side as follows:

which is equivalent to

where

Let , , , and , then

If , it is easy to check that . Thus, infers that

for some constant independent of and . This then proves .

By with and we get

Hence, holds.

Finally, it follows from the standard elliptic regularity theory for Poisson equation and the trace inequality (cf. [11]) that

Hence holds. The proof is complete.

As a non-trivial byproduct, the above stability estimates can be used conveniently to establish the existence and uniqueness of solutions to problem – as defined in Definition ?.

The proof is based on the well known Fredholm Alternative Principle (cf. [11]). First, it is easy to check that the sesquilinear form on the right-hand side of satisfies a Gärding’s inequality on the space . Second, to apply the Fredholm Alternative Principle we need to prove that solutions to the adjoint problem of – is unique. It is easy to verify that the adjoint problem is associated with the sesquilinear form

which differs from only in the sign of the last term. As a result, all the stability estimates for problem – still hold for its adjoint problem. Since the adjoint problem is a linear problem (so is problem –), the stability estimates immediately infers the uniqueness. Finally, the Fredholm Alternative Principle then implies that problem – has a unique solution . The proof is complete.

## 3Multi-modes representation of the solution and its finite modes approximations

The first goal of this section is to develop a multi-modes representation for the solution to problem – in terms of powers of the parameter . We first postulate such a representation and then prove its validity by establishing some energy estimates for all the mode functions. The second goal of this section is to establish an error estimate for finite modes approximations of the solution. Both the multi-modes representation and its finite modes approximations play a pivotal role in our overall solution procedure for solving problem – as they provide the theoretical foundation for the solution procedure. Throughout this section, we use to denote the solution to problem – which is proved in Theorem ?.

We start by postulating that the solution has the following multi-modes expansion:

whose validity will be justified later. Without loss of the generality, we assume that and . Otherwise, the problem can be rescaled to this regime by a suitable change of variable. We note that the normalization implies that .

Substituting the above expansion into the Helmholtz equation and matching the coefficients of order terms for , we obtain

Similarly, the boundary condition translates to each mode function as follows:

A remarkable feature of the above multi-modes expansion is that all the mode functions satisfy the same type (nearly deterministic) Helmholtz equation and the same boundary condition. The only difference is that the Helmholtz equations have different right-hand side source terms (all of them except one are random variables), and each pair of consecutive mode functions supply the source term for the Helmholtz equation satisfied by the next mode function. This remarkable feature will be crucially utilized in Section 5 to construct our overall numerical methodology for solving problem –.

Next, we address the existence and uniqueness of each mode function .

For each , the PDE problem associated with is the same type Helmholtz problem as the original problem – (with in the left-hand side of the PDE). Hence, all a priori estimates of Theorem ? hold for each (with its respective right-hand source side function). First, we have

Thus, and hold for .

Next, we use induction to prove that and hold for all . Assume that and hold for all , then

where and denotes the Kronecker delta, and we have used the fact that and

Similarly, for we have

Hence, and hold for . So the induction argument is complete.

We now use and the elliptic theory for Possion problems directly to verify estimate .

Hence, holds for all .

With a priori estimates and in hand, the proof of existence and uniqueness of each follows verbatim the proof of Theorem ?, which we leave to the interested reader to verify. The proof is complete.

Now we are ready to justify the multi-modes representation for the solution of problem –.

The proof consists of two parts: (i) the infinite series on the right-hand side of converges in ; (ii) the limit coincides with the solution . To prove (i), we define the partial sum

Then for any fixed positive integer we have

It follows from Schwarz inequality and that for

Thus, if we have

Therefore, is a Cauchy sequence in . Since is a Banach space, then there exists a function such that

To show (ii), we first notice that by the definitions of and , it is easy to check that satisfies

for all . Where . In other words, solves the following Helmholtz problem:

By and Schwarz inequality we have

Similarly we get

Setting in immediately yields

Thus, is a solution to problem –. By the uniqueness of the solution, we conclude that . Therefore, holds in . The proof is complete.

The above proof also infers an upper bound for the error as stated in the next theorem.

Let , subtracting from we get

In other words, solves the following Helmholtz problem:

By Theorem ? and we obtain for

The proof is complete.

## 4Monte Carlo discontinuous Galerkin approximations of the mode functions

In the previous section, we present a multi-modes representation of the solution and a convergence rate estimate for its finite approximations. These results will serve as the theoretical foundation for our overall numerical methodology for approximating the solution of problem –. To compute following this approach, we need to compute the expectations of the first mode functions . This requires the construction of an accurate and robust numerical (discretization) method to compute the expectations of the solutions to the “nearly” deterministic Helmholtz problems , and , satisfied by the mode functions . The construction of such a numerical method is exactly our focus in this section. We note that due to the multiplicative structure of the right-hand side of , can not be computed directly for . On the other hand, can be computed directly because it satisfies the deterministic Helmholtz equation with the source term and homogeneous boundary condition.

The goal of this section is to develop some *Monte Carlo interior penalty discontinuous Galerkin* (MCIP-DG) methods for the above mentioned Helmholtz problems. Our MCIP-DG methods are the direct generalizations of the deterministic IP-DG methods proposed in [7] for the related deterministic Helmholtz problems. It should be noted that although various numerical methods (such as finite difference, finite element and spectral methods) can be used for the job, the IP-DG methods to be presented below are the only general purpose discretization methods which are unconditionally stable (i.e., stable without mesh constraint) and optimally convergent. This is indeed the primary reason why we choose the IP-DG methods as our spatial discretization methods.

### 4.1DG notations

Let be a quasi-uniform partition of such that . Let denote the diameter of and . denotes the standard broken Sobolev space and denotes the DG finite element space which are defined as

where is the set of all polynomials whose degrees do not exceed a given positive integer . Let denote the set of all interior faces/edges of , denote the set of all boundary faces/edges of , and . The -inner product for piecewise functions over the mesh is naturally defined by

and for any set , the -inner product over is defined by

Let and and assume global labeling number of is smaller than that of . We choose as the unit normal on outward to and define the following standard jump and average notations across the face/edge :

for . We also define the following semi-norms on :

### 4.2IP-DG method for deterministic Helmholtz problem

In this subsection we consider following deterministic Helmholtz problem and its IP-DG approximations proposed in [7].

We note that satisfies the above equations with and . As an interesting byproduct, all the results to be presented in this subsection apply to .

The IP-DG weak formulation for – is defined by (cf. [7]) seeking such that

where

and are piecewise constant nonnegative functions defined on . denotes an orthonormal basis of the edge and denotes the tangential derivative in the direction of .

Following [7] and based on the DG weak formulation , our IP-DG method for problem – is defined by seeking such that

For the above IP-DG method, it was proved in [7] that the method is unconditionally stable and its solutions satisfy some wave-number-explicit stability estimates. Its solutions also satisfy optimal order (in ) error estimates, which are described below.

An immediate consequence of is the following unconditional solvability and uniqueness result.

### 4.3MCIP-DG method for approximating for

We recall that each mode function satisfies the following Helmholtz problem:

where

Clearly, and are random variables for , and . We remark again that due to its multiplicative structure and can not be computed directly for . Otherwise, and would be easily converted into deterministic equations for , as we did early for . In other words, – is a genuine random PDE problem. On the other hand, since all the coefficients of the equations are constants, then the problem is nearly deterministic. Such a remarkable property will be fully exploited in our overall numerical methodology which will be described in the next section.

Several numerical methodologies are well known in the literature for discretizing random PDEs, Monte Carlo Galerkin and stochastic Galerkin (or polynomial chaos) methods and stochastic collocation methods are three of well-known methods (cf. [2] and the references therein). Due to the nearly deterministic structure of –, we propose to discretize it using the Monte Carlo IP-DG approach which combines the classical Monte Carlo method for stochastic variable and the IP-DG method, which is presented in the proceeding subsection, for the spatial variable.

Following the standard formulation of the Monte Carlo method (cf. [2]), let be a (large) positive integer which will be used to denote the number of realizations and be the DG space defined in Section 4.1. For each , we sample i.i.d. realizations of the source term and random medium coefficient , and recursively find corresponding approximation such that

for . Where

We point out that in order for to be computable, and , not and , are used on the right-hand side of . This (small) perturbation on the right-hand side will result in an additional discretization error which must be accounted later, see Section 5.

Next, we approximate by the following sample average

The following lemma is well known (cf. [2]).

To bound , we once again use the induction argument. To avoid some technicalities, we only provide a proof for the case when the mesh size is in pre-asymptotic regime, i.e., .

By and estimate we immediately get

which verifies and for . Suppose and hold for all , we now prove that they also hold for .

Again, by with and estimate we have