An efficient Monte Carlo interior penalty discontinuous Galerkin method for elastic wave scattering in random media

# An efficient Monte Carlo interior penalty discontinuous Galerkin method for elastic wave scattering in random media

Xiaobing Feng Department of Mathematics, The University of Tennessee, Knoxville, TN 37996, U.S.A. (xfeng@math.utk.edu). The work of this author was partially supported by the NSF grant DMS-1318486.    Cody Lorton Department of Mathematics and Statistics, University of West Florida, Pensacola, FL 32514, U.S.A. (clorton@uwf.edu).
###### Abstract

This paper develops and analyzes an efficient Monte Carlo interior penalty discontinuous Galerkin (MCIP-DG) method for elastic wave scattering in random media. The method is constructed based on a multi-modes expansion of the solution of the governing random partial differential equations. It is proved that the mode functions satisfy a three-term recurrence system of partial differential equations (PDEs) which are nearly deterministic in the sense that the randomness only appears in the right-hand side source terms, not in the coefficients of the PDEs. Moreover, the same differential operator applies to all mode functions. A proven unconditionally stable and optimally convergent IP-DG method is used to discretize the deterministic PDE operator, an efficient numerical algorithm is proposed based on combining the Monte Carlo method and the IP-DG method with the direct linear solver. It is shown that the algorithm converges optimally with respect to both the mesh size and the sampling number , and practically its total computational complexity only amounts to solving a few deterministic elastic Helmholtz equations using a Guassian elimination direct linear solver. Numerical experiments are also presented to demonstrate the performance and key features of the proposed MCIP-DG method.

E

lastic Helmholtz equations, random media, Rellich identity, discontinuous Galerkin methods, error estimates, LU decomposition, Monte Carlo method.

{AMS}

65N12, 65N15, 65N30, 78A40

Elastic wave scattering problems arise from applications in a variety of fields including geoscience, image science, the petroleum industry, and the defense industry, to name a few. Such problems have been extensively studied both analytically and numerically in the past several decades (cf. [16, 19] and the references therein). The material properties of the elastic media in which the wave propagates play a principle role in the methods used to solve the elastic wave scattering problem. Common medium characterizations include homogeneous and isotropic media, inhomogeneous and anisotropic media, and random media. As the characterization of the media becomes more complicated, so do the computations of the solutions of the associated wave equations. In the case of random media, wave forms may vary significantly for different samplings and as a result, stochastic quantities of interests such as the mean, variance, and/or higher order moments must often be sought.

In this paper we are concerned with developing efficient numerical methods for solving the elastic Helmholtz equations with random coefficients, which models the propagation in random media of elastic waves with a fixed frequency. Specifically, we consider the following random elastic Helmholtz problem:

 (1) −\rm div\,(σ(u(ω,⋅)))−k2α2(ω,⋅)u(ω,⋅) =f(ω,⋅) in D, (2) σ(u(ω,⋅))ν+ikAu(ω,⋅) =0 on ∂D,

for a.e. . Here is the stress tensor defined by

 σ(u(ω,⋅)) :=2μ∇su(ω,⋅)+λ\rm div\,u(ω,⋅)I, ∇su(ω,⋅) :=12(∇u(ω,⋅)+∇u(ω,⋅)T),

and is a constant SPD matrix. denotes the frequency of the wave. denotes the imaginary unit. () is a bounded domain with boundary , denotes the outward normal to . For each , is a real-valued random variable defined over a probability space , where denotes the density of the random media which is the main source of randomness in the above PDEs. Thus, characterizes a random wave number for the elastic medium . We also note that the notation is often called the strain tensor and is denoted by in the literature.

In this paper we mostly focus on the case of weakly random media in the sense that the elastic medium is a small random perturbation of a homogeneous background medium, that is, . Here represents the magnitude of the random fluctuation and is some random field which has a compact support on and satisfies

 P{ω∈Ω;∥η(ω,⋅)∥L∞(D)≤1}=1.

