A Boundary Integral Formulation and a Topological Energy-Based Method for an Inverse 3D Multiple Scattering Problem with Sound-Soft, Sound-Hard, Penetrable, and Absorbing Objects

Frédérique Le Louër; María-Luisa Rapún

doi:10.1515/cmam-2021-0234

Article Publicly Available

A Boundary Integral Formulation and a Topological Energy-Based Method for an Inverse 3D Multiple Scattering Problem with Sound-Soft, Sound-Hard, Penetrable, and Absorbing Objects

Frédérique Le Louër and María-Luisa Rapún

Published/Copyright: August 30, 2022

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Computational Methods in Applied Mathematics Volume 22 Issue 4

Abstract

In this paper, we study numerical methods for simulating acoustic scattering by multiple three-dimensional objects of different nature (penetrable, sound-soft, sound-hard and absorbing targets) simultaneously present in the background media. We derive and analyze a boundary integral system of equations that arises when the solution of the problem is represented via single-layer potentials. We give abstract necessary and sufficient conditions for convergence of Petrov–Galerkin discretizations and show that spectral methods satisfy these conditions. Superalgebraic convergence order of the discrete method for smooth objects is illustrated in some test cases. After that, we tackle the inverse problem of finding the shape of objects of different unknown nature from measurements of the total field at a set of receptors. We propose a numerical algorithm based on the computation of the topological energy of a weighted multifrequency least squares cost functional and present some numerical examples to illustrate its capabilities.

Keywords: Boundary Integral Methods; Spectral Methods; Inverse Scattering; Shape Reconstruction; Topological Energy; Multifrequency

MSC 2010: 65N38; 65N20; 65N21

1 Introduction

The study of the scattering of time-harmonic acoustic waves in unbounded media is a topic of great relevance in the boundary element community. There is a huge number of papers where different boundary integral formulations have been proposed and analyzed, aiming at developing fast and accurate solvers for simulating the scattering of acoustic waves (namely, for solving the so-called forward problem). Among them, it is worth highlighting the pioneering works [8, 17, 38, 40, 72, 73] and the monographs [13, 34, 39]. Conversely, in the inverse problems community, a recurrent subject is the retrieval of the support of unknown scatterers and/or their physical properties from measurements of time-harmonic acoustic waves at a few observation points. These problems arise in a wide variety of physical and engineering applications like in bio-medical imaging, non-destructive testing of materials, geophysical explorations and sonar and radar applications, to mention a few. In some of these applications, not only the support of the scatterers is unknown, but also their physical properties might be unknown. For instance, when searching for anti-personnel landmines, in principle, it is unknown whether they are made of metal or plastic. Similarly, when aiming to detect chemical deposits, it is unknown what material they are made of or if they are rusted, which is the degree of corrosion (modeled by the impedance parameter in the boundary condition). As general references for inverse problems (for theory and/or applications), we refer to the monographs [16, 35, 37, 54, 63].

In this paper, we deal with the scattering of time-harmonic waves from both the direct and the inverse points of view. In the first part of the paper, we focus on the direct problem, where a finite number of scatterers of different physical nature will be immersed in R 3 . As already mentioned, there exists a large number of papers where different boundary integral formulations have been proposed for exterior problems governed by the Helmholtz equation, most of them devoted to objects of the same nature. Among them, it is unavoidable to mention the seminal papers by T. von Petersdorff [73] and R. Kress & G. F. Roach [40], where the simultaneous appearance of sound-soft, sound-hard and penetrable objects is considered in the first one via a direct formulation, while sound-soft and sound-hard objects are studied in the second one by an indirect formulation. The recent paper [64] proposes an indirect boundary integral formulation based on single-layer potentials where sound-soft, sound-hard and penetrable objects are immersed in R 2 . In this paper, we extend this last approach to the three-dimensional case considering also the possible presence of absorbing objects. Assuming that the boundaries of the scatterers are smooth, we study invertibility of the resulting matrix operator and give necessary and sufficient conditions for the convergence of abstract Petrov–Galerkin methods in the full range of the Sobolev spaces defined on the boundaries of the scatterers. This is done by following the systematic analysis techniques developed by F. J. Sayas and collaborators (see [22, 21, 42, 43, 65, 67, 68]) which combine the use of suitable factorizations of the matrix operators with Fredholm theory. We show that the proposed formulation is equivalent to the original problem if - λ 2 is not a Dirichlet eigenvalue of the Laplace operator inside any of the objects and - μ i 2 is a not Dirichlet eigenvalue for the Laplace operator inside the penetrable object Ω i (𝑖 being each of the indexes for the penetrable objects), where 𝜆 and μ i are the wavenumbers in the exterior domain and inside the 𝑖-th penetrable object, respectively. After that, we show that, when choosing the spaces of spherical harmonics to approximate the unknowns, the method has superalgebraic convergence order. We numerically corroborate this fact for two representative configurations of scatterers when considering different excitation frequencies. It is mandatory here to highlight the papers [26, 29] by I. G. Graham and collaborators dealing with the full discretization of all the boundary integral operators involved in our formulation when spherical harmonics are considered. We point out that a possibility to circumvent the appearance of spurious eigenvalues is using a Brakhage–Werner potential [7, 49, 57] or a Burton–Miller formulation [8]. Techniques in the present paper can be adapted to these types of formulations with the drawback of having to deal with hypersingular operators. For this extension, we refer to the results in the paper [66] which deals with the use of Brakhage–Werner potentials for a transmission problem in 2D. We end this part devoted to the forward problem by giving some hints about the numerical implementation of the spectral method for the 3D hypersingular operator in case one of these options would be preferred. This alternative combines the spectrally accurate algorithm developed for solving vector boundary integral equations in electromagnetism [27] with differential geometry tools on tangent and cotangent planes [55, pages 71–75].

In this paper, boundary conditions characterizing penetrable and absorbing obstacles are described by constant parameters. However, at least in two dimensions, it is possible to include piecewise constant homogeneous obstacles in the physical framework leading to multitrace boundary integral formulations [33]. Increasing the number of transmission interfaces would be better tackled by a volume FEM-BEM formulation as recently proposed by Domínguez, Ganesh and Sayas [20] or Caudron, Antoine and Geuzaine [12].

The second part of the paper is focused on the inverse problem. Our aim is to determine the number, location and shape of a set of unknown scatterers from measurements of the total wave at some observation points when the medium is acoustically excited by planar incident waves. We make no assumptions on the physical nature of the unknown objects. That is, we neither have a priori information about the geometry of the scatterers, nor do we know whether the waves will penetrate them or not. In case they do not penetrate, we also do not know what boundary condition the total wave will satisfy. We only deal in this work with the shape identification problem. The extension to not only identify shapes but also their nature and/or the physical parameters (interior wavenumber in penetrable objects, coefficients in penetrable and impedance boundaries) is out of the scope of the paper, but we plan to deal with this more ambitious problem in future. Shapes will be identified by an inversion non-iterative method based on the evaluation of an indicator function called the topological energy. This function is related with the sensitivity of a least squares functional that measures the discrepancy between the available data (the measured data at the observation points) and the numerically calculated responses of the system to time-harmonic excitations. It is closely related to the topological derivative [23, 56, 70] of such least squares functional but has the advantage that it does not need to know in advance the nature of the unknown scatterers. The topological energy was first introduced by N. Dominguez, V. Gibiat and Y. Esquerre in [19] to invert time-dependent ultrasound data and successfully used in [10] for visualizing micro and nanostructures in a 3D holography problem where all objects had the same nature. It has also been tested against experimental ultrasound [18] and microwave [11] data in laboratory controlled experiments where multiple objects also had the same nature. As far as we know, it has only been applied for shape identification of simultaneous objects of different nature in [64], where multiple penetrable, sound-soft and sound-hard objects were considered in 2D configurations. Here we will extend this approach to the three-dimensional case including absorbing objects as well. We remark that only few papers deal with simultaneous occurrence of objects of different properties. Among them, it is worth mentioning the qualitative methods in [15, 62, 9, 1]. In comparison with these direct methods (linear sampling, factorization, probe and monotonicity methods), our algorithm has the advantage of having a much smaller computational cost. Some numerical experiments showing the capabilities and limitations of the method will be shown at the end of this part.

Notation

Throughout this work, 𝐶 (also C ′ ) will denote a general positive constant independent of the discretization parameter 𝑛. We will use the standard Sobolev spaces H m ⁢ ( Ω ) with m ≥ 0 , where Ω is a general bounded or unbounded domain, and the Sobolev spaces H s ⁢ ( Γ ) for smooth closed surfaces. For a detailed explanation of these spaces in the setting we will be using them, we refer to [53]. We will also deal with the Sobolev spaces defined on the unit sphere 𝐵 for which we refer to [55]. The notation ( ⋅ , ⋅ ) will be used for sesquilinear forms (linear and conjugate linear in the left and the right components, respectively), both for inner products and antiduality products.

2 Forward Problem

In this section, we first mathematically model the scattering of acoustic waves in R 3 produced by objects of different nature. Secondly, we propose and analyze an indirect boundary integral equation method where the unknowns are some densities defined on the boundaries of the objects. After that, convergence of abstract Petrov–Galerkin methods is studied. Finally, we consider a spectral method, giving numerical evidences of its superalgebraic convergence order.

2.1 Mathematical Formulation

Let us consider a finite set of simply connected open scatterers Ω i ⊂ R 3 , i = 1 , … , N scat , with non-intersecting closures, where N scat is the total number of scatterers. Let us denote by N T , N D , N N and N I the numbers of penetrable (transmission boundary conditions), sound-soft (Dirichlet condition), sound-hard (Neumann condition) and absorbing scatterers (impedance boundary condition), satisfying N scat = N T + N D + N N + N I . To refer to each kind of object, we will use the index sets I T = { 1 , … , N T } , I D = { N T + 1 , … , N T + N D } , I N = { N T + N D + 1 , … , N T + N D + N N } and I I = { N T + N D + N N + 1 , … , N T + N D + N N + N I } . In other words, the considered N scat objects will be ordered such that the first N T ones will be penetrable, the subsequent N D ones will be sound-soft, the next N N ones are sound-hard and the last N I ones will be absorbing scatterers.

For simplicity, we assume that the boundaries Γ i = ∂ ⁡ Ω i can be globally parameterized by spherical coordinates as C ∞ surfaces. The extension to more general geometries would need several coordinate charts to cover the surface.

In this work, we deal with the scattering of time-harmonic acoustic waves, namely, with a time dependence of the form e - i ⁢ ω ⁢ t , where ω > 0 is the frequency. Objects of different nature will be immersed in R 3 which will scatter a known incident wave u inc , assumed to solve the Helmholtz equation Δ ⁢ u inc + λ 2 ⁢ u inc = p inc in R 3 for the wavenumber λ = ω c (where c > 0 is the wave velocity in absence of scatterers). Here, p inc is a source with compact support in R 3 ∖ ⋃ i = 1 N scat Ω i ¯ , or else p inc = 0 in the case of plane waves scattering. The wavenumbers for the penetrable objects are defined as μ i = ω c i for the wave speeds c i > 0 , i ∈ I T . We also introduce the transmission parameters ν i > 0 for i ∈ I T , and the impedance coefficients α i for i ∈ I I , which are assumed to satisfy Re ⁡ ( α i ) ≥ 0 . The problem is then modeled by the following set of equations:

(2.1) Δ ⁢ u + μ i 2 ⁢ u = 0 in ⁢ Ω i , i ∈ I T ,

(2.2) Δ ⁢ u + λ 2 ⁢ u = p inc in ⁢ R 3 ∖ ⋃ i = 1 N scat Ω i ¯ ,

(2.3) u | Γ i int - u | Γ i ext = 0 on ⁢ Γ i , i ∈ I T ,

(2.4) ν i ⁢ ∂ n i ⁡ u | Γ i int - ∂ n i ⁡ u | Γ i ext = 0 on ⁢ Γ i , i ∈ I T ,

(2.5) u | Γ i ext = 0 on ⁢ Γ i , i ∈ I D ,

(2.6) ∂ n i ⁡ u | Γ i ext = 0 on ⁢ Γ i , i ∈ I N ,

(2.7) ∂ n i ⁡ u | Γ i ext + ı ⁢ α i ⁢ λ ⁢ u | Γ i ext = 0 on ⁢ Γ i , i ∈ I I ,

(2.8) lim r → ∞ ⁡ r ⁢ ( ∂ r ⁡ ( u - u inc ) - ı ⁢ λ ⁢ ( u - u inc ) ) = 0 ,

where n i is the exterior normal vector to Γ i , i = 1 , … , N scat , and the superscripts ^ext and ^int refer to exterior and interior traces, respectively. Condition (2.8) is the well-known Sommerfeld radiation condition at infinity, which has to be satisfied uniformly for all unitary directions x | x | ∈ R 3 , r := | x | . Existence and uniqueness of a solution to (2.1)–(2.8) can be derived by combining the results in [73] and [50]. Although we focus on acoustic problems where wavenumbers are real, our work can be extended to more general settings where wavenumbers and/or transmission parameters are complex. Conditions on these parameters for existence and uniqueness of a solution can be derived from [17, 38, 41, 72, 73]. In particular, the results extend to the case μ i , λ ∈ { ( 1 + ı ) ⁢ ξ , ξ > 0 } that arises in the study of time-harmonic solutions to the heat equation [65, 67].

2.2 Indirect Boundary Integral Formulation

In this section, we propose a boundary integral formulation based on the use of single-layer potential operators which is an extension to that proposed in [64] for the counterpart problem in R 2 when I I = ∅ . Its analysis will be based on the techniques previously used in the papers [65, 67, 68].

We quote elementary results from [53, Chapters 6 and 7] to introduce the classical potential theory in Sobolev spaces. Let us introduce first the fundamental solution to the Helmholtz equation for a generic wavenumber 𝜌,

(2.9) ϕ ρ ⁢ ( x , y ) = exp ⁡ ( ı ⁢ ρ ⁢ | x - y | ) 4 ⁢ π ⁢ | x - y | .

Notice that, for ρ = 0 , we have that ϕ 0 is the fundamental solution to the Laplace equation. Then, for each smooth boundary Γ j and any given density φ j ∈ H s - 1 / 2 ⁢ ( Γ j ) , s > - 1 2 , we define the classical single-layer potential as

