A shape calculus based method for a transmission problem with a random interface

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A shape calculus based method for a transmission problem with a random interface✩ Alexey Chernov a,∗,1 , Duong Pham a,2 , Thanh Tran b a

Hausdorff Center for Mathematics and Institute for Numerical Simulation, University of Bonn, Endenicher Allee 64, Bonn 53115, Germany b

School of Mathematics and Statistics, The University of New South Wales, Sydney 2052, Australia

article

info

Article history: Available online xxxx Dedicated to E. Rank on the occasion of his 60th anniversary Keywords: Shape derivative Transmission problem Random domain Boundary integral equations Spectral discretization Moment equations

abstract The present work is devoted to an approximation of the statistical moments of the solution of a class of elliptic transmission problems in R3 with uncertainly located transmission interfaces. In this model, the diffusion coefficient has a jump discontinuity across the random transmission interface which models linear diffusion in two different media separated by an uncertain surface. We apply the shape calculus approach to approximate the solution perturbation by the so-called shape derivative. Correspondingly, statistical moments of the solution are approximated by the moments of the shape derivative. We characterize the shape derivative as a solution of a related homogeneous transmission problem with nonzero jump conditions, which is solved by the boundary integral equation method. A rigorous theoretical framework is developed, and the theoretical findings are supported by and illustrated in two particular examples. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Elliptic transmission or interface problems arise in many fields in science and engineering, such as tomography, deformation of an elastic body with inclusions, stationary groundwater flow in heterogeneous medium, fluid–structure interaction, scattering of an elastic body and many others. Combined with the state-of-the-art hardware, advanced numerical schemes are capable of producing a highly accurate and efficient deterministic numerical simulation, provided that the problem data are known exactly. However, in real applications, a complete knowledge of the problem parameters is not realistic for many reasons. First, the simulation parameters are often estimated from measurements which can be inexact, e.g. due to imperfect measurement devices. Second, the parameters are estimated based on a large but finite number of system samples (snapshots); this information can be incomplete or stochastic. Finally, parameters of the system originate from a mathematical model which is itself only an approximation of the actual process. Under such circumstances, highly accurate results of a single deterministic simulation for one particular set of problem parameters are of limited use. An important paradigm, becoming rapidly popular

✩ This work has been funded by BMBF and the Group of Eight Australia within the DAAD-Go8 Project ‘‘Numerical methods for elliptic transmission problems on uncertain interfaces’’, Project ID 56266715 and RG123838. ∗ Corresponding author. E-mail addresses: [email protected] (A. Chernov), [email protected] (D. Pham), [email protected] (T. Tran). 1 Current address: Institute for Mathematics, Carl von Ossietzky University, 26111 Oldenburg, Germany.

2 Current address: Vietnamese German University, Le Lai street, Binh Duong New City, Binh Duong Province, Vietnam. http://dx.doi.org/10.1016/j.camwa.2015.06.021 0898-1221/© 2015 Elsevier Ltd. All rights reserved.

2

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

over the past years, see e.g. [1–11] and the references therein, is to treat the lack of knowledge via modelling uncertain parameters as random fields. If the forward solution operator is continuous, the solution of the forward problem with random parameters becomes a well-defined random field (since a composition of a continuous and a measurable function is measurable). Efficient numerical approximation of the random (or stochastic) solution and its probabilistic characteristics, e.g. statistical moments, is a highly non-trivial task representing numerous new interdisciplinary challenges: from regularity analysis and numerical analysis to modelling and efficient parallel large scale computing. As already mentioned in the beginning, the case of uncertain but sufficiently regular interfaces appears in numerous applications, i.e. in microbiology and medicine. Interior domains may represent the interior of biological cells (or, say, a volume occupied by a human organ, tumour, etc.) and exterior domains the outer medium (the surrounding part of the human body). Cells of the same type (as well as human organs) have a similar shape which, however, is somewhat different between two particular species, so that the family of all possible shapes can be viewed as a perturbation of a nominal one. The interface between the interior and the exterior domain can be reasonably understood as the location of the jump of the material parameters, e.g. as the jump of the diffusivity. In this article we develop a deterministic method for numerical solution for a class of transmission problems with randomly perturbed interfaces. The equation to be solved is of the form

−∇ · (α∇ u) = f

in D± ,

where D− is a random bounded domain in R3 and D+ = R3 \D− is its complement. The domains share a common random surface Γ , and the coefficient function α takes distinct constant values in D− and D+ , respectively. The solution u is subject to jump conditions across Γ . A precise description of the model problem is deferred until Section 2.3, where a probabilistic perturbation model for the surface Γ (and thus D± ) will be rigorously introduced. Within this model, the transmission interface depends on the ‘‘random event’’ ω and the parameter ϵ ≥ 0 controlling the amplitude of the perturbation. Therefore, the solution u depends on ω and ϵ , and will be denoted by uϵ (ω). The case ϵ = 0 corresponds to the zero perturbation. In the present paper we aim at estimating probabilistic properties of the solution perturbation uϵ (ω) − u0 when the perturbation parameter is small, namely ϵ ≪ 1. More precisely, we exploit the ideas from the recent publications [3,9,12–14] and propose to approximate the statistical moments of the solution perturbation by the moments of the linearized solution, i.e. for a fixed (small) value of the perturbation parameter ϵ , the kth order statistical moments of the solution perturbation are approximated by

Mk [uϵ − u0 ] ≈ ϵ k M k [u′ ]

(1.1)

and similarly

Mk [uϵ − E[uϵ ]] ≈ ϵ k Mk [u′ ]. ′

(1.2)

ϵ

Here u is the shape derivative of u formally understood as the linear order term in the asymptotic expansion uϵ (x, ω) = u0 (x) + ϵ u′ (x, ω) + · · · ,

ϵ → 0,

(1.3)

for almost all random events ω ∈ Ω at a certain fixed point x in the Euclidean space R . The notion of the shape derivative has been introduced in the context of the shape optimization (see e.g. the monograph [15] and the references therein) and allows to quantify sensitivity of the solution of a PDE to small perturbations of the boundary. It is worth mentioning that similar concepts have been developed (in the deterministic framework) and termed domain derivatives within the inverse problem community, see e.g. [16,17] for treatment of related elliptic and parabolic transmission problems in the case of the bounded outer subdomain D+ . Although very intuitive, (1.3) cannot be used as a rigorous definition of u′ (x, ω). In particular, the existence of the shape derivative and the convergence of the asymptotic expansion (1.3) are unclear without further assumptions. In the first part of this article (Section 3) we develop a rigorous mathematical theory of existence of the shape derivative for the class of elliptic transmission problems under consideration. Similarly to [14, Lemma 1], we obtain a characterization of the shape derivative u′ (x, ω) as a solution of a deterministic transmission problem on a fixed interface. Our contribution in this section is two-fold: (i) we extend the notion of shape derivatives for interface problems developed in [14] to the case of unbounded domains and higher order moments, and (ii) we fill the gaps in the existing literature where only a limited rigorous discussion on existence of shape derivatives is presented. We point out the rigorous results in the deterministic setting in [16,17] and also in [18,19] where a similar concept of the Fréchet domain derivative has been presented and rigorously analysed. This derivative in fact coincides with the notion of the material derivative rather than the shape derivative in the terminology of [15] and Definition 3.5 below. We also refer to the related work [20] and references therein. As mentioned above, for almost all ω ∈ Ω the shape derivative u′ (·, ω) is a solution of a deterministic problem in R3 with (in general) nonhomogeneous jump conditions but with vanishing volume source term. The second contribution of this article is the derivation and analysis of boundary integral equation methods [21–23] which are used to solve this transmission problem on deterministic domains with deterministic interfaces. A tensorization argument is then used to obtain the approximation (1.1) for the statistical moments. 3

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

3

Finally, we illustrate the accuracy of the linearization approach by considering several examples setting on the unit sphere

Γ := {|x| = 1} with radial perturbation. The first two examples involve a pre-determined solution with radial symmetry, so that the exact and the linearized solutions as well as their second moments are available explicitly. We observe that in this particular case the linearization error for the second order statistical moments is of the order O (ϵ 3 ) rather than o(ϵ 2 ) as confirmed by the theory. The slightly higher order can be explained by the smooth dependence of the solution on the perturbation parameter. The third example involves non-symmetric data so that the linearized solution is not available explicitly. To solve this problem numerically we use the sparse spectral tensor product BEM developed in [24]. This method exploits the underlying geometry of the formulation and uses the basis of spherical harmonics being the eigenfunctions of the integral operator governing the problem. In summary, the main results of this paper are:

• Theorem 3.12 which shows the existence of the shape derivative (of the solution uϵ (ω) to the given transmission problem) as the solution of a transmission problem on a fixed domain;

• Theorem 3.13 which shows that the k-moment of the solution to the given transmission problem can be approximated by that of the shape derivative;

• Theorem 4.3 which shows that the k-moment of the shape derivative can be computed as the unique solution of a boundary integral equation. Hence we propose a whole process to compute the k-moment of the solution to the given transmission problem. The paper is organized as follows. Section 2 contains the description of the random surface perturbation model and the rigorous formulation of the model transmission problem, preceded by the details on the function spaces involved in the analysis. Section 3 contains a generalization of shape calculus to the case of unbounded domains, the definition and characterization of the material and shape derivatives for the underlying model transmission problem, and a rigorous proof and error bounds for the approximation (1.1). Section 4 contains the details of the boundary reduction for the linearized problem. Section 5 contains two examples, one analytical and one numerical, illustrating the accuracy of the method. 2. Model elliptic transmission problem on a random interface We start with some preliminary definitions and notations in Section 2.1. Section 2.2 contains the description of a model for the random surface perturbation. We introduce the randomized model problem in the strong form in Section 2.3. The Sobolev spaces to be used are defined in Section 2.4. 2.1. Bochner spaces and statistical moments Throughout this paper we denote by (Ω , Σ , P) a generic complete probability space and let X be a Banach space. For any 1 ≤ k ≤ ∞, the Bochner space Lk (Ω , X ) is defined as usual by Lk (Ω , X ) := v : Ω → X , measurable : ∥v∥Lk (Ω ,X ) < ∞





(2.1)

with the norm

∥v∥Lk (Ω ,X ) :=

 1/k   ∥v(ω)∥kX dP(ω) , Ω

 

esssup ∥v(ω)∥X , ω∈Ω

1 ≤ k < ∞,

(2.2)

k = ∞.

The elements of L (Ω , X ) are called random fields. We remark that for a part of the subsequent analysis we may restrict to the special case when X is a Hilbert space. In particular, when X1 and X2 are two separable Hilbert spaces, their tensor product X1 ⊗ X2 is a separable Hilbert space with the natural inner product extended by linearity from ⟨v ⊗ a, w ⊗ b⟩X1 ⊗X2 = ⟨v, w⟩X1 ⟨a, b⟩X2 , cf. e.g. [25, p. 20], [26, Definition 12.3.2, p. 298]. In this paper we work with k-fold tensor products of Hilbert spaces k

X (k) := X ⊗ · · · ⊗ X

(2.3)

with the natural inner product satisfying ⟨v1 ⊗ · · · ⊗ vk , w1 ⊗ · · · ⊗ wk ⟩X (k) = ⟨v1 , w1 ⟩X . . . ⟨vk , wk ⟩X . Definition 2.1. For a random field v ∈ Lk (Ω , X ), its k-order moment M k [v] is an element of X (k) defined by k