At the end of the paper, we will also present an idea on how to extend the numerical method and algorithm of this paper to a more general media cases. We note that the boundary condition given in (2) is known as the first order absorbing boundary condition (ABC) and this boundary condition simulates an unbounded domain by absorbing plane waves that come into the boundary in a normal direction (cf. [9]). We also note that since is compactly supported on that . This choice was made to ensure that (2) was indeed a first order ABC for every .

Numerical approximations of random and stochastic partial differential equations (SPDEs) have gained a lot of interests in recent years because of ever increasing needs for modeling the uncertainties or noises that arise in industrial and engineering applications [1, 2, 4, 16, 19, 25]. Two main numerical methods for random SPDEs are the Monte Carlo (finite element) method and the stochastic Galerkin method. The Monte Carlo method obtains a set of independent identically distributed (i.i.d.) solutions by sampling the PDE coefficients, and calculates the mean of the solution via a statistical average over all the sampling in the probability space [4]. The stochastic Galerkin method, on the other hand, reduces the SPDE into a high dimensional deterministic equation by expanding the random coefficients in the equation using the Karhunen-Loève or Wiener Chaos expansions [1, 2, 3, 6, 8, 21, 25, 27, 26]. In general, these two methods become computationally expensive when a large number of degrees of freedom is involved in the spatial discretization, especially for Helmholtz-type equations. Indeed both methods become computationally prohibitive in the case that the frequency is large, because solving a deterministic Helmholtz-type problem with large frequency is equivalent to solving a large indefinite linear system of equations. Furthermore, it is well-known that standard iterative methods perform poorly when applied to linear systems arising from Helmholtz-type problems [10]. The Monte Carlo method requires solving the boundary value problem many times with different sampling coefficients, while the stochastic Galerkin method usually leads to a high dimensional deterministic equation that may be too expensive to solve.

Recently, we have developed a new efficient multi-modes Monte Carlo method for modeling acoustic wave propagation in weakly random media [11]. To solve the governing random Helmholtz equation, the solution is first represented by a sum of mode functions, where each mode satisfies a Helmholtz equation with deterministic coefficients and a random source. The expectation of each mode function is then computed using a Monte Carlo interior penalty discontinuous Galerkin (MCIP-DG) method. We take advantage that the deterministic Helmholtz operators for all the modes are identical, and employ an solver for obtaining the numerical solutions. Since the discretized equations for all the modes have the same constant coefficient matrix, by using the decomposition matrices repeatedly, the solutions for all samplings of mode functions are obtained in an efficient way by performing simple forward and backward substitutions. This leads to a tremendous saving in the computational costs. Indeed, as discussed in [11], the computational complexity of the proposed algorithm is comparable to that of solving a few deterministic Helmholtz problem using the a Gaussian elimination direct solver. Due to the similarities between the scalar and elastic Helmholtz operators, it is natural to extend the multi-modes MCIP-DG method of [11] to the elastic case for solving (1)–(2). On the other hand the scalar and elastic Helmholtz operators have different behaviors and kernel spaces so a separate study must be carried out to construct and analyze the multi-modes MCIP-DG method for the elastic Helmholtz problem. This is exactly the primary goal of this paper.

The rest of the paper is organized as follows. In Section 1 we present a complete PDE analysis of problem (1)–(2), including frequency-explicit solution estimates along with existence and uniqueness of solutions. In Section 2 the multi-modes expansion of the solution is defined and the convergence of the expansion is also demonstrated. Moreover, error estimates are derived for its finite term approximations. In Section 3 we formulate our MCIP-DG method and derive error estimates for the method. Section 4 lays out the overall multi-modes MCIP-DG algorithm. Computational complexity and convergence rate analysis are carried out for the algorithm. In Section 5 we present several numerical experiments to demonstrate the performance and key features of the proposed multi-modes MCIP-DG method and the overall algorithm. Finally, in Section 6 we describe an idea on how to extend the proposed MCIP-DG method and algorithm to the cases where more general (i.e., non-weak) random media must be considered.