S j ρ ⁢ φ j := ∫ Γ j ϕ ρ ⁢ ( ⋅ , y ) ⁢ φ j ⁢ ( y ) ⁢ d σ ⁢ ( y ) : R 3 ∖ Γ j → C ,

The single-layer potential is bounded from H s - 1 / 2 ⁢ ( Γ j ) to H loc s + 1 / 2 ⁢ ( R 3 ∖ Γ j ) . Evaluations of the Dirichlet and Neumann traces of S j ρ ⁢ φ j on a boundary Γ i , where 𝑖 might be equal or different from 𝑗, involves the following two boundary integral operators defined for any density φ j ∈ H s ⁢ ( Γ j ) , s ∈ R , by

V i ⁢ j ρ ⁢ φ j := ∫ Γ j ϕ ρ ⁢ ( ⋅ , y ) ⁢ φ j ⁢ ( y ) ⁢ d σ ⁢ ( y ) : Γ i → C , J i ⁢ j ρ ⁢ φ j := ∫ Γ j ∂ n i ⁢ ( ⋅ ) ⁡ ϕ ρ ⁢ ( ⋅ , y ) ⁢ φ j ⁢ ( y ) ⁢ d σ ⁢ ( y ) : Γ i → C .

The jump relations of the single-layer potential operator read as follows (see [13, Chapter 7]): for all i , j ,

(2.10) ( S j ρ ⁢ φ j ) | Γ i int = V i ⁢ j ρ ⁢ φ j , ( S j ρ ⁢ φ j ) | Γ i ext = V i ⁢ j ρ ⁢ φ j , ∂ n i ⁡ ( S j ρ ⁢ φ j ) | Γ i int = δ i ⁢ j 2 ⁢ φ j + J i ⁢ j ρ ⁢ φ j , ∂ n i ⁡ ( S j ρ ⁢ φ j ) | Γ i ext = - δ i ⁢ j 2 ⁢ φ j + J i ⁢ j ρ ⁢ φ j ,

where δ i ⁢ j is the Kronecker symbol defined as δ i ⁢ j = 1 if i = j and δ i ⁢ j = 0 if i ≠ j .

For simplicity, we introduce the following notation for the traces of the incident wave:

(2.11) f i D := u inc | Γ i , f i N := ∂ n i ⁡ u inc | Γ i , i = 1 , … , N scat ,

where the superscripts 𝐷 and 𝑁 stand for Dirichlet and Neumann.

We are ready now to introduce our indirect boundary integral formulation of problem (2.1)–(2.8), which will be based on the ansatz