M [v] := Ω

 v(ω) ⊗ · · · ⊗ v(ω) dP(ω).    k-times

(2.4)

4

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

In the case k = 1, the statistical moment M 1 [v] coincides with the mean value of v and is denoted by E[v]. If k ≥ 2, the statistical moment M k [v] is the k-point autocorrelation function of v . The quantity M k [v − E[v]] is termed the kth central moment of v . We distinguish in particular second order moments: the correlation and covariance defined by Cor[v] := M 2 [v]

and

Cov[v] := M 2 [v − E[v]].

(2.5)

In this paper we work with X being Sobolev spaces of real-valued functions defined on a domain U ⊆ R3 yielding, in particular, the representation Cor[v](x, y ) :=

 Ω

v(x, ω)v(y , ω) dP(ω),

x, y ∈ U .

(2.6)

We observe that Cor[v] is defined on the Cartesian product U × U. Similarly, M k [v] is defined on the k-fold Cartesian product U × · · · × U. Here, the dimension of the underlying domain grows rapidly with increasing moment order k. 2.2. Representation of random interfaces Consider a fixed bounded domain D0− ⊂ R3 and let D0+ := R3 \D0− be its complement. Then the interface Γ 0 = D0− ∩ D0+ is a closed manifold in R3 . For the subsequent analysis we assume that Γ 0 is at least of the class C 1,1 . This implies that the

outward normal vector n0 to Γ 0 is Lipschitz continuous: n0 ∈ C 0,1 (Γ 0 ). The partition R3 = D0+ ∪ D0− and the interface Γ 0 will be fixed throughout the paper and will be called the nominal partition and nominal interface, respectively. In the present paper we utilize the domain perturbation model based on the speed method (see e.g. the monograph [15] and the references therein) and random domain perturbation model from [3,9,12–14]. Suppose κ ∈ Lk (Ω , C 0,1 (Γ 0 )) is a random field, i.e. for almost any realization ω ∈ Ω , we have κ(·, ω) ∈ C 0,1 (Γ 0 ). For some sufficiently small, nonnegative ϵ we consider a family of random interfaces of the form

Γ ϵ (ω) = {x + ϵκ(x, ω)n0 (x) : x ∈ Γ 0 },

ω ∈ Ω.

ϵ

(2.7) ϵ

Here, the uncertainty of the surfaces Γ (ω) is represented by the uncertainty in κ(·, ω). Notice that the interface Γ (ω)|ϵ=0 is identical with Γ 0 and therefore is a deterministic closed manifold. Moreover, the limit Γ ϵ (ω) → Γ 0 as ϵ → 0 is well defined in Lk (Ω , C 0,1 ). If we identify Γ ϵ and Γ 0 with the functions defining their graphs, then ϵ

∥Γ − Γ ∥Lk (Ω ,C 0,1 ) = ϵ 0

 Ω

∥κ(·, ω) ∥

n0 kC 0,1 (Γ 0 )

 1k ≤ 2ϵ∥κ∥Lk (Ω ,C 0,1 (Γ 0 )) ∥n0 ∥C 0,1 (Γ 0 ) . dP(ω)

(2.8)

This implies that for almost all ω ∈ Ω and a sufficiently small ϵ ≥ 0 the surface Γ ϵ (ω) is a Lipschitz continuous closed manifold separating the interior domain Dϵ− (ω) and its complement Dϵ+ (ω) := R3 \Dϵ− . The shape calculus in Section 3 requires a somewhat stronger smoothness assumption on κ , namely that the realizations of κ belong to C 1 (Γ 0 ). From (2.7) we observe that the mean random interface is represented by

E[Γ ϵ ] = x + ϵ E[κ(x, ·)]n0 (x), x ∈ Γ 0 .





Without loss of generality, we may assume that the random perturbation amplitude κ(x, ω) is centred, i.e.,

E[κ(x, ·)] = 0 ∀x ∈ Γ 0 .

(2.9)

In this case

E[Γ ϵ ] = Γ 0

and

Cov[κ](x, y ) = Cor[κ](x, y ).

Concerning stochastic regularity of κ , we impose later on the condition that κ ∈ Lk (Ω , C 1 (Γ 0 )) for some integer k, which is sufficient for existence of the kth order moment (cf. Theorem 3.13). No additional restrictions on κ (a particular covariance structure, etc.) will be imposed. 2.3. The model problem As shown above, for a sufficiently small value ϵ ≥ 0 the surface perturbation model (2.7) generates a well defined partition of R3 into a bounded Lipschitz domain Dϵ− (ω) and its complement Dϵ+ (ω) = R3 \Dϵ− separated by the closed Lipschitz manifold Γ ϵ (ω) = Dϵ− (ω)∩ Dϵ+ (ω). We consider a piecewise constant diffusion function subjected to this partition:

α ϵ (x, ω) =



α− , α+ ,

x ∈ Dϵ− (ω), x ∈ Dϵ+ (ω),

(2.10)

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

5

where α− and α+ are two positive constants independent of x, ϵ , and ω. Having this we introduce the model elliptic transmission problem as a problem of finding uϵ satisfying

  −∇ · α ϵ (x, ω)∇ uϵ (x, ω) = f (x) in Dϵ± (ω), ϵ

(2.11a)

ϵ

[u (x, ω)] = 0 on Γ (ω),   ∂ uϵ α ϵ (x, ω) ϵ (x, ω) = 0 on Γ ϵ (ω), ∂n

(2.11b) (2.11c)

uϵ (x, ω) = O(|x|−1 ) as |x| → +∞. ϵ

(2.11d) ϵ

ϵ

ϵ

ϵ

Here, ∂/∂ n denotes the normal derivative on Γ (ω), i.e. ∂/∂ n = n (x, ω) · ∇ , where n (x, ω) is the unit normal vector to the interface Γ ϵ (ω) pointing into the interior of Dϵ+ (ω). Let uϵ− (ω) and uϵ+ (ω) be the restrictions of uϵ (ω) on Dϵ− (ω) and Dϵ+ (ω), respectively. Then the jump [uϵ (ω)] is understood to be uϵ− (ω) − uϵ+ (ω) on Γ ϵ (ω) in the sense of trace for each sample ω. Similarly



 ϵ ϵ ∂ uϵ ϵ ∂ u− ϵ ∂ u+ α (x, ω) ϵ (x, ω) = α− ( x , ω) − α (x, ω), + ∂n ∂ nϵ ∂ nϵ ϵ

x ∈ Γ ϵ (ω).

The function f ∈ H 1 (R3 ) is assumed to be independent of ω and thereby it represents a deterministic source function in R3 . We remark that the case of a sufficiently smooth but nonconstant diffusion α± = α± (x, ω) can be treated with similar arguments. Here we consider the aforementioned setting to keep the presentation reasonably simple. The model problem (2.11) represents a stationary diffusion in R3 with piecewise constant diffusivity in the interior and exterior domain. The uncertainty in the random solution uϵ (x, ω) is implied by the uncertain location of the transmission interface Γ ϵ (ω). The solution depends nonlinearly on the interface and a linearization process will first be used to linearize the initial problem. The tool in this process is shape calculus which will be presented in Section 3. In what follows we address the problem of approximation of the statistical moments

E[uϵ ],

Mk [uϵ − u0 ],

and M k [uϵ − E[uϵ ]],

k ≥ 2,

(2.12)

with this strategy as well as the rigorous control of the approximation error. 2.4. Sobolev spaces In this section we introduce the function spaces needed for the forthcoming analysis. Let U be a bounded domain in R3 . The Sobolev space H 1 (U ) is defined, as usual, as the space of all distributions which together with their first order partial derivatives are square integrable. The Sobolev space H 1/2 (∂ U ) is defined by H 1/2 (∂ U ) = {g : ∂ U → R | g = v on ∂ U (in the trace sense) for some v ∈ H 1 (U )} and equipped with the following norm

∥g ∥H 1/2 (∂ U ) := inf{∥v∥H 1 (U ) : v ∈ H 1 (U )and g = v|∂ U }. The dual of H 1/2 (∂ U ) is denoted by H −1/2 (∂ U ). 1/2 The tensor product of Sobolev spaces on the k-fold Cartesian product ∂ U k = ∂ U × · · · × ∂ U is denoted by Hmix (∂ U k ). s k These spaces will be used later on for characterization of statistical moments. We also use the notation Hmix (K ) for the tensor product H s (K )(k) where K is a compact subset of R3 . Proper treatment of the transmission problem (2.11) in unbounded domains in R3 requires a special care. Following [22], for an unbounded domain U ⊂ R3 we introduce the space 1 Hw (U ) :=



v ∈ D ′ (U ) : ∥v∥Hw1 (U ) =

  1/2  |v(x)|2  |∇v|2 + dx < +∞ . 1 + |x|2 U

(2.13)

Specifically, for a given partition R3 = Dϵ− ∪ Dϵ+ we define the space 1 Wϵ := v = (v− , v+ ) ∈ H 1 (Dϵ− ) × Hw (Dϵ+ ) : [v ]Γ ϵ = 0





(2.14)

which is a weighted Sobolev space on Dϵ− ∪ Dϵ+ with corresponding norm and seminorm

∥v∥Wϵ

 1/2 := ∥v− ∥2H 1 (Dϵ ) + ∥v+ ∥2H 1 (Dϵ ) , −

w

+

 |v|Wϵ :=

Dϵ−

2

|∇v− | dx +

1/2



2

Dϵ+

|∇v+ | dx

.

(2.15)

The following lemma which will be frequently used in the rest of the paper states the equivalence between the norm

∥·∥Wϵ and seminorm |·|Wϵ . The proof of this result follows by the Friedrichs inequality and the technique in the proof of [22, Theorem 2.10.10].

Lemma 2.2. The seminorm |·|Wϵ is also a norm in Wϵ which is equivalent to ∥·∥Wϵ .

6

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

3. Shape calculus The aim of the present section is the systematic development of the linearization theory for the solution uϵ of the model problem (2.11) with respect to the shape of the perturbed interface Γ ϵ . This techniques is also known as shape calculus and originates from shape optimization; see [15] and the references therein. For this purpose, in the first three subsections that follow, we temporarily stay away from randomness and consider only deterministic perturbed interfaces. 3.1. Perturbation of deterministic interfaces In this subsection we collect several properties of perturbed interfaces which are important for the subsequent analysis. Assume that the perturbation function κ is a fixed deterministic function in W 1,∞ (Γ 0 ), in particular κ is independent of ω. Then Γ ϵ is defined by

Γ ϵ := {x + ϵκ(x)n0 (x) : x ∈ Γ 0 },

ϵ > 0.

(3.1)

As already noticed in Section 2.2, Γ ϵ is a closed Lipschitz manifold in R3 provided 0 ≤ ϵ ≤ ϵ0 and ϵ0 is sufficiently small. In this case Γ ϵ introduces a decomposition of R3 into the interior and exterior subdomains Dϵ− and Dϵ+ , respectively. Following [15], we define a mapping T ϵ : R3 → R3 which transforms Γ 0 into Γ ϵ and D0± into Dϵ± , respectively, by T ϵ (x) := x + ϵ κ( ˜ x)n˜ 0 (x),

x ∈ R3 ,

(3.2) 1,∞

˜ are any smoothness-preserving extensions of κ and n into R . We require in particular that κ˜ ∈ W (R3 ). where κ˜ and n Without loss of generality we assume that the extension κ˜ vanishes outside a sufficiently large ball BR := {x ∈ R3 : |x| < R} containing Γ ϵ for any 0 ≤ ϵ ≤ ϵ0 . This implies that the perturbation mapping T ϵ (x) is an identity in the complement BcR := R3 \BR , i.e. 0

