# 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.

elmholtz equation, random media, Rellich identity, discontinuous Galerkin method, error estimate, LU decomposition, Monte Carlo method.

65N12, 65N15, 65N30,

## 1 Introduction

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, 2, 3, 15, 16] 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:

(1) | |||||

(2) |

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

(3) |

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 (1)–(2) arises in the modeling of the wave propagation in complex environments, such as composite materials, oil reservoir and geological basins [9, 12]. 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, 2, 5, 15, 16] 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 (1)–(2) when the medium is weakly random defined by (3). 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 4, where the error estimates for the approximation of each mode function is also obtained. In Section 5, a numerical procedure for solving (1)–(2) 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.

## 2 PDE analysis

### 2.1 Preliminaries

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

(4) | ||||

(5) |

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 (1) with the following slightly more general nonhomogeneous boundary condition:

(6) |

### 2.2 Wave-number-explicit solution estimates

In this subsection we shall derive stability estimates for the solution of problem (1),(6) which is defined in Definition 2.1. 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, 8]). As a byproduct of the stability estimates, the existence and uniqueness of solutions to problem (1),(6) can be conveniently established.

Setting in (7) yields

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

(11) | ||||

(12) |

Applying the Cauchy-Schwarz inequality to (12) produces

Thus, (10) holds. Applying Cauchy-Schwarz to (11) yields

To this one can apply (10) with and obtain (9). The proof is complete.

Let , then there hold

(13) | ||||

(14) | ||||

(13) follows immediately from applying the divergence theorem to

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

Then (14) 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.

###### Remark 2.2

(14) could be called a stochastic Rellich identity for the Laplacian.

We are now ready to state and prove our wave-number-explicit estimate for solutions of problem (1),(6) defined in Definition 2.1.

Let be a solution of (1),(6) and be the smallest number such that contains the domain . Then there hold the following estimates:

(15) | ||||

(16) |

provided that . Where is some positive constant independent of and , and

(17) |

Moreover, if and , there also holds

(18) |

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 (7) yields

(19) |

Using (13) and (14) after taking the real part of (19) 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 2.2 we can bound right-hand side as follows:

which is equivalent to

(20) |

where

Let , , , and , then

If , it is easy to check that . Thus, (20) infers that

(21) |

for some constant independent of and . This then proves (15).

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

Hence (18) holds. The proof is complete.

###### Remark 2.3

By the definition of , we see that . In practice, this is not a restrictive condition because is often taken to be proportional to the wave length. Hence, .

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

Let and . For each fixed pair of positive number and such that , there exists a unique solution to problem (7)–(8).

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 (7) 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 (7)–(8) 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 (7)–(8) still hold for its adjoint problem. Since the adjoint problem is a linear problem (so is problem (7)–(8)), the stability estimates immediately infers the uniqueness. Finally, the Fredholm Alternative Principle then implies that problem (7)–(8) has a unique solution . The proof is complete.

###### Remark 2.4

The uniqueness of the adjoint problem can also be proved using the classical unique continuation argument (cf. [13]).

## 3 Multi-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 (1)–(2) 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 (1)–(2) as they provide the theoretical foundation for the solution procedure. Throughout this section, we use to denote the solution to problem (1)–(2) which is proved in Theorem 2.2.

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

(22) |

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 (1) and matching the coefficients of order terms for , we obtain

(23) | ||||

(24) | ||||

(25) |

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

(26) |

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 (1)–(2).

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

Let . Then for each , there exists a unique solution (understood in the sense of Definition 2.1) to problem (24),(26) for and problem (25),(26) for . Moreover, for , satisfies

(27) | ||||

(28) |

where

(29) |

Moreover, if , there also holds

(30) |

where .

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

(31) | ||||

(32) |

Next, we use induction to prove that (27) and (28) hold for all . Assume that (27) and (28) hold for all , then

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

Similarly, for we have

Hence, (27) and (28) hold for . So the induction argument is complete.

We now use (27) and the elliptic theory for Possion problems directly to verify estimate (30).

Hence, (30) holds for all .

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

Now we are ready to justify the multi-modes representation (22) for the solution of problem (1)–(2).