## 1 PDE analysis

### 1.1 Preliminaries

Standard function and space notations are adopted in this paper. denotes the space of all complex vector-valued square-integrable functions on , and denotes the standard complex vector-valued Sobolev space. For any and we let and denote the standard complex-valued -inner products on and , respectively. We also define the special function spaces

 H1+(D) :={v∈H1(D);∇u∣∣∂D∈L2(∂D)}, V :={v∈H1+(D);\rm div\,(σ(u))∈L2(D)}.

Without loss of generality, we assume that the domain . Throughout this paper we also assume that is a convex polygonal or a smooth domain that satisfies a star-shape condition with respect to the origin, i.e. there exists a positive constant such that

 x⋅ν≥c0 for all x∈∂D.

will denote the space of vector-valued square integrable functions on the probability space . will denote the expectation operator given by

 (3) E(v):=∫ΩvdP % for all v∈L2(Ω),

and the abbreviation will stand for almost surely.

With these conventions in place, we introduce the following definition. {definition} Let . A function is called a weak solution to problem (1)–(2) if it satisfies the following identity:

 (4) ∫Ωa(u,v)dP=∫Ω(f,v)DdP∀v∈L2(Ω,H1(D)),

where

 (5) a(w,v) :=2μ(∇sw,∇sv)D+λ(\rm div\,w,\rm div\,v)D−k2(α2w,v)D +ik⟨Aw,v⟩∂D.

To simplify the analysis throughout the rest of this paper we introduce the following special semi-norm on :

 (6) |v|1,D:=λ∥\rm div\,v∥L2(D)+2μ∥∇sv∥L2(D) for all v∈H1(D).
###### Remark \thetheorem

a) Korn’s second inequality ensures that the semi-norm defined above is equivalent to the standard semi-norm.

b) By using Lemma 1.2 below, it is easy to show that any solution of (4)– (5) satisfies .

In [5], it was shown that estimates for solutions of the deterministic elastic Helmholtz problem are optimal in when the solution satisfies a Korn-type inequality on the boundary of the form

 ∥∇u∥2L2(∂D)≤~K[∥∇u∥2L2(D)+∥∇su∥2L2(∂D)+∥\rm div\,u∥2L2(∂D)],

where is a positive constant independent of . In the next subsection, a similar result is shown for solutions satisfying a stochastic Korn-type inequality on the boundary given by

 (7) E(∥∇u∥2L2(∂D)) ≤~K[E(∥∇u∥2L2(D))+E(∥∇su∥2L2(∂D))+E(∥\rm div\,% u∥2L2(∂D))].

The Korn-type inequality above is just a conjecture at this point. For this reason we introduce the special function space

 V~K :={v∈H1+(D);vsatisfies (???)},