0

3

T ϵ (x) = x ∀x ∈ BcR .

(3.3)

For the ease of notation we abbreviate V (x) := κ( ˜ x)n˜ 0 (x),

x ∈ R3 .

(3.4) ϵ

In [15], V is called the velocity field of the mapping T . The following result is straightforward.

3

Lemma 3.1. Assuming κ˜ ∈ W 1,∞ (R3 ) and κ( ˜ x) = 0 for x ∈ BcR , there hold V ∈ H 1 (R3 )



and

∂ V (x) = 0 ∀x ∈ BcR , l = 1, 2, 3, m = 0, 1. ∂ xm l m

Recall the definition (2.14) of the weighted space Wϵ associated to the splitting R3 = Dϵ− ∪ Dϵ+ . It can be proved that a function v belongs to Wϵ if and only if the composition v ◦ T ϵ belongs to W0 , and there hold

∥(v ϵ )− ∥H 1 (Dϵ− ) ≃ ∥(v ϵ ◦ T ϵ )− ∥H 1 (D0 ) −

∥(v ϵ )+ ∥Hw1 (Dϵ ) ≃ ∥(v ϵ ◦ T ϵ )+ ∥Hw1 (D0 ) +

(3.5)

+

ϵ

ϵ

ϵ

∥v ∥Wϵ ≃ ∥v ◦ T ∥W0 . In the subsequent analysis, for any 3 by 3 matrix A(x) whose entries are functionals of x ∈ U ⊂ R3 , we denote

∥A(·)∥Lp (U ) := max

i,j=1,2,3

  Ai,j (·) p



L (U )

,

1 ≤ p ≤ ∞,

where Aij are components of A. The following three lemmas state some important properties of the mapping T ϵ which will be used later in this section. Until the end of this section we assume that T ϵ is defined by (3.2) and (3.3) with κ˜ ∈ C 1 (R3 ), and denote its Jacobian matrix and Jacobian determinant by JT ϵ and γ (ϵ, ·), respectively. ⊤ Lemma 3.2. Consider A(ϵ, ·) := γ (ϵ, ·)JT−ϵ1 JT−⊤ ϵ , where JT ϵ is the transpose of JT ϵ . Then there hold

lim ∥A(ϵ, ·) − I ∥L∞ (R3 ) = 0

(3.6)

ϵ→0

and

  A(ϵ, ·) − I

lim  

ϵ→0

ϵ

  − A′ (0, ·) 

L∞ (R3 )

= 0.

Here, A′ (0, ·) is the Gâteaux derivative of A(ϵ, ·) at ϵ = 0, namely A′ (0, x) = lim

ϵ→0

A(ϵ, x) − I (x)

ϵ

,

x ∈ R3 .

(3.7)

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

7

Proof. Denoting V (x) := (V1 (x), V2 (x), V3 (x))⊤ , the Jacobian matrix and the Jacobian of T ϵ are given by

∂ V1 (x) ∂ x2 ∂ V2 (x) 1+ϵ ∂ x2 ∂ V3 (x) ϵ ∂ x2

∂ V1 (x) ∂ x1 ∂ V2 (x) ϵ ∂ x1 ∂ V3 (x) ϵ ∂ x1

    JT ϵ (x) =    

 ∂ V1 (x)  ∂ x3  ∂ V2 (x)   ϵ  ∂ x3  ∂ V3 (x)  1+ϵ ∂ x3

ϵ

1+ϵ

ϵ

(3.8)

and 3 3    ∂ Vk (x) ∂ Vl (x) ∂ Vl (x) ∂ Vk (x)  ∂ Vk (x)   + ϵ2 − γ (ϵ, x) = 1 + ϵ ∂ xk ∂ xk ∂ xl ∂ xk ∂ xl k,l=1 k=1 k̸=l

∂ Vi (x) ∂ Vj (x) ∂ Vk (x)   ∂ x1 ∂ x2 ∂ x3 i,j,k=1   =: 1 + ϵγ1 (x) + ϵ 2 γ2 (x) + ϵ 3 γ3 (x). + ϵ3

3 

sign(i, j, k)

(3.9)

Here sign(i, j, k) denotes the sign of the permutation (i, j, k). The entries Aij (ϵ, x), i, j = 1, 2, 3, of the matrix A(ϵ, x) are given by

 Aij (ϵ, x) = γ (ϵ, x)−1

δij +

4 

 ϵ n hijn (x) ,

(3.10)

n =1

where hijn is a polynomial of partial derivatives of V and δij is the Kronecker delta. Using Lemma 3.1, we deduce

γn , hijn ∈ L∞ (R3 ) ∩ L2 (R3 ), lim ∥γ (ϵ, ·)∥L∞ (R3 ) > 0,

i, j = 1, 2, 3 and n = 1, . . . , 4,

(3.11)

ϵ→0

where γ1 , γ2 , γ3 are defined by (3.9) and γ4 := 0 for notational convenience later. In particular, for sufficiently small ϵ > 0, there holds

γ (ϵ, x) = 1 + ϵγ1 (x) + ϵ 2 γ2 (x) + ϵ 3 γ3 (x) ≥ c > 0 ∀x ∈ R3 .

(3.12)

Consider from now on sufficiently small ϵ > 0. It follows from (3.10) and (3.12) that the ij-entry of the matrix A(ϵ, x) − I is Aij (ϵ, ·) − δij = ϵ γ (ϵ, ·)−1

4 

  ϵ n−1 hijn − δij γn .

(3.13)

n=1

Hence, (3.11) yields

  Aij (ϵ, ·) − δij  ∞ L

(R3 )

→ 0 as ϵ → 0,

proving (3.6). From (3.13), we have Aij (ϵ, ·) − δij

ϵ

= γ (ϵ, ·)−1

4 

  ϵ n−1 hijn − δij γn .

(3.14)

n=1

Taking the limit when ϵ goes to 0, noting that γ (ϵ, ·) → 1, we obtain A′ij (0, ·) = hij1 − δij γ1 ,

i, j = 1 , 2 , 3 .

(3.15)

Subtracting (3.15) from (3.14) side by side, we obtain Aij (ϵ, ·) − δij

ϵ

− A′ij (0, ·) = γ (ϵ, ·)−1

4  n =2

 ϵ n−1 (hijn − δij γn ) − (hij1 − δij γ1 )(γ (ϵ, ·) − 1) .

(3.16)

8

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

Noting (3.11), we infer

  Aij (ϵ, ·) − δij

lim  

ϵ→0

  − A′ij (0, ·) 

ϵ

= 0,

L∞ (R3 )

proving (3.7), completing the proof of the lemma.



Lemma 3.3. For any function v ∈ L2 (R3 ), there holds

 

 



ϵ 2  lim   1 + |·| γ (ϵ, ·) v ◦ T − v 

ϵ→0

L2 (R3 )

= 0.

Proof. Since T ϵ (x) = x for any x ∈ BcR , see (3.3), the Jacobian satisfies

γ (ϵ, x) = 1 for any x ∈ BcR .

(3.17)

Therefore,

       1 + |·|2 γ (ϵ, ·) − 1 (v ◦ T ϵ ) 2

L (R3 )

      ϵ  2 |·| = 1 + γ (ϵ, ·) − 1 (v ◦ T )  2 L (BR )  ≤ 1 + R2 ∥γ (ϵ, ·) − 1∥L∞ (R3 ) ∥v ◦ T ϵ ∥L2 (R3 ) ≤ C ϵ ∥v ◦ T ϵ ∥L2 (R3 ) .

Using the change of variables y = T ϵ (x) and noting (3.11), we have

∥v ◦ T ϵ ∥2L2 (R3 ) =

 R3

 −1 |v(y )|2 γ (ϵ, (T ϵ )−1 (y )) dy ≤ C ∥v∥L2 (R3 ) .

Therefore,

 

 

lim  1 + |·|2 γ (ϵ, x) − 1 (v ◦ T ϵ )

ϵ→0





L2 (R3 )

= 0.

(3.18)

Furthermore, (3.3) also gives

      1 + |·|2 v ◦ T ϵ − v   2

L (R3 )

     ϵ 2  |·| = v ◦ T − v 1 +  

≤ L2 (BR )



1 + R2 ∥v ◦ T ϵ − v∥L2 (BR ) .

Note that limϵ→0 ∥v ◦ T ϵ − v∥L2 (BR ) = 0 if v is continuous. By using a density argument we deduce that limϵ→0 ∥v ◦ T ϵ − v∥L2 (BR ) = 0 for v ∈ L2 (BR ). Hence,

     ϵ 2  lim 1 + |·| v ◦ T − v  2 ϵ→0 

L (R3 )

= 0.

The above identity and (3.18) together with the triangle inequality give the required result.



Lemma 3.4. For any function v ∈ H 1 (R3 ), there holds

 

2 lim   1 + |·|

ϵ→0



    γ (ϵ, ·)(v ◦ T ϵ ) − v − div v V   2 3 = 0. ϵ L (R )

Proof. Noting (3.3), Lemma 3.1, (3.17) and the triangle inequality, we obtain

  γ (ϵ, ·)(v ◦ T ϵ ) − v    − div(v V )  2 3  1 + |·|2 L (R ) ϵ    γ (ϵ, ·)(v ◦ T ϵ ) − v    2  |·| = 1 + − div v V 2  ϵ L (BR )    γ (ϵ, ·)(v ◦ T ϵ ) − v   . − div v V  2  ϵ L (BR )      γ (ϵ, ·) − 1  v ◦ Tϵ − v  ϵ    ≤ (v ◦ T ) − v divV  + − V · ∇v  2 . ϵ ϵ L (BR ) L2 (BR )

(3.19)

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

9

Recall from (3.9) that γ1 = divV . It follows from (3.12) that

γ (ϵ, ·) − 1 (v ◦ T ϵ ) − v divV = γ1 (v ◦ T ϵ − v) + ϵ(γ2 + ϵγ3 )(v ◦ T ϵ ). ϵ Employing the density argument as in proof of Lemma 3.3 and noting (3.11), we obtain lim ∥γ1 (v ◦ T ϵ − v)∥L2 (BR ) = 0.

ϵ→0

Noting (3.11), we deduce lim ∥ϵ(γ2 + ϵγ3 )(v ◦ T ϵ )∥L2 (BR ) = 0.

ϵ→0

Hence,

    γ (ϵ, ·) − 1 ϵ  (v ◦ T ) − v divV  = 0. lim 2 ϵ→0  ϵ L (BR ) The second term on the right hand side of (3.19) also tends to zero by a density argument, noting that V = ∂ T ϵ /∂ϵ at ϵ = 0. This completes the proof of the lemma.  3.2. Material and shape derivatives 1 In this subsection, for notational convenience we use the notation Dϵ for Dϵ− or Dϵ+ , and H 1 (Dϵ ) for H 1 (Dϵ− ) or Hw (Dϵ+ ), when a result holds for both spaces.

Definition 3.5. For any sufficiently small ϵ , let v ϵ be an element in H 1 (Dϵ ) or H 1/2 (Γ ϵ ). The material derivative of v ϵ , denoted by v˙ , is defined by

vϵ ◦ T ϵ − v0 , ϵ→0 ϵ

v˙ := lim

(3.20)

