Computing the Sound of the Sea in a Seashell
CCOMPUTING THE SOUND OF THE SEA IN A SEASHELL
JONATHAN BEN-ARTZI, MARCO MARLETTA, AND FRANK R ¨OSLER
Abstract.
The question of whether there exists an approximation procedure to compute theresonances of any Helmholtz resonator, regardless of its particular shape, is addressed. A positiveanswer is given, and it is shown that all that one has to assume is that the resonator chamber isbounded and that its boundary is C . The proof is constructive, providing a universal algorithmwhich only needs to access the values of the characteristic function of the chamber at any requestedpoint. Introduction
This paper provides an affirmative answer to the following question:
Does there exist a universal algorithm for computing the resonances of the Laplacianin R \ U for any open bounded set U ⊂ R ? Any domain U ⊂ R d gives rise to resonances, i.e. special frequencies that are ‘nearly’ eigenvalues ofthe Laplacian in R d \ U (with appropriate boundary conditions). The most famous example is the“sound of the sea” in a seashell: when we hold a seashell against our ear we hear frequencies thatare nearly the eigenvalues of the Laplacian in the closed cavity with Neumann boundary conditions.This is also known as a Helmholtz resonator [18]. We therefore refer to U as a “resonator”.We are interested in computing (Dirichlet) resonances in a way that is independent of the domain U itself; i.e. U is the input of the problem, and the computation returns the associated resonances.We work in the plane, i.e. d = 2, but our results can be adapted to higher dimensions. To ourbest knowledge this is the first time this question is addressed. Furthermore, the proof of existenceprovides an actual algorithm (that is, the proof is constructive). We test this algorithm on somestandard examples, and compare to known results.The framework required for this analysis is furnished by the Solvability Complexity Index (SCI),which is an abstract theory for the classification of the computational complexity of problems thatare infinite-dimensional. This framework has been developed over the last decade by Hansen andcollaborators (cf. [15, 4]) and draws inspiration from the seminal result [10] on solving quinticequations via a tower of algorithms . We therefore emphasize that ours is an abstract resultin analysis, not in numerical analysis.
Resonances.
Let us first define what we mean by a resonance. Let U ⊂ R d be an open set andassume that ∂U ∈ C . Let H be the Laplacian in L ( R d \ U ) with homogeneous Dirichlet boundaryconditions on ∂U . Resonances of H can be defined via analytic continuation of the associatedDirichlet-to-Neumann (DtN) operators. Indeed, we can start from a characterization of eigenvaluesand perform analytic continuation as follows. Denote the spectral parameter by k (with the branchcut of the complex square root running along the positive real line). Then k is an eigenvalue of H if and only if there exists a function u ∈ L ( R d \ U ), such that ( − ∆ − k ) u = 0 in R d \ U and u = 0on ∂U . The existence of such a u is equivalent to the following. Let R > U ⊂ B R (the open ball of radius R around 0). Denote by M inner ( k ) and M outer ( k ), respectively, the innerand outer DtN maps on ∂B R associated with − ∆ − k . M inner ( k ) and M outer ( k ) can be shown to be Date : September 9, 2020.2010
Mathematics Subject Classification.
Key words and phrases.
Helmholtz resonator, Dirichlet-to-Neumann map, Solvability Complexity Index, Com-putational complexity.JBA and FR acknowledge support from an Engineering and Physical Sciences Research Council FellowshipEP/N020154/1 and MM acknowledges support from an Engineering and Physical Sciences Research Council GrantEP/T000902/1. a r X i v : . [ m a t h . SP ] S e p JONATHAN BEN-ARTZI, MARCO MARLETTA, AND FRANK R ¨OSLER analytic operator-valued functions of k in the upper half plane C + := { z ∈ C | Im( z ) > } . If thereexists an eigenfunction u of H with eigenvalue k , then u | ∂B R ∈ ker( M inner ( k ) + M outer ( k )) . Conversely, if φ ∈ ker( M inner ( k ) + M outer ( k )), then there exist solutions u inner and u outer of ( − ∆ − k ) u = 0 in B R \ U and R d \ B R respectively, such that u inner | ∂B R = u outer | ∂B R = φ and ∂ ν u inner | ∂B R = − ∂ ν u outer | ∂B R (note that the normal vector ν is always taken to point to the exterior of the domainconsidered, that is, away from 0 for B R and towards 0 for R d \ B R ). Hence, the function u := (cid:40) u inner on B R \ Uu outer on R d \ B R (1.1)is in H ( R d \ U ). This shows that a complex number k ∈ C + in the upper half plane is a (squareroot of an) eigenvalue if and only if ker( M inner ( k ) + M outer ( k )) (cid:54) = ∅ . A resonance can now be definedas follows. Definition 1.1.
Let us use the same symbols M inner ( · ) , M outer ( · ) to denote the meromorphic con-tinuations of the DtN maps to all of C . A number k ∈ C − := { z ∈ C | Im( z ) < } in the lower halfplane is called a resonance of H , ifker( M inner ( k ) + M outer ( k )) (cid:54) = { } . (1.2) Remark . Resonances in the sense of Definition 1.1 are well-defined. Indeed, M outer ( k ) is well-defined for all k ∈ C and the only poles of M inner ( k ) are equal to the Dirichlet eigenvalues of B R \ U .But these all lie on the real axis, which is excluded in Definition 1.1.Moreover, an argument similar to the one surrounding eq. (1.1) shows that the nontriviality ofker( M inner ( k ) + M outer ( k )) is independent of the value of R as long as U ⊂ B R . Remark . In two dimensions, the DtN maps M inner , M outer can actually be analytically continuedto the logarithmic cover of C , rather than just C − (cf. [12, Sec. 3.1.4]). However, since the vastmajority of the literature is concerned with resonances near the real axis, we decided to simplify ourpresentation by considering only resonances in C − . A strategy for dealing with the general case hasbeen outlined in [6, Sec. 4.2.2].We show that resonances can be computed as the limit of a sequence of approximations, each ofwhich can be computed precisely using finitely many arithmetic operations and accessing finitelymany values of the characteristic function χ U of U . The proof is constructive: we define an algorithmand prove its convergence. We emphasize that this single algorithm is valid for any open U ⊂ R with ∂U ∈ C . We implement this algorithm and compare its output to known results.1.2. The Solvability Complexity Index.
The Solvability Complexity Index (SCI) addressesquestions which are at the nexus of pure and applied mathematics, as well as computer science:
How do we compute objects that are “infinite” in nature if we can only handle afinite amount of information and perform finitely many mathematical operations?Indeed, what do we even mean by “computing” such an object?
These broad topics are addressed in the sequence of papers [15, 4, 5]. Let us summarize the maindefinitions and discuss how these relate to our problem of finding resonances:
Definition 1.4 (Computational problem) . A computational problem is a quadruple (Ω , Λ , Ξ , M ),where(i) Ω is a set, called the primary set ,(ii) Λ is a set of complex-valued functions on Ω, called the evaluation set ,(iii) M is a metric space,(iv) Ξ : Ω → M is a map, called the problem function . Definition 1.5 (Arithmetic algorithm) . Let (Ω , Λ , Ξ , M ) be a computational problem. An arith-metic algorithm is a map Γ : Ω → M such that for each T ∈ Ω there exists a finite subset Λ Γ ( T ) ⊂ Λsuch that(i) the action of Γ on T depends only on { f ( T ) } f ∈ Λ Γ ( T ) , OMPUTING THE SOUND OF THE SEA IN A SEASHELL 3 (ii) for every S ∈ Ω with f ( T ) = f ( S ) for all f ∈ Λ Γ ( T ) one has Λ Γ ( S ) = Λ Γ ( T ),(iii) the action of Γ on T consists of performing only finitely many arithmetic operations on { f ( T ) } f ∈ Λ Γ ( T ) . Definition 1.6 (Tower of arithmetic algorithms) . Let (Ω , Λ , Ξ , M ) be a computational problem. A tower of algorithms of height k for Ξ is a family Γ n ,n ,...,n k : Ω → M of arithmetic algorithms suchthat for all T ∈ Ω Ξ( T ) = lim n k →∞ · · · lim n →∞ Γ n ,n ,...,n k ( T ) . Definition 1.7 (SCI) . A computational problem (Ω , Λ , Ξ , M ) is said to have a Solvability Com-plexity Index (SCI) of k if k is the smallest integer for which there exists a tower of algorithms ofheight k for Ξ. If a computational problem has solvability complexity index k , we writeSCI(Ω , Λ , Ξ , M ) = k. Setting of the problem and main result.
Let us describe the elements comprising ourcomputational problem, followed by our main theorem. (i) The primary set Ω . We defineΩ = (cid:8) U ⊂ R open | U is bounded and ∂U ∈ C (cid:9) . (ii) The evaluation set Λ . The evaluation set, which describes the data at our algorithm’s dis-posal, is comprised of (all points in) the set U , as well as some basic special functions which weassume we can access at no additional computational cost and any (nonnegative) real number (suchas e , π or √ ∪ Λ ∪ Λ ∪ Λ ∪ Λ , (1.3)where Λ = { U (cid:55)→ χ U ( x ) | x ∈ R } Λ = { U (cid:55)→ e inθ | θ ∈ [0 , π ) n ∈ N } Λ = { U (cid:55)→ J n ( z ) | z ∈ C − , n ∈ N } Λ = { U (cid:55)→ H (1) n ( z ) | z ∈ C − , n ∈ N } Λ = { U (cid:55)→ x | x ∈ [0 , ∞ ) } , where J n , H (1) n denote the Bessel and Hankel functions of the first kind, respectively. Including thecharacteristic functions U (cid:55)→ χ U ( x ) in Λ means that for every x ∈ R we can test whether x isincluded in U or not. (iii) The metric space M . M is the space cl( C ) of all closed subsets of C equipped with theAttouch-Wets metric, generated by the following distance function: Definition 1.8 (Attouch-Wets distance) . Let
A, B be closed sets in C . The Attouch-Wets distance between them is defined as d AW ( A, B ) = ∞ (cid:88) i =1 − i min (cid:40) , sup | x |
A, B ⊂ C are bounded, then d AW is equivalent to the Hausdorff distance d H . Fur-thermore, it can be shown (cf. [3, Ch. 3]) that(1.4) d H ( A n ∩ B, A ∩ B ) → B ⊂ C compact ⇒ d AW ( A n , A ) → . JONATHAN BEN-ARTZI, MARCO MARLETTA, AND FRANK R ¨OSLER (iv) The problem function Ξ . The problem function Ξ( U ) = Res( U ) is the map that associatesto each U ∈ Ω the set of resonances (as defined in Definition 1.1) of H .With these at hand, we can state our main theorem: Theorem 1.9.
There exists a family of arithmetic algorithms { Γ n } n ∈ N such that Γ n ( U ) → Res( U ) as n → + ∞ for any U ∈ Ω , where the convergence is in the sense of the Attouch-Wets metric. Thatis, SCI(Ω , Λ , Res( · ) , (cl( C ) , d AW )) = 1 .Remark . In Subsection 5.2 we outline a strategy for adapting this result to Neumann boundaryconditions, and implement a numerical model consisting of four circular holes: this is relevant forunderstanding resonances caused by oil-drilling platforms [13], for instance.Theorem 1.9 is proved by identifying an operator of the form Id L + K ( k ), where K ( k ) is in someSchatten class C p , with the property that ker(Id L + K ( k )) (cid:54) = { } if and only if ker( M inner ( k ) + M outer ( k )) (cid:54) = { } . The former is equivalent to det (cid:100) p (cid:101) (Id L + K ( k )) = 0, where det (cid:100) p (cid:101) denotes the (cid:100) p (cid:101) -perturbation determinant. This determinant is approximated via a finite element procedure on B R \ U with explicit error bounds. A thresholding procedure then yields an approximation for thezero set of det (cid:100) p (cid:101) (Id L + K ( k )).1.4. Discussion.
The study of cavity resonances is much older than the study of quantum mechan-ical resonances. The foundational work is generally ascribed to Helmholtz [18], who in the 1850s hadconstructed devices which were designed to identify special frequencies from within a sound wave.These
Helmholtz resonators consisted of a small aperture at one end to admit sound, and a largeropening at the other to emit it.The operator theoretic foundations of the study of cavity resonances are usually based on variousapproaches to the theory of scattering by obstacles. A good early treatment of both quantum-mechanical and obstacle scattering may be found in the seminal text of Newton [25]. The bookof Lax and Phillips [20] is perhaps the most famous work on scattering theory in a semi-abstractsetting, though the whole subject was extensively researched in the Soviet school by mathematicianssuch as Faddeev and Pavlov, see, e.g. [14], and also the monograph of Yafaev [31].There is extensive work in the applied mathematics literature too, devoted to estimating reso-nances in various asymptotic regimes, either geometrical or semiclassical: see [16, 17, 23, 30, 28] asstarting points. Assemblies of resonators are proposed as cloaking devices in a variety of contexts,see, e.g., [1]; the mathematical treatment of acoustic waveguides is very similar to that of optical orquantum waveguides. A very up-to-date overview of the current state of the art in the subject maybe found in [12].Numerical approaches use a variety of techniques of which complex scaling (rediscovered as ‘per-fectly matched layers’), boundary integral techniques and combinations of these with special finiteelement methods are the most common.
Organization of the paper . In Section 2 we analyze the inner and outer DtN maps and obtain a newexpression equivalent to (1.2) in which M inner ( k ) + M outer ( k ) is replaced by an expression of theform Id L +perturbation. Section 3 is dedicated to the construction of a finite element method forthe approximation of this perturbation, and in Section 4 the algorithm for approximating Res( U ) isdefined and shown to converge. Finally, in Section 5 we provide some numerical examples.2. Formulas for the inner and outer DtN maps
As a first step, we prove a weaker version of Theorem 1.9. We show that under the additionalassumption that the diameter of the resonator is known a priori , there exists an algorithm thatcomputes the set of resonances in one limit. To this end we defineΩ R := (cid:8) U ⊂ R open | ∂U ∈ C , U ⊂ B R − (cid:9) where B ρ is the open ball about the origin of radius ρ >
0. Then one has:
Theorem 2.1.
There exists a family of arithmetic algorithms { Γ Rn } n ∈ N depending on R , such that Γ Rn ( U ) → Res( U ) as n → + ∞ for any U ∈ Ω R , where the convergence is in the sense of theAttouch-Wets metric. That is, SCI(Ω R , Λ , Res( · ) , (cl( C ) , d AW )) = 1 . OMPUTING THE SOUND OF THE SEA IN A SEASHELL 5
U B R − B R Figure 1.
A sketch of the resonator U Henceforth, the radius
R > ∂B R e n ( θ ) := e i nθ √ πR , θ ∈ [0 , π ) , n ∈ Z , which will be used frequently throughout the paper.2.1. The outer map.
The outer DtN map acts on a function φ ∈ H ( ∂B R ) as M outer ( k ) φ = ∂ ν u, where u ∈ L ( R \ B R ) solves the problem (cid:40) ( − ∆ − k ) u = 0 in R \ B R u = φ on ∂B R . In the orthonormal basis { e n } n ∈ Z the map M outer ( k ) has the explicit representation M outer ( k ) = diag − k H (1) | n | (cid:48) ( kR ) H (1) | n | ( kR ) , n ∈ Z , where H (1) ν denote the Hankel functions of the first kind. Using well-known identities for Besselfunctions (cf. [21]), this can be rewritten as(2.1) M outer ( k ) = diag | n | R − k H (1) | n |− ( kR ) H (1) | n | ( kR ) , n ∈ Z . Moreover, it can be seen that H (1) | n |− ( z ) H (1) | n | ( z ) ∼ z | n | as | n | → + ∞ , so that the matrix elements of M outer grow linearly in | n | .2.2. The inner map.
We first show that the inner DtN map can be decomposed into a part thatdoes not depend on U and a bounded part: JONATHAN BEN-ARTZI, MARCO MARLETTA, AND FRANK R ¨OSLER
Lemma 2.2.
Let M inner , denote the free inner DtN map, i.e. the one with U = ∅ . Then thereexists a meromorphic operator-valued function C (cid:51) k (cid:55)→ K ( k ) with values in the bounded operatorson L ( ∂B R ) such that (2.2) M inner ( k ) = M inner , ( k ) + K ( k ) . Proof.
For the inner DtN map, we first isolate the contribution of the resonator using the followingcalculation. Let φ ∈ H ( ∂B R ). Then M inner ( k ) φ = ∂ ν u | ∂B R , where u solves ( − ∆ − k ) u = 0 in B R \ U ,u = φ on ∂B R ,u = 0 on ∂U. (2.3)Define a new map S ( k ) : H ( ∂B R ) → H ( B R ) by S ( k ) φ = w , where w solves (cid:40) ( − ∆ − k ) w = 0 in B R ,w = φ on ∂B R . That is, S ( k ) φ is the harmonic extension of φ into B R . The regularity properties of S ( k ) follow from[24, Th. 4.21]. Next choose any smooth radial function ρ satisfying (cid:40) ρ ≡ ∂B R ρ ≡ B R − (2.4)(e.g. ρ ( r ) can be chosen as a piecewise polynomial in r ; then the values ρ ( r ) can be computed from r in finitely many algebraic steps) and introduce the function v := u − ρS ( k ) φ on B R \ U . Then by construction we have ( − ∆ − k ) v = (2 ∇ ρ · ∇ + ∆ ρ ) S ( k ) φ in B R \ U ,v = 0 on ∂B R ,v = 0 on ∂U. In operator terms, this means that v = ( H D − k ) − T ρ S ( k ) φ, (2.5)where H D denotes the Laplacian on L ( B R \ U ) with homogeneous Dirichlet boundary condition on ∂ ( B R \ U ) and where we have denoted T ρ := 2 ∇ ρ · ∇ + ∆ ρ. Substituting the new representation for v into u = ρS ( k ) φ + v , we immediately obtain M inner ( k ) φ = ∂ ν u | ∂B R = ∂ ν ( ρS ( k ) φ ) + ∂ ν ( H D − k ) − T ρ S ( k ) φ = ∂ ν ρ (cid:124)(cid:123)(cid:122)(cid:125) =0 S ( k ) φ + ρ (cid:124)(cid:123)(cid:122)(cid:125) =1 ∂ ν S ( k ) φ + ∂ ν ( H D − k ) − T ρ S ( k ) φ = M inner , ( k ) φ + ∂ ν ( H D − k ) − T ρ S ( k ) φ. (2.6)Let us next study the operator ∂ ν ( H D − k ) − T ρ S ( k ) appearing in the last term in (2.6). Theharmonic extension operator S ( k ) extends to a bounded operator S ( k ) : L ( ∂B R ) → H ( B R )(cf. [24, Th. 6.12]). Moreover, we have that T ρ : H ( B R ) → H − ( B R ) is bounded. Under theassumption that ∂U is C , we may conclude from elliptic regularity that( H D − k ) − T ρ S ( k ) : L ( ∂B R ) → H ( B R ) OMPUTING THE SOUND OF THE SEA IN A SEASHELL 7 as a bounded operator. Finally, by the trace theorem, the normal derivative operator ∂ ν is continuousfrom H ( B R ) to L ( ∂B R ), so we conclude that ∂ ν ( H D − k ) − T ρ S ( k ) : L ( ∂B R ) → L ( ∂B R )is bounded. This completes the proof. (cid:3) The decomposition in Lemma 2.2 allows us to reduce the study of M inner to that of M inner , .Next, we note that M inner , can in fact be represented as a diagonal operator in the basis { e n } n ∈ Z in a very similar manner to M outer . Indeed, one has M inner , ( k ) = diag (cid:32) k J (cid:48)| n | ( kR ) J | n | ( kR ) , n ∈ Z (cid:33) (2.7) = diag (cid:18) | n | R − k J | n | +1 ( kR ) J | n | ( kR ) , n ∈ Z (cid:19) , where J ν denote the Bessel functions of the first kind. It can be seen that J | n | +1 ( z ) J | n | ( z ) ∼ z | n | as | n | → + ∞ , so that the matrix elements of M inner , grow linearly in | n | .2.3. Approximation procedure.
From equations (2.1), (2.2) and (2.7) we have M inner ( k ) + M outer ( k ) = 2 R N + H ( k ) + J ( k ) + K ( k ) , where J ( k ) = diag (cid:18) − k J | n | +1 ( kR ) J | n | ( kR ) , n ∈ Z (cid:19) H ( k ) = diag (cid:18) − k H (1) | n |− ( kR ) H (1) | n | ( kR ) , n ∈ Z (cid:19) N = diag (cid:0) | n | , n ∈ Z (cid:1) K ( k ) = ∂ ν ( H D − k ) − (2 ∇ ρ · ∇ + ∆ ρ ) S ( k ) . Note that H , J , K are bounded and analytic in k . In order to determine whether ker( M inner + M outer ) (cid:54) = { } we want to transform M inner + M outer into an operator of the form I + (compact). Tothis end, we introduce the bijective operator N := N + P = diag( ..., , , , , , , , ... ), where P denotes the projection onto the zeroth component. This new notation brings M inner + M outer intothe form M inner ( k ) + M outer ( k ) = 2 R N + H ( k ) + J ( k ) + K ( k ) − R P = 2 R N (cid:18) I + R N − (cid:0) H ( k ) + J ( k ) + K ( k ) (cid:1) N − − N − P N − (cid:19) N = 2 R N (cid:18) I + R N − (cid:0) H ( k ) + J ( k ) + K ( k ) (cid:1) N − − P (cid:19) N . Because N = diag( ..., √ , √ , , , , √ , √ , ... ) is bijective from its domain to L ( ∂B R ), we haveker( M inner + M outer ) = { } ⇔ ker (cid:18) I + R N − (cid:0) H ( k ) + J ( k ) + K ( k ) (cid:1) N − − P (cid:19) = { } . Now, observe that for all p > R N − (cid:0) H ( k ) + J ( k ) + K ( k ) (cid:1) N − − P ∈ C p , the p -Schatten class(because H , J , K are bounded and N − is in C p ). Therefore, its perturbation determinant det (cid:100) p (cid:101) is defined, where (cid:100)·(cid:101) denotes the ceiling function (cf. [11, Sec. XI.9]). Our task is reduced to findingthe zeros k ∈ C − of the analytic functiondet (cid:100) p (cid:101) (cid:18) I + R N − (cid:0) H ( k ) + J ( k ) + K ( k ) (cid:1) N − − P (cid:19) . In order to approximate this determinant by something computable, we first prove that a squaretruncation of the operator matrix will converge in Schatten norm.
JONATHAN BEN-ARTZI, MARCO MARLETTA, AND FRANK R ¨OSLER
Lemma 2.3.
Let k ∈ C − and for n ∈ N let P n : L ( ∂B R ) → span { e − n , . . . e n } be the orthogonalprojection. Then there exists a constant C > depending only on the set U such that (cid:13)(cid:13)(cid:13) N − ( H + J + K ) N − − P n N − ( H + J + K ) N − P n (cid:13)(cid:13)(cid:13) C p ≤ Cn − + p . Proof.
The proof is given by a direct calculation: (cid:13)(cid:13)(cid:13) N − ( H + J + K ) N − − P n N − ( H + J + K ) N − P n (cid:13)(cid:13)(cid:13) C p ≤ (cid:13)(cid:13)(cid:13) ( I − P n ) N − ( H + J + K ) N − (cid:13)(cid:13)(cid:13) C p + (cid:13)(cid:13)(cid:13) P n N − ( H + J + K ) N − ( I − P n ) (cid:13)(cid:13)(cid:13) C p ≤ (cid:107)H + J + K(cid:107) L ( L ( ∂B R )) (cid:13)(cid:13)(cid:13) ( I − P n ) N − (cid:13)(cid:13)(cid:13) C p ≤ C (cid:107)H + J + K(cid:107) L ( L ( ∂B R )) n − p p . (cid:3) Our problem is therefore reduced to computing the matrix elements of the truncated operator P n N − ( H + J + K ) N − P n − P . This is done by performing a finite element procedure on B R \ U .This is the objective of the next section. Remark . We remind the reader that the algorithm is assumed to have access to the values ofBessel and Hankel functions (recall (1.3)) so that the values of P n N − ( H + J ) N − P n − P can becomputed in finitely arithmetic operations. Hence we only need to compute P n N − KN − P n .3. Matrix elements of the inner DtN map
The goal of this section is to compute the matrix elements K αβ := (cid:90) ∂B R e β K ( k ) e α dσ where { e α } α ∈ Z is the orthonormal basis of L ( ∂B R ) defined in the previous section and where werecall that K ( k ) = ∂ ν ( H D − k ) − T ρ S ( k ). This is done by converting the integral on the boundary ∂B R to an integral on the interior, using Green’s theorem. To this end, we extend the basis functionsto the interior, by defining (in polar coordinates) E β ( r, θ ) := ρ ( r ) e β ( θ ) , where ρ was defined in (2.4), and denote v α := ( H D − k ) − T ρ S ( k ) e α f α := T ρ S ( k ) e α (3.1)so that v α = ( H D − k ) − f α . Now note that by Green’s theorem K αβ = (cid:90) ∂B R e β ∂ ν v α dσ = (cid:90) B R \ U E β ∆ v α dx + (cid:90) B R \ U ∇ E β · ∇ v α dx = (cid:90) B R \ U E β ( − f α − k v α ) dx + (cid:90) B R \ U ∇ E β · ∇ v α dx = (cid:90) B R \ U ∇ E β · ∇ v α dx − k (cid:90) B R \ U E β v α dx − (cid:90) B R \ U E β f α dx. (3.2)Next, note that the integrand in the third term on the right hand side is explicitly given. Indeed,one can easily see that S ( k ) e α ( r, θ ) = J α ( kr ) J α ( kR ) e iαθ √ πR , (3.3) OMPUTING THE SOUND OF THE SEA IN A SEASHELL 9 where J α are the Bessel functions of the first kind. Hence, the last integral in (3.2) can be explicitlyapproximated by quadrature, with a known bound for the error. In order to compute the first andsecond terms in the right hand side of (3.2), we need to approximate the function v α , which is thesolution of the boundary value problem (cid:40) ( − ∆ − k ) v α = (2 ∇ ρ · ∇ + ∆ ρ ) S ( k ) e α in B R \ U ,v α = 0 on ∂ ( B R \ U ) . (3.4)Note again that the right hand side of this equation is explicitly given in terms of (3.3).3.1. Finite element approximation of v α . For notational convenience we write O := B R \ U .
In the current subsection we prove:
Proposition 3.1.
Let v α be the solution of (3.4) . For small discretization parameter h > thereexists a piecewise linear function v hα which is computable from a finite subset of (1.3) in finitelymany algebraic steps, which satisfies the error estimate (cid:107) v α − v hα (cid:107) H ( O ) ≤ Ch (cid:107) f α (cid:107) H ( O ) , where C is independent of h and α , and where f α is defined in (3.1) . The reader should bear in mind that all constants obtained in our estimates will actually dependon k . However, these constants are uniformly bounded, as long as k varies in a compact setseparated from the spectrum of H D (recall that H D is the Laplacian on L ( O ) with homogeneousDirichlet boundary conditions on ∂ O ). Note that, since σ ( H D ) is strictly positive and discrete, anycompact subset of C − has this property.Let h > neighbors of a grid point j ∈ h Z . Definition 3.2.
Let j = ( j , j ) ∈ h Z . The set of neighbors of j is the set of all grid pointscontained in the cell [ j − h, j + 2 h ] × [ j − h, j + 2 h ] (cf. Figure 2). Remark . The apparent asymmetry in the definition above is merely due to the fact that a point j defines a cell that is “northeast” of it, i.e. [ j , j + h ] × [ j , j + h ].Define a uniform rectangular grid by L h := (cid:8) j ∈ h Z ∩ B R (cid:12)(cid:12) all neighbors of j are contained in O (cid:9) and note that L h can be completely constructed from the knowledge of χ U ( j ) , j ∈ h Z ∩ B R infinitely many steps. Next, define the open polygonal domain O h generated by L h as follows. O h := int (cid:18) (cid:91) j ∈ L h [0 , h ] + j (cid:19) . Note that as soon as h < ∂ O )) = ∂U )) one will have O h ⊂ O and for all x ∈ ∂ O h (cf. Figure 2) √ h ≤ dist( x, ∂ O ) ≤ √ h. (3.5) ∂ O hjh O h ≥ √ h Figure 2.
Illustration of eq. (3.5). The point j is included in L h , if and only if all thesolid black points are in O . These solid black points are the neighbors of j . Finally, we choose the obvious triangulation T h for O h obtained by splitting each square [0 , h ] + j along its diagonal and denote by W h the corresponding P1 finite element space with zero boundaryconditions, i.e. W h := (cid:8) w ∈ H ( O ) (cid:12)(cid:12) w ≡ O \ O h and w piecewise linear in O h w.r.t. T h (cid:9) Next, we will construct a finite element approximation for the solution of problem (3.4) startingfrom C´ea’s classical lemma:
Lemma 3.4 (C´ea’s Lemma, [9, Th. 2.4.1]) . Let W h ⊂ H ( O ) be a family of subspaces and assumethe bilinear form a : H ( O ) × H ( O ) → C and linear functional f ∈ H − ( O ) satisfy the hypothesesof the Lax-Milgram theorem. Let v be the solution of the problem a ( v, w ) = (cid:104) f, w (cid:105) for all w ∈ H ( O ) and v h be the solution of a ( v h , w h ) = (cid:104) f, w h (cid:105) for all w h ∈ W h . (3.6) Then there exists
C > s.t. (cid:107) v − v h (cid:107) H ( O ) ≤ C inf w h ∈W h (cid:107) v − w h (cid:107) H ( O ) (3.7) for all h > . We apply this lemma to our present situation by choosing W h = W h a ( v, w ) = (cid:82) O ∇ v · ∇ w dxf = f α v = v α as in (3.4) v hα = the solution of (3.6) . Remark . Note that the finite element approximation v hα is not computable from the informationin (1.3) in finitely many algebraic operations. An additional step of numerical integration is necessaryto achieve this. We will revisit this issue at the end of the section.Lemma 3.4 tells us that in order to estimate the approximation error (cid:107) v α − v hα (cid:107) H ( O ) it is enoughto estimate the projection error inf w h ∈ W h (cid:107) v α − w h (cid:107) H ( O ) . For the latter, we have (cid:107) v α − w h (cid:107) H ( O ) ≤ (cid:107) v α (cid:107) H ( O\O h ) + (cid:107) v α − w h (cid:107) H ( O h )OMPUTING THE SOUND OF THE SEA IN A SEASHELL 11 for any w h ∈ W h . To estimate the right hand side, we further decompose O h into a boundary layerand an interior part. More precisely, we write O h = O int h ∪ O bdry h , where O bdry h is the union of all cells having nonempty intersection with ∂ O h , and O int h = O h \ O bdry h (cf. Figure 3). OO bdry h O int h K Figure 3.
Left: Sketch of O , O h and O bdry h . Right: Sketch of the chosen triangulation. Hence, for any w h ∈ W h we have (cid:107) v α − w h (cid:107) H ( O ) ≤ (cid:107) v α − w h (cid:107) H ( O\O int h ) + (cid:107) v α − w h (cid:107) H ( O int h ) (3.8)Next, estimate the all three terms on the r.h.s. of (3.8) individually. To estimate the two last termswe will choose a particular “test function” u h ∈ W h . We first recall a classical interpolation estimatefor linear finite elements [9, Th. 3.1.6]: Theorem 3.6 (Classical interpolation estimate) . Let K be an element of T h and let Π hK denote thelinear interpolation operator on H ( K ) , i.e. Π hK w is the affine function with (Π hK w )( j ) = w ( j ) forall corners j of K . Then there exists C > (independent of h, w, K ) such that (cid:107) w − Π hK w (cid:107) H ( K ) ≤ Ch (cid:107) D w (cid:107) L ( K ) , (cid:107) w − Π hK w (cid:107) L ( K ) ≤ Ch (cid:107) D w (cid:107) L ( K ) . Now, we choose an explicit function u h ∈ W h as follows. • For each finite element K ⊂ O int h , let u h | K = Π hK ( v α ) be the linear interpolation of v α , i.e. u h ( j ) = v α ( j ) for all nodes j ∈ ∂K ∩ L h . • For boundary elements K ⊂ O bdry h , we set u h ( i ) = 0 , for boundary nodes i,u h ( j ) = v α ( j ) , for interior nodes j. Here a boundary node is any i ∈ L h ∩ ∂ O h . Note that testing whether a node is a boundarynode can be done in finitely many steps, because boundary nodes are precisely those nodeswhose neighbors are not all contained in O .This defines a continuous, piecewise linear function in W h . We have the following error estimates. Lemma 3.7 (third term in (3.8)) . One has the interior error estimate (cid:107) v α − u h (cid:107) H ( O int h ) ≤ Ch (cid:107) v α (cid:107) H ( O int h ) . Proof.
This follows immediately from Theorem 3.6. (cid:3)
Lemma 3.8 (first term in (3.8)) . One has the boundary layer error estimate (cid:107) v α − u h (cid:107) H ( O\O int h ) ≤ Ch (cid:107) v α (cid:107) H ( O\O int h ) . Proof.
We begin with (cid:107) v α − u h (cid:107) L ( O\O int h ) . Let x ∈ O \ O int h and note that similarly to eq. (3.5)there exists y ∈ ∂ O with | x − y | ≤ √ h . By [8, Ch. 9, Remark 7] this implies | v α ( x ) | = | v α ( x ) − v α ( y ) | ≤ | x − y |(cid:107)∇ v α (cid:107) L ∞ ( O\O int h ) . (3.9)Hence (cid:107) v α (cid:107) L ∞ ( O\O int h ) ≤ √ h (cid:107)∇ v α (cid:107) L ∞ ( O\O int h ) ≤ Ch (cid:107) v α (cid:107) H ( O\O int h ) , (3.10)where the last line follows from the Gagliardo-Sobolev-Nirenberg inequality. Moreover, we note thatby the definition of u h , one has (cid:107) u h (cid:107) L ∞ ( O\O int h ) ≤ (cid:107) v α (cid:107) L ∞ ( O\O int h ) . These facts allow us to estimate (cid:107) v α − u h (cid:107) L ( O\O int h ) ≤ (cid:107) v α (cid:107) L ( O\O int h ) + (cid:107) u h (cid:107) L ( O\O int h ) ≤ |O \ O int h | (cid:0) (cid:107) v α (cid:107) L ∞ ( O\O int h ) + (cid:107) u h (cid:107) L ∞ ( O\O int h ) (cid:1) ≤ Ch (cid:107) v α (cid:107) L ∞ ( O\O int h ) ≤ Ch (cid:107) v α (cid:107) H ( O\O int h ) (3.11)where we have used (3.10) in the last line.Turning to the gradient, a similar estimate is obtained as follows. By definition of u h , one has |∇ u h ( x ) | ≤ (cid:40) x ∈ O \ O hCh (cid:107) v α (cid:107) L ∞ ( K ) for x ∈ K ⊂ O bdry h In particular, for every K ⊂ O bdry h one has (cid:107)∇ u h (cid:107) L ( K ) ≤ C (cid:107) v α (cid:107) L ∞ ( K ) ≤ Ch (cid:107) v α (cid:107) H ( K ) by (3.9).(3.12)Thus (cid:107)∇ v α − ∇ u h (cid:107) L ( O\O int h ) ≤ (cid:107)∇ v α (cid:107) L ( O\O int h ) + (cid:107)∇ u h (cid:107) L ( O\O int h ) = (cid:107)∇ v α (cid:107) L ( O\O int h ) + (cid:107)∇ u h (cid:107) L ( O bdry h ) = (cid:107)∇ v α (cid:107) L ( O\O int h ) + (cid:88) K ⊂O bdry h (cid:107)∇ u h (cid:107) L ( K ) ≤ |O \ O int h |(cid:107)∇ v α (cid:107) L ∞ ( O\O int h ) + Ch (cid:88) K ⊂O bdry h (cid:107) v α (cid:107) H ( K ) ≤ Ch (cid:107) v α (cid:107) H ( O\O int h ) + Ch (cid:107) v α (cid:107) H ( O bdry h ) ≤ Ch (cid:107) v α (cid:107) H ( O\O int h ) (where (3.12) was used in the fourth line and (3.10) in the last line). Hence the gradient estimatebecomes (cid:107)∇ v α − ∇ u h (cid:107) L ( O\O int h ) ≤ Ch (cid:107) v α (cid:107) H ( O\O int h ) (3.13)Combining estimates (3.11) and (3.13) completes the proof. (cid:3) Lemmas 3.7 and 3.8, together with (3.7) and (3.8), finally yield the estimate (cid:107) v α − v hα (cid:107) H ( O ) ≤ C inf w h ∈ W h (cid:107) v α − w h (cid:107) H ( O ) ≤ C (cid:107) v α − u h (cid:107) H ( O ) ≤ C (cid:104) (cid:107) v α − u h (cid:107) H ( O\O int h ) + (cid:107) v α − u h (cid:107) H ( O int h ) (cid:105) ≤ C (cid:104) h (cid:107) v α (cid:107) H ( O\O int h ) + h (cid:107) v α (cid:107) H ( O int h ) (cid:105) ≤ Ch (cid:107) v α (cid:107) H ( O )OMPUTING THE SOUND OF THE SEA IN A SEASHELL 13 (for h small enough). Finally, by elliptic regularity, we obtain the estimate (cid:107) v α − v hα (cid:107) H ( O ) ≤ Ch (cid:107) f α (cid:107) H ( O ) (3.14)where we recall that f α was defined in (3.1). Remark . As we noted in Remark 3.5, the approximate solution v hα isnot computable from a finite subset of (1.3) in finitely many steps. However, the classical theory offinite elements shows that numerical integration does not change the error estimate (3.14) (see e.g.[9, Sec. 4.1]).This concludes the proof of Proposition 3.1.3.2. Matrix elements, revisited.
Going back to (3.2) (and taking into account (3.3)), in orderto compute the matrix elements K αβ , we have to approximate the integrals I := (cid:90) O E β v α dx,I := (cid:90) O ∇ E β · ∇ v α dx,I := (cid:90) O E β f α dx. (3.15)This will be a simple consequence of Prop. 3.1 and Th. 3.6. Let Π h denote the linear interpolationoperator w.r.t. the triangulation T h , i.e. Π h u ( j ) = u ( j ) for all j ∈ L h and Π h u | K is an affinefunction for all elements K . Proposition 3.10.
We have (cid:12)(cid:12)(cid:12)(cid:12) I − (cid:90) O Π h E β v hα dx (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ch β (cid:107) f α (cid:107) H ( O ) , (3.16) (cid:12)(cid:12)(cid:12)(cid:12) I − (cid:90) O Π h ∇ E β · ∇ v hα dx (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ch β (cid:107) f α (cid:107) H ( O ) , (3.17) (cid:12)(cid:12)(cid:12)(cid:12) I − (cid:90) O Π h E β Π h f α dx (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ch β (cid:107) f α (cid:107) H ( O ) . (3.18) Proof.
We begin with (3.16) and assume w.l.o.g. h <
1. Constants independent of h, α, β will bedenoted C and their value may change from line to line. By the triangle and H¨older inequalities, wehave (cid:12)(cid:12)(cid:12)(cid:12) I − (cid:90) O Π h E β v hα dx (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) E β (cid:107) L (cid:107) v α − v hα (cid:107) L + (cid:107) v hα (cid:107) L (cid:107) E β − Π h E β (cid:107) L ≤ Ch (cid:107) E β (cid:107) L (cid:107) f α (cid:107) H + Ch (cid:107) E β (cid:107) H (cid:0) (cid:107) v α (cid:107) L + (cid:107) v hα − v α (cid:107) L (cid:1) ≤ Ch (cid:107) E β (cid:107) L (cid:107) f α (cid:107) H + Ch (cid:107) E β (cid:107) H (cid:0) (cid:107) f α (cid:107) L + h (cid:107) f α (cid:107) H (cid:1) ≤ Ch (cid:107) E β (cid:107) H (cid:107) f α (cid:107) H ≤ Ch β (cid:107) f α (cid:107) H , where we have used Prop. 3.1 & Th. 3.6 in the second line, and Prop. 3.1 & elliptic regularity inthe third line. The estimate (cid:107) E β (cid:107) H ≤ Cβ used in the last line is evident from the definition of E β . We omit the calculations leading to (3.17) and (3.18), since they are entirely analogous. (cid:3) Remark . We note that (by a straightforward calculation in polar coordinates) the integral(3.15) vanishes unless α = β . In fact, one has (cid:90) O E β f α dx = 1 RJ α ( kR ) δ αβ (cid:90) RR − ρ ( r ) (cid:2) ρ (cid:48) ( r )( krJ α − ( kr ) − αJ α ( kr )) + r ∆ ρ ( r ) J α ( kr ) (cid:3) dr, where δ αβ denotes the Kronecker symbol. This one-dimensional integral can easily be approximatedby a Riemann sum. While this fact is irrelevant to our theoretical treatment, it may improveperformance and precision of numerical implementations. Let us now go back to (3.2). Prop 3.10 allows us to define an approximate n × n matrix K h by( K h ) αβ := (cid:90) O (Π h ∇ E β ) · ∇ v hα dx − k (cid:90) O (Π h E β ) v hα dx − (cid:90) O (Π h E β )(Π h f α ) dx, (3.19)whose matrix elements satisfy the error estimate |K αβ − ( K h ) αβ | ≤ C ( k ) β (cid:16) h (cid:107) f α (cid:107) H ( O ) + h (cid:107) f α (cid:107) H ( O ) (cid:17) . (3.20)The norms on the right hand side can be estimated by explicit calculation. One has (cid:107) f α (cid:107) H ( O ) ≤ Cα (cid:107) f α (cid:107) H ( O ) ≤ Cα Plugging these bounds into (3.20), we finally arrive at |K αβ − ( K h ) αβ | ≤ C ( k ) (cid:16) h β α + h β α (cid:17) . (3.21)We summarize these results in the following Theorem 3.12.
For any n ∈ N one has the operator norm error estimate ∗ (cid:107) P n K P n − K h (cid:107) L ( H ) ≤ C ( k )( h n + h n ) , where C ( k ) is independent of n, h and bounded for k in a compact set disjoint from the real line.Proof. This follows from the generalized Young inequality: (cid:107) ( K − K h ) u (cid:107) L ≤ max (cid:40) sup α ≤ n n (cid:88) β =1 |K αβ − ( K h ) αβ | , sup β ≤ n n (cid:88) α =1 |K αβ − ( K h ) αβ | (cid:41) · (cid:107) u (cid:107) L . Since |K αβ − ( K h ) αβ | ≤ C ( h n + h n ), by (3.21), this immediately implies the assertion. (cid:3) Remark . We emphasize that all terms on the right hand side of (3.19) can be computed fromthe data provided in (1.3) in finitely many steps. Indeed, the mass and stiffness matrices m , s of theP1 finite element method can be computed exactly and satisfy (cid:90) O φψ dx = (cid:88) i,j ∈ L h φ ( i ) m ij ψ ( j ) (cid:90) O ∇ φ · ∇ ψ dx = (cid:88) i,j ∈ L h φ ( i ) s ij ψ ( j )for all φ, ψ ∈ W h . 4. The algorithm
Proof of Theorem 2.1.
From Theorem 3.12 (by choosing h ∼ n − ) and Lemma 2.3, weobtain a computable approximation K h ( n ) such that (cid:13)(cid:13)(cid:13) N − ( H + J + K ) N − − P n N − ( H + J + K h ( n ) ) N − P n (cid:13)(cid:13)(cid:13) C p ≤ Cn − + p . (4.1)Using the Lipschitz continuity of perturbation determinants (cf. [27, Th. 6.5]) we conclude Lemma 4.1.
Let k ∈ C − and define A Rn ( k ) := R P n N − (cid:0) H ( k ) + J ( k ) + K h ( n ) ( k ) (cid:1) N − P n − P .Then there exists C > , which is independent of k for k in a compact subset of C − , such that (cid:12)(cid:12)(cid:12)(cid:12) det (cid:100) p (cid:101) (cid:18) I + R N − (cid:0) H ( k ) + J ( k ) + K ( k ) (cid:1) N − − P (cid:19) − det (cid:100) p (cid:101) (cid:0) I + A Rn ( k ) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) ≤ Cn − + (cid:100) p (cid:101) . ∗ We are abusing notation here, by identifying K h with K h ⊕ span { e ,...,e n } ⊥ OMPUTING THE SOUND OF THE SEA IN A SEASHELL 15
Fix a compact set Q in C − and define the grid G n = n ( Z + i Z ). Because det (cid:100) p (cid:101) ( · · · ) is analyticin k , it can only have finitely many zeros in Q and all are of finite order. That is, near a zero k onehas det (cid:100) p (cid:101) (cid:18) I + R N − (cid:0) H ( k ) + J ( k ) + K ( k ) (cid:1) N − − P (cid:19) ≤ ˜ C ( k − k ) ν (4.2)for some ˜ C, ν > k in a sufficiently small neighborhood of k . Next, we define the R and Q dependent algorithm byΓ Q,Rn : Ω R → cl( C )Γ Q,Rn ( U ) := (cid:26) k ∈ G n ∩ Q (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12) det (cid:100) p (cid:101) (cid:0) I + A Rn ( k ) (cid:1)(cid:12)(cid:12) ≤ n ) (cid:27) . Theorem 4.2.
For any U ∈ Ω R , we have that Γ Q,Rn ( U ) → Res( U ) ∩ Q in Hausdorff distance as n → + ∞ .Proof. Let us write A R ( k ) := R N − (cid:0) H ( k )+ J ( k )+ K ( k ) (cid:1) N − − P for notational convenience. Theproof consists of two steps. First we prove that any convergent sequence k n ∈ Γ Q,Rn ( U ) necessarilyconverges to a zero of det (cid:100) p (cid:101) ( I + A R ), and second we prove that for every zero k of det (cid:100) p (cid:101) ( I + A R )there exists a sequence k n ∈ Γ Q,Rn ( U ) converging to k . Together, these two facts imply Hausdorffconvergence.Let us begin with step 1. Assume that k n ∈ Γ Q,Rn ( U ) converges to some ˜ k ∈ C . By the definitionof Γ Q,Rn , we have | det( I + A Rn ( k n )) | ≤ n ) . Hence, by Lemma 4.1 we have that | det (cid:100) p (cid:101) ( I + A R ( k n )) | ≤ n ) + Cn − + (cid:100) p (cid:101) → n → + ∞ . Hence we havedet (cid:100) p (cid:101) ( I + A R (˜ k )) = lim n → + ∞ det( I + A R ( k n )) = 0 . Conversely, suppose that det( I + A R ( k )) = 0 for some k ∈ Q . Then there exists a sequence k n ∈ G n ∩ Q such that | k − k n | ≤ n . Hence, by (4.2) we have | det (cid:100) p (cid:101) ( I + A R ( k n )) | ≤ ˜ C n ν and employing Lemma 4.1 again, we obtain | det (cid:100) p (cid:101) ( I + A Rn ( k n )) | ≤ ˜ C n ν + Cn − + (cid:100) p (cid:101) . (4.3)Clearly, for large enough n , the right hand side of (4.3) will be less than n ) , hence k n ∈ Γ Q,Rn ( U )for n large enough. This completes the proof. (cid:3) It remains to extend our argument from a single compact set Q to the entire complex plane. Thisis done via a diagonal-type argument. We choose a tiling of C − , where we start with a rectangle Q = (cid:8) z ∈ C − (cid:12)(cid:12) | Re( z ) | ≤ , | Im( z ) + | ≤ (cid:9) and then add rectangles in a counterclockwise spiralmanner as shown in Figure 4. Note again that each individual rectangle Q i is well separated from { k | k ∈ σ ( H D ) } ⊂ R . Next, we define our algorithm as follows. We letΓ R ( U ) := Γ Q ,R ( U )Γ R ( U ) := Γ Q ,R ( U ) ∪ Γ Q ,R ( U )Γ R ( U ) := Γ Q ,R ( U ) ∪ Γ Q ,R ( U ) ∪ Γ Q ,R ( U )...Γ Rn ( U ) := n (cid:91) j =1 Γ Q j ,Rn ( U ) . Re z Im zQ Q Q Q Q Q Q Q Q · · ·· · · Figure 4.
Tiling of the lower half plane
Theorem 4.2 ensures that Γ Rn ( U ) ∩ B → Res( U ) ∩ B as n → + ∞ in Hausdorff sense for each compactsubset B of C − . The proof of Theorem 2.1 now follows from the characterization of the Attouch-Wetsmetric appearing in (1.4).4.2. Proof of Theorem 1.9.
It remains to extend the above algorithm from Ω R to Ω. In thissection we will define a new algorithm, Γ n , based on the algorithm Γ Rn constructed above. The maindifficulty consists in “locating U ” by testing only finitely many points.The algorithm below uses two radii, R and r . R is used for running Γ Rn , while r is successivelyincreased to search for hidden components of U . The following pseudocode gives the definition ofΓ n . n = 0 : Let U ∈ Ω and initialise R = r = 1. n > : Consider the lattice L n = n − Z ∩ B r (recall the choice h ∼ n − before eq. (4.1)). Forevery j ∈ L n , test whether j ∈ U . For all j ∈ L n ∩ U , compute | j | . • If | j | ≤ R − j ∈ L n ∩ U , define Γ n ( U ) := Γ Rn ( U ∩ B R − ), increment r by 1 andproceed to n + 1. • Else, increment r by 1, set R := r and repeat the current step.This process generates increasing sequences, that we denote by R n and r n , and defines an algorithmΓ n : Ω → cl( C ). Note that r n diverges to + ∞ , because it gets incremented by at least 1 in everystep. Lemma 4.3.
The sequence { R n } n ∈ N is eventually constant and if N > is such that R n = R N forall n > N , one has U ⊂ B R N − Proof.
The fact that { R n } n ∈ N is eventually constant follows immediately from the boundedness of U and the above pseudocode. Now let N ∈ N be as in the assertion. Assume for contradiction that U (cid:42) B R N − . Then, there exists x ∈ U with | x | > R N −
1, and, since U \ B R N − is open, thereexists ε > B ε ( x ) ⊂ U \ B R N − . But then, as soon as n > N is large enough such that n − < ε and r n > | x | , there would exist j ∈ L n with j ∈ B ε ( x ). According to the pseudocode, then, R n would be incremented by 1, contradicting the fact that R n = R N for all n > N . (cid:3) By definition, the output of the algorithm Γ n is Γ R n n ( U ∩ B R n − ). Since R n is eventually constant,there is N ∈ N such that R n ≡ R N for all n ≥ N . Hence we havelim n → + ∞ Γ n ( U ) = lim n → + ∞ Γ R n n ( U ∩ B R n − )= lim n → + ∞ Γ R N n ( U ∩ B R N − )= lim n → + ∞ Γ R N n ( U )= Res( U ) , OMPUTING THE SOUND OF THE SEA IN A SEASHELL 17 where the second and third lines follow from Lemma 4.3 and the last line follows from the proof ofTheorem 2.1. This completes the proof of Theorem 1.9.5.
Numerical Results
Although the algorithm described in Sections 3 and 4 was never designed for computationalefficiency, we show in this section that it is possible to obtain surprisingly good numerical resultswith a MATLAB implementation which differs only in two minor respects:(a) The decomposition M inner = M inner , + K obtained in Lemma 2.2 introduces artificial singularitiesat the spectrum of H D (the Dirichlet Laplacian introduced in eq. (2.5)), which lead to numericalinstabilities. For this reason, rather than approximating the solution v α of (3.4), we directlyimplemented a finite element approximation of the solution u of (2.3). The matrix elements of M inner can then be calculated via Green’s formula, similarly to (3.2).(b) For reasons of performance, instead of using the triangulation procedure outlined in Section 3.1,we use the meshing tool Distmesh , cf. [26], to triangulate our domain.5.1.
A circular resonator chamber.
The first geometry we consider is that of a circular resonatorchamber, connected to open space by a narrow channel. Figure 5 shows the triangulation of thedomain B R \ U . The specific parameter values in this example are: outer radius of the DtN sphere R = 3, outer radius of the resonator chamber r outer = 2, inner radius of the resonator chamber r inner = 1 .
8, width of the opening d = 1 .
3, meshing parameter h = 0 . d r inner Figure 5.
Example triangulation of a circular Helmholtz resonator used in our implementation.
Figure 6 shows a logarithmic contour plot of the computed determinant | det( I + A n ( k )) | in thecomplex plane for n = 20. For the sake of comparison we included the Dirichlet eigenvalues of a closed resonator chamber in the image (red dots). One can see that next to each of the Dirichleteigenvalues there is zero of the determinant, as we would expect. . . . . . . . . − . − . Re( k ) I m ( k ) − − l og ( | d e t( I + A n ( k )) | ) Figure 6.
Logarithmic contour plot of | det( I + A n ( k )) | for n = 20 and opening width d = 1 . The resonance effects can in fact be seen on the level of the associated PDE. Figure 7 showsthe finite element approximation of the solution u of (2.3) with right hand side φ = e α , once for k far away from a resonance, and once for k near a resonance. As can be seen, when k is near theresonance the wave penetrates into the chamber, whereas when k is far from a resonance the solutionappears to be nearly 0 inside. − − − . . − − − . . Figure 7.
Solution of (2.3) with boundary values φ = e α for α = 5. Left: k = 1 . k = 2 . − . Finally, our results show what happens to the locations of the resonances when the opening width d is varied. Figure 8 shows the contour plot of | det( I + A n ( k )) | for an opening width of d = 1 .
0. Onecan clearly see that the zeros of det( I + A n ( k )) move closer to the real axis and become narrower,as is expected from the abstract theory (cf. [19, Ch. 23.3]). OMPUTING THE SOUND OF THE SEA IN A SEASHELL 191 .
32 1 .
33 1 . − − · − Re( k ) I m ( k ) .
08 2 . . − − . . · − Re( k ) I m ( k ) . . . . . . . . − . − . Re( k ) I m ( k ) − − − − l og ( | d e t( I + A n ( k )) | ) Figure 8.
Top: Logarithmic contour plot of | det( I + A n ( k )) | for n = 20 and openingwidth d = 1 . Table 1 summarizes the values of the first three eigenvalues and associated resonances in bothcases (that is, with openings d = 1 . d = 1 . Eigenvalue Resonance d = 1 . d = 1 . first 1 . . − . . − . . . − . . − . . . − . . − . Table 1.
Approximate numerical values of resonances and eigenvalues corresponding tothe contour plots in Figures 6 and 8 respectively.
Four circular Neumann holes.
In many applications it is necessary to consider not onlyDirichlet, but also Neumann boundary conditions. An important real life situation in which this isthe case is presented by resonance effects of water waves caused by submerged objects. A concretequestion which has entailed a vast body of scientific literature is this: do the pillars of offshorestructures, such as oil drilling platforms, introduce a positive interference of water waves that coulddamage the structure itself? At what frequencies are such positive interferences expected to occur?(cf. [29, 22, 13].) This section contains a short study of the latter question, together with somenumerical experiments.
Figure 9.
Triangulation of B R \ U for U consisting of four circular Neumann holes. The approach outlined in Sections 2, 3 and 4 can be adapted to Neumann boundary conditionson the boundary of the obstacle U , given a set of stronger assumptions. In this section, we willbriefly outline the ideas and present some numerical results.Let now Ω N denote the class of bounded domains U (cid:98) B R − with C boundary, on which weconsider the Neumann Laplacians H := − ∆ N on L ( R \ U ) . In addition, let us assume that each obstacle U comes with an associated regular parametrization γ U : [ a U , b U ] ∪ [ a U , b U ] · · · ∪ [ a Un U , b Un U ] → R (parametrized by its arc length) of its boundary. DefineΛ N := Λ ∪ (cid:40) U (cid:55)→ γ U ( x ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x ∈ n U (cid:91) k =1 [ a Uk , b Uk ] (cid:41) , where Λ was defined in (1.3). Under this additional hypotheses, Sections 2, 3 and 4 (with modifica-tions) show that SCI(Ω N , Λ N , Res( · ) , (cl( C ) , d AW )) = 1 . Indeed, the discussions in Sections 2 and 4 remain true with trivial changes for Neumann conditionson ∂U . The only major difference is in Section 3.1, where the finite element approximation of v α isconstructed. A version of Prop. 3.1 in the Neumann case can be obtained from our new assumptionsas follows.(i) For each U , choose a uniform discretization of the intervals [ a Uk , b Uk ]. This induces an orienteddiscretization of ∂U via γ U .(ii) Using well known meshing algorithms (e.g. [7, Ch. 5]), construct a fitted mesh of B R \ U basedon the boundary discretization from (i).(iii) Apply [2, Th. 4.1] to obtain the desired FEM error estimate.We implemented the algorithm thus obtained (with the caveats (a), (b) mentioned at the beginningof the section) for an obstacle consisting of four circular holes, as shown in Figure 9.This geometry had previously been studied in the context of water waves interacting with acircular array of cylinders [13]. Figure 10 shows the output of our algorithm for four circular holesof radius 0.6, situated at the corners of a square of edge length 2. OMPUTING THE SOUND OF THE SEA IN A SEASHELL 210 0 . . . . . − . − . − . Re( k ) I m ( k ) − l og ( | d e t( I + A n ( k )) | ) Figure 10.
Logarithmic contour plot of | det( I + A n ( k )) | for n = 30 in a strip below thereal axis The heat map in Figure 10 suggests 3 resonances located approximately at the following values:
Resonance first 0 . − . . − . . − . Table 2.
Approximate numerical values of resonances corresponding to the contour plotin Figure 10.
A comparison with [13, Fig. 8] shows good agreement.
References [1] H. Ammari, K. Imeri, and W. Wu. A mathematical framework for tunable metasurfaces. Part I.
Asymptot. Anal. ,114(3-4):129–179, 2019.[2] J. W. Barrett and C. M. Elliott. A finite-element method for solving elliptic equations with Neumann data on acurved boundary using unfitted meshes.
IMA J. Numer. Anal. , 4(3):309–325, 1984.[3] G. Beer.
Topologies on closed and closed convex sets , volume 268 of
Mathematics and its Applications . KluwerAcademic Publishers Group, Dordrecht, 1993.[4] J. Ben-Artzi, M. J. Colbrook, A. C. Hansen, O. Nevanlinna, and M. Seidel. Computing spectra – On the solvabilitycomplexity index hierarchy and towers of algorithms. arXiv e-prints, 1508.03280 , Aug 2015.[5] J. Ben-Artzi, A. C. Hansen, O. Nevanlinna, and M. Seidel. New barriers in complexity theory: On the solvabilitycomplexity index and the towers of algorithms.
Comptes Rendus Math. , 353(10):931–936, Oct 2015.[6] J. Ben-Artzi, M. Marletta, and F. R¨osler. Computing Scattering Resonances. arXiv e-prints, 2006.03368 , June2020.[7] G. Borouchaki. Delaunay triangulation and meshing: Application to finite elements.
Hermes, Paris , 1998.[8] H. Brezis.
Functional analysis, Sobolev spaces and partial differential equations . Universitext. Springer, NewYork, 2011.[9] P. G. Ciarlet.
The finite element method for elliptic problems , volume 40 of
Classics in Applied Mathematics .Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002. Reprint of the 1978 original[North-Holland, Amsterdam].[10] P. Doyle and C. McMullen. Solving the quintic by iteration.
Acta Math. , 163(3-4):151–180, 1989.[11] N. Dunford and J. T. Schwartz.
Linear operators. Part II . Wiley Classics Library. John Wiley & Sons, Inc., NewYork, 1988. Spectral theory. Selfadjoint operators in Hilbert space, With the assistance of William G. Bade andRobert G. Bartle, Reprint of the 1963 original, A Wiley-Interscience Publication.[12] S. Dyatlov and M. Zworski.
Mathematical theory of scattering resonances , volume 200 of
Graduate Studies inMathematics . American Mathematical Society, Providence, RI, 2019.[13] D. Evans and R. Porter. Near-trapping of waves by circular arrays of vertical cylinders.
Applied Ocean Research ,19(2):83–99, 1997.[14] L. D. Faddeev and B. S. Pavlov. Scattering theory and automorphic functions. In
Fifty years of mathematicalphysics , volume 2 of
World Sci. Ser. 21st Century Math. , pages 75–101. World Sci. Publ., Hackensack, NJ, 2016.[15] A. C. Hansen. On the solvability complexity index, the n -pseudospectrum and approximations of spectra ofoperators. J. Amer. Math. Soc. , 24(1):81–124, 2011.[16] E. Harrell and B. Simon. The mathematical theory of resonances whose widths are exponentially small.
DukeMath. J. , 47(4):845–902, 1980. [17] B. Helffer and J. Sj¨ostrand. R´esonances en limite semi-classique.
M´em. Soc. Math. France (N.S.) , (24-25):iv+228,1986.[18] H. Helmholtz. Die Lehre von den Tonempfindungen als physiologische Grundlage f¨ur die Theorie der Musik (Onthe sensations of tone as a physiological basis for the theory of music).
Braunschweig, Germany: Friedrich Viewegund Sohn , 1863.[19] P. Hislop and I. Sigal.
Introduction to spectral theory , volume 113 of
Applied Mathematical Sciences . Springer-Verlag, New York, 1996.[20] P. D. Lax and R. S. Phillips.
Scattering theory , volume 26 of
Pure and Applied Mathematics . Academic Press,Inc., Boston, MA, second edition, 1989.[21] N. N. Lebedev.
Special functions and their applications . Revised English edition. Translated and edited by RichardA. Silverman. Prentice-Hall, Inc., Englewood Cliffs, N.J., 1965.[22] C. M. Linton and D. V. Evans. The interaction of waves with arrays of vertical circular cylinders.
J. Fluid Mech. ,215:549–569, 1990.[23] A. Martinez and L. N´ed´elec. Optimal lower bound of the resonance widths for the Helmholtz resonator.
Ann.Henri Poincar´e , 17(3):645–672, 2016.[24] W. McLean.
Strongly elliptic systems and boundary integral equations . Cambridge University Press, Cambridge,2000.[25] R. G. Newton.
Scattering theory of waves and particles . McGraw-Hill Book Co., New York-Toronto, Ont.-London,1966.[26] P.-O. Persson and G. Strang. A simple mesh generator in Matlab.
SIAM Rev. , 46(2):329–345, 2004.[27] B. Simon. Notes on infinite determinants of Hilbert space operators.
Advances in Math. , 24(3):244–273, 1977.[28] J. Sj¨ostrand and M. Zworski. Complex scaling and the distribution of scattering poles.
J. Amer. Math. Soc. ,4(4):729–769, 1991.[29] B. H. Spring and P. L. Monkmeyer. Interaction of plane waves with vertical cylinders.
Coastal EngineeringProceedings , 1(14):107, Jan. 1974.[30] P. Stefanov. Approximating resonances with the complex absorbing potential method.
Comm. Partial DifferentialEquations , 30(10-12):1843–1862, 2005.[31] D. R. Yafaev.
Mathematical scattering theory , volume 158 of
Mathematical Surveys and Monographs . AmericanMathematical Society, Providence, RI, 2010. Analytic theory.
E-mail address : [email protected]
E-mail address : [email protected]