for some independent of . We also introduce a special parameter , that will be used in the solution estimates presented in the next section. is defined as

 (8) ~α:={1 if u∈V~K2 otherwise ,

where is the solution to (1)–(2).

### 1.2 Frequency-explicit solution estimates

In this subsection, we derive stability estimates for the solution of problem (4)–(5). Since the main concern of this paper is the case when is large, we make the assumption that to simplify some of the estimates. The goal of this section is to derive solution estimates that are explicitly dependent on the frequency . These frequency-explicit estimates play a pivotal role in the development of numerical methods for deterministic wave equations (cf. [13], [14],[12]). We will obtain existence and uniqueness of solutions to (4)–(5) as a direct consequence of the estimates established in this subsection.

We begin with a number of technical lemmas which will be used in the proof of our solution estimates. Our analysis follows the analysis for the deterministic elastic Helmholtz equations carried out in [5] with many changes made to accommodate the inclusion of the random field .

For many of the estimates derived in this paper, it is important to note that the matrix in (2) is a real symmetric positive definite matrix. Thus, there exist positive constants and such that

 (9) cA∥u∥2L2(∂D)≤⟨Au,u⟩∂D≤CA∥u∥2L2(∂D),

for all .

{lemma}

Suppose solves (4)–(5). Then for any , satisfies the following estimates:

 (10) E(|u|21,D) ≤((1+ε)2k2+δ1)E(∥u∥2L2(D))+14δ1E(∥f∥2L2(D)), (11) E(∥u∥2L2(∂D)) ≤δ2cAkE(∥u∥2L2(D))+14δ2cAkE(∥f∥2L2(D)).
{proof}

Setting in (4) and taking the real and imaginary parts separately yields

 (12) Re∫Ω(f,u)DdP=∫Ω|u|21,D−k2∥∥(1+εη)u∥∥2L2(D)dP, (13) Im∫Ω(f,u)DdP=k∫Ω⟨Au,u⟩∂DdP.

(10) is obtained by rearranging the terms in (12) and applying the Cauchy-Schwarz inequality. Applying the Cauchy-Schwarz inequality along with (9) to (13) produces

 cAkE(∥u∥2L2(∂D))

We divide both sides of the above inequality by . This yields (11). The proof is complete.

Korn’s second inequality is essential to establishing solution estimates for the deterministic elastic Helmholtz equation (cf. [5]). Here we state a stochastic version of this inequality. {lemma}[Stochastic Korn’s Inequality] Let , then there exists a positive constant such that

 (14)
{proof}

There exists such that for fixed , satisfies Korn’s second inequality

 ∥∥∇su(ω,⋅)∥∥2L2(D)+∥∥u(ω,⋅)∥∥2L2(D)≥K∥∥u(ω,⋅)∥∥2H1(D).

For a proof of this inequality see [24]. Integrating over yields (14). The proof is complete.

Rellich identities for the elastic Helmholtz operator are also essential in the proof of solution estimates for the deterministic elastic Helmholtz equation (cf. [5]). Here we state stochastic versions of these Rellich identities. {lemma}[Stochastic Rellich Identities] Suppose that . Then for , we have the following stochastic Rellich identities:

 (15) (16) Re∫Ω(2μ(∇su,∇sv)D+λ(\rm div\,u,\rm div\,v)D)dP =2−d2∫Ω(2μ∥∇su∥2L2(D)+λ∥\rm div\,u∥2L2(D))dP +12∫Ω(2μ⟨x⋅ν,|\rm div\,u|2⟩∂Ω+λ⟨x⋅ν,|∇su|2⟩∂Ω)dP.
{proof}

For any we obtain the following identities from Proposition 2 and Lemma 5 of [5]:

 Re(u(ω,⋅),v(ω,⋅))D=−d2∥∥u(ω,⋅)∥∥2L2(D)+12⟨x⋅ν,|u(ω,⋅)|2⟩∂Ω, =2−d2(2μ∥∥∇su(ω,⋅)∥∥2L2(D)+λ∥∥\rm div\,% u(ω,⋅)∥∥2L2(D)) +12(2μ⟨x⋅ν,|\rm div\,u(ω,⋅)|2⟩∂Ω+λ⟨x⋅ν,|∇su(ω,⋅)|2⟩∂Ω).

(15) and (16) are obtained by integrating the above identities over . The proof is complete.

The following two lemmas relate higher order norms of the solution of (4)–(5) to the -norm of . {lemma} Suppose that solves (4)–(5). Then for all ,

 (17) 2μKE(∥u∥2H1(D))
{proof}

We obtain (17) by combining (10) and (14).

{lemma}

Suppose that solves (4)–(5) with and . Then there exists a positive constant , independent of and such that

 (18) E(∥u∥2H2(D))≤C(1+(1+ε)2k2)2E(∥u∥2L2(D))+CE(∥f∥2L2(D)), (19) E(∥∇u∥2L2(∂D))≤C(k3E(∥u∥2L2(D))+1kE(∥f∥2L2(D))).
{proof}

Regularity theory for elliptic problems [18, 20] implies for a.e.

 ∥u(ω,⋅)∥2H2(D) ≤C(∥f(ω,⋅)∥2L2(D)+((1+ε)2k2)2∥u(ω,⋅)∥2L2(D) +CAk2∥u(ω,⋅)∥2H12(∂D)+∥u(ω,⋅)∥2L2(D)) ≤C(∥f(ω,⋅)∥2L2(D)+(1+(1+ε)2k2)2∥u(ω,⋅)∥2L2(D)

Taking the expectation on both sides and using Lemma 1.2 yield

 E(∥u∥2H2(D)) +Ck2μK((k2(1+ε)2+δ+2μ)E(∥u∥2L2(D))+14δ(∥f∥2L2(D))).

Hence (18) holds with .

To prove (19), we note that is piecewise smooth. Thus, by the trace inequality, (17), and (18) the following inequalities hold

 E(∥∇u∥2L2(∂D)) ≤C∫Ω∥∇u∥L2(D)∥u∥H2(D)dP ≤CkE(∥∇u∥2L2(D))+CkE(∥u∥2H2(D)) ≤CkμK(((1+ε)2k2+δ+2μ)E(∥u∥2L2(D))+14δE(∥f∥2L2(D))) +Ck((1+(1+ε)2k2)2E(∥u∥2L2(D))+E(∥f∥2L2(D))).

Letting in the above inequality yields

 E(∥∇u∥2L2(∂D)) ≤C(k(1+(1+ε)2k2)+1k(1+(1+ε)2k2)2)E(∥u∥2L2(D)) +C(1(1+ε)2k+1k)E(∥f∥2L2(D)).

By the assumptions and we obtain (19). The proof is complete.

With these technical lemmas in place we are now ready to prove the main result for this section.

{theorem}

Suppose that solves (4)–(5) with . Further, let be a convex polygonal or a smooth domain and be the smallest number such that . Then the following estimates hold

 (20) E(∥u∥2L2(D)+∥u∥2L2(∂D)+c0k2|u|21,∂D) ≤C0(k~α−2+1k2)2E(∥f∥2L2(D)), (21) E(∥u∥2H1(D))≤C0(k~α−1+1k2)2E(∥f∥2L2(D)),

provided that . Here is a positive constant independent of and , and is defined by (8). Moreover, if the following estimate also holds

 (22) E(∥u∥2H2(D))≤C0(k~α+1k2)2E(∥f∥2L2(D)).
{proof}

Step 1: We begin by proving (20). Let in (4). By taking the real part and rearranging terms we find

To this identity we apply the Rellich identities in (15) and (16), after another rearrangement of terms we get

 dk22E(∥u∥2L2(D)) ≤12∫Ω(k2⟨x⋅ν,|u|2⟩∂Ω−2μ⟨x⋅ν,|∇su|2⟩∂Ω−λ⟨x⋅ν,|\rm div\,u|2⟩∂Ω)dP +d−22E(|u|21,D)+Re∫Ω(k2ε(η(2+εη)u,v)D+(f,v)D)dP +Im∫Ωk⟨Au,v⟩∂DdP.

Using the fact that is star-shaped and along with the Cauchy-Schwarz inequality we obtain

 (23) dk22E(∥u∥2L2(D)) ≤k2R2E(∥u∥2L2(∂D))−c02E(|u|21,∂D)+d−22E(|u|21,D) +R2δ2E(∥f∥2L2(D))+Rδ22E(∥∇u∥2L2(D))

From (10) we find

 (24) E(|u|21,D)−k2E(∥u∥2L2(D))≤(k2ε(2+ε)+δ4)∥u∥2L2(D)+14δ4∥f∥2L2(D).

By adding times (24) to (23), grouping like terms and letting , we get

 (25) ≤k2γR2δ1E(∥u∥2L2(D))+(k2γRδ12+Rδ22)E(∥∇u∥2L2(D)) +(k2R2+kCAR2δ3)E(∥u∥2L2(∂D))+kCAδ32E(∥∇u∥2L2(∂D)) +(k2γ+δ4)∥u∥2L2(D)+14δ4∥f∥2L2(D)+R2δ2E(∥f∥2L2(D)).

Step 2: The source of the different values of in (20) comes from different treatments for the term. In particular, if we apply (7) to control this term. Otherwise, we apply (19) to control this term. We first prove (20) with . Applying (11), (17) and (19) to (25) yields.

 ≤1μK(k2γRδ12+Rδ22)((k2(1+γ)+δ5+2μ)E(∥u∥2L2(D))+12δ5E(∥f∥2L2(D))) +(k2R2+kCAR2δ3)(δ6cAkE(∥u∥2L2(D))+14δ6cAkE(∥f∥2L2(D)))) +CkCARδ32(k3E(∥u∥2L2(D))+1kE(∥f∥2L2(D))) +d−12((γk2+δ4)E(∥u∥2L2(D))+14δ4E(∥f∥2L2(D))) +k2γR2δ1E(∥u∥2L2(D))+R2δ2E(∥f∥2L2(D)).

In order to apply (19) we require . Since , it can be easily shown that . Thus,

 (26) c1E(∥u∥2L2(D))+c02E(|u|21,∂D)+12E(|u|21,D)≤c2E(∥f∥2L2(D)),

where

 c1 −(k2R2+kCAR2δ3)δ6cAk−CCAk4Rδ32−k2γR2δ1,

and

 c2 :=12μKδ5(k2γRδ12+Rδ22)+14δ6cAk(k2R2+kCAR2δ3) +CkCARδ32k+d−18δ4+R2δ2.

In the third term of , we have used the fact that to include a coefficient of to the term. This has been done to simplify constants later. Setting

 δ1=12k,δ2=μK16R(3+2μ),δ3=18CCAk2R,δ4=116(d−1),δ5=k22,δ6=cAk28(kR+8CC2Ak2R2),

and using the fact that yields

 c1=k22−k22(d−1+3kRμK+2kRK+kR)γ−k24.

Also using the fact that yields . It is easy to check that . Therefore, (26) becomes

Multiplying both sides by and applying (11) with implies (20) with .

Step 3: If , we apply (14) and obtain.

 kCA2δ3E(∥∇u∥2L2(∂D)−14(k2E(∥u∥2L2(D))+c0E(|u|21,∂D)+E(|u|21,D)) ≤kCA2δ3E(∥∇u∥2L2(∂D)−14min{k2K,2μK,c0λ,2c0μ} ⋅(E(∥∇u∥2L2(D))+E(∥∇su∥2L2(∂D))+E(∥\rm div\,% u∥2L2(∂D))).

By choosing , using the fact that , and applying (7) we find

 kCA2δ3E(∥∇u∥2L2(∂D)−14(k2E(∥u∥2L2(D))+c0E(|u|21,∂D)+E(|u|21,D))≤0.

We apply this inequality to (25) as well as (11) and (17) in order to find

 ≤1μK(k2γRδ12+Rδ22)((k2(1+γ)+δ5+2μ)E(∥u∥2L2(D))+12δ5E(∥f∥2L2(D))) +(k2R2+kCAR2δ3)(δ6cAkE(∥u∥2L2(D))+14δ6cAkE(∥f∥2L2(D)))) +d−12((γk2+δ4)E(∥u∥2L2(D))+14δ4E(∥f∥2L2(D))) +k2γR2δ1E(∥u∥2L2(D))+R2δ2E(∥f∥2L2(D)).

Thus,

 (27) c1E(∥u∥2L2(D))+c04E(|u|21,∂D)+14E(|u|21,D)≤c2E(∥f∥2L2(D)),

where

 c1 −(k2R2+kCAR2δ3)δ6cAk−k2γR2δ1, c2 :=12μKδ5(k2γRδ12+Rδ22)+14δ6cAk(k2R2+kCAR2δ3) +d−18δ4+R2δ2.

Setting

 δ1=1k,δ2=μK16R(3+2