if the limit exists in the corresponding space H 1 (D0 ) or H 1/2 (Γ 0 ). The shape derivative of v ϵ is defined by

 v˙ − ∇v 0 · V v = v˙ − ∇Γ 0 v 0 · V ′

if v ϵ ∈ H 1 (Dϵ ), if v ϵ ∈ H 1/2 (Γ ϵ ),

(3.21)

where ∇Γ 0 denotes the surface gradient. Lemma 3.6. If v ′ is a shape derivative of v ϵ ∈ H 1 (Dϵ ), then for any compact set K ⊂⊂ D0 we have

vϵ − v0 ϵ→0 ϵ

v ′ = lim

in H 1 (K ).

(3.22)

Proof. Given K ⊂⊂ D0 , there exists an ϵ0 > 0 such that K ⊂⊂ Dϵ for all 0 ≤ ϵ ≤ ϵ0 . We denote by T : [0, ϵ0 ] × R3 → R3 the mapping given by

T (ϵ, x) := T ϵ (x),

∀(ϵ, x) ∈ [0, ϵ0 ] × R3 .

We also denote by v˜ (ϵ, x) := v ϵ (x) for any 0 ≤ ϵ ≤ ϵ0 and x ∈ Dϵ . By the definition of material derivative, we have

  ∂ v˙ = v˜ (ϵ, T (ϵ, ·)) , ∂ϵ ϵ=0

in H 1 (K ).

Applying the chain rule, we obtain

∂ v˜ ∂ T (0, ·) (0, T (0, ·)) + ∇ v˜ (0, T (0, ·)) · ∂ϵ ∂ϵ ∂ v˜ (0, ·) 0 1 = + ∇v · V , in H (K ). ∂ϵ

v˙ =

This implies

v′ =

∂ v˜ (0, ·) vϵ − v0 = lim ϵ→0 ∂ϵ ϵ

in H 1 (K ).



10

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

Remark 3.7. The limit in the above lemma does not hold in H 1 (D0 ) since, in general, v ϵ is not defined in D0 . Similar definitions can be introduced for vector functions v. The following lemmas state some useful properties of material and shape derivatives which will be used frequently in the remainder of the paper. Lemma 3.8. Let v˙ , w ˙ be material derivatives, and v ′ , w′ be shape derivatives of v ϵ , wϵ in H 1 (Dϵ ), ϵ ≥ 0, respectively. Then the following statements are true. (i) The material and shape derivatives of the product v ϵ w ϵ are v˙ w 0 + v 0 w ˙ and v ′ w 0 + v 0 w ′ , respectively. ϵ ϵ 0 0 2 (ii) The material and shape derivatives of the quotient v /w are (˙v w − v 0 w)/(w ˙ ) and (v ′ w 0 − v 0 w ′ )/(w 0 )2 , respectively, provided that all the fractions are well-defined. (iii) If v ϵ = v for all ϵ ≥ 0, then v˙ = ∇v 0 · V = ∇v · V and v ′ = 0. (iv) If ϵ

J1 (D ) :=



ϵ

Dϵ

ϵ

v dx,

J2 (D ) :=

 Γϵ

v ϵ dσ ,

and dJi (Dϵ )|ϵ=0 := lim

Ji (Dϵ ) − Ji (D0 )

ϵ

ϵ→0

,

i = 1, 2,

then ϵ

dJ1 (D )|ϵ=0 =

 D0

v dx + ′

 Γ0

  v 0 V , n0 dσ

and dJ2 (Dϵ )|ϵ=0 =

 Γ0

v ′ dσ +



 Γ0

   ∂v 0 + divΓ 0 (n0 ) v 0 V , n0 dσ . ∂n

Proof. Statements (i)–(iii) can be obtained by using elementary calculations. Statement (iv) is proved in [15, pages 113, 116].  Lemma 3.9. The material and shape derivatives of the normal field nϵ are given by

˙ = n′ = −∇Γ 0 κ. n Proof. We start by recalling that the material and the shape derivative of surface fields are identical in the case of normal surface perturbation (3.4). Particularly, from (3.4) and (3.21) we find

˙ − n′ = ∇Γ 0 n0 · κ n0 = 0. n Recall that the unit normal vector field nϵ of the perturbed interface Γ ϵ is related to that of the reference interface Γ 0 by ϵ 0 J −⊤ ϵ (T (x)) n (x) . nϵ ◦ T ϵ (x) =  T−⊤ J ϵ (T ϵ (x)) n0 (x) T

Therefore,

˙ = lim n

nϵ ◦ T ϵ (x) − n0 (x)

ϵ

ϵ→0

 =

lim

ϵ

ϵ→0

 =

ϵ JT−⊤ ϵ (T (x)) − I

− lim

  −⊤ ϵ  J ϵ (T (x))n0 (x) − 1 T

ϵ  −⊤ ϵ  d JT ϵ (T (x))n0 (x)   −  dϵ ϵ→0

 ϵ dJT−⊤ ϵ (T (x))    dϵ

ϵ=0

ϵ=0

n0 (x)  lim  −⊤ ϵ→0 J ϵ (T ϵ (x)) n0 (x) T

 n0 (x),

(3.23)

noting from (3.8) that lim JT−⊤ ϵ = lim JT ϵ = I .

ϵ→0

ϵ→0

  Since I = JT−ϵ1 (T ϵ (x)) JT ϵ (x) for all x ∈ R3 , we have 0 = ddϵ JT−ϵ1 JT ϵ |ϵ=0 , which together with the product rule and (3.8) yields 





d     d  −⊤ ϵ d JT ϵ (T (x))  = −(JT 0 )−⊤ (JT⊤ϵ ) (JT 0 )−1 = − (JT ϵ ) = −JV⊤ . dϵ dϵ dϵ ϵ=0 ϵ=0 ϵ=0

(3.24)

A. Chernov et al. / Computers and Mathematics with Applications (

   0 We also have, using the fact that JT−⊤ 0 n  = 1,      1 d  −⊤ 0 2  d  −⊤ 0    −⊤ 0  d  −⊤ 0     J ϵ n    JT ϵ n J ϵ n  = = JT 0 n    dϵ T dϵ T 2 dϵ ϵ=0 ϵ=0 ϵ=0    1 d  −1 −⊤  0 0 1 ⊤ = JT ϵ JT ϵ n , n = − (JV + JV ) n0 , n0 . 2 dϵ 2

)

–

11

(3.25)

A simple calculation reveals that and (JV⊤ + JV )n0 = ∇κ + ∇κ, n0 n0 .

JV⊤ = ∇κ (n0 )⊤





(3.26)

Inserting (3.24)–(3.26) into (3.23), we obtain

˙ = −JV⊤ n0 + n

   1 ⊤ (JV + JV ) n0 , n0 n0 = −∇κ + ∇κ, n0 n0 = −∇Γ 0 κ,

2

finishing the proof of the lemma.



3.3. Shape derivative of solutions of transmission problems In this subsection, we shall discuss the existence of material and shape derivatives of the solutions of transmission problems on perturbed interfaces. Consider a deterministic problem with respect to the reference interface Γ 0 :

−α 1u0 = f in D0− ∪ D0+ ,  0 u = 0 on Γ 0 ,   ∂ u0 = 0 on Γ 0 , α ∂n u0 (x) = O (|x|−1 )

(3.27a) (3.27b) (3.27c)

when |x| → ∞.

(3.27d)

The perturbed problem corresponding to the perturbed interface Γ ϵ is given by

−α ϵ 1uϵ = f ϵ

αϵ

ϵ

∂u ∂n

(3.28a)

ϵ

[u ] = 0



in Dϵ− ∪ Dϵ+ ,



on Γ ,

(3.28b)

= 0 on Γ ϵ ,

(3.28c)

uϵ (x) = O (|x|−1 ) when |x| → ∞,

(3.28d)

where (cf. (2.10)) ϵ

α (x) =



α− , α+ ,

x ∈ Dϵ− x ∈ Dϵ+ .

Lemma 3.10. Let u0 and uϵ be solutions of (3.27) and (3.28), respectively. Suppose f ∈ L2 (R3 ) ∩ W0∗ and κ ∈ C 1 (Γ 0 ), then lim uϵ ◦ T ϵ − u0 W = 0.





ϵ→0

(3.29)

0

Here, W0∗ denotes the dual space of W0 with respect to the L2 -inner product. Proof. By multiplying both sides of (3.28a) with an arbitrary function v ∈ C0∞ (R3 ) and integrating over Dϵ− ∪ Dϵ+ , we obtain

 R3

f v dx = −α−

 Dϵ−

1uϵ (x) v(x) dx − α+

 Dϵ+

1uϵ (x) v(x) dx.

(3.30)

Applying Green’s identity and noting (3.28c), we obtain

 Dϵ+ ∪Dϵ−

α ϵ (x) ∇ uϵ (x) · ∇v(x) = ⟨f , v⟩L2 (R3 )

∀v ∈ C0∞ (R3 ).

(3.31)

12

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

Since the space C0∞ (R3 ) is dense in Wϵ (see [22, Remark 2.9.3]), there holds

 Dϵ+ ∪Dϵ−

α ϵ (x) ∇ uϵ (x) · ∇v ϵ (x) = ⟨f , v ϵ ⟩L2 (R3 )

∀v ϵ ∈ Wϵ .

(3.32)

Choosing v ϵ = uϵ gives

|uϵ |2Wϵ ≃ ⟨f , uϵ ⟩L2 (R3 ) ≤ ∥f ∥Wϵ∗ ∥uϵ ∥Wϵ . It follows from Lemma 2.2 that

∥uϵ ∥Wϵ . ∥f ∥Wϵ∗ ≃ ∥f ∥W0∗ .

(3.33)

On the other hand, using the change of variables x = T ϵ (y ) in (3.32), we have (noting that α ϵ (T ϵ (y )) = α(y ))

 D0+ ∪D0−

α(y ) (∇w(y ))⊤ A(ϵ, y ) ∇(uϵ ◦ T ϵ )(y ) dy =

 D0+ ∪D0−

f (T ϵ (y )) w(y )γ (ϵ, y ) dy ,

(3.34)

for any w ∈ W0 . We also obtain from problem (3.27)

 D0+ ∪D0−

α(y ) (∇w(y ))⊤ ∇ u0 (y ) dy =

 D0+ ∪D0−

f (y ) w(y ) dy ,

(3.35)

for any w ∈ W0 . Subtracting (3.35) from (3.34) we deduce

 D0+ ∪D0−

  α(y )∇w(y )⊤ ∇ (uϵ ◦ T ϵ )(y ) − u0 (y ) dy

 =− D0+ ∪D0−

 + D0+ ∪D0−

  ⊤  α(y ) ∇w(y ) A(ϵ, y ) − I ∇(uϵ ◦ T ϵ )(y ) dy   γ (ϵ, y )f (T ϵ (y )) − f (y ) w(y ) dy ∀w ∈ W0 .

Choosing in (3.36) w = uϵ ◦ T ϵ − u0 gives

 D0+ ∪D0−

  2   α(y ) ∇ (uϵ ◦ T ϵ )(y ) − u0 (y )  dy 

=− D0+ ∪D0−

  ⊤   α(y ) ∇ (uϵ ◦ T ϵ )(y ) − u0 (y ) A(ϵ, y ) − I ∇(uϵ ◦ T ϵ )(y ) dy

 (uϵ ◦ T ϵ )(y ) − u0 (y )  dy D0+ ∪D0− 1 + |y |2      .  A(ϵ, ·) − I L∞ (R3 ) ∥∇(uϵ ◦ T ϵ )∥L2 (R3 ) ∇ uϵ ◦ T ϵ − u0 L2 (R3 )    uϵ ◦ T ϵ − u0       ϵ 2  |·| + 1 + γ (ϵ, ·) f ◦ T − f 2 3  2 3   L (R ) 2 L (R ) 1 + |·| 