(2.12) u = { S j μ j ⁢ φ j in ⁢ Ω j , j ∈ I T , u inc + ∑ j = 1 N scat S j λ ⁢ ψ j in ⁢ R 3 ∖ ⋃ j = 1 N scat Ω j ¯ .

Notice that, by definition, 𝑢 solves (2.1), (2.2) and (2.8). For the remaining equations, using the jump relations (2.10), we conclude that the field 𝑢 defined in (2.12) satisfies the transmission conditions (2.3) if and only if

(2.13) V i ⁢ i μ i ⁢ φ i - ∑ j = 1 N scat V i ⁢ j λ ⁢ ψ j = f i D , i ∈ I T ,

while the transmission conditions (2.4) read as

(2.14) ν i ⁢ ( 1 2 ⁢ φ i + J i ⁢ i μ i ⁢ φ i ) + ( 1 2 ⁢ ψ i - ∑ j = 1 N scat J i ⁢ j λ ⁢ ψ j ) = f i N , i ∈ I T .

For the sound-soft objects, we find that conditions (2.5) are equivalent to

(2.15) ∑ j = 1 N scat V i ⁢ j λ ⁢ ψ j = - f i D , i ∈ I D ,

while for the sound-hard objects, conditions (2.6) are equivalent to

(2.16) - 1 2 ⁢ ψ i + ∑ j = 1 N scat J i ⁢ j λ ⁢ ψ j = - f i N , i ∈ I N .

Finally, for the absorbing scatterers, conditions (2.7) are equivalent to

(2.17) - 1 2 ⁢ ψ i + ∑ j = 1 N scat J i ⁢ j λ ⁢ ψ j + ı ⁢ α i ⁢ λ ⁢ ∑ j = 1 N scat V i ⁢ j λ ⁢ ψ j = - f i N - ı ⁢ α i ⁢ λ ⁢ f i D , i ∈ I I .

To simplify notation, we group the densities in the column vectors

(2.18) φ = ( φ i ) , i ∈ I T , ψ P = ( ψ i ) , i ∈ I P , where ⁢ P ∈ { T , D , N , I } ,

and the functions appearing in the right-hand sides in equations (2.13)–(2.17) by defining the column vectors

(2.19) f P D := ( f i D ) , f P N := ( f i N ) , i ∈ I P , where ⁢ P ∈ { T , D , N , I } .

We also introduce the diagonal operators

(2.20) V μ := diag ⁡ ( V i ⁢ i μ i ) , J μ := diag ⁡ ( J i ⁢ i μ i ) , N := diag ⁡ ( ν i ⁢ I ) , i ∈ I T , A = diag ⁡ ( ı ⁢ λ ⁢ α i ⁢ I ) , i ∈ I I ,

and the full operators

V P ⁢ Q λ := ( V i ⁢ j λ ) , J P ⁢ Q λ := ( J i ⁢ j λ ) , i ∈ I P , j ∈ I Q , where ⁢ P , Q ∈ { T , D , N , I } .

Using this notation, equations (2.13)–(2.17) can be collected as

(2.21) [ V μ - V T ⁢ T λ - V T ⁢ D λ - V T ⁢ N λ - V T ⁢ I λ N ⁢ ( 1 2 ⁢ I + J μ ) 1 2 ⁢ I - J T ⁢ T λ - J T ⁢ D λ - J T ⁢ N λ - J T ⁢ I λ 0 V D ⁢ T λ V D ⁢ D λ V D ⁢ N λ V D ⁢ I λ 0 J N ⁢ T λ J N ⁢ D λ - 1 2 ⁢ I + J N ⁢ N λ J N ⁢ I λ 0 J I ⁢ T λ + AV I ⁢ T λ J I ⁢ D λ + AV I ⁢ D λ J I ⁢ N λ + AV I ⁢ N λ - 1 2 ⁢ I + J I ⁢ I λ + AV I ⁢ I λ ] ⏟ = ⁣ : M ⁢ [ φ ψ T ψ D ψ N ψ I ] ⏟ = ⁣ : Ψ = [ f T D f T N - f D D - f N N - f I N - Af I D ] ⏟ = ⁣ : f .

To prove invertibility, we will use the Fredholm alternative. With this aim, we introduce first the notation

H P s = ∏ i ∈ I P H s ⁢ ( Γ i ) , where ⁢ P ∈ { T , D , N , I } ,

for the product of Sobolev spaces of the same order 𝑠 related to each type of boundary (namely, H P s is the product of the spaces H s ⁢ ( Γ i ) for all the boundaries of the type 𝑃, where P = T , D , N , I stands for transmission, Dirichlet, Neumann and impedance, respectively).

We also consider the diagonal operators corresponding to the fundamental solution to the Laplace equation related with Dirichlet and transmission boundaries,

V P 0 := diag ⁡ ( V i ⁢ i 0 ) , where ⁢ P ∈ { T , D } .

We will see now that

M : H T s × H T s × H D s × H N s × H I s → H T s + 1 × H T s × H D s + 1 × H N s × H I s

is an isomorphism for all s ∈ R . To do so, we start by gathering in the next lemma some regularity properties of the boundary integral operators involved in the operator 𝐌. They are either well-known results or easy to derive from standard results of pseudodifferential operators (see [13, 53, 71]).

Lemma 1

For all s ∈ R ,

the bounded operator V i ⁢ i ρ : H s ⁢ ( Γ i ) → H s + 1 ⁢ ( Γ i ) is invertible if and only if - ρ 2 is not a Dirichlet eigenvalue for the Laplace operator in Ω i ,
the operators V i ⁢ i ρ - V i ⁢ i 0 : H s ⁢ ( Γ i ) → H s + 1 ⁢ ( Γ i ) are compact,
the operators V i ⁢ j ρ : H s ⁢ ( Γ j ) → H s + 1 ⁢ ( Γ i ) are compact for all i ≠ j ,
the operators J i ⁢ j ρ : H s ⁢ ( Γ j ) → H s ⁢ ( Γ i ) are compact.

The next two results are straightforward adaptations of [64, Theorems 2.2 and 2.3]. However, proofs are included for the sake of completeness.

Theorem 1

Assume that - λ 2 is not a Dirichlet eigenvalue for the Laplace operator in Ω i for all i ∈ I T . Then the operator V T ⁢ T λ : H T s → H T s + 1 is an isomorphism for all s ∈ R .

Proof

By Lemma 1, the operator V T ⁢ T λ - V T 0 : H T s → H T s + 1 is compact, and therefore, V T ⁢ T λ is Fredholm of index zero. The proof is reduced now to showing that V T ⁢ T λ is one-to-one. Indeed, it is sufficient to prove injectivity for a single value of 𝑠. Let us assume then that V T ⁢ T λ ⁢ ψ T = 0 with ψ T ∈ H T - 1 / 2 and define u = ∑ j ∈ I T S j λ ⁢ ψ j . Then u i := u | Ω i ∈ H 1 ⁢ ( Ω i ) satisfies the Helmholtz equation Δ ⁢ u i + λ 2 ⁢ u i = 0 in Ω i , and by the jump relations of the single-layer potential (2.10), it follows that u | Γ i int = ∑ j ∈ I T V i ⁢ j λ ⁢ ψ j = 0 for all i ∈ I T . Then, by uniqueness of the solution of the interior Dirichlet problems (recall that - λ 2 is not a Dirichlet eigenvalue), it follows that u = 0 in Ω i for all i ∈ I T . In the same way, u * := u | R 3 ∖ ⋃ i ∈ I T Ω i ¯ ∈ H loc 1 ⁢ ( R 3 ∖ ⋃ i ∈ I T Ω i ¯ ) satisfies the Sommerfeld radiation condition at infinity and satisfies Δ ⁢ u * + λ 2 ⁢ u * = 0 in R 3 ∖ ⋃ i ∈ I T Ω i ¯ with u * | Γ i ext = 0 for all i ∈ I T . We find then by uniqueness of the exterior Dirichlet problem that u = 0 in R 3 ∖ ⋃ i ∈ I T Ω i ¯ . Finally, by the jump relations of the single-layer potential,

ψ i = ( 1 2 ⁢ ψ i + J i ⁢ i λ ⁢ ψ i ) - ( - 1 2 ⁢ ψ i + J i ⁢ i λ ⁢ ψ i ) = ∂ n i ⁡ u | Γ i int - ∂ n i ⁡ u | Γ i ext = 0 for all ⁢ i ∈ I T ,

which finishes the proof. ∎

Theorem 2

Assume that, for all i ∈ I T , the numbers - μ i 2 are not Dirichlet eigenvalues for the Laplace operator in Ω i and that - λ 2 is not a Dirichlet eigenvalue in Ω i for i = 1 , … , N scat . Then the operator

M : H T s × H T s × H D s × H N s × H I s → H T s + 1 × H T s × H D s + 1 × H N s × H I s

is an isomorphism for all s ∈ R .

Proof

We follow the same steps as in the proof of Theorem 1. We introduce first the principal part of the operator 𝐌 which is defined as

(2.22) M 0 := [ V T 0 - V T 0 0 0 0 1 2 ⁢ N 1 2 ⁢ I 0 0 0 0 0 V D 0 0 0 0 0 0 - 1 2 ⁢ I 0 0 0 0 0 - 1 2 ⁢ I ]

and can be decomposed as

M 0 = [ I 0 0 0 0 1 2 ⁢ N ⁢ ( V T 0 ) - 1 I 0 0 0 0 0 I 0 0 0 0 0 I 0 0 0 0 0 I ] ⁢ [ V T 0 - V T 0 0 0 0 0 1 2 ⁢ ( N + I ) 0 0 0 0 0 V D 0 0 0 0 0 0 - 1 2 ⁢ I 0 0 0 0 0 - 1 2 ⁢ I ] .

The two triangular-matrix operators involved in the decomposition above are invertible (notice that N + I is invertible since it is a diagonal matrix with elements ( ν i + 1 ) ⁢ I and ν i + 1 > 0 ), which implies invertibility of M 0 . Using Lemma 1, we find now that the operator

M - M 0 : H T s × H T s × H D s × H N s × H I s → H T s + 1 × H T s × H D s + 1 × H N s × H I s

is compact, and therefore, 𝐌 is Fredholm of index zero, and it is enough now to show that 𝐌 is one-to-one for s = - 1 2 . Let us assume then that

(2.23) M ⁢ [ φ ψ T ψ D ψ N ψ I ] = [ 0 0 0 0 0 ] ,

where φ , ψ T ∈ H T - 1 / 2 , ψ D ∈ H D - 1 / 2 , ψ N ∈ H N - 1 / 2 and ψ I ∈ H I - 1 / 2 , and define

u = { S i μ i ⁢ φ i in ⁢ Ω i , i ∈ I T , ∑ i = 1 N scat S i λ ⁢ ψ i in ⁢ R 3 ∖ ⋃ i ∈ I T Ω i ¯ .

By definition, u | Ω i ∈ H 1 ⁢ ( Ω i ) for all i ∈ I T and u | R 3 ∖ ⋃ i ∈ I T Ω i ¯ ∈ H loc 1 ⁢ ( R 3 ∖ ⋃ i ∈ I T Ω i ¯ ) , and 𝑢 satisfies

(2.24) Δ ⁢ u + μ i 2 ⁢ u = 0 in ⁢ Ω i , i ∈ I T ,

(2.25) Δ ⁢ u + λ 2 ⁢ u = 0 in ⁢ R 3 ∖ ⋃ i ∈ I T Ω i ¯ ,

(2.26) u | Γ i int - u | Γ i ext = 0 on ⁢ Γ i , i ∈ I T ,

(2.27) ν i ⁢ ∂ n i ⁡ u | Γ i int - ∂ n i ⁡ u | Γ i ext = 0 on ⁢ Γ i , i ∈ I T ,

(2.28) u | Γ i ext = 0 on ⁢ Γ i , i ∈ I D ,

(2.29) ∂ n i ⁡ u | Γ i ext = 0 on ⁢ Γ i , i ∈ I N ,

(2.30) ∂ n i ⁡ u | Γ i ext + ı ⁢ α i ⁢ λ ⁢ u | Γ i ext = 0 on ⁢ Γ i , i ∈ I I ,

(2.31) lim r → ∞ ⁡ r ⁢ ( ∂ r ⁡ u - ı ⁢ λ ⁢ u ) = 0 .

Let us introduce the index set I D ⁢ N ⁢ I := I D ∪ I N ∪ I I . Condition (2.25) implies that

(2.32) Δ ⁢ u + λ 2 ⁢ u = 0 in ⁢ R 3 ∖ ⋃ i = 1 N scat Ω i ¯ ,

and therefore, 𝑢 solves the homogeneous problem (2.24), (2.26)–(2.32), whose unique solution is u = 0 in R 3 ∖ ⋃ i ∈ I D ⁢ N ⁢ I Ω i ¯ . Taking exterior Dirichlet and Neumann traces on Γ i for i ∈ I D ⁢ N ⁢ I , we deduce that

u | Γ i ext = 0 , ∂ n i ⁡ u | Γ i ext = 0 on ⁢ Γ i , i ∈ I D ⁢ N ⁢ I .

For the interior Dirichlet traces, we use the jump relations in (2.10) to derive now that

u | Γ i int = ∑ j = 1 N scat V i ⁢ j λ ⁢ ψ j = u | Γ i ext = 0 , i ∈ I D ⁢ N ⁢ I .

By (2.25), 𝑢 solves the interior Dirichlet problem for i ∈ I D ⁢ N ⁢ I , and since by hypothesis - λ 2 is not a Dirichlet eigenvalue of the Laplace operator in Ω i , we find that u = 0 in ⋃ i ∈ I D ⁢ N ⁢ I Ω i . Therefore, ∂ n i ⁡ u | Γ i int = 0 on Γ i for i ∈ I D ⁢ N ⁢ I , and from the jump relations, we find now that

ψ i = ∂ n i ⁡ u | Γ i int - ∂ n i ⁡ u | Γ i ext = 0 for all ⁢ i ∈ I D ⁢ N ⁢ I .

Namely,

(2.33) ψ D = 0 , ψ N = 0 , ψ I = 0 .

It only remains now to prove that φ = ψ T = 0 . To do that, notice first that, since u = 0 in Ω i for i ∈ I T , we can take the interior Dirichlet trace to deduce that

0 = u | Γ i int = V i ⁢ i μ i ⁢ φ i , i ∈ I T ,

and use that - μ i 2 is not an eigenvalue of the Laplace operator in Ω i and Lemma 1 to derive that φ i = 0 for i ∈ I T , i.e.

(2.34) φ = 0 .

Finally, the first group in equation (2.23) and (2.33)–(2.34) imply that V T ⁢ T λ ⁢ ψ T = 0 , and by Theorem 1, we also get that ψ T = 0 . ∎

Remark

The use of single-layer potentials introduces spurious eigenmodes in the problem. A possible alternative to circumvent this problem is using mixed or Brakhage–Werner potentials (first introduced independently in the papers [7, 49, 57]). The formulation has the drawback of the appearance of weakly singular and hypersingular boundary integral operators. Results in the present paper can be adjusted to the Brakhage–Werner formulation taking into account the results in [66]. Alternatively, we could keep single-layer potentials but use a Burton–Miller formulation of the boundary conditions (see [8]), which again promotes the appearance of hypersingular operators.

2.3 Numerical Approximation

We consider two families of finite-dimensional spaces depending on a parameter n → ∞ ,

X i n ⊂ H - 1 / 2 ⁢ ( Γ i ) , Y i n ⊂ H 1 / 2 ⁢ ( Γ i ) , dim ⁡ X i n = dim ⁡ Y i n , i = 1 , … , N scat ,

and the product spaces

X P n := ∏ i ∈ I P X i n , Y P n := ∏ i ∈ I P Y i n , where ⁢ P ∈ { T , D , N , I } .

To simplify notation, we have assumed that the dimensions of all components are the same, but we could have created families X P n := ∏ i ∈ I P X i n i , Y P n := ∏ i ∈ I P Y i n i with dim ⁡ X i n i = dim ⁡ Y i n i directed on the parameter n = ( n 1 , … , n N scat ) which has to diverge componentwise.

Our aim is to analyze Petrov–Galerkin discretizations of the system of boundary integral equations (2.21) of the form

(2.35) { find ⁢ φ n , ψ T n ∈ X T n , ψ D n ∈ X D n , ψ N n ∈ X N n , ψ I n ∈ X I n ⁢ such that ( M ⁢ [ φ n ψ T n ψ D n ψ N n ψ I n ] , [ η n p n ξ n r n s n ] ) = ( f , [ η n p n ξ n r n s n ] ) for all η n ∈ X T n , p n ∈ Y T n , ξ n ∈ X D n , r n ∈ Y N n , s n ∈ Y I n .

The approximate solution to problem (2.1)–(2.8) is then defined by introducing the discrete densities in the definition of the single-layer potential in (2.12), namely,

(2.36) u n := { S j μ j ⁢ φ j n in ⁢ Ω j , j ∈ I T , u inc + ∑ j = 1 N scat S j λ ⁢ ψ j n in ⁢ R 3 ∖ ⋃ j = 1 N scat Ω j ¯ .

In order to give necessary and sufficient conditions for convergence of schemes of the form (2.35), we recall first some well-known general properties on Petrov–Galerkin schemes (see [39, 75]). For a general isomorphism A : V → W ′ between two Hilbert spaces 𝑉 and 𝑊, a Petrov–Galerkin scheme consists of two families of finite-dimensional subspaces V n ⊂ V and W n ⊂ W with dim ⁡ V n = dim ⁡ W n and a discretization scheme

{ find ⁢ v n ∈ V n ⁢ such that ( A ⁢ v n , w n ) = ( A ⁢ v , w n ) for all ⁢ w n ∈ W n .

For brevity, we call this method the Petrov–Galerkin { V n ; W n } scheme for A : V → W ′ . The method is said to be stable if there exists a constant γ > 0 independent of 𝑛 such that, for 𝑛 large enough, the inf-sup condition

sup 0 ≠ w n ∈ W n ⁡ | ( A ⁢ v n , w n ) | ∥ w n ∥ ≥ γ ⁢ ∥ v n ∥ for all ⁢ v n ∈ V n

holds. Stability is equivalent to unique solvability of the discrete equations plus the inequality

∥ v n ∥ ≤ ( ∥ A ∥ γ ) ⁢ ∥ v ∥ ,

and to the Céa estimate

(2.37) ∥ v - v n ∥ ≤ ( ∥ A ∥ γ ) ⁢ inf u n ∈ V n ⁡ ∥ v - u n ∥ .

Convergence is equivalent to stability plus the approximation property

inf u n ∈ V n ⁡ ∥ v - u n ∥ → 0 for all ⁢ v ∈ V .

We are ready now to characterize convergent methods of the form (2.35).

Theorem 3

Let s ∈ R , and assume that

X P n ⊂ H P s ∩ H P - s - 1 ⁢ for ⁢ P ∈ { T , D , N , I } and Y P n ⊂ H P - s ⁢ for ⁢ P ∈ { T , N , I } .

Then convergence of the method (2.35) is equivalent to the simultaneous convergence of the { X i n ; Y i n } scheme for I : H s ⁢ ( Γ i ) → H s ⁢ ( Γ i ) for all i ∈ I T ∪ I N ∪ I I and of the { X i n ; X i n } scheme for V i ⁢ i 0 : H s ⁢ ( Γ i ) → H s + 1 ⁢ ( Γ i ) for all i ∈ I T ∪ I D .

Proof

Convergence of Petrov–Galerkin methods is preserved by compact perturbations (see [39, Theorem 13.7]), and therefore, we can equivalently show stability for the principal part of the operator 𝐌, which is the operator M 0 : H T s × H T s × H D s × H N s × H I s → H T s + 1 × H T s × H D s + 1 × H N s × H I s introduced in (2.22). This operator can be decomposed as follows:

M 0 = [ V T 0 0 0 0 0 0 I 0 0 0 0 0 V D 0 0 0 0 0 0 I 0 0 0 0 0 I ] ⏟ = ⁣ : M 1 0 ⁢ [ I - I 0 0 0 1 2 ⁢ N 1 2 ⁢ I 0 0 0 0 0 I 0 0 0 0 0 - 1 2 ⁢ I 0 0 0 0 0 - 1 2 ⁢ I ] ⏟ = ⁣ : M 2 0 .

Notice now that M 2 0 is a bounded isomorphism in X T n × X T n × X D n × X N n × X I n since ν i > 0 . Therefore, stability of the Petrov–Galerkin method { X T n × X T n × X D n × X N n × X I n ; X T n × Y T n × X D n × Y N n × Y I n } for the operator M 0 is equivalent to stability of the same method for the operator M 1 0 . Furthermore, M 1 0 is diagonal, and stability reduces now to show the stability of the Galerkin methods { X i n ; X i n } for the operators V i ⁢ i 0 : H s ⁢ ( Γ i ) → H s + 1 ⁢ ( Γ i ) for i ∈ I T ∪ I D and of the Petrov–Galerkin methods { X i n ; Y i n } for the identity operators I : H s ⁢ ( Γ i ) → H s ⁢ ( Γ i ) for i ∈ I T ∪ I N ∪ I I . ∎

2.4 Boundary Integral Formulation on the Unit Sphere

As already mentioned, we assume that each domain Ω j can be globally described in spherical coordinates. This means that there exists a bijective parameterization x j : Γ → Γ j , where Γ = ∂ ⁡ B is the boundary of the unit ball 𝐵 (namely, Γ is the unit sphere), such that

(2.38) ∫ Γ j ψ ⁢ d s ⁢ ( x ) = ∫ Γ ψ ⁢ ( x j ⁢ ( ξ ) ) ⁢ | x j ′ ⁢ ( ξ ) | ⁢ d s ⁢ ( ξ ) ,

where | x j ′ | is the Jacobian of x j . Using these parameterizations, an equivalent parameterized version of the system of boundary integral equations (2.21) can be obtained, where now the unknown densities are defined on Γ and integral operators are also defined on the unit sphere. The right-hand side contains nothing but parameterized versions of the Dirichlet and Neumann traces of the incident wave on the unit sphere.

Let us consider now the ( n + 1 ) 2 -dimensional space S n of spherical harmonics of degree less than or equal to 𝑛, which can be defined in terms of the Legendre functions P ℓ j as the span of the following orthonormal basis system in L 2 ⁢ ( Γ ) (see [16, Section 2.3] or [55, Chapter 2]):

(2.39) Y ℓ , j ⁢ ( ξ ) = ( - 1 ) ( j + | j | ) / 2 ⁢ 2 ⁢ ℓ + 1 4 ⁢ π ⁢ ( ℓ - | j | ) ! ( ℓ + | j | ) ! , P ℓ | j | ⁢ ( cos ⁡ θ ) , e ı ⁢ j ⁢ ϕ , 0 ≤ ℓ ≤ n , | j | ≤ ℓ ,

for ξ = ( sin ⁡ θ ⁢ cos ⁡ ϕ , sin ⁡ θ ⁢ sin ⁡ ϕ , cos ⁡ θ ) ⊤ ∈ Γ . It is well known that Sobolev spaces on Γ = ∂ ⁡ B can be equivalently characterized as

H s ⁢ ( Γ ) = { u = ∑ ℓ = 0 ∞ ∑ j = - ℓ ℓ c ℓ , j ⁢ Y ℓ , j | ∑ ℓ = 0 ∞ ∑ j = - ℓ ℓ ( 1 + ℓ 2 ) s ⁢ | c ℓ , j | 2 < ∞ } ,

endowed with the norm ∥ u ∥ s = ( ∑ ℓ = 0 ∞ ∑ j = - ℓ ℓ ( 1 + ℓ 2 ) s ⁢ | ( u , Y ℓ , j ) | 2 ) 1 / 2 . The space S n satisfies the approximation property in H s ⁢ ( Γ ) since the truncation operator

(2.40) π n ⁢ u = ∑ ℓ = 0 n ∑ j = - ℓ ℓ ( u , Y ℓ , j ) ⁢ Y ℓ , j

verifies that π n ⁢ u → u in H s ⁢ ( Γ ) for all u ∈ H s ⁢ ( Γ ) . Furthermore, for all t ≥ s , we have

(2.41) ∥ u - π n ⁢ u ∥ s ≤ C s , t ⁢ ( 1 + n 2 ) ( s - t ) / 2 ⁢ ∥ u ∥ t for all ⁢ u ∈ H t ⁢ ( Γ ) ,

where C s , t > 0 depends on 𝑠 and 𝑡 but it is independent of 𝑛. We have then just seen that the { S n ; S n } method for I : H s ⁢ ( Γ ) → H s ⁢ ( Γ ) is convergent (take s = t in the equality above to show stability).

Concerning the operator V 0 defined on the sphere Γ,

(2.42) V 0 ⁢ φ 0 := ∫ Γ 1 4 π | ⋅ - y ^ | ⁢ φ 0 ⁢ ( y ^ ) ⁢ d σ ⁢ ( y ^ ) : Γ → C ,

we want to point out first that the functions defined in (2.39) are their eigenfunctions (see [16, Section 3.6]). More precisely, they satisfy

V 0 ⁢ Y ℓ , j = 1 2 ⁢ ℓ + 1 ⁢ Y ℓ , j ,

and consequently, we have that V 0 ⁢ u n ∈ S n for all u n ∈ S n . Furthermore, by definition,

( V 0 ⁢ u , v ) = ( u , V 0 ⁢ v ) for all ⁢ u ∈ H s ⁢ ( Γ ) , v ∈ H - s - 1 ⁢ ( Γ ) .

This means that the solution to

{ find ⁢ u n ∈ S n ⁢ such that ( V 0 ⁢ u n , v n ) = ( V 0 ⁢ u , v n ) for all ⁢ v n ∈ S n

is simply u n = π n ⁢ u , and therefore, also the { S n ; S n } method for V 0 : H s ⁢ ( Γ ) → H s + 1 ⁢ ( Γ ) is convergent.

We have reviewed all the required ingredients to show that a particular example of convergent method for solving (2.21) is a pure Galerkin method with X i n = Y i n = S n for all 𝑖.

Theorem 4

Let S N T + N scat n be the product of spherical harmonics spaces

S N T + N scat n := S n × ⋯ × S n ⏟ N T + N scat times .

Then the method

(2.43) { find ⁢ Ψ n := [ φ n , ψ T n , ψ D n , ψ N n , ψ I n ] ⊤ ∈ S N T + N scat n ⁢ such that ( M ⁢ Ψ n , Ξ n ) = ( f , Ξ n ) for all ⁢ Ξ n ∈ S N T + N scat n

is convergent in H N T + N scat s ⁢ ( Γ ) := H T s ⁢ ( Γ ) × H T s ⁢ ( Γ ) × H D s ⁢ ( Γ ) × H N s ⁢ ( Γ ) × H I s ⁢ ( Γ ) for all s ∈ R . Furthermore, if the solution Ψ to the parameterized version of system (2.21) belongs to H N T + N scat t ⁢ ( Γ ) , then

(2.44) ∥ Ψ - Ψ n ∥ s ≤ C s , t ⁢ ( 1 + n 2 ) ( s - t ) / 2 ⁢ ∥ Ψ ∥ t

for arbitrary s < t with C s , t > 0 depending only on 𝑠 and 𝑡. In addition, the approximate solution u n to the solution 𝑢 of problem (2.1)–(2.8) defined by (2.36) satisfies

(2.45) | u ⁢ ( x ) - u n ⁢ ( x ) | ≤ C s , t , x ′ ⁢ ( 1 + n 2 ) ( s - t ) / 2 ⁢ ∥ Ψ ∥ t for all ⁢ x ∈ R 3 ∖ ⋃ i ∈ I D ⁢ N ⁢ I Ω i ¯ ,

where C s , t , x ′ depends only on s , t and 𝐱.

Proof

By Theorem 3, since for all s ∈ R the { S n ; S n } method is convergent both for V 0 : H s ⁢ ( Γ ) → H s + 1 ⁢ ( Γ ) and I : H s ⁢ ( Γ ) → H s ⁢ ( Γ ) , we already have that the { S N T + N scat n ; S N T + N scat n } method converges for

M : H N T + N scat s ⁢ ( Γ ) → H T s + 1 ⁢ ( Γ ) × H T s ⁢ ( Γ ) × H D s + 1 ⁢ ( Γ ) × H N s ⁢ ( Γ ) × H I s ⁢ ( Γ ) .

Furthermore, the bound (2.44) is immediate from the Céa estimate (2.37) and relation (2.41). Finally, to show estimate (2.45), we just need to note that if x ∈ Ω j for j ∈ I T , then from the definitions of 𝑢 and u n in terms of the single-layer potentials, we can write

where C j ⁢ ( x , s ) is a positive constant that only depends on 𝐱 and 𝑠. Analogously, for x ∈ R 3 ∖ ⋃ j = 1 N scat Ω j ¯ ,

| u ⁢ ( x ) - u n ⁢ ( x ) | = | ∑ j = 1 N scat ( S j λ ⁢ ψ j ⁢ ( x ) - S j λ ⁢ ψ j n ⁢ ( x ) ) | ≤ C ′ ⁢ ( x , s ) ⁢ ∥ Ψ - Ψ n ∥ s ,

C ′ ⁢ ( x , s ) being a positive constant that only depends on 𝐱 and 𝑠. Estimate (2.45) is straightforward now from inequality (2.44). ∎

Notice that, when the boundaries of the objects are C ∞ and data is given by planar incident waves or more general incident waves generated by sources with a compact support away from the objects, then we can ensure that Ψ belongs to H N T + N scat t ⁢ ( Γ ) for all t ∈ R , which implies superalgebraic convergence order.

2.5 Fully Discrete Spectral Approximation

In this section, we present the main ideas of the fully discrete approximation of the weakly singular boundary integral equation system (2.21) based on the spherical parameterization introduced in equation (2.38) followed by a projection on the finite-dimensional spaces X i n and Y i n spanned by the scalar spherical harmonics functions (2.39) for which M. Ganesh and I. G. Graham [26] demonstrated superalgebraic convergence in the case of Dirichlet, Neumann or Robin problems. The idea of the method originates from K. E. Atkinson’s work [4]. Similar quadrature rules also lead to exponentially convergent Nyström type methods [74].

2.5.1 Quadrature Rule

We denote by θ , ϕ the spherical coordinates of any point ξ ∈ Γ , i.e.

ξ ⁢ ( θ , ϕ ) = ( sin ⁡ θ ⁢ cos ⁡ ϕ , sin ⁡ θ ⁢ sin ⁡ ϕ , cos ⁡ θ ) ⊤ , ( θ , ϕ ) ∈ ] 0 ; π [ × [ 0 ; 2 ⁢ π [ ∪ { ( 0 , 0 ) ; ( 0 , π ) } .

The computational implementation of integrals of continuous functions on Γ makes use of quadrature rules of the form

∫ Γ f ⁢ ( ξ ) ⁢ d s ⁢ ( ξ ) ≈ ∑ j = 1 m α j ⁢ f ⁢ ( ξ j )

for a given set of points ξ j ∈ Γ and suitable weights α j , where j = 1 , … , m . Following [26], we select the 2 ⁢ ( n + 1 ) × ( n + 1 ) -point rectangle-Gauss rule

(2.46) ∫ Γ f ⁢ ( ξ ) ⁢ d s ⁢ ( ξ ) ≈ ∑ r = 0 2 ⁢ n + 1 ∑ τ = 1 n + 1 μ r ⁢ ν τ ⁢ f ⁢ ( ξ r , τ ) ,

where the points and weights are given by

ξ r , τ = ( sin ⁡ θ τ ⁢ cos ⁡ ϕ r , sin ⁡ θ τ ⁢ sin ⁡ ϕ r , cos ⁡ θ τ ) ⊤ ,

with θ τ = cos - 1 ⁡ z τ , where z τ , τ = 1 , … , n + 1 , are the zeros of the Legendre polynomial of degree n + 1 , ν τ are the corresponding Gauss-Legendre weights, and μ r = π n + 1 and ϕ r = r ⁢ π n + 1 for r = 0 , … , 2 ⁢ n + 1 . The quadrature rule (2.46) is exact for functions defined in S n . Some other quadrature rules could be used. For a discussion on this topic, we refer to [16, 29] and references therein. Implementation of the numerical method (2.43) described in the previous section involves projections (2.40) in some spaces S n and therefore computations of the products ( ⋅ , Y ℓ , j ) that are evaluated with (2.46).

2.5.2 Spherical Parameterization and Singular Integral Approximations

We consider that, for i = 1 , … , N scat , we have Γ i = x i ⁢ ( Γ ) = c i + q i ⁢ ( Γ ) , where c i is the location of the center of the 𝑖-th object and the C 1 -diffeomorphisms q i : Γ → Γ i are shape parameterizations. The Jacobian | x i ′ | of the change of variables q i : Γ → Γ i and the parameterized normal vectors n ~ i directed outward Γ i can be computed via the formulae

| x i ′ ⁢ ( ξ ) | = | t θ ⁢ ( ξ ) × t ϕ ⁢ ( ξ ) | and n ~ i ⁢ ( ξ ) = t θ ⁢ ( ξ ) × t ϕ ⁢ ( ξ ) | x i ′ ⁢ ( ξ ) | ,

where t θ ⁢ ( ξ ) = ∂ ⁡ q i ⁢ ( ξ ) ∂ ⁡ θ and t ϕ ⁢ ( ξ ) = 1 sin ⁡ θ ⁢ ∂ ⁡ q i ⁢ ( ξ ) ∂ ⁡ ϕ .

A fully detailed explanation of the implementation of the smooth and weakly singular acoustic integral operators involved in our method is nicely specified in [26, 28]. Furthermore, by carefully following the analysis presented in [26], it can be proved that superalgebraic convergence is preserved after full discretization of the equations derived in the present paper. This will be numerically tested in the forthcoming section.

On the other hand, as previously mentioned, the boundary integral equation system derived in Section 2.2 can be adapted to Burton–Miller formulations based on single-layer representation of the exterior field too. Formulation (2.21) would involve the hypersingular operator. In this case, the full system has to be multiplied by the Jacobian of the parameterizations to handle the hypersingularity. Then, compared to the discrete schemes presented in [26, 36], the spherical integral reformulations of the operators V i ⁢ j ρ , J i ⁢ j ρ considered in this work involve the Jacobian | x i ′ | as below:

V ~ i ⁢ j ρ ⁢ φ ~ ⁢ ( ξ ) = | x i ′ ⁢ ( ξ ) | ⁢ ∫ Γ ϕ ρ ⁢ ( c i + q i ⁢ ( ξ ) , c j + q j ⁢ ( ζ ) ) ⁢ φ ~ ⁢ ( ζ ) ⁢ | x j ′ ⁢ ( ζ ) | ⁢ d σ ⁢ ( ζ ) , J ~ i ⁢ j ρ ⁢ φ ~ ⁢ ( ξ ) = | x i ′ ⁢ ( ξ ) | ⁢ ∫ Γ n ~ i ⁢ ( ξ ) ⋅ ∇ x ⁡ ϕ ρ ⁢ ( c i + q i ⁢ ( ξ ) , c j + q j ⁢ ( ζ ) ⁢ φ ~ ⁢ ( ζ ) ) ⁢ | x j ′ ⁢ ( ζ ) | ⁢ d σ ⁢ ( ζ ) ,

where φ ~ ∈ H - 1 / 2 ⁢ ( Γ ) . When i ≠ j , the kernel is smooth and one can use the quadrature rule (2.46). When i = j , the kernel is split in a smooth part and a weakly singular part O ⁢ ( | x i ⁢ ( ξ ) - x i ⁢ ( ζ ) | - 1 ) rewritten in the form

1 | x i ⁢ ( ξ ) - x i ⁢ ( ζ ) | ⁢ f ⁢ ( ξ , ζ ) = 1 | ξ - ζ | × | ξ - ζ | | x i ⁢ ( ξ ) - x i ⁢ ( ζ ) | ⁢ f ⁢ ( ξ , ζ ) ,

where f ⁢ ( ξ , ζ ) represents a continuous bivariable function. We point out that | ξ - ζ | - 1 is nothing but the kernel of the operator V 0 defined at equation (2.42). Next, the variable 𝝃 is rotated onto the north pole of the unit sphere, and then the smooth functions are projected as in equation (2.40) on the set S n ′ spanned by spherical harmonics of order n ′ = a ⁢ n + 1 , a > 1 , where 𝑛 is the index used in the definition of S n . The numerical treatment of the weakly singular boundary integral operators is handled analytically through the formula

∫ Γ 1 4 ⁢ π ⁢ | ξ - ζ | ⁢ π n ′ ⁢ ( u ) ⁢ ( ζ ) ⁢ d ζ = V 0 ⁢ π n ′ ⁢ ( u ) = ∑ ℓ = 0 n ′ ∑ j = - ℓ ℓ 1 2 ⁢ ℓ + 1 ⁢ ( u , Y ℓ , j ) ⁢ Y ℓ , j

for some n ′ ∈ N , in which 𝑢 represents a continuous function on the unit sphere Γ. Moreover, the identity operator on S n is discretized as the mass matrix weighted by the Jacobian of the parameterization as in formula (2.38).

For readers interested in the numerical implementation of the hypersingular operator (which would appear in Burton–Miller and Brakhage–Werner formulations avoiding spurious eigenmodes), it suffices to consider the well-known formulae

∂ n ⁢ ( x ) ⁡ ∂ n ⁢ ( y ) ⁡ ϕ ρ ⁢ ( x , y ) = - n ⁢ ( x ) ⋅ n ⁢ ( y ) ⁢ Δ ⁢ ϕ ρ ⁢ ( x , y ) + n ⁢ ( x ) ⋅ rot x ⁡ ( ∇ y ⁡ ϕ ρ ⁢ ( x , y ) ∧ n ⁢ ( y ) ) .

The operators rot Γ i := n i ⋅ rot and rot Γ j := - n j ∧ ∇ are surface derivatives. Thanks to the integration by parts formulae, the spherical parameterization of the hypersingular operator is obtained by combining V i ⁢ j ρ and the equalities

(2.47) | x j ′ | ⁢ ( rot Γ j ⁡ u ) ∘ x j = T j ⁢ rot Γ ⁡ ( u ∘ x j ) , | x i ′ | ⁢ ( rot Γ i ⁢ u ) ∘ x i = rot Γ ⁢ ( T i * ⁢ ( u ∘ x i ) ) ,

with

T j = t θ ⊗ e θ + t ϕ ⊗ e ϕ , T i * = e θ ⊗ t θ + e ϕ ⊗ t ϕ , e θ = ∂ ⁡ ξ ∂ ⁡ θ , e ϕ = 1 sin ⁡ θ ⁢ ∂ ⁡ ξ ∂ ⁡ ϕ ,

that are a tensor version of the parameterized identities derived in [55, equations (2.5.177) and (2.5.206)] in a local coordinate system. The practical formulae (2.47) enable us to get spectrally accurate discrete versions of boundary integral operators with higher order singularities as it was already demonstrated in [44, 45] for the elastodynamic and electromagnetic boundary integral operators. Indeed, being applied to the spherical harmonics functions, the surface derivatives are computed analytically on the unit sphere Γ. However, it requires the computational implementation of the spectrally accurate algorithm developed for solving vector boundary integral equations with tangential unknown fields presented in [27].

2.6 Numerical Examples

In this section, we illustrate the performance of the fully discrete formulation introduced in the previous section. We will observe superalgebraic convergence order in terms of the degree 𝑛 of the considered spherical harmonics space S n and compare the magnitude of errors in terms of the excitation frequency 𝜔. As is well known for objects with the same nature, we will check that, when dealing with simultaneous presence of objects of different nature, for achieving a given error level, the higher the frequency, the larger 𝑛.

We consider first a simple geometrical configuration where all objects have the same shape and size. For the exterior medium, we set c = 1 , which means that the associated exterior wavenumber for a given frequency 𝜔 is λ = ω . Figure 4 shows the location of the considered objects, all of them being ellipsoids of semi-axes of length 0.15, 0.25 and 0.3. Different colors have been used to identify the nature of the objects: green, yellow, magenta and blue correspond to penetrable, sound-soft, sound-hard and absorbing objects, respectively. The physical parameters for the object Ω 1 in light-green color are c 1 = 1 2 (therefore, μ 1 = 2 ⁢ ω for a given frequency 𝜔) and ν 1 = 1 2 , while for the other penetrable object Ω 2 , in dark-green color, the parameters are c 2 = 2 (i.e. μ 2 = ω 2 ) and ν 2 = 1 2 . Blue objects correspond to absorbing scatterers with impedance parameter α 7 = 1 + 3 ⁢ ı for the dark-blue object Ω 7 and α 8 = 3 + ı for the light-blue one ( Ω 8 ).

Figure 4

Geometrical configuration with eight ellipsoids of the same size. The parameter values for the penetrable objects are ( c 1 , ν 1 ) = ( 1 2 , 1 ) (light-green object) and ( c 2 , ν 2 ) = ( 2 , 1 2 ) (dark-green object). Objects in yellow ( Ω 3 and Ω 4 ) and magenta ( Ω 5 and Ω 6 ) correspond to sound-soft and sound-hard scatterers, respectively. The impedance parameters for the absorbing objects are α 7 = 1 + 3 ⁢ ı (light-blue object) and α 8 = 3 + ı (dark-blue object).

(a)

3D view

(b)

Section by the plane x = 0

(c)

Section by the plane z = - 1

We have checked the performance of the method for different parameter values 𝑛 for the space S n of spherical harmonics of degree 𝑛 and for varying values of the frequency 𝜔 ranging from 1 to 25. For this purpose, we have considered the function

u 0 ⁢ ( ⋅ ) := ϕ λ ⁢ ( ⋅ , x 0 ) ,

defined in terms of the fundamental solution (2.9) that models a source excitation from the point x 0 . We locate the source at x 0 = ( 0.05 , - 1.02 , 1.53 ) ⊤ , which is inside one of the Dirichlet objects, close to its center. Then it is rather easy to check that the function 𝑢, defined as

(2.48) u := { 0 in ⁢ Ω 1 ∪ Ω 2 , u 0 in ⁢ R 3 ∖ ⋃ i = 1 8 Ω i ¯ ,

is the unique solution to following problem:

{ Δ ⁢ u + μ i 2 ⁢ u = 0 in ⁢ Ω i , i ∈ { 1 , 2 } , Δ ⁢ u + λ 2 ⁢ u = 0 in ⁢ R 3 ∖ ⋃ i = 1 8 Ω i ¯ , u | Γ i int - u | Γ i ext = - u 0 on ⁢ Γ i , i ∈ { 1 , 2 } , ν i ⁢ ∂ n i ⁡ u | Γ i int - ∂ n i ⁡ u | Γ i ext = - ∂ n i ⁡ u 0 | Γ i on ⁢ Γ i , i ∈ { 1 , 2 } , u | Γ i ext = u 0 on ⁢ Γ i , i ∈ { 3 , 4 } , ∂ n i ⁡ u | Γ i ext = ∂ n i ⁡ u 0 | Γ i on ⁢ Γ i , i ∈ { 5 , 6 } , ∂ n i ⁡ u | Γ i ext + ı ⁢ α i ⁢ λ ⁢ u | Γ i ext = ∂ n i ⁡ u 0 | Γ i + ı ⁢ α i ⁢ λ ⁢ u 0 | Γ i on ⁢ Γ i , i ∈ { 7 , 8 } , lim r → ∞ ⁡ r ⁢ ( ∂ r ⁡ u - ı ⁢ λ ⁢ u ) = 0 .

On account of the previous sections, it is clear that we can obtain the function 𝑢 as the ansatz

u = { S i μ i ⁢ φ i in ⁢ Ω i , i ∈ { 1 , 2 } , ∑ i = 1 8 S i λ ⁢ ψ i in ⁢ R 3 ∖ ⋃ i = 1 8 Ω i ¯ ,

where the densities vector Ψ = [ φ , ψ T , ψ D , ψ N , ψ I ] ⊤ with φ , ψ T , ψ D , ψ N , ψ I defined as in (2.18) solves

M ⁢ Ψ = [ - u 0 , T D , - u 0 , T N , u 0 , D D , u 0 , N N , u 0 , I N + Au 0 , I D ] ⊤ ,

where vectors u 0 , T D , u 0 , T N , u 0 , D D , u 0 , N N , u 0 , I N and u 0 , I D are defined as in (2.19) with u 0 playing the role of u inc in (2.11), while the operators 𝐀 and 𝐌 are defined in (2.20) and (2.21), respectively.

After full discretization as explained in Section 2.5, we compare the numerically approximated function u n with the closed-form expression 𝑢 in (2.48) at a given set of points. Table 1 gathers the relative root mean square error

e n := ∑ k = 1 400 | u ⁢ ( y k ) - u n ⁢ ( y k ) | 2 ∑ k = 1 400 | u ⁢ ( y k ) | 2 ,

where the observation points { y 1 , … , y 100 } are randomly located on the sphere of radius 0.1 and centered at the center of the penetrable object Ω 1 (therefore, they are located inside Ω 1 ), { y 101 , … , y 200 } are the counterpart points on the sphere of radius 0.1 for the penetrable object Ω 2 , { y 201 , … , y 300 } are random points located on the sphere of radius 3 centered at the origin and { y 301 , … , y 400 } are random points located on the sphere of radius 6 and centered at the origin (therefore, the last 200 points are exterior to all the objects).

Table 1

Relative RMS errors e n for several values of 𝑛 and for different values of 𝜔 for the geometrical configuration in Figure 4.

𝑛	ω = 1	ω = 5	ω = 10	ω = 15	ω = 20	ω = 25
5	2.31 (−04)	6.44 (−04)	2.64 (−03)	1.57 (−02)	1.19 (−01)	1.40 (−00)
10	5.27 (−07)	2.06 (−06)	7.13 (−06)	2.77 (−05)	7.28 (−05)	2.75 (−04)
15	9.88 (−09)	2.16 (−08)	5.90 (−08)	2.69 (−07)	4.92 (−07)	1.81 (−06)
20	3.76 (−14)	9.09 (−11)	3.02 (−10)	1.41 (−09)	2.80 (−09)	9.77 (−09)

Data in Table 1 corroborate the theoretically predicted superalgebraic order of convergence (notice that, from (2.45), we expect superalgebraic order in the RMS norm for C ∞ surfaces). It is also observed that, as expected, it is necessary to slightly increase the value of 𝑛 if the frequency is increased. It is also noticeable that, by simply selecting n = 10 , we can guarantee a relative error of less than one thousandth for all the considered frequencies.

Figure 8

Geometrical configuration with six objects of different shapes and sizes. The parameter values for the penetrable objects are ( c 1 , ν 1 ) = ( 1 2 , 1 ) (light-green object) and ( c 2 , ν 2 ) = ( 2 , 1 2 ) (dark-green object). Objects in yellow ( Ω 3 ) and magenta ( Ω 4 ) correspond to sound-soft and sound-hard scatterers, respectively. The impedance parameters for the absorbing objects are α 5 = 1 + 3 ⁢ ı (light-blue object) and α 6 = 3 + ı (dark-blue object).

(a)

3D view

(b)

Section by the plane x = 0

(c)

Section by the plane y = 0

To end this section, we repeat the experiment for the configuration in Figure 8, where objects of different shapes and sizes are considered, and the minimum distance among objects has been reduced. The color code in Figure 8 is the same as in Figure 4, where green, yellow, magenta and blue stand for penetrable, sound-soft, sound-hard and absorbing objects, respectively. The parameters for the blue and green objects are those for the objects represented with the same colors in Figure 4. The source point is located now at x 0 = ( 0.02 , - 1.03 , 0.71 ) ⊤ , almost at the center of the absorbing object with a shape that resembles a bean. Observation points are selected as before, one hundred of them inside each penetrable object and the remaining two hundred ones are located in random positions at spheres of radius 3 and 6. Despite having different sizes and more complicated shapes, and objects being closer to each other than in the previous example, superalgebraic convergence order can be observed in Table 2. However, the magnitude of errors is now bigger and convergence is slower (compare with Table 1). We remark that, although it is not shown here for brevity, we have performed the counterpart examples when locating the source point close to the center of the sound-soft object (in yellow in Figure 8), the sound-hard (magenta ball) and the absorbing ellipsoid in dark blue. We have also kept the same geometry but interchanged the physical nature of the objects. Our conclusion is that, in all cases, we achieve superalgebraic order. Nevertheless, depending on the location of the source point, working precision is achieved for different values of 𝑛. When the source is located inside the ellipsoid, independently of its nature, errors are almost the same as in Table 1. For the ball, errors are at least one order of magnitude smaller and convergence is much faster, while for the rounded tetrahedron, errors are similar to those in Table 2, being slightly smaller. Errors in Table 2 present therefore the worst test case scenario for this geometrical configuration.

Table 2

Relative RMS errors e n for several values of 𝑛 and for different values of 𝜔 for the geometrical configuration in Figure 8.

𝑛	ω = 1	ω = 5	ω = 10	ω = 15	ω = 20	ω = 25
5	2.95 (−02)	8.48 (−02)	1.15 (−00)	1.18 (−00)	4.72 (−00)	3.75 (−00)
10	4.36 (−03)	8.33 (−03)	7.34 (−02)	1.24 (−01)	2.54 (−01)	9.16 (−01)
15	4.32 (−04)	6.39 (−04)	2.07 (−03)	1.83 (−03)	5.60 (−03)	1.19 (−01)
20	5.93 (−05)	9.62 (−05)	2.42 (−04)	2.61 (−04)	2.67 (−04)	1.33 (−02)
25	6.36 (−06)	9.89 (−06)	2.17 (−05)	1.80 (−05)	2.47 (−05)	1.42 (−04)

Results in Tables 1 and 2 are in accordance with those presented in [26], where the interested reader can find an extensive gallery of simulations where one sound-soft, sound-hard or penetrable object of varying shape and size is immersed in R 3 . For the performance in case of having several objects of different shapes but the same physical properties at high frequency regimes, we refer to [28]. In case of dealing with smooth but not C ∞ scatterers, the convergence rate would be slower. For a related study in the two-dimensional case, where the performance of a spectral discretization and a piecewise polynomial approximation are compared for C ∞ , C 3 and C 2 objects, we refer to [65, 67].

3 Inverse Problem

This section is devoted to the study of an inverse problem related to the forward problem considered in the previous section. The aim is to recover the shape of an unknown set of objects of different nature from near-field data measured at a set of observation points after actively exciting the medium by incident planar waves at different frequencies.

3.1 Description of the Problem

Assume that the medium is excited via a known planar incident wave of the form u inc ⁢ ( x ) = exp ⁡ ( ı ⁢ λ ⁢ d ⋅ x ) , where 𝜆 is the exterior wavenumber (recall that λ = ω c , where 𝑐 is the wave speed and 𝜔 the excitation frequency) and 𝐝 is the incident direction satisfying | d | = 1 . Measurements of the total wave u meas , 1 , … , u meas , N obs are taken at a set of N obs observation points y 1 , … , y N obs ∈ R 3 .

We are interested in finding the number of objects, N scat , and the regions Ω 1 , … , Ω N scat occupied by the scatterers without knowing their nature or any a priori information about them (physical parameters, location, approximated size, shape, …). Our aim is therefore to find the set of objects Ω := ⋃ i = 1 N scat Ω i such that the solution u Ω to forward problem (2.1)–(2.8) (with p inc = 0 ) coincides at the observation points with our measurements,

(3.1) u Ω ⁢ ( y j ) = u meas , j , j = 1 , … , N obs ,

without knowing either the numbers N T , N D , N N and N I of penetrable, sound-soft, sound-hard and absorbing objects, nor the interior wavenumbers μ i and transmission parameters ν i , i = 1 , … , N T , and impedance parameters α i for i = 1 , … , N I . We want to emphasize here that this paper is devoted to recover the scatterers. The problem of recovering not only the shape of the objects but also identifying their nature and the physical parameters μ i , ν i and α i is more involved and will be the subject of future research.

To alleviate the intrinsic ill-posedness of the problem, as commonly done in the literature, we will impose conditions (3.1) in a weaker form, replacing the original inverse problem by a minimization problem where the objective is to find the objects Ω minimizing the least squares functional

(3.2) J ⁢ ( R 3 ∖ Ω ¯ ) := 1 2 ⁢ ∑ j = 1 N obs | u Ω ⁢ ( y j ) - u meas , j | 2 .

When measurements for N inc incident waves for incident unitary directions d 1 , … , d N inc , produced for a set of N freq different frequencies ω 1 , … , ω N freq are available, we propose to consider a weighted least squares functional of the form

(3.3) J ⁢ ( R 3 ∖ Ω ¯ ) := 1 2 ⁢ ∑ k = 1 N freq r k ⁢ ( ∑ ℓ = 1 N inc ∑ j = 1 N obs | u Ω k , ℓ ⁢ ( y j ) - u meas , j k , ℓ | 2 ) ,

where u meas , j k , ℓ and u Ω k , ℓ are the measured total wave and the solution to forward problem (2.1)–(2.8) with objects Ω for an incident wave at the frequency ω k and incident direction d ℓ . Different weights r k will be considered for different frequencies ω k to balance the magnitude of the measurements depending on the considered frequency. This will be clarified in the next section.

3.2 Topological Derivative and Topological Energy-Based Methods

As already mentioned in the introduction, there is an increasing interest in the development of efficient and fast non-iterative methods for solving shape identification problems. One of them is based on the computation of the topological gradient of the involved shape functional.

The topological derivative or topological gradient [23, 56, 70] of a shape functional 𝐽 defined on a region ℛ is a scalar function that measures the sensitivity of 𝐽 to produce an infinitesimal perturbation at each point x ∈ R . In our setting, where no initial guess is available, we should set R = R 3 . When the considered infinitesimal perturbations are balls of radius ε > 0 , it provides an asymptotic expansion of the form

J ⁢ ( R ∖ B ε ⁢ ( x ) ¯ ) = J ⁢ ( R ) + f ⁢ ( ε ) ⁢ G ⁢ ( x ) + o ⁢ ( f ⁢ ( ε ) ) as ⁢ ε → 0 ,

where 𝐺 is the topological gradient at the point x ∈ R and the positive decreasing function 𝑓 satisfies f ⁢ ( ε ) → 0 as ε → 0 . The topological gradient 𝐺 can be then used as a damage indicator function because, since 𝑓 is positive, we observe that, at the points where 𝐺 attains the largest negative values, we have J ⁢ ( R ∖ B ε ⁢ ( x ) ¯ ) < J ⁢ ( R ) when ε → 0 , namely, the cost functional is expected to decrease. The method consists then in identifying the points where the largest negative values of 𝐺 are attained as belonging to the objects. In consequence, we approximate the true set of defects by the set

(3.4) Ω 0 := { x ∈ R 3 , G ⁢ ( x ) < ( 1 - C ) ⁢ min y ∈ R 3 ⁡ G ⁢ ( y ) } ,

where 0 < C < 1 is a constant threshold that in principle could be calibrated (notice that the closer to 1, the bigger Ω 0 ).

The topological gradient 𝐺 and the function 𝑓 depend on the boundary conditions (BC) on the boundary of the ball B ε ⁢ ( x ) . Closed-form expressions for these functions are well known for cost functionals of the form (3.2) when all the objects are assumed to be of the same nature. For the pure transmission conditions ( I D = I N = I I = ∅ ) with μ i = μ 0 and ν i = ν 0 for all i ∈ I T , Dirichlet ( I T = I N = I I = ∅ ), Neumann ( I T = I D = I I = ∅ ) or impedance conditions ( I T = I D = I N = ∅ ) with α i = α 0 for all i ∈ I I , the functions 𝑓 and 𝐺 are given by the following expressions (see [24, 30, 48, 69] or the recent paper [47]):

f ⁢ ( ε ) = { 4 3 ⁢ π ⁢ ε 3 for transmission or Neumann BC , 4 ⁢ π ⁢ ε for Dirichlet BC , 4 ⁢ π ⁢ ε 2 for impedance BC ,

and

(3.5) G ⁢ ( x ) = { Re ⁡ ( ( ν 0 ⁢ μ 0 2 - λ 2 ) ⁢ u ∅ ⁢ ( x ) ⁢ p ∅ ⁢ ( x ) + 3 ⁢ ( 1 - ν 0 ) 2 + ν 0 ⁢ ∇ ⁡ u ∅ ⁢ ( x ) ⋅ ∇ ⁡ p ∅ ⁢ ( x ) ) for transmission BC , Re ⁡ ( - λ 2 ⁢ u ∅ ⁢ ( x ) ⁢ p ∅ ⁢ ( x ) + 3 2 ⁢ ∇ ⁡ u ∅ ⁢ ( x ) ⋅ ∇ ⁡ p ∅ ⁢ ( x ) ) for Neumann BC , Re ⁡ ( - u ∅ ⁢ ( x ) ⁢ p ∅ ⁢ ( x ) ) for Dirichlet BC , Re ⁡ ( ı ⁢ α 0 ⁢ λ ⁢ u ∅ ⁢ ( x ) ⁢ p ∅ ⁢ ( x ) ) for impedance BC ,

where u ∅ and p ¯ ∅ solve direct and adjoint problems in the background medium without defects (namely, with Ω = ∅ and therefore I T = I D = I N = I I = ∅ ) and 𝐺 is defined for all x ∈ R 3 ∖ { y 1 , … , y N obs } . Specifically, u ∅ is nothing but the incident wave and p ∅ can be computed in closed form using the fundamental solution ϕ λ defined in (2.9),

(3.6) u ∅ ⁢ ( x ) = exp ⁡ ( ı ⁢ λ ⁢ d ⋅ x ) , p ∅ ⁢ ( x ) = ∑ j = 1 N obs ϕ λ ⁢ ( x , y j ) ⁢ ( u ∅ ⁢ ( y j ) - u meas , j ) ¯ .

We observe that, despite the dependence of 𝑓 and 𝐺 on the boundary conditions, these formulae involve the same functions u ∅ and p ∅ . What is more, in all cases, the evaluation of 𝐺 is extremely cheap and fast from the computational point of view as closed-form expressions for u ∅ , p ∅ , ∇ ⁡ u ∅ and ∇ ⁡ p ∅ are known. As a result, obtaining the set Ω 0 defined by (3.4) is very simple because one only needs to evaluate 𝐺 at a grid of points covering the region where objects are searched.

Nevertheless, in our setting, it is not only the number of objects that is unknown, it is also possible that there are objects of different nature in the same media. Therefore, it is not clear which function 𝐺 or which combination of those functions could work. Even more, there could be several penetrable objects with different unknown parameters μ i and ν i and several absorbing objects corresponding to different unknown impedance parameters α i , in which case it is not clear how to extend the formula given in (3.5).

A way to overcome all these difficulties is to replace the topological gradient by a closely related indicator function: the topological energy. This function was firstly introduced in [19] for object detection from time-dependent ultrasound measurements and used in [10] for the imaging of micro and nanostructures in a 3D holography problem. Its application to process experimental data obtained in laboratories has produced very satisfactory results, both in ultrasound [18] and in microwave [11] imaging. In the mentioned papers, all objects were of the same nature, but the topological energy was used for being more suitable to overcome oscillations than the topological derivative. In [64], the method was tested against synthetic data for multiple penetrable, sound-soft and sound-hard objects in a 2D setting. Here we will extend this approach to the three-dimensional case including absorbing objects as well. We will use it as a qualitative imaging tool. The rigorous mathematical analysis of this method is out of the scope of the present paper. We refer to [3, 5, 6, 31, 51, 59, 61] for the analysis of the topological derivative method in some specific cases that might be extensible to topological energy methods. In case of time-dependent measurements (which is not the subject of the present paper), the topological energy indicator function can be studied as a time-reversal method, as done in [19].

For our model problem, the topological energy associated to the cost functional (3.2) is given by the function

E ⁢ ( x ) = - | u ∅ ⁢ ( x ) | 2 ⁢ | p ∅ ⁢ ( x ) | 2 for all ⁢ x ∈ R 3 ∖ { y 1 , … , y N obs } ,

where the functions u ∅ and p ∅ are given in (3.6). In analogy to the topological derivative case, the points at which the function 𝐸 attains its largest negative values will be used to characterize the searched scatterers by considering the approximate region

(3.7) Ω 0 := { x ∈ R 3 , E ⁢ ( x ) < ( 1 - C ) ⁢ min y ∈ R 3 ⁡ E ⁢ ( y ) } ,

where 0 < C < 1 is a constant threshold. We remark that, in the original paper [19], the topological energy is defined with the opposite sign, and in consequence, objects are described from the points where the largest positive values are attained.

When measurements u meas , j k , ℓ of the total wave for different incident directions d ℓ , ℓ = 1 , … , N inc , and different frequencies ω k , k = 1 , … , N freq , are available at the observation points y j , j = 1 , … , N obs , we consider the weighted cost functional defined in (3.3). Its topological energy is given by the function

(3.8) E ⁢ ( x ) = ∑ k = 1 N freq r k ⁢ E k ⁢ ( x ) , x ∈ R 3 ∖ { y 1 , … , y N obs } ,

defined in terms of the monofrequency topological energy functions

(3.9) E k ⁢ ( x ) = - ∑ ℓ = 1 N inc | u ∅ k , ℓ ⁢ ( x ) | 2 ⁢ | p ∅ k , ℓ ⁢ ( x ) | 2 ,

with

u ∅ k , ℓ ⁢ ( x ) = exp ⁡ ( ı ⁢ λ k ⁢ d ℓ ⋅ x ) , p ∅ k , ℓ ⁢ ( x ) = ∑ j = 1 N obs ϕ λ k ⁢ ( x , y j ) ⁢ ( u ∅ k , ℓ ⁢ ( y j ) - u meas , j k , ℓ ) ¯ ,

where λ k = ω k c , 𝑐 being the wave velocity in the pristine medium.

The weights r k that appear in functional (3.3) are intended to weight the contribution of each frequency since each individual topological energy E k might have a different order of magnitude. As in [11, 64], we propose to select r k after evaluating E k as follows. First, in order to numerically obtain the domain Ω 0 defined in (3.7), we need to consider a bounded region R test ⊂ R 3 ∖ { y 1 , … , y N obs } where the scatterers are assumed to be located (in our examples, it will be a cube), and a discrete set of points R test N ⊂ R test consisting of 𝑁 points where the topological energy 𝐸 will be evaluated. Then our approximated objects are numerically computed as the set

(3.10) Ω 0 := { x ∈ R test N , E ⁢ ( x ) < ( 1 - C ) ⁢ min y ∈ R test N ⁡ E ⁢ ( y ) } ,

where 𝐸 is defined in (3.8) with weights

(3.11) r k = - min y ∈ R test N ⁡ E k ⁢ ( y ) ,

obtained after evaluation of the 𝑘-th frequency indicator function E k introduced in (3.9). This choice guarantees that each summand in (3.8) satisfies min x ∈ R test N ⁡ r k ⁢ E k ⁢ ( x ) = - 1 , making the contribution of each frequency somehow equally important. This kind of weight was first suggested in [58] for weighting topological derivative data in an electromagnetic setting and afterwards successfully used for different applications related with acoustic problems [48, 47], thermal imaging [32, 60] and structural health monitoring [25, 52].

We will illustrate in the next section the performance of this indicator function. We want to point out here that some other choices for r k could be considered. Notice also that Ω 0 depends on the selection of the discrete points, R test N , in three different places: (i) the evaluation of the topological energy itself, (ii) the computation of r k and (iii) the computation of the minimum that appears in the threshold criteria. For some insights about the selection of optimal meshes for the evaluation of indicator functions, we refer to [2]. Finally, we remark that the selection of 𝐶 will be heuristic, as it is commonly done in most of the papers devoted to direct sampling methods (topological derivative-based, linear sampling, MUSIC algorithms, etc.). Performing a calibration for some tests with known geometries as done in [14] for the selection of the threshold constant for the linear sampling method could be explored. However, that paper highlights that optimal constants depend on the working frequency and on the nature of the scatterers, which means that calibration is expected to become more involved for our problem where simultaneous presence of objects of different nature will occur. In the forthcoming section, we will just fix a test region R test N , define the weights as in (3.11) and show results for representative values of 𝐶 to illustrate the robustness of the method.

3.3 Numerical Examples

Let us explore now the performance of the topological energy-based method for the two configurations already introduced in Section 2.6. In the first test case, all objects have identical shape and size but different physical properties; see Figure 4. We chose this example to illustrate the behavior of the method in terms of the nature of the objects. Secondly, the configuration in Figure 8 is selected to visualize the ability of the method to recover more complicated shapes in more demanding situations where objects are closer to each other and have different sizes. We emphasize that our aim is to recover shapes, not to identify the nature of the objects or the actual physical parameters. For both configurations, the test region is R test = [ - 3 , 3 ] 3 .

Figure 12

(a) Test region R test = [ - 3 , 3 ] 3 (gray cube), incidence directions (arrows) and observation points (blue crosses). (b) True ellipsoids for the configuration with eight objects with the same size and shape. (c) True objects in the configuration with six objects of different size and shape.

(a)

(b)

(c)

Synthetically generated measured data for N inc = 25 incident waves corresponding to 25 random incident directions on the unitary sphere have been created by solving the corresponding forward problem by using the fully discrete boundary integral method described in previous sections. The degree of spherical harmonics was selected as n = 12 for the configuration with eight ellipsoids and n = 20 for the configuration with six objects of different sizes. According to the experiments performed in Section 2.6, these choices should guarantee a relative error below 5 % in both cases for all the frequencies 1 ≤ ω ≤ 25 . In addition, to simulate experimental measurement errors, a random relative noise of a 5 % size was added to the synthetic data. Measurements were computed at N obs = 80 observation points randomly located on the sphere of radius 6. Figure 12 resumes schematically the setups.

First we will consider monofrequency data, namely, for a given frequency ω k , we consider the indicator function E k defined in (3.9). Afterwards, data corresponding to several frequencies will be simultaneously processed via the indicator function 𝐸 defined in (3.8).

Figure 19

Topological energy for different frequencies 𝜔 at the planes x = 0 and z = - 1 . First column: ω = 1 . Second column: ω = 5 . Third column: ω = 10 .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 26

Topological energy for different frequencies 𝜔 at the planes x = 0 and z = - 1 . First column: ω = 15 . Second column: ω = 20 . Third column: ω = 25 .

(a)

(b)

(c)

(d)

(e)

(f)

Let us start with the example where true objects are eight ellipsoids. Figure 19 shows the color maps corresponding to the evaluation of the topological energy function for the frequencies ω = 1 , 5 and 10. To facilitate the visualization of results, we have selected the x = 0 and z = - 1 planes because four out of the eight ellipsoids are intersected by the plane x = 0 and the remaining four are intersected by the plane z = - 1 (Figure 4). We observe that, depending on the selected frequency, results are rather different and complementary. We observe that, for the lowest frequency ω = 1 , the largest negative values are attained in the region occupied by the two sound-soft objects and only one of the absorbing objects is vaguely identified (compare blue regions in Figure 19 (a), (d) with Figure 4 (b), (c)). The remaining objects are completely disregarded even though all have exactly the same shape and size. For ω = 5 , the behavior is different: not only the two sound-soft objects are identified, but also the largest negative values concentrate now inside a penetrable object and in a sound-soft one. A further increase to ω = 10 (see panels (c), (f)) promotes a more blurred color map where the largest negative values concentrate now both inside some of the objects, on the boundaries of some others and in some spurious small regions.

Figure 26 is the counterpart of Figure 19 for the frequencies ω = 15 , 20 and 25. The evaluation of the topological energy indicator function at high frequencies produces color maps that are more difficult to interpret at plain sight. Now, the number of oscillations and of spurious regions increase. However, they seem to contain relevant information about details concerning the shape of the objects. In particular, high frequencies are precisely the ones able to identify the location of some objects missed by the lower ones (for instance, sound-hard objects are clearly recovered when ω = 15 ). They also highlight some small regions around the boundaries of the objects.

In view of the results in Figures 19 and 26, it is clear that different frequencies entail complementary information, so we can anticipate that their combined use will provide better results than their individual usages. It is also evident that the order of magnitude of the minimum value of the topological energy function depends on the selected frequency. This led to the need to introduce the weights r k in the shape functional (3.3), and their subsequent definition according to formula (3.11). Owing to the use of these weights, we are somehow making all the identified points (blue regions the figures) equally important in the weighted indicator function 𝐸. The values of 𝐸 at the planes x = 0 and z = - 1 corresponding to weighting the monofrequency results for ω = 1 , … , 25 are shown in Figure 29. It is remarkable that, combining the monofrequency indicator functions, all objects seem to be distinguishable, although depending on their nature, some of them are more clearly identified. We also observe that this combination produces less spurious regions and smoother plots.

Figure 29

Weighted topological energy for ω = 1 , … , 25 at the planes x = 0 and z = - 1 .

(a)

(b)

Figure 33

Approximated domain Ω 0 defined in (3.10) for different values of 𝐶.

(a)

C = 0.3

(b)

C = 0.4

(c)

C = 0.5

The previous figures only show the evaluation of the topological energy functions at two planes and produce qualitative plots, in which the points where the largest negative values are attained (regions in blue colors) are identified with damaged regions. For quantitative reconstructions in 3D, we compute the set Ω 0 defined in (3.10) from the evaluation of the topological energy functions at a discrete set of points consisting in a grid of 50 × 50 × 50 points uniformly distributed in the test region R test = [ - 3 , 3 ] 3 . To do that, we identify first the connected components of Ω 0 by using the Matlab function bwconncomp and approximate each boundary with the help of the function boundary. For further details on the practical implementation of Ω 0 in Matlab, we refer to [46, Section 4.2].

As can be extrapolated in view of Figure 29, depending on the selected constant 𝐶, we could obtain approximated regions where not all objects are identified. Our reconstructed objects for the cut-off values C = 0.3 , 0.4 and 0.5 are shown in Figure 33. By definition, the larger 𝐶, the bigger Ω 0 . This promotes that, for low values of 𝐶, some objects are missed, while when all objects are detected, some of them are overestimated in size and some others are underestimated. However, the location is sharp in all cases, and the approximations are rather robust in terms of 𝐶, which means that 𝐶 should be only roughly calibrated. It is important to remark that reconstructions are not perfect, but we are not using any a priori information about the shape, size, number or physical properties of the objects, which makes therefore this problem a highly demanding one.

Figure 37

(a) Weighted topological energy for ω = 1 , … , 15 at the plane x = 0 . (b) Weighted topological energy for ω = 1 , … , 15 at the plane z = - 1 . (c) Approximated domain Ω 0 defined in (3.10) for C = 0.6 .

(a)

(b)

(c)

In view of Figure 29, it might be thought that the highest frequencies are less useful than the lowest ones. However, if we just consider the frequencies in the range 1 ≤ ω ≤ 15 , then the identification of all the objects becomes more complicated, and even increasing the value of the constant 𝐶, some objects are missed, while some spurious regions appear. This can be observed in Figure 37. Although it is not shown here for brevity, we have checked that, when considering frequencies in the range 1 ≤ ω ≤ 20 , reconstructions are similar to those for ω ∈ [ 1 , 25 ] or even slightly better. However, since we do not know a priori the most suitable set of frequencies, we recommend to use an ample set since spurious regions are somehow compensated and high frequencies enrich details.

Let us examine now the performance of the method for a more complicated setting where, as anticipated at the beginning of this section, six objects of different nature and different sizes and shapes are considered (Figure 8). In this case, three of the scatterers are intersected by the plane x = 0 and the remaining ones are intersected by the plane y = 0 . Figures 44 and 51 show the values of the monofrequency indicator functions E k for ω k = 1 , 5, 10, 15, 20 and 25. The behavior of the method is qualitatively as in the case when all objects were ellipsoids. In spite of having objects of different sizes and that the higher frequencies seem to contain very few information about the objects (Figure 51), when combining the 25 monofrequency indicator functions, we are able to satisfactorily capture the true geometrical configuration, as can be seen in Figures 54 and 58, which are the counterparts to Figures 29 and 33. It is remarkable now that not only the location and size but also shapes can be accurately recovered. Notice that three of the objects obtained for C = 0.3 sharply resemble their associated true objects, while only a small part of the bean-shaped object is devised. The increase of the constant to C = 0.4 permits the reconstruction of five out of the six true objects with reasonable accuracy in their shapes. A further increase to C = 0.5 allows us to identify the six true objects at the correct location but at the price of loosing some accuracy in the size and shape of some of the objects that were already identified for the lowest constants.

Figure 44

Topological energy for different frequencies 𝜔 at the planes x = 0 and y = 0 . First column: ω = 1 . Second column: ω = 5 . Third column: ω = 10 .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 51

Topological energy for different frequencies 𝜔 at the planes x = 0 and y = 0 . First column: ω = 15 . Second column: ω = 20 . Third column: ω = 25 .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 54

Weighted topological energy for ω = 1 , … , 25 at the planes x = 0 and y = 0 .

(a)

(b)

Figure 58

Approximated domain Ω 0 defined in (3.10) for different values of 𝐶.

(a)

C = 0.3

(b)

C = 0.4

(c)

C = 0.5

Figure 62

(a) Test region R test = [ - 3 , 3 ] 3 (gray cube), incidence directions (arrows) and observation points (blue crosses). (b), (c) Weighted topological energy for 15 random frequencies 1 ≤ ω ≤ 25 at the planes x = 0 and y = 0 .

(a)

(b)

(c)

Figure 66

Approximated domain Ω 0 defined in (3.10) for different values of 𝐶.

(a)

C = 0.3

(b)

C = 0.4

(c)

C = 0.5

To end this section, we want to show how the method behaves when drastically reducing the number and the quality of the available data. To do that, we have reduced the number of incident waves from N inc = 25 to N inc = 15 and the number of observation points from N obs = 80 to N obs = 50 . Figure 62 (a) shows the location of the observation points and the incidence directions considered in this situation (compare with Figure 12 (a)). In addition, the level of relative noise has been increased from a 5 % to a 15 %. Furthermore, only 15 random frequencies between 1 and 25 have been considered. Certainly, the method looses accuracy. However, we are able to obtain reasonable reconstructions even in such extreme case, as can be observed in Figures 62 (b), (c) and 66, where we appreciate the recovery of all but the penetrable ellipsoid, which was also the one that presented more difficulties to be identified in the more favorable previous case. Probably, this is due to the fact that it is smaller and surrounded by other bigger objects which do not allow the incident waves to strike so directly, nor to be directly scattered afterwards.

4 Conclusions

A boundary element method has been proposed and analyzed for solving direct three-dimensional multiple acoustic scattering problems in which a finite number of scatterers that respond in a different way to incident excitations (penetrable, sound-soft, sound-hard and absorbing scatterers) are present in the medium. Our approach makes use of single-layer potentials. Due to the fact that double-layer potentials are not involved, the resulting system of equations is simpler than those that would be obtained when using direct formulations or mixed potentials. However it might fail when either - λ 2 is a Dirichlet eigenvalue of the Laplacian inside any of the scatterers, 𝜆 being the exterior wavenumber, or when - μ i 2 is a Dirichlet eigenvalue for the Laplacian in the penetrable object Ω i , where μ i is the corresponding wavenumber in Ω i . Spurious eigenvalues or proximity to them would produce a huge increase in the condition number of the associated discrete matrix. This could be then avoided by using Brakhage–Werner potentials instead of single-layer potentials, but the price to pay is having to work with a more complicated system involving hypersingular operators.

We have given necessary and sufficient conditions for the convergence of abstract Petrov–Galerkin methods, showing that, when all the boundaries of the scatterers can be globally described in spherical coordinates and all densities are approximated in the space of spherical harmonics of degree 𝑛, the method converges superalgebraically. Our numerical experiments corroborate this convergence order and show that only a very reduced number 𝑛 is needed to achieve reasonably small relative errors over a fairly large frequency range. We want to emphasize that the theoretical convergence analysis is not limited to spectral methods, and convergence of other type of discretizations, i.e. based on piecewise polynomials, could be studied by verifying our equivalent conditions.

In the second part, we have used the proposed formulation to numerically generate measurements of the total wave at some observation points to be processed with the aim of recovering the true objects using no a priori information other than the knowledge of the measured data, the incident waves, the employed frequencies and the properties of the background medium. Namely, the number, size, shape, location and physical nature of the objects were unknown. Our method is based on the evaluation of suitable linear combinations of topological energy indicator functions associated to monofrequency data. In view of our experiments, we conclude that the ability of the monofrequency indicators to identify objects depends not only on the frequency itself but also on the size and the nature of each object. Processing multifrequency data overcomes this uncertainty, and the resultant method is shown to be very effective for shape reconstruction of multiple objects in different demanding situations. Although the numerical study focuses on the detection of a few disjoint homogeneous particles, we are confident that the topological energy indicator function can produce accurate initial guesses for recovering the structure of the more complex heterogeneous medium studied in [12, 33, 20]. In future, we plan to extend the method to the highly challenging problem of recovering both, the shape of the objects and their nature and/or the physical parameters involved (i.e. the interior wavenumber in penetrable objects and the coefficients in penetrable and impedance boundaries).

Dedicated to the memory of Francisco-Javier Sayas (1968–2019) who will be deeply missed as a scientist, and as a friend.

Funding source: Agencia Estatal de Investigación

Award Identifier / Grant number: PID2020-114173RB-I00

Funding statement: Grant PID2020-114173RB-I00 funded by MCIN/AEI/10.13039/501100011033.

Acknowledgements

The helpful suggestions of both referees on an earlier version of the paper are gratefully acknowledged.

References

[1] A. Albicker and R. Griesmaier, Monotonicity in inverse obstacle scattering on unbounded domains, Inverse Problems 36 (2020), no. 8, Article ID 085014. 10.1088/1361-6420/ab98a3Search in Google Scholar

[2] H. Ammari, Y. T. Chow and K. Liu, Optimal mesh size for inverse medium scattering problems, SIAM J. Numer. Anal. 58 (2020), no. 1, 733–756. 10.1137/18M122159XSearch in Google Scholar

[3] H. Ammari, J. Garnier, V. Jugnon and H. Kang, Stability and resolution analysis for a topological derivative based imaging functional, SIAM J. Control Optim. 50 (2012), no. 1, 48–76. 10.1137/100812501Search in Google Scholar

[4] K. E. Atkinson, The numerical solution of Laplace’s equation in three dimensions, SIAM J. Numer. Anal. 19 (1982), no. 2, 263–274. 10.1007/978-3-0348-6314-8_1Search in Google Scholar

[5] C. Bellis, M. Bonnet and F. Cakoni, Acoustic inverse scattering using topological derivative of far-field measurements-based L 2 cost functionals, Inverse Problems 29 (2013), no. 7, Article ID 075012. 10.1088/0266-5611/29/7/075012Search in Google Scholar

[6] M. Bonnet and F. Cakoni, Analysis of topological derivative as a tool for qualitative identification, Inverse Problems 35 (2019), no. 10, Article ID 104007. 10.1088/1361-6420/ab0b67Search in Google Scholar

[7] H. Brakhage and P. Werner, Über das Dirichletsche Aussenraumproblem für die Helmholtzsche Schwingungsgleichung, Arch. Math. 16 (1965), 325–329. 10.1007/BF01220037Search in Google Scholar

[8] A. J. Burton and G. F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary-value problems, Proc. Roy. Soc. Lond. Ser. A 323 (1971), 201–210. 10.1098/rspa.1971.0097Search in Google Scholar

[9] F. Cakoni and D. Colton, Qualitative Methods in Inverse Scattering Theory: An Introduction, Interact. Mech. Math., Springer, Berlin, 2006. Search in Google Scholar

[10] A. Carpio, T. G. Dimiduk, M.-L. Rapún and V. Selgas, Noninvasive imaging of three-dimensional micro and nanostructures by topological methods, SIAM J. Imaging Sci. 9 (2016), no. 3, 1324–1354. 10.1137/16M1068669Search in Google Scholar

[11] A. Carpio, M. Pena and M.-L. Rapún, Processing the 2D and 3D Fresnel experimental databases via topological derivative methods, Inverse Problems 37 (2021), no. 10, Article ID 105012. 10.1088/1361-6420/ac21c8Search in Google Scholar

[12] B. Caudron, X. Antoine and C. Geuzaine, Optimized weak coupling of boundary element and finite element methods for acoustic scattering, J. Comput. Phys. 421 (2020), Article ID 109737. 10.1016/j.jcp.2020.109737Search in Google Scholar

[13] G. Chen and J. Zhou, Boundary Element Methods. Computational Mathematics and Applications, Academic Press, London, 1992. Search in Google Scholar

[14] D. Colton, K. Giebermann and P. Monk, A regularized sampling method for solving three-dimensional inverse scattering problems, SIAM J. Sci. Comput. 21 (2000), no. 6, 2316–2330. 10.1137/S1064827598340159Search in Google Scholar

[15] D. Colton and A. Kirsch, A simple method for solving inverse scattering problems in the resonance region, Inverse Problems 12 (1996), no. 4, 383–393. 10.1088/0266-5611/12/4/003Search in Google Scholar

[16] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, 2nd ed., Appl. Math. Sci. 93, Springer, Berlin, 1998. 10.1007/978-3-662-03537-5Search in Google Scholar

[17] M. Costabel and E. Stephan, A direct boundary integral equation method for transmission problems, J. Math. Anal. Appl. 106 (1985), no. 2, 367–413. 10.1016/0022-247X(85)90118-0Search in Google Scholar

[18] N. Dominguez and V. Gibiat, Non-destructive imaging using the time domain topological energy method, Ultrasonics 50 (2010), 367–372. 10.1016/j.ultras.2009.08.014Search in Google Scholar PubMed

[19] N. Dominguez, V. Gibiat and Y. Esquerre, Time domain topological gradient and time reversal analogy: An inverse method for ultrasonic target detection, Wave Motion 42 (2005), no. 1, 31–52. 10.1016/j.wavemoti.2004.09.005Search in Google Scholar

[20] V. Domínguez, M. Ganesh and F. J. Sayas, An overlapping decomposition framework for wave propagation in heterogeneous and unbounded media: Formulation, analysis, algorithm, and simulation, J. Comput. Phys. 403 (2020), Article ID 109052. 10.1016/j.jcp.2019.109052Search in Google Scholar

[21] V. Domínguez, M. Lyon and C. Turc, Well-posed boundary integral equation formulations and Nyström discretizations for the solution of Helmholtz transmission problems in two-dimensional Lipschitz domains, J. Integral Equations Appl. 28 (2016), no. 3, 395–440. 10.1216/JIE-2016-28-3-395Search in Google Scholar

[22] V. Domínguez and F.-J. Sayas, Overlapped BEM-FEM for some Helmholtz transmission problems, Appl. Numer. Math. 57 (2007), no. 2, 131–146. 10.1016/j.apnum.2006.02.001Search in Google Scholar

[23] H. Eschenauer, V. Kobelev and A. Schumacher, Bubble method for topology and shape optimization of structures, Struct. Optim. 8 (1994), 42–51. 10.1007/BF01742933Search in Google Scholar

[24] G. R. Feijoo, A new method in inverse scattering based on the topological derivative, Inverse Problems 20 (2004), no. 6, 1819–1840. 10.1088/0266-5611/20/6/008Search in Google Scholar

[25] J. F. Funes, J. M. Perales, M.-L. Rapún and J. M. Vega, Defect detection from multi-frequency limited data via topological sensitivity, J. Math. Imaging Vision 55 (2016), no. 1, 19–35. 10.1007/s10851-015-0611-ySearch in Google Scholar

[26] M. Ganesh and I. G. Graham, A high-order algorithm for obstacle scattering in three dimensions, J. Comput. Phys. 198 (2004), no. 1, 211–242. 10.1016/j.jcp.2004.01.007Search in Google Scholar

[27] M. Ganesh and S. C. Hawkins, A spectrally accurate algorithm for electromagnetic scattering in three dimensions, Numer. Algorithms 43 (2006), no. 1, 25–60. 10.1007/s11075-006-9033-7Search in Google Scholar

[28] M. Ganesh and S. C. Hawkins, Simulation of acoustic scattering by multiple obstacles in three dimensions, ANZIAM J. 50 (2008), C31–C45. 10.21914/anziamj.v50i0.1451Search in Google Scholar

[29] I. G. Graham and I. H. Sloan, Fully discrete spectral boundary integral methods for Helmholtz problems on smooth closed surfaces in R 3 , Numer. Math. 92 (2002), no. 2, 289–323. 10.1007/s002110100343Search in Google Scholar

[30] B. B. Guzina and M. Bonnet, Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics, Inverse Problems 22 (2006), no. 5, 1761–1785. 10.1088/0266-5611/22/5/014Search in Google Scholar

[31] B. B. Guzina and F. Pourahmadian, Why the high-frequency inverse scattering by topological sensitivity may work, Proc. A. 471 (2015), no. 2179, Article ID 20150187. 10.1098/rspa.2015.0187Search in Google Scholar PubMed PubMed Central

[32] M. Higuera, J. M. Perales, M. L. Rapún and J. M. Vega, Solving inverse geometry heat conduction problems by postprocessing steady thermograms, Int. J. Heat Mass Transfer 149 (2019), Article ID 118490. 10.1016/j.ijheatmasstransfer.2019.118490Search in Google Scholar

[33] R. Hiptmair and C. Jerez-Hanckes, Multiple traces boundary integral formulation for Helmholtz transmission problems, Adv. Comput. Math. 37 (2012), no. 1, 39–91. 10.1007/s10444-011-9194-3Search in Google Scholar

[34] G. C. Hsiao and W. L. Wendland, Boundary Integral Equations, Appl. Math. Sci. 164, Springer, Berlin, 2008. 10.1007/978-3-540-68545-6Search in Google Scholar

[35] V. Isakov, Inverse Problems for Partial Differential Equations, 2nd ed., Appl. Math. Sci. 127, Springer, New York, 2006. Search in Google Scholar

[36] O. Ivanyshyn and R. Kress, Identification of sound-soft 3D obstacles from phaseless data, Inverse Probl. Imaging 4 (2010), no. 1, 131–149. 10.3934/ipi.2010.4.131Search in Google Scholar

[37] A. Kirsch, An Introduction to the Mathematical Theory of Inverse Problems, 3rd ed., Appl. Math. Sci. 120, Springer, Cham, 2021. 10.1007/978-3-030-63343-1Search in Google Scholar

[38] R. E. Kleinman and P. A. Martin, On single integral equations for the transmission problem of acoustics, SIAM J. Appl. Math. 48 (1988), no. 2, 307–325. 10.1137/0148016Search in Google Scholar

[39] R. Kress, Linear Integral Equations, 2nd ed., Springer, Berlin, 1996. Search in Google Scholar

[40] R. Kress and G. F. Roach, On mixed boundary value problems for the Helmholtz equation, Proc. Roy. Soc. Edinburgh Sect. A 77 (1977), no. 1–2, 65–77. 10.1017/S0308210500018047Search in Google Scholar

[41] R. Kress and G. F. Roach, Transmission problems for the Helmholtz equation, J. Math. Phys. 19 (1978), no. 6, 1433–1437. 10.1063/1.523808Search in Google Scholar

[42] A. R. Laliena, M.-L. Rapún and F.-J. Sayas, Symmetric boundary integral formulations for Helmholtz transmission problems, Appl. Numer. Math. 59 (2009), no. 11, 2814–2823. 10.1016/j.apnum.2008.12.030Search in Google Scholar

[43] A. R. Laliena and F.-J. Sayas, Theoretical aspects of the application of convolution quadrature to scattering of acoustic waves, Numer. Math. 112 (2009), no. 4, 637–678. 10.1007/s00211-009-0220-zSearch in Google Scholar

[44] F. Le Louër, A high order spectral algorithm for elastic obstacle scattering in three dimensions, J. Comput. Phys. 279 (2014), 1–17. 10.1016/j.jcp.2014.08.047Search in Google Scholar

[45] F. Le Louër, Spectrally accurate numerical solution of hypersingular boundary integral equations for three-dimensional electromagnetic wave scattering problems, J. Comput. Phys. 275 (2014), 662–666. 10.1016/j.jcp.2014.07.022Search in Google Scholar

[46] F. Le Louër and M.-L. Rapún, Topological sensitivity for solving inverse multiple scattering problems in three-dimensional electromagnetism. Part I: One step method, SIAM J. Imaging Sci. 10 (2017), no. 3, 1291–1321. 10.1137/17M1113850Search in Google Scholar

[47] F. Le Louër and M.-L. Rapún, Topological sensitivity analysis revisited for time-harmonic wave scattering problems. Part I: The free space case, Eng. Comput. 39 (2022), no. 1, 232–271. 10.1108/EC-06-2021-0327Search in Google Scholar

[48] F. Le Louër and M.-L. Rapún, Detection of multiple impedance obstacles by non-iterative topological gradient based methods, J. Comput. Phys. 388 (2019), 534–560. 10.1016/j.jcp.2019.03.023Search in Google Scholar

[49] R. Leis, Zur Dirichletschen Randwertaufgabe des Aussenraumes der Schwingungsgleichung, Math. Z. 90 (1965), 205–211. 10.1007/BF01119203Search in Google Scholar

[50] T.-C. Lin, On an integral equation approach for the exterior Robin problem for the Helmholtz equation, J. Math. Anal. Appl. 126 (1987), no. 2, 547–555. 10.1016/0022-247X(87)90061-8Search in Google Scholar

[51] Y. K. Ma, P. S. Kin and W. K. Park, Analysis of topological derivative function for a fast electromagnetic imaging of perfectly conducting cracks, Prog. Electro-Magn. Res. 122 (2012), 311–325. 10.2528/PIER11092901Search in Google Scholar

[52] A. Martínez, J. A. Güemes, J. M. Perales and J. M. Vega, SHM via topological derivative, Smart Mat. Struct. 27 (2018), Article ID 085002. 10.20868/UPM.thesis.65446Search in Google Scholar

[53] W. McLean, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge University, Cambridge, 2000. Search in Google Scholar

[54] J. L. Mueller and S. Siltanen, Linear and Nonlinear Inverse Problems with Practical Applications, Comput. Sci. Eng. 10, Society for Industrial and Applied Mathematics, Philadelphia, 2012. 10.1137/1.9781611972344Search in Google Scholar

[55] J.-C. Nédélec, Acoustic and Electromagnetic Equations. Integral Representations for Harmonic Problems, Appl. Math. Sci. 144, Springer, New York, 2001. 10.1007/978-1-4757-4393-7Search in Google Scholar

[56] A. A. Novotny and J. Sokoł owski, Topological Derivatives in Shape Optimization, Interact. Mech. Math., Springer, Heidelberg, 2013. 10.1007/978-3-642-35245-4Search in Google Scholar

[57] O. I. Panič, On the solubility of exterior boundary-value problems for the wave equation and for a system of Maxwell’s equations, Russian Math. Surveys 20 (1965), 221–226. Search in Google Scholar

[58] W.-K. Park, Multi-frequency topological derivative for approximate shape acquisition of curve-like thin electromagnetic inhomogeneities, J. Math. Anal. Appl. 404 (2013), no. 2, 501–518. 10.1016/j.jmaa.2013.03.040Search in Google Scholar

[59] W.-K. Park, Performance analysis of multi-frequency topological derivative for reconstructing perfectly conducting cracks, J. Comput. Phys. 335 (2017), 865–884. 10.1016/j.jcp.2017.02.007Search in Google Scholar

[60] M. Pena and M.-L. Rapún, Application of the topological derivative to post-processing infrared time-harmonic thermograms for defect detection, J. Math. Ind. 10 (2020), Paper No. 4. 10.1186/s13362-020-0072-9Search in Google Scholar

[61] J. Pommier and B. Samet, The topological asymptotic for the Helmholtz equation with Dirichlet condition on the boundary of an arbitrarily shaped hole, SIAM J. Control Optim. 43 (2004), no. 3, 899–921. 10.1137/S036301290241616XSearch in Google Scholar

[62] R. Potthast, A survey on sampling and probe methods for inverse problems, Inverse Problems 22 (2006), no. 2, R1–R47. 10.1088/0266-5611/22/2/R01Search in Google Scholar

[63] A. G. Ramm, Inverse Problems. Mathematical and Analytical Techniques with Applications to Engineering, Springer, New York, 2005. Search in Google Scholar

[64] M.-L. Rapún, On the solution of direct and inverse multiple scattering problems for mixed sound-soft, sound-hard and penetrable objects, Inverse Problems 36 (2020), no. 9, Article ID 095014. 10.1088/1361-6420/ab98a2Search in Google Scholar

[65] M.-L. Rapún and F.-J. Sayas, Boundary integral approximation of a heat-diffusion problem in time-harmonic regime, Numer. Algorithms 41 (2006), no. 2, 127–160. 10.1007/s11075-005-9002-6Search in Google Scholar

[66] M.-L. Rapún and F.-J. Sayas, Indirect methods with Brakhage–Werner potentials for Helmholtz transmission problems, Numerical Mathematics and Advanced Applications, Springer, Berlin (2006), 1146–1154. 10.1007/978-3-540-34288-5_115Search in Google Scholar

[67] M.-L. Rapún and F.-J. Sayas, Boundary element simulation of thermal waves, Arch. Comput. Methods Eng. 14 (2007), no. 1, 3–46. 10.1007/s11831-006-9000-4Search in Google Scholar

[68] M.-L. Rapún and F.-J. Sayas, Mixed boundary integral methods for Helmholtz transmission problems, J. Comput. Appl. Math. 214 (2008), no. 1, 238–258. 10.1016/j.cam.2007.02.028Search in Google Scholar

[69] B. Samet, S. Amstutz and M. Masmoudi, The topological asymptotic for the Helmholtz equation, SIAM J. Control Optim. 42 (2003), no. 5, 1523–1544. 10.1137/S0363012902406801Search in Google Scholar

[70] J. Sokoł owski and J.-P. Zolésio, Introduction to Shape Optimization. Shape Sensitivity Analysis, Springer Ser. Comput. Math. 16, Springe, Berlin, 1992. 10.1007/978-3-642-58106-9Search in Google Scholar

[71] M. E. Taylor, Partial Differential Equations II. Qualitative Studies of Linear Equations, 2nd ed., Appl. Math. Sci. 116, Springer, New York, 2011. 10.1007/978-1-4419-7052-7Search in Google Scholar

[72] R. H. Torres and G. V. Welland, The Helmholtz equation and transmission problems with Lipschitz interfaces, Indiana Univ. Math. J. 42 (1993), no. 4, 1457–1485. 10.1512/iumj.1993.42.42067Search in Google Scholar

[73] T. von Petersdorff, Boundary integral equations for mixed Dirichlet, Neumann and transmission problems, Math. Methods Appl. Sci. 11 (1989), no. 2, 185–213. 10.1002/mma.1670110203Search in Google Scholar

[74] L. Wienert, Die numerische Approximation von Randintegraloperatoren für die Helmholtzgleichung R 3 , PhD thesis, University of Linz, 1999. Search in Google Scholar

[75] J. Xu and L. Zikatanov, Some observations on Babuška and Brezzi theories, Numer. Math. 94 (2003), no. 1, 195–202. 10.1007/s002110100308Search in Google Scholar

Received: 2021-12-10

Revised: 2022-06-25

Accepted: 2022-07-09

Published Online: 2022-08-30

Published in Print: 2022-10-01

Articles in the same Issue

https://doi.org/10.1515/cmam-2021-0234

Keywords for this article

Boundary Integral Methods; Spectral Methods; Inverse Scattering; Shape Reconstruction; Topological Energy; Multifrequency