+





1 + |y |2 γ (ϵ, y )f (T ϵ (y )) − f (y )

implying

 ϵ  u ◦ T ϵ − u0  . ∥A(ϵ, ·) − I ∥L∞ (R3 ) ∥∇(uϵ ◦ T ϵ )∥L2 (R3 ) W0      ϵ 2  +  1 + |·| γ (ϵ, ·)f ◦ T − f  .  L2 (R3 )

Hence, applying Lemma 3.2, noting (3.33) and Lemma 3.3, we obtain lim uϵ ◦ T ϵ − u0 W = 0,



ϵ→0



0

finishing the proof of this lemma.



(3.36)

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

13

Lemma 3.11. Let u0 and uϵ be solutions of (3.27) and (3.28), respectively. Assume that f ∈ H 1 (R3 ) ∩ W0∗ and κ ∈ C 1 (Γ 0 ). Then, uϵ has a material derivative belonging to W0 which is the solution to the following equation with unknown z:

 D0+ ∪D0−

α(y ) ∇ z (y ) · ∇w(y ) dy = −



 ⊤ α(y )∇ u0 (y )A′ (0, y ) ∇w(y ) dy

D0+ ∪D0−

 + D0+ ∪D0−

div (V (y )f ) w(y ) dy

∀w ∈ W0 .

(3.37)

Proof. The uniqueness and existence of the solution z ∈ W0 to the above equation is confirmed by [22, Theorem 2.10.14]. Let z ϵ := (uϵ ◦ T ϵ − u0 )/ϵ . Our task is to prove that limϵ→0 ∥z ϵ − z ∥W0 = 0. Dividing (3.36) by ϵ we obtain

 0

α(y ) ∇ z ϵ (y ) · ∇w(y ) dy = −

0

 0

0

α(y )∇(uϵ ◦ T ϵ )(y )

A(ϵ, y ) − I 

D+ ∪D−

D+ ∪D−

⊤ ∇w(y ) dy

ϵ

ϵ

 + 0

0

D+ ∪D−

γ (ϵ, y )f (T (y )) − f (y ) w(y ) dy ∀w ∈ W0 . ϵ

(3.38)

Subtracting (3.37) from (3.38) yields

 D0+ ∪D0−

α(y ) ∇ (z ϵ (y ) − z (y )) · ∇w(y ) dy

  A(ϵ, y ) − I α(y ) ∇(uϵ ◦ T ϵ )(y ) − ∇ u0 (y )A′ (0, y ) · ∇w(y ) dy ϵ D0+ ∪D0−     ϵ γ (ϵ, y )f (T (y )) − f (y ) − div V (y )f (y ) w(y ) dy + ϵ D0+ ∪D0− 

=−

=: I1 (w) + I2 (w).

(3.39)

The first integral on the right hand side of (3.39) can be written as I1 (w) =

 D0+ ∪D0−

α(y ) ∇(uϵ ◦ T ϵ )(y )

 + D0+ ∪D0−



A(ϵ, y ) − I

 − A′ (0, y ) · ∇w(y ) dy

ϵ  ϵ  α(y ) ∇ (u ◦ T ϵ )(y ) − u0 (y ) A′ (0, y ) · ∇w(y ) dy ,

which converges to 0 due to (3.29) and Lemma 3.2. The second integral on the right hand side of (3.39) also converges to 0 due to Lemma 3.4. Therefore, we have

 lim

ϵ→0 D0 ∪D0 + −

α(y ) ∇ (z ϵ (y ) − z (y )) · ∇w(y ) dy = 0 ∀w ∈ W0 .

(3.40)

We choose in (3.39) w = z ϵ − z. Then the absolute value of the first integral on the right hand side of (3.39) can be estimated as

  |I1 (z − z )| =  ϵ

ϵ

D0+ ∪D0−

α(y ) ∇ u (y )

 + D0+ ∪D0−



A(ϵ, y ) − I

ϵ

   − A (0, y ) · ∇ z ϵ (y ) − z (y ) dy ′

     α(y ) ∇ uϵ (y ) − u0 (y ) A′ (0, y ) · ∇ z ϵ (y ) − z (y ) dy ,   A(ϵ, ·) − I

. ∥∇ uϵ ∥L2 (R3 )  

ϵ

  − A′ (0, ·) 

L∞ (R3 )

∥∇(z ϵ − z )∥L2 (R3 )

     + ∇ uϵ − u0 L2 (R3 ) A′ (0, ·)L∞ (R3 ) ∥∇(z ϵ − z )∥L2 (R3 ) .

(3.41)

The absolute value of the second integral in (3.39) when w = z ϵ − z is bounded by

       γ (ϵ, y )f (T ϵ (y )) − f (y ) 2  |I2 (z − z )| ≤  |·| 1 + − div V ( y ) f ( y )  2 ϵ ϵ

L (R3 )

∥z ϵ − z ∥W0 .

(3.42)

14

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

Inequalities (3.41) and (3.42) give

   A(ϵ, ·) − I  ′  ∥z ϵ − z ∥W0 ≤ ∥∇ uϵ ∥L2 (R3 )  − A ( 0 , ·)  ∞ 3 ϵ L (R )     ϵ  ′ 0     + ∇ u − u L2 (R3 ) A (0, ·) L∞ (R3 )        γ (ϵ, y )f (T ϵ (y )) − f (y ) 2  |·| + 1 + − div V ( y ) f ( y )  2 3 . ϵ L (R )

(3.43)

Using this together with (3.29) and Lemma 3.2, we can deduce from (3.43) lim ∥z ϵ − z ∥W0 = 0. 

(3.44)

ϵ→0

Hence, we have shown that the solution of the transmission problem (3.28) has a material derivative, and thus a shape derivative. The latter turns out to be the solution of a transmission problem on the nominal interface Γ 0 . We state and prove that result in the following theorem. Theorem 3.12. Under the assumption of Lemma 3.11, the shape derivative u′ of uϵ exists and is the solution of the transmission problem

 ′ 0 1 D0+   u′  = 0 in D− ∪  0   u = gD on Γ ∂ u′ = gN on Γ 0 α   ∂n    ′    u (x) = O |x|−1 as |x| → ∞,

(3.45)

where

    ∂ u0 κ and gN := ∇Γ 0 · κ α∇Γ 0 u0 . gD := − ∂n 

Proof. The existence of u′ is confirmed by Lemma 3.11. In this proof only, for notational convenience, we use nϵ± to indicate the normal vector to Γ ϵ pointing outwards Dϵ± , respectively. Note here that nϵ = nϵ− = −nϵ+ . From (3.32) we deduce

α−

 Dϵ−

∇ uϵ−

· ∇v dx + α+

 Dϵ+

∇ uϵ+ · ∇v dx = ⟨f , v⟩L2 (R3 )

∀v ∈ C0∞ (R3 ).

(3.46)

Denoting J (Dϵ± ) := α±

 Dϵ±

∇ uϵ± (x) · ∇v(x) dx

and using Green’s formula, we obtain J (Dϵ± ) = −α±

 Dϵ±

uϵ± (x)1v(x) dx + α±

 Γϵ

uϵ± (x)

∂v dσ =: J1 (Dϵ± ) + J2 (Dϵ± ). ∂ nϵ±

By Lemma 3.8, u′ 1v is the shape derivative of uϵ 1v . On the other hand, by Lemmas 3.8–3.9, the shape derivative of           ∂v  ϵ 0 ϵ ∂v  ′ ∂v  0 0 = ∇v · n is −∇ is u − u ∇ . 0 v · ∇Γ 0 V , n , so that the shape derivative of u 0 v · ∇Γ 0 V , n ϵ ϵ 0 Γ Γ ∂n  ∂n  ∂n  Γϵ

Γϵ

Using Lemma 3.8, we deduce dJ1 (Dϵ± )|ϵ=0

 = −α± 0

u′± (x)1v dx

D±

− α±



u0 1v V , n0± dσ



Γ0

Γ0



and

      ∂v ∂  0 ∂v   0 0 − u ∇ v · ∇ V , n d σ + α u V , n0± dσ 0 0 ± Γ Γ 0 0 0 ∂ n± ∂ n± Γ0 Γ 0 ∂ n±    ∂v + α± divΓ 0 (n0± )u0 0 V , n0± dσ , ∂ n± Γ0

dJ2 (Dϵ± )|ϵ=0 = α±





u′±

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

15

since u0− = u0+ on the interface Γ 0 by (3.27b). Therefore, differentiating by ϵ both sides of (3.46), using Green’s formula, the jump condition (3.27c) and noting that 1v = ∆Γ 0 v + divΓ 0 (n0 )∂v/∂ n0 + ∂ 2 v/(∂ n0 )2 , we obtain 0 = α−

 0

∇ u′ · ∇v dx + α+

 0





u0 ∇Γ 0 v · ∇Γ 0 V , n0− − α+



Γ0





u0 V , n0− ∆Γ 0 v dσ − α+



Γ0

D+

D−

− α−

∇ u′ · ∇v dx − α−



u0 V , n0+ ∆Γ 0 v dσ



Γ0



u0 ∇Γ 0 v · ∇Γ 0 V , n0+ .



Γ0





(3.47)

Applying the tangential Green formula on the third and the fourth integrals on the right hand side of the above identity and the product rule, the above identity can be written as 0 = α−

 0

∇ u′ · ∇v dx + α+



∇ u′ · ∇v dx + 0

Γ0

D+

D−



  (α− ∇Γ 0 u0− − α+ ∇Γ 0 u0+ ) · ∇Γ 0 v V , n0− dσ .

(3.48)

We choose in (3.48) v ∈ C0∞ (D± ) to obtain

α 1u′ (x) = 0,

x ∈ D0± .

(3.49)

We now choose v ∈ C0∞ (R3 ) and applying the Green’s identity to the first two integrals on the right hand side of (3.48), noting (3.49), to obtain 0 = α−

 Γ0

v

∂ u′− dσ + α+ ∂ n0−

 Γ0

v

∂ u′+ + ∂ n0+

 Γ0

  (α− ∇Γ 0 u0− − α+ ∇Γ 0 u0+ ) · ∇Γ 0 v V , n0− dσ .

(3.50)

Applying the tangential Green formula on the surface Γ 0 to the last term on the right hand side of the above identity, we deduce



      ∂ u′ v∇Γ 0 · V , n0− α∇Γ 0 u0 dσ , v α 0 dσ = ∂n Γ0 Γ0

yielding



α

    ∂ u′ = ∇Γ 0 · V , n0− α∇Γ 0 u0 on Γ 0 . 0 ∂n

(3.51)

Recalling the transmission conditions (3.28b), we have for any smooth function v

 Γϵ

[uϵ ] v dσ = 0.

Differentiating by ϵ both sides, applying Lemma 3.8 we have

 0 = d

ϵ

Γϵ

[u ] v dσ



 =d

Γ−ϵ

uϵ− v dσ



 −

Γ+ϵ

uϵ+ v dσ

   ∂(u0− v) 0 0 + div ( n )( u v) V , n0− dσ 0 Γ − − 0 ∂ n− Γ0 Γ0       ∂(u0+ v) 0 ′ 0 0 − + divΓ 0 (n+ )(u+ v) V , n0+ dσ (u+ v) dσ − 0 0 0 ∂ n Γ Γ +      0   ′    0  ∂v   ∂u 0 0 = u v dσ + v V , n− dσ + u + divΓ 0 (n− ) v V , n0− dσ 0 0 0 0 0 ∂ n ∂ n− Γ Γ Γ    0  ′   ∂u v V , n0− dσ , u v dσ + = 0 Γ0 Γ 0 ∂ n−   noting that u0 = 0. Hence, there holds  0  ′  ∂u  u =− V , n0− =: gD . 0 ∂n 

=

(u0− v)′ dσ +





(3.52)

1 Hence, from (3.49), (3.51) and (3.52), the shape derivative u′ ∈ H 1 (D0− ) × Hw (D0+ ) is the weak solution of the transmission problem (3.45). 

16

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

3.4. Random interfaces In Section 3.2, we have defined material and shape derivatives in which the quantity κ(x) does not contain uncertainty. Since the transmission problem (2.11) is posed on a domain with a random interface (see (2.7)), the shape derivative also depends on ω, and it is necessary to compute the mean and the covariance fields of the random solutions. The result is given in the following theorem. Theorem 3.13. Let uϵ (ω) be the solution of the transmission problem (2.11) with the random interface Γ ϵ (ω) given by (2.7), and let u0 denote the solution of the transmission problem with the reference interface Γ 0 . Assume that the perturbation function κ belongs to Lk (Ω , C 1 (Γ 0 )) for an integer k and f ∈ H 1 (R3 ) ∩ W0∗ . Then, for any compact subset K ⊂⊂ D0± , the expectation and the kth order central moments of the solution uϵ (ω) can be approximated, respectively, by

E[uϵ ] = u0 + o(ϵ) in H 1 (K )

(3.53)

1 Mk [uϵ − E[uϵ ]] = ϵ k Mk [u′ ] + o(ϵ k ) in Hmix (K k ).

(3.54)

and

Moreover 1 Mk [uϵ − u0 ] = ϵ k M k [u′ ] + o(ϵ k ) in Hmix (K k ).

(3.55)

Proof. It follows from Lemma 3.6 and Theorem 3.12 that uϵ (x, ω) = u0 (x) + ϵ u′ (x, ω) + ϵ h(ϵ, x, ω) in H 1 (K ),

(3.56)

where h satisfies limϵ→0 ∥h(ϵ, ·, ·)∥Lk (Ω ,H 1 (K )) = 0. This implies

E[uϵ (x, ·)] = u0 (x) + ϵ E[u′ (x, ·)] + ϵ E[h(ϵ, x, ·)] in H 1 (K ). Here, u′ is the solution of (3.45) in which the function κ defining gD and gN depends on ω and satisfies E[κ] = 0; see (2.9). Since u′ depends linearly on κ , it can be seen that E[u′ ] satisfies (3.45) with zero boundary data, and hence E[u′ ] = 0, yielding (3.53). The definition of the statistical moments (2.4) and simple calculations reveal

Mk [uϵ − u0 ] − ϵ k Mk [u′ ] = ϵ k Mk [u′ + h] − Mk [u′ ] .





It follows from [27, Corollary V.5.1] that

∥Mk [u′ + h] − Mk [u′ ]∥H 1

mix

(K k )

 ≤ E ∥(u′ + h) ⊗ · · · ⊗ (u′ + h) − u′ ⊗ · · · ⊗ u′ ∥H 1

mix

Then by the triangle inequality, binomial formula and Hölder’s inequality with p =

 



E =E 

(K k )

E ∥v1 ⊗ · · · ⊗ vk ∥H 1

(K k )

vi =u′ or h, (v1 ,...,vk )̸=(u′ ,...,u′ )



≤

=

mix



mix

vi =u′ or h, (v1 ,...,vk )̸=(u′ ,...,u′ )





 v1 ⊗ · · · ⊗ vk H 1



  E ∥v1 ∥H 1 (K ) . . . ∥vk ∥H 1 (K )

vi =u′ or h, (v1 ,...,vk )̸=(u′ ,...,u′ )

=

 k     k j k−j E ∥h∥H 1 (K ) ∥u′ ∥H 1 (K ) j =1

j

 1q  1p  k     k jp ′ (k−j)q ≤ E ∥h∥H 1 (K ) E ∥u ∥H 1 (K ) j =1

=

j

 k−k j  kj  k     k E ∥h∥kH 1 (K ) E ∥u′ ∥kH 1 (K ) j =1

j

k j

 (K k )

=: E .

k and q = k− j

A. Chernov et al. / Computers and Mathematics with Applications (

=

k    k j =1

j

)

–

17

∥h∥jLk (Ω ,H 1 (K )) ∥u′ ∥kLk−(jΩ ,H 1 (K ))

= o(1) and (3.55) follows. An analogous estimate holds for

M k [uϵ − E[uϵ ]] − ϵ k Mk [u′ ] = ϵ k Mk [u′ + (h − E[h])] − Mk [u′ ] ,



proving (3.54).





The above lemma states in particular that M k [uϵ − u0 ], M k [uϵ − E[uϵ ]] and ϵ k M k [u′ ] coincide in the limit ϵ → 0, indicating that ϵ k M k [u′ ] may be a good approximation for M k [uϵ − u0 ] and M k [uϵ − E[uϵ ]] if ϵ is small. On the other hand, the task of approximation of ϵ k M k [u′ ] is significantly simpler than approximation of M k [uϵ − u0 ] or M k [uϵ − E[uϵ ]] and reduces to solving the homogeneous transmission problem (3.45). 4. Boundary reduction 4.1. Computation of u′ In this section we briefly recall boundary integral equation methods to solve (3.45). We rewrite here this problem for convenience, with the dependency of the solution on the random event ω explicitly shown. 1 For a.a. ω ∈ Ω , find u′ (ω) ∈ H 1 (D0− ) × Hw (D0+ ) satisfying

 ′ 0 1   u′ (ω) = 0 in D±   0  u (ω) = gD (ω) on Γ ′ ∂ u (ω)  = gN (ω) on Γ 0  α ∂ n0      ′ |u (x, ω)| = O |x|−1 as |x| → ∞.

(4.1)

It is well-known that (see e.g. [21]) the solution of (4.1) can be represented as

  ′   ∂ u− (·, ω)   V (x) − W u′− (·, ω) (x), 0 ′ ∂ n u (x, ω) =   ′    W u′ (·, ω) (x) − V ∂ u+ (·, ω) (x), + ∂ n0

x ∈ D0− ,

(4.2)

x ∈ D0+ ,

where V and W are, respectively, the single and double layer potentials defined by V w(x) =



1 4π

1 Γ0

w(y ) dσy ,

|x − y |

W v(x) =

1 4π

 Γ0

∂ 1 v(y ) dσy , ∂ n0y |x − y |

x ∈ D0± ,

for w ∈ H −1/2 (Γ 0 ) and v ∈ H 1/2 (Γ 0 ). Recalling the definition of the interior and exterior Dirichlet-to-Neumann operators:



1



S± v := V −1 ∓ I + K v, 2

where (see e.g., [21, Sections 1.2, 1.4])

V w(x) := lim V w(y ) ≡



1

|x − y | ±  1 K v(x) := lim W v(y ) ∓ v(x) ≡ y →x y →x y ∈D0

Γ 0 \{x}

w(y )dsy for x ∈ Γ 0 ,

∂ 1 v(y )dsy for x ∈ Γ 0 , 0 |x − y | ∂ n y y ∈D0 ±  1 ∂ 1 K ′ w(x) := lim n0x · ∇y V w(y ) ± w(x) ≡ w(y )dsy for x ∈ Γ 0 , 0 |x − y | y →x 0 2 ∂ n Γ \{x} x 0 2

Γ 0 \{x}

y ∈D±

D v(x) := − lim n0x · ∇y W v(y ) ≡ −(n0x × ∇x )V (n0y × ∇y v)(x) for x ∈ Γ 0 , y →x y ∈D0 ±

we obtain

∂ u′± = S± u′± . ∂ n0

(4.3)

18

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

Hence it follows from (4.2) that u′ (·, ω) =



(V S− − W )(u′− (·, ω)) =: E− (u′− (·, ω)) (W − V S+ )(u′+ (·, ω)) =: E+ (u′+ (·, ω))

in D0− , in D0+ .

(4.4)

Therefore, it suffices to compute u′± (·, ω) on Γ 0 . Denoting

[α S ] := α− S− − α+ S+ , we obtain from the jump conditions in (4.1) that u′− (ω) = u′+ (ω) + gD (ω) on Γ 0

(4.5)

[α S ] u′± (ω) = gN (ω) − (α∓ S∓ )gD (ω) on Γ 0 .

(4.6)

and

Hence it suffices to solve (4.6) for either u′+ or u′− on Γ 0 , and then use (4.5) to compute the other. In the following we show for example how u′+ is computed by (4.6). We note that for a fixed ω ∈ Ω , the right hand side gN (ω) − (α− S− )gD (ω) belongs to H −1/2 (Γ 0 ) and thus the solution u′+ (ω) of (4.6) belongs to H 1/2 (Γ 0 ). The variational form for (4.6) is: Find u′+ (ω) ∈ H 1/2 (Γ 0 ) satisfying B(u′+ (ω), v) = ⟨gN (ω) − (α− S− )gD (ω), v⟩

∀v ∈ H 1/2 (Γ 0 ),

(4.7)

with the bilinear form B(·, ·) and the duality pairing ⟨·, ·⟩ given by B(v, w) :=

 Γ0



([α S ] v)w dσ and ⟨g , v⟩ :=

Γ0

g v dσ

∀v, w ∈ H 1/2 (Γ 0 ), g ∈ H −1/2 (Γ 0 ).

(4.8)

We next show the continuity and ellipticity of the operator [α S ] which confirms existence of the unique solution of Eq. (4.6) for a fixed arbitrary ω. Lemma 4.1. The bilinear form B(·, ·) : H 1/2 (Γ 0 ) × H 1/2 (Γ 0 ) → R is bounded, i.e.

|B(v, w)| ≤ C1 ∥v∥H 1/2 (Γ 0 ) ∥w∥H 1/2 (Γ 0 )

∀v, w ∈ H 1/2 (Γ 0 ),

(4.9)

and H 1/2 (Γ )-elliptic, i.e. B(v, v) ≥ C2 ∥v∥2H 1/2 (Γ 0 )

∀v ∈ H 1/2 (Γ 0 ),

(4.10)

where the positive constants C1 and C2 are independent of v . Proof. The boundedness of the bilinear form B is derived directly from the boundedness of V −1 and K . To prove ellipticity we first note that the hypersingular operator D is H 1/2 (Γ 0 )-semi-elliptic for all closed interface Γ 0 , i.e.,

∀v ∈ H 1/2 (Γ 0 );

⟨D v, v⟩L2 (Γ 0 ) ≥ C |v|H 1/2 (Γ 0 )

(4.11)

∂u

see e.g. [23, Corollary 6.25]. The Cauchy data (u− , ∂ n−0 ) on Γ 0 satisfy



u−



1

I −K

2  ∂ u−  =   0 D ∂n



 u−   ∂ u−  . 1 ′ I +K ∂ n0 2 V

(4.12)

Substituting (4.3) into the second equation of (4.12) gives

∂ u− = D u− + ∂ n0



1 2

I + K′



V −1



1 2

 I +K

u−

on Γ 0 .

This equation and (4.3) yield

 S− = D +

1 2

I + K′



V −1



1 2

 I +K

.

Noting that K ′ is the adjoint operator of K , we have

      1 1 ⟨S− v, v⟩ = ⟨D v, v⟩ + V −1 I + K v, I +K v 2

2

∀v ∈ H 1/2 (Γ 0 ).

(4.13)

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

19

Similarly, the exterior Dirichlet-to-Neumann operator S+ satisfies

 S+ = −D −

1 2

I − K′



V −1



1 2

 I −K

and

      1 1 ⟨S+ v, v⟩ = − ⟨D v, v⟩ − V −1 I − K v, I −K v 2

∀v ∈ H 1/2 (Γ 0 ).

2

(4.14)

From (4.13), (4.14), (4.11) and noting the H 1/2 -ellipticity of the inverse operator of V , we derive

      1 1 ⟨[α S ] v, v⟩ = (α− + α+ ) ⟨D v, v⟩ + α− V −1 I + K v, I +K v 2 2 Γ0       1 1 I − K v, I −K v + α + V −1 2

2

Γ0

  2   1 I + K v & (α− + α+ ) |v|2H 1/2 (Γ 0 ) + α−   1/2  2 H

  2  1  & |v|2H 1/2 (Γ 0 ) +  I + K v  2  1/2 H

(Γ 0 )

(Γ 0 )

  2  1  + α+  I − K v   2

  2  1  + I − K v  2  1/2 H

(Γ 0 )

H 1/2 (Γ 0 )

.

(4.15)

Applying the triangle inequality to the last two terms on the right hand side of the inequality above, we obtain

⟨[α S ] v, v⟩ & |v|2H 1/2 (Γ 0 ) + ∥v∥2H 1/2 (Γ 0 ) & ∥v∥2H 1/2 (Γ 0 ) completing the proof of the lemma.

∀v ∈ H 1/2 (Γ 0 ),



4.2. Computation of the k-moments of u′ In practice we need to compute the k-moment of u′ instead of u′ itself. In this subsection, we show how to compute these moments from u′± (ω) on the interface Γ 0 . Tensorizing and integrating both sides of (4.4) we deduce (noting that E[u′ ] = 0) Cov[u′ ](x1 , x2 ) =



(E−,x1 ⊗ E−,x2 )Cor[u′− ](x1 , x2 ), (E+,x1 ⊗ E+,x2 )Cor[u′+ ](x1 , x2 ),

x1 , x2 ∈ D0− , x1 , x2 ∈ D0+ ,

(4.16)

and in general

Mk [u′ ](x1 , . . . , xk ) =



(E−,x1 ⊗ · · · ⊗ E−,xk )Mk [u′− ](x1 , . . . , xk ), (E+,x1 ⊗ · · · ⊗ E+,xk )Mk [u′+ ](x1 , . . . , xk ),

x1 , . . . , xk ∈ D0− , x1 , . . . , xk ∈ D0+ ,

(4.17)

where E±,xi denotes the acting of E± on u′± (xi ). Eq. (4.17) suggests that the k-moment of the solution u′ in D0± can be computed from the k-moment of the Dirichlet data u′± |Γ 0 on the transmission interface. It is noted that the tensorization of both jump conditions (4.5) and (4.6) results in the cross products of the terms on the right hand sides of these two equations. In particular, the tensorization of the right hand side of (4.5) is in general not easy to compute. Therefore, differently from the case of computing u′ , it is not convenient now to solve the tensorization of (4.6) for only one of M k [u′+ ] and M k [u′− ], and obtain the other one from the tensorization of (4.5). We compute both M k [u′+ ] and M k [u′− ] from (4.6). Consider the tensor product operator [α S ](k) := [α S ] ⊗ · · · ⊗ [α S ] which is a linear mapping 1/2

−1/2

[α S ](k) : Hmix (Γ 0 × · · · × Γ 0 ) → Hmix (Γ 0 × · · · × Γ 0 ), see [28, Proposition 2.4] for more details. Tensorization of equation (4.6) yields for almost all ω ∈ Ω [α S ](k) u′± (ω) ⊗ · · · ⊗ u′± (ω) = ⊗ki=1 gN (ω) − (α∓ S∓ )gD (ω)









−1/2

in Hmix (Γ 0 × · · · × Γ 0 ).

Taking the mean of both sides yields a deterministic k-moment problem, namely (k)

[α S ](k) M k [u′± ] = g± , (k)

where g± := E[⊗ki=1

 gN (ω) − (α∓ S∓ )gD (ω) ].



(4.18)

20

A. Chernov et al. / Computers and Mathematics with Applications (

(k)

)

–

−1/2

Recalling (4.7), the variational formulation for finding M k [u′± ] reads: Given g± ∈ Hmix (Γ 0 × · · · × Γ 0 ), find M k [u′± ] ∈

1/2 Hmix

(Γ 0 × · · · × Γ 0 ) satisfying   1/2 (k) B (M k [u′± ], v) = g± , v ∀v ∈ Hmix (Γ 0 × · · · × Γ 0 ), (4.19)   where B (·, ·) = [α S ](k) ·, · and ⟨·, ·⟩, obtained by tensorization of B(·, ·) and ⟨·, ·⟩ defined in (4.8), are (respectively) a 1/2 −1/2 1/2 bilinear form defined in Hmix (Γ 0 ×· · ·× Γ 0 ), and the duality pairing between Hmix (Γ 0 ×· · ·× Γ 0 ) and Hmix (Γ 0 ×· · ·× Γ 0 ). Proposition 2.4 in [28] implies. 1/2

1/2

1/2

Lemma 4.2. The bilinear form B (·, ·) : Hmix (Γ 0 ×· · ·× Γ 0 )× Hmix (Γ 0 ×· · ·× Γ 0 ) → R is bounded and Hmix (Γ 0 ×· · ·× Γ 0 )elliptic, i.e.,

B (v, w) ≤ C1 ∥v∥H 1/2 (Γ 0 ×···×Γ 0 ) ∥w∥H 1/2 (Γ 0 ×···×Γ 0 ) ,

(4.20)

≤ B (v, v)

(4.21)

mix

mix

and C2 ∥v∥2 1/2

Hmix (Γ 0 ×···×Γ 0 ) 1/2

for all v, w ∈ Hmix (Γ 0 × · · · × Γ 0 ). We are now ready to state the second main result of the paper. Theorem 4.3. The k-moment M k [u′± ] exists as the unique solution of (4.19). In particular, for k = 2 the covariance of [u′± ] can be computed as follows. 1/2

Corollary 4.4. The covariance Cov[u′± ](x, y ) ∈ Hmix (Γ 0 × Γ 0 ) is the unique solution to the following equation:

    ([α S ] ⊗ [α S ]) Cov[u′± ](x, y ) = (∇Γ ,x ⊗ ∇Γ ,y ) · Cov[κ](x, y ) α∇Γ ,x u0 (x) α∇Γ ,y u0 (y )   0   0    ∂ u (y ) ∂ u (x) + (α∓ S∓ ) ⊗ (α∓ S∓ ) Cov[κ](x, y ) ∂ nx ∂ ny   0      ∂ u (y ) − ∇Γ ,x · ⊗(α∓ S∓ ) Cov[κ](x, y ) α∇Γ ,x u0 (x) ∂ ny   0      ∂ u (x) . − (α∓ S∓ ) ⊗ ∇Γ ,y · Cov[κ](x, y ) α∇Γ ,y u0 (y ) ∂ nx

(4.22)

(2)

Proof. A simple calculation reveals that the right hand side of (4.22) is indeed g± . The result is a consequence of Theorem 4.3.  5. Examples In this section, we consider the transmission problem (2.11) where the random interface Γ ϵ (ω) is given by

Γ ϵ (ω) = {x + ϵκ(x, ω)n(x) : x ∈ S}. Here, the reference interface Γ 0 is the unit sphere S. The perturbation parameter κ(x, ω) = a(ω) has a constant (but random) value over the whole Γ 0 . The random variable a(ω) takes values in [−1, 1] and is centred so that E[κ] ≡ E[a] = 0. We consider two models for the interface perturbations:

• a(ω) is a uniformly distributed random variable with values in [−1, 1] and probability density function (PDF) ρ1 (t ) = 1/2, so that  1 1 Cov[κ](x, y ) ≡ E[a2 ] = t 2 ρ1 (t ) dt = ; (5.1) −1

3

• a(ω) is a random variable with values in [−1, 1] and PDF ρ2 (t ) =  1 1 Cov[κ](x, y ) ≡ E[a2 ] = t 2 ρ2 (t ) dt = . −1

35 32

(1 − t )3 (1 + t )3 , so that

9

The interface Γ ϵ (ω) is the sphere centred at the origin and having radius Rϵ (ω) = 1 + ϵ a(ω).

(5.2)

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

21

5.1. Analytic examples In this section, we choose the right hand side f to be f (x) =

 (4 |x|2 − 1)2 0

if |x| ≤ 1/2, if |x| ≥ 1/2,

where |x| is the Euclidean norm of the vector x. Then solution of the transmission problem with respect to the random interface Γ ϵ (ω) can be analytically computed as follows:

   8 β 2 1 1 1 6 4 2   |x| − |x| + |x| + + , −   α 21 5 6 24 α R  − − ϵ (ω)   1 β uϵ (x, ω) = |x|−1 + ,  105α− Rϵ (ω)     1   |x|−1 , 105α+ where we denote β :=

α− −α+ 105α− α+ 0

if |x| ≤ if

1 2

1 2

,

≤ |x| ≤ Rϵ (ω),

(5.3)

if |x| ≥ Rϵ (ω)

to simplify the expressions. In particular, the exact solution u0 of the transmission problem

on the reference interface Γ is given by (5.3) where Rϵ (ω) = R0 (ω) = 1, i.e.,

   8 2 1 1 1 6 4 2   | | | | | | x − x + x + + β, −   α− 21 5 6 24α−    1 u0 (x) = |x|−1 + β,  105α−     1   |x|−1 , 105α+

if |x| ≤ if

1 2

1 2

,

≤ |x| ≤ 1,

(5.4)

if |x| ≥ 1.

5.1.1. Uniform distribution a(ω) ∼ ρ1 (x) = 12 Eq. (5.3) and simple calculations yield for ϵ small enough

   u0 (x) + β 1 ln(1 + ϵ) − ln(1 − ϵ) − 1 , if |x| < 1, 2 2ϵ E[uϵ (x, ·)] = (5.5)  0 u (x), if |x| > 1.  ln(1+ϵ)−ln(1−ϵ) ϵ 2n Elementary calculus reveals that −1 = ∞ n=1 2n+1 . Therefore, the mean value E[u] in (5.9) agrees with our 2ϵ result (3.53) in Theorem 3.13. In this example, the linearized error is O (ϵ 2 ). By using (5.3) and a simple calculation, we obtain uϵ (x, ω) − E[uϵ (x, ·)] = β



1

−

1 ln(1 + ϵ) − ln(1 − ϵ)

1 + ϵ a(ω) 2  ∞ ∞  =β (−ϵ a(ω))n − n =1

n=1



2ϵ

ϵ

2n



2n + 1

in the domain |x| < 1 and uϵ (x, ω) − E[uϵ (x, ·)] = 0 when |x| > 1 for any fixed x and a sufficiently small ϵ . Recall that a(ω) has finite moments of arbitrary order (in particular E[a2 ] = 31 ) and therefore Covuϵ (x, y ) =

 2 β 3 0,

ϵ 2 + O (ϵ 3 ),

if |x| < 1 and |y | < 1,

(5.6)

if |x| > 1 or |y | > 1.

We test the accuracy of our shape calculus method by computing the covariance of uϵ (ω) via the covariance of the shape derivative. From (5.4), we obtain

  0    ∂ u− ∂ u0+  ∂ u0 1 1 gD = −κ(·, ω) = −κ(·, ω) − = a(ω) − = β a(ω), ∂n ∂n ∂ n |x|=1 105α− 105α+     gN = ∇Γ 0 · κ α∇Γ 0 u0 = 0. 

(5.7)

22

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

Hence, the unique solution of (3.45), and equivalently of (4.22), is the piecewise constant function u′ (x, ω) =



β a(ω), 0,

if |x| < 1, if |x| > 1.

This and the relation E[a2 ] = Cov[u′ ](x, y ) =

 2 β 3 0,

,

1 3

(5.8)

imply if |x| < 1 and |y | < 1, if |x| > 1 or |y | > 1.

This and (5.6) agree with the theoretical result (3.54) and the linearized error in this example is O (ϵ 3 ).

(1 − t )3 (1 + t )3 5.1.2. Polynomial distribution a(ω) ∼ ρ2 (t ) = 35 32 For this distribution of interface perturbation equation (5.3) yields for sufficiently small ϵ    7 (ln(1 + ϵ) − ln(1 − ϵ))(15ϵ 6 − 45ϵ 4 + 45ϵ 2 − 15) + 66ϵ 5 − 80ϵ 3 + 30ϵ  0  ( x ) + β − 1 u ,  96 ϵ7 E[uϵ (x, ·)] = if |x| < 1,    0 u (x), if |x| > 1.

(5.9)

(ln(1+ϵ)−ln(1−ϵ))(15ϵ 6 −45ϵ 4 +45ϵ 2 −15)+66ϵ 5 −80ϵ 3 +30ϵ

7 − 1 = 19 ϵ 2 + O (ϵ 4 ). Therefore, the mean value The Taylor expansion gives 96 ϵ7 E[u] in (5.9) agrees with our result (3.53) in Theorem 3.13. Moreover, Eq. (5.3) implies

   − 1 + 91 ϵ 2 + O (ϵ 4 ) 1 + ϵ a(ω)     = β −ϵ a(ω) + ϵ 2 a(ω)2 − 19 + O (ϵ 3 )

uϵ (x, ω) − E[uϵ (x, ·)] = β



1

in the domain |x| < 1 and uϵ (x, ω) − E[uϵ (x, ·)] = 0 when |x| > 1 for any fixed x and a sufficiently small ϵ . Therefore Covuϵ (x, y ) =

 2 β 9 0,

ϵ 2 + O (ϵ 3 ),

if |x| < 1 and |y | < 1,

(5.10)

if |x| > 1 or |y | > 1.

On the other hand, the shape derivative u′ is explicitly given by (5.8) as the unique solution of (3.45) with transmission conditions (5.7). Recalling that E[a2 ] = 19 we find Cov[u′ ](x, y ) =

 2 β 9 0,

,

if |x| < 1 and |y | < 1, if |x| > 1 or |y | > 1.

This and (5.10) agree with our theoretical result (3.54) and the linearized error in this example is O (ϵ 3 ). Comparing the obtained results with (3.53) and (3.54) we observe that the mean field and the covariance are approximated at a slightly better rate in ϵ : O (ϵ 2 ) instead of o(ϵ) for the mean and O (ϵ 3 ) instead of o(ϵ 2 ) for the covariance. This improvement can be explained by the smooth dependence of the solution uϵ on the perturbation parameter ϵ . 5.2. Numerical example In this example, we solve the problem (2.11) where the right hand side f is given by f (x) = 2 [x21 + x22 + (x3 − 1)2 ]−1/2 (1 − |x|2 ) − 4 [x21 + x22 + (x3 − 1)2 ]−1/2 (|x|2 − x3 )

− 6[x21 + x22 + (x3 − 1)2 ]1/2 .

(5.11)

The deterministic solution of the transmission problem with the reference interface Γ u0± (x) =

1

α±

[x21 + x22 + (x3 − 1)2 ]1/2 (1 − |x|2 ),

x ∈ D0±

0

= S is then (5.12)

Following the method discussed in Section 3, the covariance of the solution is approximated by the covariance of the shape derivative (see Theorem 3.13), which can be obtained by solving Eqs. (4.22) posed on the reference interface Γ 0 = S. In this section we consider the case of the uniform perturbation of the interface a(ω) ∼ ρ1 (t ) = 21 implying Cov[κ] = 31 . Notice that in the case of the polynomial interface perturbation a(ω) ∼ ρ2 (t ) =

35 32

(1 − t )3 (1 + t )3 with the covariance Cov[κ] =

1 9

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

23

Fig. 1. Convergence of the absolute error |Var[u′ ](x) − Var[u′p ](x)| for three points x inside and outside the unit sphere with respect to the order of the hyperbolic cross p.

(and also of any other sufficiently regular perturbation of the interface with constant covariance) provide the same solution as for the uniform perturbation up to a scaling with a constant factor. 2−σ The right hand sides and the solutions of these equations belong to the space Hmix (Γ 0 × Γ 0 ) for any σ > 0. To solve these equations numerically we use the hyperbolic cross tensor approximation spaces of spherical harmonics which are defined by Spδ := span Yℓ,m : ℓ ∈ δp , mi = −ℓi , . . . , ℓi for i = 1, 2 ,





(5.13)

where

  2  ( 1 + ℓi ) ≤ 1 + p . δp := ℓ = (ℓ1 , ℓ1 ) ∈ N2 :

(5.14)

i=1

The Galerkin method was used to find the approximate solutions Cov[u′p,± ] ∈ Spδ of (4.22). It is shown in [24] that the use

of the space Spδ yields a convergence rate of p−(2−σ −t ) and demands only O p2 log p unknowns, where t is the order of





the Sobolev norm in which the errors are computed. The same convergence rate of p −(2−σ −t ) is achieved when using the standard full tensor product approximation of degree p which meanwhile requires O p4 unknowns. We computed the variance of u′ (x, ω) at three points x = (0, 0, 0.2), (0, 0, 0.5) and (0, 0, 5) inside and outside the unit sphere. The convergence curves for the absolute error

|Var[u′ ](x) − Var[u′p ](x)| with respect to the order of the hyperbolic cross p are presented in Fig. 1. References [1] I. Babuška, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal. 45 (2007) 1005–1034. [2] A. Barth, C. Schwab, N. Zollinger, Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients, Numer. Math. 119 (2011) 123–161. [3] A. Chernov, C. Schwab, First order kth moment finite element analysis of nonlinear operator equations with stochastic data, Math. Comp. 82 (2013) 1859–1888. [4] A. Cohen, R. De Vore, C. Schwab, Convergence rates of best N-term Galerkin approximations for a class of elliptic sPDEs, Found. Comput. Math. 10 (2010) 615–646. [5] A. Cohen, R. Devore, C. Schwab, Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE’s, Anal. Appl. (Singap.) 9 (2011) 11–47. [6] R. Forster, R. Kornhuber, A polynomial chaos approach to stochastic variational inequalities, J. Numer. Math. 18 (2010) 235–255. [7] C.J. Gittelson, An adaptive stochastic Galerkin method for random elliptic operators, Math. Comp. 82 (2013) 1515–1541. [8] I.G. Graham, F.Y. Kuo, D. Nuyens, R. Scheichl, I.H. Sloan, Quasi-Monte Carlo methods for elliptic PDEs with random coefficients and applications, J. Comput. Phys. 230 (2011) 3668–3694. [9] H. Harbrecht, R. Schneider, C. Schwab, Sparse second moment analysis for elliptic problems in stochastic domains, Numer. Math. 109 (2008) 385–414. [10] C. Schwab, C.J. Gittelson, Sparse tensor discretizations of high-dimensional parametric and stochastic PDEs, Acta Numer. 20 (2011) 291–467. [11] C. Schwab, R.A. Todor, Karhunen–Loève approximation of random fields by generalized fast multipole methods, J. Comput. Phys. 217 (2006) 100–122. [12] A. Chernov, Abstract sensitivity analysis for nonlinear equations and applications (Numerical Mathematics and Advanced Applications), in: K. Kunisch, G. Of, O. Steinbach (Eds.), Proceedings of ENUMATH 2007, Graz, Austria, Springer, 2008, pp. 407–414. [13] H. Harbrecht, On output functionals of boundary value problems on stochastic domains, Math. Methods Appl. Sci. 33 (2010) 91–102.

24

A. Chernov et al. / Computers and Mathematics with Applications (

)

–

[14] H. Harbrecht, J. Li, First order second moment analysis for stochastic interface problems based on low-rank approximation, ESAIM Math. Model. Numer. Anal. 47 (2013) 1533–1552. [15] J. Sokołowski, J.-P. Zolésio, Introduction to Shape Optimization. Shape sensitivity analysis, in: Springer Series in Computational Mathematics, vol. 16, Springer-Verlag, Berlin, 1992. [16] F. Hettlich, W. Rundell, The determination of a discontinuity in a conductivity from a single boundary measurement, Inverse Problems 14 (1998) 67–82. [17] F. Hettlich, W. Rundell, Identification of a discontinuous source in the heat equation, Inverse Problems 17 (2001) 1465–1482. [18] A. Kirsch, The domain derivative and two applications in inverse scattering theory, Inverse Problems 9 (1993) 81–96. [19] R. Potthast, Fréchet Differenzierbarkeit von Randintegraloperatoren und Randwertproblemen zur Helmholtzgleichung und den zeitharmonischen Maxwellgleichungen (Ph.D. thesis), Georg-August-Universität Göttingen, 1994. [20] K. Eppler, A shape calculus analysis for tracking type formulations in electrical impedance tomography, J. Inverse Ill-Posed Probl. 17 (2009) 733–751. [21] G.C. Hsiao, W.L. Wendland, Boundary Integral Equations, in: Applied Mathematical Sciences, vol. 164, Springer-Verlag, Berlin, 2008. [22] S.A. Sauter, C. Schwab, Boundary Element Methods, in: Springer Series in Computational Mathematics, vol. 39, Springer-Verlag, Berlin, 2011, Translated and expanded from the 2004 German original. [23] O. Steinbach, Numerical Approximation Methods for Elliptic Boundary Value Problems. Finite and boundary elements, Springer, New York, 2008, Translated from the 2003 German original. [24] A. Chernov, T.D. Pham, Sparse spectral BEM for elliptic problems with random input data on a spheroid, Adv. Comput. Math. 41 (2015) 77–104. [25] W.A. Light, E.W. Cheney, Approximation Theory in Tensor Product Spaces, in: Lecture Notes in Mathematics, vol. 1169, Springer-Verlag, Berlin, 1985, p. vii+157. [26] J.-P. Aubin, Applied Functional Analysis, second ed., in: Pure and Applied Mathematics (New York), Wiley-Interscience, New York, 2000, With exercises by Bernard Cornet and Jean-Michel Lasry, Translated from the French by Carole Labrousse. [27] K. Yosida, Functional Analysis, in: Die Grundlehren der Mathematischen Wissenschaften, Band 123, Academic Press Inc., New York, 1965. [28] T. von Petersdorff, C. Schwab, Sparse finite element methods for operator equations with stochastic data, Appl. Math. 51 (2006) 145–180.