Efficient resonance computations for Helmholtz problems based on a Dirichlet-to-Neumann map
Juan Carlos Araujo-Cabarcas, Christian Engstrom, Elias Jarlebring
EEfficient resonance computations for Helmholtz problems based on aDirichlet-to-Neumann map
Juan Carlos Araujo-Cabarcas , Christian Engstr¨om , Elias Jarlebring Abstract
We present an efficient procedure for computing resonances and resonant modes of Helmholtzproblems posed in exterior domains. The problem is formulated as a nonlinear eigenvalue problem(NEP), where the nonlinearity arises from the use of a Dirichlet-to-Neumann map, which accountsfor modeling unbounded domains. We consider a variational formulation and show that the spec-trum consists of isolated eigenvalues of finite multiplicity that only can accumulate at infinity. Theproposed method is based on a high order finite element discretization combined with a special-ization of the Tensor Infinite Arnoldi method. Using Toeplitz matrices, we show how to specializethis method to our specific structure. In particular we introduce a pole cancellation technique inorder to increase the radius of convergence for computation of eigenvalues that lie close to the polesof the matrix-valued function. The solution scheme can be applied to multiple resonators with avarying refractive index that is not necessarily piecewise constant. We present two test cases toshow stability, performance and numerical accuracy of the method. In particular the use of a highorder finite element discretization together with TIAR results in an efficient and reliable methodto compute resonances.
Keywords:
Nonlinear eigenvalue problems, Helmholtz problem, scattering resonances, Dirichlet-to-Neumannmap, Arnoldi’s method, Bessel functions, Matrix functions
1. Introduction
In this paper we consider the numerical approximation of resonances in an open system, wherethe solutions satisfy the Helmholtz equation for a given refractive index η ( x ). In general, resonancesof an operator are defined as poles of the resolvent operator taken in a particular generalized sense[53, 38]. For Helmholtz equation Lenoir et al. [34] have shown that resonances are solutions to anonlinear eigenvalue problem (NEP) with a Dirichlet-to-Neumann (DtN) map G ( λ ) on an artificialboundary Γ. The pair ( u, λ ) is a scattering resonance pair if∆ u + λ η u = 0 in Ω ,∂u∂n = G ( λ ) u on Γ , (1) Department of Mathematics and Mathematical Statistics, Ume˚a University, MIT-Huset, 90187 Ume˚a, Sweden Department of Mathematics, Royal Institute of Technology (KTH), Stockholm, SeRC Swedish e-Science ResearchCenter
Preprint submitted to arXiv September 5, 2018 a r X i v : . [ m a t h . NA ] J un here ∂u/∂n is the normal derivative and the non-negative function η − λ theDtN operator G ( λ ) (which we formalize in Section 2) depends in a nonlinear way on the eigenvalue λ . In this work we present a new computational approach for approximating resonances of (1) inan efficient and accurate way.We present in Section 2.1 a variational formulation of the PDE (1) with the DtN map and showthat the spectrum of the operator function consists of isolated eigenvalues of finite multiplicity,which accumulate only at infinity. This variational formulation is the base for the finite elementmethod (FEM) in Section 3. In the FEM-implementation we discretize the operator function usinghigh order Gauss-Lobatto shape functions and apply the p -version of the finite element method.The lower part of the spectrum of the operator function is then well approximated by the matrixfunction. Then very accurate approximations of the scattering resonance pairs are obtained if thenonlinear matrix eigenvalue problem can be solved accurately.The considered NEP is of the type: find λ ∈ C in an open subset of the complex plane and anon-zero u ∈ C n such that T ( λ ) u = 0 . (2)In our case T is meromorphic in C , with poles in the region of interest defined as scaled roots ofHankel functions. Many numerical methods for the NEP (2) have been developed in the numericallinear algebra community, in particular when T is holomorphic in a large domain. Since the NEPwith an arbitrary normalization is a system of nonlinear equations, Newton’s method can be appliedand considerably improved, e.g., by using block variations that can compute several eigenpairssimultaneously [33]. However, Newton-type methods bear the danger that some eigenvalues closeto a given target could be missed.There are also generalizations of successful methods for linear eigenvalue problems, e.g., thenonlinear Arnoldi method [50], the Jacobi-Davidson method [11] and LOBPCG [47]. Numericalmethods for NEP can be based on contour integrals of the generalized resolvent T − ( λ ) [2, 12, 18,8, 25]. The DtN-map is the source of the complicated nonlinearity in the eigenvalue parameter λ .Several other problems have been approached with artificial boundary conditions and NEPs; see,e.g., the NEP arising in the modeling of an electromagnetic cavity [35], the model of an optical fiberin [31] and bi-periodic slabs [48]. There are also several approaches leading to NEPs developed inthe context of photonic crystals [46, 19, 32, 16]. To our knowledge, none of the methods developedin those papers have been adapted to resonance problems of our type.Our methods belong to a class of methods which can be interpreted as Krylov methods, eitherfor an infinite-dimensional problem, or as a dynamically increasing companion linearization of anapproximation of the problem [28, 9, 21]. For recent developments and problems see [51] and[10]. Our approach is based on the infinite Arnoldi (IAR) method [28]. In particular, we adaptthe variant with a tensor representation of the basis presented in [27], called the tensor infiniteArnoldi method (TIAR). This method is designed to find all eigenvalues close to a given shift,where the radius of convergence depends on properties of the matrix-valued function (2). Theinfinite Arnoldi method (both IAR and TIAR) requires access a procedure to compute a quantityinvolving the derivatives of the problem. In Section 4 we derive an algorithm to compute thenecessary quantities for our NEP which contains nonlinearities expressed as quotients of Hankelfunctions. In order to improve the convergence of the method, we introduce a pole cancellationtechnique which transforms the problem by removing poles.In Section 5 we provide a characterization of the performance of our approach. In particular,2e show that the new infinite Arnoldi method together with a p -FEM strategy is an efficient andreliable tool for resonance calculations with the DtN map. a Ω a Γ a Ω Ω Ω Figure 1: Example geometry of resonators Ω i , i = 1 , ,
2. Background and preliminaries
Results for the problem (1) can be found in a considerable amount of literature; see [42, 34, 15]and references therein. For Im λ > λ < a ⊂ R be an open disk of radius a and boundary Γ a . Assume η ∈ L ∞ (Ω a ) and thatthe non-negative function η − a . A schematic setup of anexample is illustrated in Figure 1. The resonance problem restricted to Ω a is formally to findnon-trivial solutions ( u, λ ) such that (1) holds, wherethe DtN operator on the circle Γ a has the explicit form G ( λ ) u := 12 π ∞ (cid:88) ν = −∞ g ν ( λ ) e iνθ (cid:90) π u ( a, θ (cid:48) ) e − iνθ (cid:48) dθ (cid:48) , (3)where g ν ( λ ) := λ H (cid:48) ν ( λa ) H ν ( λa ) (4)and G ( λ ) : H / (Γ a ) → H − / (Γ a ) is bounded [34]. In the following, we identify the dual pairing (cid:104)· , ·(cid:105) H − / (Γ a ) × H / (Γ a ) with the L -inner product ( · , · ) Γ a over Γ a . The theory presented in [34] canwith minor changes be used in the present case to derive properties of a variational formulation ofthe problem. Let S denote the union of the set of zeros of H ν ( λa ) , ν ∈ Z . The eigenvalues of (1) aredetermined by the following variational problem: Find u ∈ H (Ω a ) \ { } and λ ∈ D := C \ { R − ∪ S } such that for all v ∈ H (Ω a )( ∇ u, ∇ v ) Ω a − λ ( η u, v ) Ω a − ( G ( λ ) u, v ) Γ a = 0 , (5)3here ( u, v ) Ω a := (cid:82) Ω a u ¯ v dx , ( ∇ u, ∇ v ) Ω a := (cid:82) Ω a ∇ u · ∇ ¯ v dx , and( G u, v ) Γ a = ∞ (cid:88) ν = −∞ λa H (cid:48) ν ( λa ) H ν ( λa ) ˆ u ν ¯ˆ v ν , ˆ ϕ ν = 1 √ π (cid:90) π ϕ ( a, θ (cid:48) ) e − iνθ (cid:48) dθ (cid:48) . (6)Let ( u, v ) := ( ∇ u, ∇ v ) Ω a + ( η u, v ) Ω a . Then the norm || u || := (cid:112) ( u, u ) is equivalent to thestandard norm on H (Ω a ). Define the operator F ( λ ) : H (Ω a ) → H (Ω a ) by( F ( λ ) u, v ) := ( λ + 1)( η u, v ) Ω a + ( G ( λ ) u, v ) Γ a . (7)An operator formulation of (5) is to find u ∈ H (Ω a ) \ { } and λ ∈ D such that( I − F ( λ )) u = 0 . (8)Let L ( H (Ω a )) denote the space of bounded linear operators on H (Ω a ). An operator A ∈L ( H (Ω a )) is Fredholm if it has finite-dimensional kernel ker A and cokernal coker A . The indexof a Fredholm operator is ind A := dim ker A − dim coker A . The essential spectrum of A denoted σ ess ( A ) is the set of of complex numbers such that A − λI is not Fredholm. Let σ ∞ ( A ) ⊂ σ ess ( A )denote the set of eigenvalues of infinite multiplicity.Let G ν max ( λ ) denote the operator (3) truncated after | ν | = ν max and let I − F ν max ( λ ), λ ∈ D denote the operator defined by(( I − F ν max ( λ )) u, v ) := ( u, v ) − ( λ + 1)( η u, v ) Ω a − ( G ν max ( λ ) u, v ) Γ a . (9)The operator function with a truncated DtN map (9) is the base for the numerical method andthe following proposition shows that the spectrum is discrete with no finite accumulation point inthe right-half plane. Proposition 2.1.
Let ν max ∈ { , , . . . , ∞} . Then the H (Ω a ) -spectrum of the operator valuedfunction I −F ν max : D → L ( H (Ω a )) consists of isolated eigenvalues of finite multiplicity. Moreover,these eigenvalues can not have a finite accumulation point in the right-half plane.Proof. The function F ν max : D → L ( H (Ω a )) is holomorphic [34, prop 4] and I − F ν max ( λ ) is aFredholm operator of index zero [34, prop 5]. Moreover, we have ia H (cid:48) ν ( ia ) H ν ( ia ) = a K (cid:48) ν ( a ) K ν ( a ) , (10)where K ν ( a ) = π i ν +1 H ν ( ia ), ν ∈ Z are the modified Bessel functions [52]. The expression (10) isnegative since K ν ( a ) > K (cid:48) ν ( a ) < a > G ν max ( i ) u, u ) Γ a < u, u ) − ( F ν max ( i ) u, u ) = ( u, u ) − ( G ν max ( i ) u, u ) Γ a ≥ (cid:107) u (cid:107) , which shows that the resolvent set of I − F is non-empty. Hence, from the analytical Fredholmtheorem follows that the spectrum is discrete and all eigenvalues are of finite multiplicity [20,Theorem 5.1]. In the right half-plane the eigenvalues can therefore only accumulate at the poles.The location of the poles is a function of a , but the spectrum is for large enough a independent of a , which implies that we have no accumulation in the right half-plane.4n Section 4.4, we improve convergence of the infinite Arnoldi method by multiplying theoriginal matrix-valued function with a polynomial. In Proposition 2.2, we prove basic propertiesof the underlaying operator-valued function. This proposition shows that the canceled poles willbe in the essential spectrum of the modified operator. Proposition 2.2.
Take
B ⊂ D in the right-half plane such that F has exactly one pole z in B .Then is ˜ T ( λ ) : H (Ω a ) → H (Ω a ) , ˜ T ( λ ) := ( λ − z ) I − ˜ F ( λ ) , ˜ F ( λ ) := ( λ − z ) F ( λ ) (11) holomorphic in B , { z } = σ ess ( ˜ T ) = σ ∞ ( ˜ T ) , and ˜ T ( z ) is of rank one.Proof. Take without loss of generality a = 1. Since all zeros of H ν ( λ ) are simple [1, p 370], wehave by definition H ν ( λ ) = ( λ − z ) f ( λ ) with f holomorphic and f ( z ) (cid:54) = 0. Hence, H (cid:48) ν ( λ ) H ν ( λ ) = 1 λ − z + f (cid:48) ( λ ) f ( λ ) , which implies that ˜ T is holomorphic in B . Since H (Ω a ) is infinite dimensional, it follows that˜ T ( z ) is not Fredholm. Assume H ν z ( z ) = 0, H ν ( z ) (cid:54) = 0, ν (cid:54) = ν z . Then ( ˜ T ( z ) u, v ) = − z ˆ u ν z ¯ˆ v ν z ,which implies that ˜ T ( z ) and ˜ T ( z ) ∗ are of rank one and hence z ∈ σ ∞ ( ˜ T ). Remark Proposition 2.2 shows that a pole z after multiplication by λ − z is an eigenvalueof infinite multiplicity of the new operator function. This change of the spectral properties alsoapply more generally. Let H denote an infinite dimensional Hilbert space and assume that F is a L ( H ) -valued finitely meromorphic function at z . Hence, F ( λ ) = ∞ (cid:88) m = − s ( λ − z ) m F m , where F − s , F − s +1 , . . . , F − are finite rank. Define the bounded operator function ˜ T ( λ ) := ( λ − z ) s I − ˜ F ( λ ) , ˜ F ( λ ) := ( λ − z ) s F ( λ ) . (12) Then z ∈ σ ∞ ( ˜ T ) , since ˜ T ( z ) = − F − s is finite rank. However if F − s is not finite rank, z will ingeneral not be an eigenvalue of infinite multiplicity; see for example [18].
3. Discretization with the finite element method
The disk Ω a , depicted in Figure 1, is covered with a regular and quasi uniform finite elementmesh T consisting of quadrilateral elements { K i } Ni =1 . Let ρ i be the diameter of the largest ballcontained in K i and denote by h the maximum mesh size h := max ρ i . Let P p denote the space ofpolynomials on R of degree ≤ p and set γ := { h, p } [44, Ch 4]. We define the finite element space S γ (Ω a ) := { u ∈ H (Ω a ) : u | K i ∈ P p ( K i ) for K i ∈ T } , and N := dim( S γ (Ω a )) [5]. Furthermore, allour computations are done in the approximated domain Ω γa by using curvilinear elements followingstandard procedures [5].From (5) we state the finite element problem: Find u γ ∈ S γ (Ω γa ) \ { } and λ γ ∈ D such that( ∇ u γ , ∇ v ) Ω γa − λ γ ( η u γ , v ) Ω γa − ( G ν max ( λ γ ) u γ , v ) Γ γa = 0 for all v ∈ S γ (Ω γa ) . (13)5e showed in Section 2.1 that the analytic operator function I − F ν max with a truncated DtNmap is Fredholm-valued on D and that the resolvent set is non-empty. Hence, all eigenvaluesare isolated and of finite multiplicity. Moreover, from Karma [29, 30] follows that any sequence { λ γ } , dim S γ (Ω γa ) → ∞ of approximative eigenvalues of (13) converges to an eigenvalue λ of(5). The eigenfunctions are in H (Ω a ) and are piecewise analytic if the interfaces are analyticcurves. Optimal convergence is expected under the assumptions that all interfaces are resolved bycurvilinear cells and the eigenvalues are semi-simple. Hence, exponential convergence is expectedwith p -FEM and optimal converge rates are expected with h -FEM [4, 40]. Let { ϕ , . . . , ϕ N } be a basis of S γ (Ω γa ). Then u γ ∈ S γ (Ω γa ) can be written in the form u γ = N (cid:88) j =1 ξ j ϕ j (14)and the entries in the finite element matrices are A ij = ( ∇ ϕ j , ∇ ϕ i ) Ω γa , M ij = ( η ϕ j , ϕ i ) Ω γa , G ij ( λ ) = l (cid:88) ν = − l λa H (cid:48) ν ( λa ) H ν ( λa ) Q νij (15)with Q νij = ˆ ϕ νj ¯ˆ ϕ νi , ˆ ϕ νj = 1 √ π (cid:90) π ϕ j ( a, θ (cid:48) ) e − iνθ (cid:48) dθ (cid:48) , i, j = 1 , . . . , N. (16)The nonlinear matrix eigenvalue problem is then: Find the eigenpars ( λ, ξ ) ∈ D × C N \ { } suchthat T ( λ ) ξ := (cid:0) A − λ M − G ( λ ) (cid:1) ξ = 0 . (17)In the assembly routine, we only store the vectors q νj = ˆ ϕ νj and then compute the matrices as Q ν = q ν ⊗ ¯ q ν . Notice that q νj = 0 for nodes x j such that supp( ϕ j ) ∩ Γ a = ∅ . If the ordering of nodesis such that the N a boundary nodes are placed first, then the Q ν matrices have a dense upper-leftblock of size N a × N a . This lost of sparsity is taken into account in the assembly routine, wherewe allow extra entries in the sparsitty pattern of the FE matrices.It is a standard technique to use Gauss-Legendre quadratures for the evaluation of integrals inFE. However, the trace integral (16) requires that the number of quadrature points is increasedlinearly with ν , because the integrand contains e iνθ ( x , x ) .All numerical experiments have been carried out using the finite element library deal.II [7] withGauss-Lobatto shape functions [45, Sec. 1.2.3]. For fast assembly and computations with complexnumbers the package PETSc [6] is used.The computational platform was provided by the High-Performance Computing Center North(HPC2N) at Ume˚a University, and all experiments were run on the distributed memory systemAbisko. The jobs were run in serial on an exclusive node: during the process, no other jobs wererunning on the same node. Node specifications: four AMD Opteron 6238 processors with a totalof 48 cores per node. 6 . Specialization of the infinite Arnoldi method Our algorithm for solving the NEP (17) is based on the tensor infinite Arnoldi method (TIAR)[27, Algorithm 2], which is an improvement of the infinite Arnoldi method (IAR) [28]. The versionof the infinite Arnoldi method considered here can be viewed as the standard Arnoldi methodapplied to the companion linearization arising from a Taylor expansion of an analytic matrix-valued function T [28]. More precisely, it can be derived from a particular companion linearizationof the Taylor expansion of T and Arnoldi’s method for eigenvalue problems. The particular choiceof companion linearization provides a structure such that the truncation parameter in a certainsense can be increased to infinity and the algorithm is therefore equivalent to Arnoldi’s methodapplied to an infinite matrix.Hence, the algorithm generates a Hessenberg matrix and the eigenvalues of this matrix corre-spond to eigenvalue approximations of the NEP.TIAR is a variant of IAR where a tensor representation of the basis of IAR is used, whichresults in an algorithm of the same complexity as IAR but it requires less memory. Moreover,the usage of tensors in TIAR makes it considerably more efficient than IAR for certain types ofproblems [27] on modern computer architectures. We derive below a new version of TIAR adaptedto the special structure of (17). All of the variants of the infinite Arnoldi method require (in some way) access to derivatives.In TIAR we need a procedure to compute x = − T ( µ ) − (cid:32) k (cid:88) i =1 T ( i ) ( µ ) x i (cid:33) , (18)where k is the iteration count, and we compute eigenvalues close to a target µ .The matrix T ( µ ) is independent of k and we compute an LU-factorization of T ( µ ) beforestarting the iterations in TIAR. Without loss of generality, we write T of the form T ( λ ) = N (cid:88) j =1 A j f j ( λ ) . (19)Then, we express the NEP (17) in the form (19) with A = A, f ( λ ) = 1 (20a) A = M, f ( λ ) = − λ (20b)and for j = 3 , . . . , ν max + 3, A j = − aq j − ν max − q Tj − ν max − (21a) f j ( λ ) = g j − ν max − ( λ ) = H (cid:48) j − ν max − ( aλ ) H j − ν max − ( aλ ) λ. (21b)Hence, (18) for the problem stated in (17) can be written as x = − T ( µ ) − ν max +3 (cid:88) j =1 A j k (cid:88) i =1 x i f ( i ) j ( µ ) . (22)7n order to specialize TIAR to the considered NEP a procedure to compute the derivatives in (22)is required. The nonlinearities (21b) are products and quotients of analytic functions of the form g ( λ ) = h ( λ ) h ( λ ) p ( λ ) , (23)where h and h are analytic functions and p a polynomial. Formula (22) includes (high-order)derivatives of such scalar functions. There are high-order product and quotient rules for differenti-ation, which are explicit and lead to formulas for the high-order derivatives of (23). However, thedirect application of the high-order product and quotient rules, i.e., the general Leibniz rule, aresomewhat unsatisfactory in terms of computation time. Here we propose to use a formula involvingToeplitz matrices, and compute directly a given number of derivatives using only matrix-vectoroperations.The derivation of our Toeplitz matrix computational formula for the derivatives, can be doneby comparing terms in a Taylor expansion of g ( λ ). For our purposes it is more natural to usematrix functions in the sense of [23]. Let J µ denote a Jordan matrix with eigenvalue µ , J µ = µ
1. . . . . .. . . 1 µ . (24) The matrix function of a Jordan matrix is a triangular Toeplitz matrix [23, Definition 1.2] con-taining scaled derivatives of the scalar function, f ( J Tµ ) = f ( J µ ) T = f ( µ )0! ... . . . f ( p − ( µ )( p − · · · f ( µ )0! . (25)Therefore, the derivatives of a function for which the corresponding matrix function is availablecan be computed with the formula f ( µ ) f (cid:48) ( µ )... f ( p − ( µ ) = diag( u ) f ( J Tµ ) e , (26)where u T = [0! , . . . , ( p − g ( µ ) g (cid:48) ( µ )... g ( p − ( µ ) = diag( u ) h ( J Tµ ) h ( J Tµ ) − p ( J Tµ ) e . (27)8he TIAR algorithm can be applied to any NEP expressed in the form (19) for which efficientand reliable computation of the corresponding matrix function is available. Such matrix functionrepresentation is analogous to the input for several other methods and software packages, e.g.,the block Newton method [33] and NLEIGS [21]. This representation is suitable for many NEPsincluding the problem considered in this paper (due to the derivative computation specializationprovided here), but not in general for NEPs where the number of terms is large. Completelyanalogous to scalar-valued functions, products and quotients of matrix functions are expressed asmatrix multiplies and inverses (which is formally a consequence of [23, Theorem 1.15]). For the purpose of using formula (27), we require the matrix functions corresponding to (20)and (21b). These matrix functions are f ( S ) = I , f ( S ) = − S , and f j ( S ) = g j − ν max − ( S ) = H (cid:48) j − ν max − ( aS ) H j − ν max − ( aS ) − S, j = 3 , . . . , ν max + 3 . (28)The matrix function (28) could be computed directly if robust and efficient methods for matrixHankel functions were available. To the best of our knowledge, there are unfortunately no suchspecialized methods.Rather than computing the matrix function of the Hankel function explicitly, we use thatfor applying (27) only the matrix function of a transposed Jordan block is required. Hence, wecompute the Toeplitz matrix (25) consisting of scaled derivatives. Then, the derivatives of theHankel function can be robustly computed as follows. For notational convenience, let for any r ∈ N , H r ( z ) := [ H ( z ) , H ( z ) , . . . , H r − ( z )] T (29a) H − r ( z ) := [ H ( z ) , H − ( z ) , . . . , H − r ( z )] T (29b)= D r H r ( z ) , (29c)where D r = diag(1 , − , , − , . . . ) ∈ R r and (29c) follows from the symmetry of Hankel functions.Using the recursion formulas for Hankel functions, the infinite vector of derivatives of the Hankelfunctions can be written as H (cid:48)∞ ( z ) = B ∞ H ∞ ( z ) (30)where B ∞ is given by the infinite extension of a matrix B r formed by the sum of a Toeplitz matrixand a rank-one matrix, B r = − / − /
2. . . . . . . . .1 / − / / ∈ R r × r . (31)The relation (30) leads to a procedure to compute the derivatives summarized in the followinglemma. 9 emma 2 (Hankel function derivative recursions) . Let the vector of Hankel functions be definedby (29) and the tridiagonal matrix B r in (31) . Then, the k th derivatives of Hankel functions aregiven by d k dz k H r ( az ) = a k [ I r , , . . . , B kr + k H r + k ( az ) (32) and H ( k ) − r ( z ) = D r H ( k ) r ( z ) .Proof. By the chain rule and repeated application of (30) we have for an arbitrary derivative k , d k dz k H ∞ ( az ) = a k H ( k ) ∞ ( az ) = a k B k ∞ H ∞ ( az ) . The matrix B ∞ is tridiagonal. Therefore,[ I r , , . . . ] B k ∞ = [ I r , , . . . , B kk + r , , . . . , . (33)Then (32) follows from (33) with H ( k ) r ( z ) = [ I r , , . . . , H ( k ) ∞ ( z ) = [ I r , , . . . , B k ∞ H ∞ ( z ) . Remark (Alternative ways to compute derivatives of Hankel functions) . There are in principle,many ways to numerically compute derivatives of Hankel functions, e.g., with various discretizationschemes. An advantage of our approach is that it is exact in exact arithmetic, and appears relativelyinsensitive to round-off errors. We describe in Algorithm 1 how the computation can be performedwith only matrix-vector products. This is easily integrated with the pole cancellation techniquedescribed in Section 4.4, and the large number of derivatives required in TIAR.4.4. Pole cancellation and derivative computation algorithm
Lemma 2 does provide a procedure to compute derivatives of the Hankel functions and it canbe directly used to specialize TIAR to (17). We show in Section 5 that the direct application of themethod numerically works well for some regions of the complex plane, but unfortunately not forthe entire complex plane. The infinite Arnoldi method is designed for problems which are analyticin a large domain, and convergence cannot be guaranteed for eigenvalues outside the convergencedisk for the power series expansion at µ . Below, use a transformation of the matrix function toincrease the convergence radius by effectively cancelling poles. A similar holomorphic extensionwas successfully applied to a matrix-valued function that is used to determine surface waves in soilmechanics [49, Section 7.3.1]. Suppose z i , i = 1 , . . . , p are zeros of Hankel functions and define˜ T ( λ ) := ( λ − z ) · · · ( λ − z p ) T ( λ ) . From Proposition 2.2 follows that ˜ T is holomorphic in a domain containing z i , i = 1 , . . . , p . Hence,we have a decomposition analogous to (19) where ˜ T ( λ ) = (cid:80) Ni A i ˜ f i ( λ ), with˜ f i ( λ ) := ( λ − z ) · · · ( λ − z p ) f i ( λ ) . Note that the nonlinear terms of ˜ T are ˜ g m ( λ ) = ( λ − z ) · · · ( λ − z p ) g m ( λ ), m = − ν max , . . . , m ,which are terms of the form (23) with p ( λ ) := ( λ − z ) · · · ( λ − z p ) λ . Therefore, we can still10se formula (27). By construction, the nonlinear (matrix) eigenvalue problem associated with ˜ T has the same solutions as the NEP associated with T , except for possibly λ = z i . Moreover, thesingularity set of ˜ T is the same as T except that z , . . . , z p are not poles of ˜ T . The values λ = z i are not solutions to the NEP T , since z i is a pole. However, Proposition 2.2 shows that the polesof the original operator function are eigenvalues of infinite multiplicity of the operator functionobtained by multiplying with a polynomial. The matrix function ˜ T will then numerically haveseveral additional eigenvalues close to the poles compared to the eigenvalues of T . It is thereforeessential that an estimate of the quality of a computed eigenpair of ˜ T is based on the originalmatrix function T . In Section 5, we estimate the quality of a computed eigenpar by the standardbackward error estimate (consistent e.g. with [36]) for NEPs (cid:107) T ( λ ) v (cid:107) /α ( λ, v ) , with α ( λ, v ) := (cid:107) A (cid:107) + | λ | (cid:107) M (cid:107) + l (cid:88) ν = − l | λ | a | H (cid:48) ν ( aλ ) | / | H ν ( aλ ) |(cid:107) Q ν (cid:107) . We summarize the combination of (27) with the Hankel function derivative computation in Lemma 2in Algorithm 1. In the algorithm we have taken advantage of the structure of the numerator anddenominator in (23), and that the derivatives corresponding to all needed indexes can be computedsimultaneously. The output of the algorithm is the matrix consisting of derivatives of ˜ g i , X = (cid:2) x · · · x ν max +1 (cid:3) = ˜ g (0)0 · · · ˜ g (0) ν max ... ...˜ g ( k max )0 · · · ˜ g ( k max ) ν max . (34)Note that ˜ g i = ˜ g − i such that we can use (34) also for negative index. Remark Efficiency improvements and memory requirements. Some further improvements inthe computation of x were achieved by also using the localization of the vectors q i . More precisely,for i > p + 2 the contribution of A and A do not need to be taken into account. Note that inorder to carry out m steps of TIAR for a NEP of size N , we mainly need to store a complex basismatrix of size N × m . Hence, in our setting the memory required by TIAR is far less demandingthan the storage required for the LU-factorization used for performing the initial step in (22) .
5. Numerical simulations
In the numerical computations, we discretize (5) with a finite element method and consider twogeometries with smooth interfaces. Then we approximate the eigenvalues of the matrix problem(17) with the new version of TIAR outlined in Section 4. Note that in order to preserve the accuracyfor basis functions of degree larger than one it is, for the considered geometries, mandatory to usecurvilinear elements [14, 13]. When the eigenvalues are semi-simple, we expect optimal convergencewith the h -version and with the p -version of the finite element method [3]. In h -FEM, we fix thepolynomial order p of the basis functions and decrease the maximum size of the cells h . Then, theoptimal convergence is algebraic: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ j − λ γj λ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ cN − p , c > , (35)11 lgorithm 1: Derivative computation for (34) input : Number of derivatives k max , largest index in modulus ν max , sequence of poles z , . . . , z p output : The matrix X in (34) consisting of derivatives of ˜ g i Compute r = H k max + ν max +1 ( µ ) ∈ C k max + ν max +1 with (29) for k = 1 , . . . , k max + 1 do Compute r k = aB k max + ν max r k − where B k max + ν max +1 is given by (31) end Set q = J T e where J = J µ ∈ C k max × k max given in (24) for i = 1 , . . . , p do Set q = J T q − z i q endfor i = 1 , . . . , ν max + 1 do Compute the tridiagonal Toeplitz matrix H ∈ C k max × k max with (25) where we set f ( j ) ( µ ) := ( r j ) i , j = 0 , . . . , k max Compute the tridiagonal Toeplitz matrix H (cid:48) ∈ C k max × k max with (25) where we set f ( j ) ( µ ) := ( r j +1 ) i , j = 0 , . . . , k max Set x i = H (cid:48) H − q end j m (cid:60) λ j (cid:61) λ j j m (cid:60) λ j (cid:61) λ j .
021 766 303 207 − .
273 829 280 623 4 0 19 .
243 876 046 899 − .
274 713 999 6012 8 8 .
936 779 164 355 − .
164 935 525 246 5 19 19 .
241 527 655 113 − .
104 420 737 3523 14 8 .
783 835 782 061 − .
000 247 588 219 6 25 19 .
156 200 970 821 − .
000 653 924 318
Table 1:
Selected reference eigenvalues for the problem 5.1, computed from (37) and ordered by |(cid:61) λ j | . where N denotes the number of degrees of freedom. In p -FEM, we fix the mesh and increase p ,which result in the exponential rate of convergence (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ j − λ γj λ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ αe − βN , α, β > . (36)It is well known that hp -FEM ( p -FEM in the case of smooth interfaces) is superior to h -FEM interms of accuracy vs number of degrees of freedom [44] but the use of higher order basis functionsresults in less sparse matrices. Therefore, it is for a given matrix size more time consuming tosolve a matrix eigenvalue problem generated by a high-order finite element method. The studiedNEPs are solved with the specialization of the infinite Arnoldi method outlined in Section 4. Weillustrate the convergence of the finite element method as well as the performance of the linearalgebra solver. In particular, we show that a discretization with the p -version of the finite elementmethod outlined in Section 3 together with the new version of the infinite Arnoldi method is anefficient tool for resonance calculations. 12 -12 -10 -8 -6 -4 -2 cN ! , e ! - N = -15 -12 -9 -6 -3 λ λ λ λ λ λ a) b)c) N Nν max R e l a t i v ee rr o r R e l a t i v ee rr o r Figure 2:
Convergence for problem 5.1: a) h -FEM with p = 2 and b) p -FEM with fitting curve (black): α =1 . × , β = 0 . . c) convergence with respect to ν max . Re λ j Im λ j j Re λ j Im λ j .
499 842 − .
400 318 9 15 19 .
256 306 448 9 − .
066 528 956 842 3 .
082 426 − .
179 510 2 .
070 609 128 9 − .
533 898 079 753 3 .
662 856 165 − .
980 016 551 .
696 799 709 9 − .
386 922 639 0034 3 .
035 882 038 − .
910 209 072 18 98 .
835 675 963 6 − .
069 143 373 7915 1 .
098 616 611 0 − .
005 745 095 69 19 98 .
825 451 610 5 − .
059 689 390 7906 2 .
165 520 379 3 − .
537 312 290 13 20 99 .
406 192 996 6 − .
000 006 880 0637 4 .
370 036 082 6 − .
526 561 686 26 21 99 .
406 192 346 6 − .
000 006 849 9598 3 .
958 068 685 7 − .
528 959 556 45 22 99 .
406 192 109 7 − .
000 009 251 5939 4 .
894 990 528 7 − .
402 810 837 17 .
406 194 108 9 − .
000 009 059 28510 20 .
001 865 223 0 − .
199 793 327 65 24 99 .
253 077 460 0 − .
000 000 004 940 .
159 017 340 2 − .
631 617 902 52 25 99 .
253 077 462 6 − .
000 000 004 46512 20 .
329 616 908 4 − .
530 465 014 38 26 99 .
253 077 459 3 − .
000 000 004 44113 21 .
059 617 919 8 − .
413 472 666 47 27 99 .
253 077 463 7 − .
000 000 004 263 .
176 565 085 7 − .
084 153 047 32 .
222 941 196 1 − .
000 000 000 000 1
Table 2:
Reference eigenvalues for problem 5.2, computed with with p = 30 , N =
317 281 . In this subsection, we consider the classical dielectric disk resonator [17]. Using separation ofvariables, we obtain equation (37). A complex Newton root finder is used to compute very accurateapproximations of the resonances. These approximations are used as a benchmark for studyingthe convergence of the used finite element methods together with the new version of the infiniteArnoldi method.Consider an open disk Ω R with radius R and refractive index η ( r ) = η s in Ω R and η ( r ) = 1 inΩ + R := R \ Ω R . Exact solutions to this problem are given in [17] and the resonance relationshipreads J m ( λη s R )( H (1) m ( λR )) (cid:48) − η s J (cid:48) m ( λη s R ) H (1) m ( λR ) = 0 . (37)For each m in equation (37), we search numerically n resonances λ m , . . . , λ mn . The Newton rootfinder presented in [37] is used with machine precision stopping criterion.In Table 1, we list a selection of benchmark values λ mn computed from (37), which are used toevaluate accuracy and convergence of the numerical solutions.The solutions λ mn are classified as external resonances or internal resonances [17]. The inter-nal resonances, also called Whispering Gallery Modes (WGM), have broad applications in optics,photonics, communications, and engineering [26]. These resonances feature negative imaginaryparts that are close to the real axis. Approximation of the internal modes with FEM require thatthe oscillatory behavior inside the resonator is resolved, but the modes outside the resonator arealmost constant and therefore less demanding to approximate.Exterior resonances are characterized by having large negative imaginary parts and the cor-responding modes grow very quickly outside the resonator. Hence, FEM approximation of thesemodes is demanding.
If the DtN map (3) is placed at a = R , then only one term is non-zero because the DtN mapcoincides with the compatibility condition for derivatives. However, by shifting the disk Ω R a14istance d away from the origin and taking a > R + d , the radial symmetry is lost and more termsare necessary in the DtN map (3). This more demanding problem formulation is used to test oursolution scheme in terms of accuracy as well as convergence with respect to the number of DtNterms. Let Γ a := ∂ Ω a denote a circle of radius a centred in the origin and let Ω dR ⊂ Ω a denote theshifted disk. We choose R = 1 with the shifting distance d = R/ x -axis, η s = 2, and a = 2.Once the geometry is set, we discretize (5) with the finite element method outlined in Section3.1 and apply the new specialization of TIAR in Section 4 to a FE discretization with p = 20, a = 2and N = 19 361. In the TIAR scheme we use the shifts µ = 9 − . i to approximate λ , λ , λ and µ = 19 − . i for λ , λ , λ , where λ , . . . , λ are the reference values in Table 1. The resonances λ and λ are classified as interior resonances while the remaining values are classified as exteriorresonances [17]. For the given shifts, we compute simultaneously approximations to both types ofresonances. Hence, the used finite element space must be able to capture oscillations in the interiorof the resonator and rapid growth in the exterior.Figure 2 illustrates the convergence of a sequence of eigenvalues λ γj approaching the referencevalues λ j in Table 1. Optimal convergence rates are reached: (35) in a) and (36) in b). Forthe analysis of the ν max -dependence we work with a fixed discretization and plot in Figure 2 (c)relative errors vs ν max . We observe a preasymptotic phase until a critical ν max = ˜ ν is reachedand for ν max > ˜ ν the convergence is very rapid. Figure 2 (c) also illustrates stability of TIAR inthe sense that having more ν terms than ˜ ν in the numerical experiments did not result in largernumerical errors.As expected, errors in the approximations λ γj drop faster when λ j in Table 1 correspond toan eigenvalue with small m . The empirical rule ν max > aλ was used in [22] as an estimate ofthe necessary number of truncation terms ν max in the DtN method applied for source problems( λ ∈ R ). Our computations indicate that ν max > a Re λ can be used for the current eigenvalueproblem. For example, from Figure 2 (c), it is clear that the relative errors of the numerical λ γj decrease below 10 − for ν max > a Re λ γj . Adjacent resonators are of particular interest as they exhibit oscillatory modes that cannot beexcited by single resonators and have interesting physical properties such as special mode symme-tries (fanoresonances) that may exhibit strong coupling, and higher Q -factors [24, 41, 39].The studied geometry consists of two disks of radius R = 1 / s = R . Each disk has constant refractive index η = 2, and are surrounded by vacuum η = 1.The setting is such that supp( η − ⊂ Ω a with a = 1. We study convergence by computingreference eigenvalues λ j with the new specialization of TIAR, outlined in section 4, applied to adiscretization with p = 30 and list them in Table 2.Since no exact solution is available, we study convergence with respect to three different parameters:degrees of freedom N , ν and TIAR iterations k . In Figure 4 we show convergence for the disk dimer problem. In (a) and (b) convergencewith respect to N and in c) convergence with respect to ν . The computed relative errors of theeigenvalues satisfy the estimates (35) and (36). Hence, the numerical computation indicate thatthe asymptotic convergence rates are optimal. As expected, the preasymptotic phase is longer for15 λ λ λ λ λ Figure 3:
Fields | u j ( x , x ) | corresponding to eigenvalues λ j given with bold font in Table 2. large Re λ j [43] and eigenvalues with smaller Re λ j converge faster. This can be seen by comparingthe model (36) for different λ γj . The fitted curves following λ γ and λ γ are plotted in dashed-blackand solid-black lines in Figure 4 (b). The exponential rate for λ is β = 0 . λ is β = 0 . β > β is in agreement with Re λ = 19 . < Re λ = 98 . ν max behaves similarly as described in Section5.1. In Figure 4 (c), the real part of λ is 19 .
159 and the relative error drops below 10 − for ν max >
20. Similarly for Re λ = 99 . − for all ν max > u corresponding to λ oscillates in a confined region aroundthe resonator resembling WGMs [26]. Furthermore, the figure shows that | u ( x ) | ≈ constantfor x ∈ Γ a . Hence, in the Fourier series (16) only the ν = 0 term is necessary to accuratelyapproximate u . In agreement, Figure 4 (c) shows that the relative error in λ is approximately10 − for ν max ≥ z j of G ( λ ) defined in (15) with a = 1. In the Figure 5 (b)we illustrate a situation with a = 2, where there are poles in between the shift µ and the closest λ γj . In this case a pole gets very close to λ as illustrated in the Figure 5 (b). We evaluate theeffectiveness of the pole cancellation technique by choosing µ = 1 . − . i and running TIARwith and without pole cancellation. In Figure 6 (b) we show the convergence vs. iterations ofthe experiment described above. Without pole cancellation we only get convergence for λ (reddashed line), until stagnation at around 10 − . When using pole cancellation, ˜ λ converges untilmachine precision (solid line) and ˜ λ also converges (dotted line), which suggests that the radiusof convergence is greater when pole cancellation is used.16 -10 -8 -6 -4 -2 cN ! , e ! - N = , e ! - N = -12 -9 -6 -3 λ λ λ λ λ a) b)c) N Nν max R e l a t i v ee rr o r R e l a t i v ee rr o r Figure 4:
Convergence for problem 5.2: a) h -FEM with p = 2 and b) p -FEM. The fitting curves in black use: α = 26 . , β = 0 . , α = 1 . × , β = 0 . . c) convergence with respect to ν max . λλ r z λλ γ z µ a) b) λ λ Re λ Re λ Im λ Figure 5: a) Eigenvalues λ computed with several shifts µ and a = 1 . Selected reference eigenvalues λ r , included inTable 2, are marked with red stars and poles z in blue bullets. b) Situation illustrating pole cancellation with a = 2 . -20 -15 -10 -5 -20 -15 -10 -5 a) b) ˜ λ ˜ λ λ µ = 20 − . i µ = 1 . − . i R e l a t i v e r e s i du a l n o r m iteration k iteration k Figure 6: a) TIAR eigenvalue convergence vs iterations. b) Case illustrating Pole cancellation. We show convergencewithout pole canceling in red dashed line and convergence using pole canceling in solid line.
Performance comparison:
Below we briefly discuss the performance of the proposed NEP solver.In the startup phase, we compute for given shift µ the LU factorization of T ( µ ) and refer to thetime spent as LU time. The LU factorization is only computed once and it is used in the inverseoperation (22).In Figure 7 we show performance plots for these routines and give in colors the relative errorof λ for each computation. By drawing a vertical line in Figure 7 (a), it becomes clear thatLU computation becomes more expensive for p -FEM than h -FEM. This is expected as matricesgenerated from p -FEM are denser than those from h -FEM. Moreover, from the plot we see thatthe TIAR performance is fairly balanced for both p -FEM and h -FEM depending mostly on N .Regarding accuracy, we can draw a horizontal line in any of the plots in Figure 7 and picktwo computations that lie close to this line, meaning that both took similar computational times.Then by comparing the errors given in colors, it is apparent that p -FEM reachs lower errors than18 | λ j − λ γj | / | λ j | − − − . × − N . × s ) 440 3 .
23 0 .
47 228 . s ) 2464 18 .
58 4 .
74 149 . Table 3:
Comparison for fixed errors on λ in problem 5.2 with shift µ = 20 − . i . -1 p -FEM h -FEM, p = 2 p -FEM h -FEM, p = 2 − − − − − − a) b) LU TIAR C P U t i m e ( s ) N N
Figure 7: a) LU and b) TIAR performance for problem 5.2. In colors we give relative errors for p -FEM (bullets) and h -FEM (squares) computed for λ that correspond to the bullets/squares in black. h -FEM for all computations. In Table 3 we list data for the numerical computation of λ andcompare the performance for h -FEM and p -FEM by keeping accuracy fixed. The columns with p = 1 , h -FEM (low polynomial orders and several mesh refinements) and thecolumn with p = 11 is p -FEM. In the last column ( p = 20) we list data for a highly accurate p -FEM discretization. The simulations suggest that for the given problems 5.1 and 5.2 p -FEMsurpasses h -FEM in terms of performance for the proposed NEP solver.
6. Conclusions and outlook
We have proposed a fast and efficient method for computing resonances and resonant modes ofHelmholtz problems posed as a NEP. The finite element method is used for discretizing the prob-lem, and the resulting NEP is solved with a new specialization of the infinite Arnoldi method. Inthe numerical experiments we observe that the performance of the TIAR iterations scale linearlywith the problem size and it is stable with respect to the number of terms in the Dirichlet-to-Neumann map. A pole cancellation technique was successfully applied to increase the radius ofconvergence when the shift is close to a pole. The h and the p version of the finite element methodwere used to discretize the Fredholm operator function and we showed that the nonlinear matrixeigenvalue solver performs well in all cases. The exponential convergence of the p -version of FE,for problems with smooth interfaces, together with the new version of the infinite Arnoldi methodis therefore an efficient tool for resonance calculations. Moreover, we expect the same performancefor hp -FEM, since the sparsity of the matrices do not critically affect the performance of the TIAR19ethod. Acknowledgments.
We gratefully acknowledge the support of the Swedish Research Councilunder Grant No. 621-2012-3863 and 621-2013-4640. J. Ara´ujo also thanks the department of Math-ematics at KTH Royal Institute of Technology very much for the kind hospitality and GiampaoloMele for interesting discussions held during the visit.
References [1] M. Abramowitz and I. A. Stegun, editors.
Handbook of Mathematical Functions . Applied Mathematics SeriesNo. 55. National Bureau of Standards, Washington D.C., 1970.[2] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami, and K. Kimura. A numerical method for nonlinear eigenvalueproblems using contour integrals.
JSIAM Letters , 1:52–55, 2009.[3] I. Babuˇska and J. Osborn. Eigenvalue problems. In
Handbook of numerical analysis, Vol. II , Handb. Numer.Anal., II, pages 641–787. North-Holland, Amsterdam, 1991.[4] I. Babuˇska, Q. Buo, and J. E. Osborn. Regularity and numerical solution of eigenvalue problems with piecewiseanalytic data.
SIAM J. Numer. Anal. , 26(6):1534–1560, 1989.[5] I. Babuˇska and B. Q. Guo. The h, p and h-p version of the finite element method: Basis theory and applications.
Adv. Eng. Softw. , 15(3-4):159–174, November 1992.[6] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith. Efficient management of paral-lelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors,
Modern Software Tools in Scientific Computing , pages 163–202. Birkh¨auser Press, 1997.[7] W. Bangerth, T. Heister, L. Heltai, G. Kanschat, M. Kronbichler, M. Maier, B. Turcksin, and T. D. Young.The deal.II library, version 8.2.
Archive of Numerical Software , 3, 2015.[8] M. Van Barel and P. Kravanja. Nonlinear eigenvalue problems and contour integrals.
J. Comput. Appl. Math. ,292:526–540, 2016.[9] R. Van Beeumen, K. Meerbergen, and W. Michiels. A rational Krylov method based on Hermite interpolationfor nonlinear eigenvalue problems.
SIAM J. Sci. Comput. , 35(1):A327–A350, 2013.[10] T. Betcke, N. J. Higham, V. Mehrmann, C. Schr¨oder, and F. Tisseur. NLEVP: A collection of nonlineareigenvalue problems.
ACM Trans. Math. Softw. , 39(2):1–28, 2013.[11] T. Betcke and H. Voss. A Jacobi-Davidson type projection method for nonlinear eigenvalue problems.
FutureGeneration Computer Systems , 20(3):363–372, 2004.[12] W.-J. Beyn. An integral method for solving nonlinear eigenvalue problems.
Linear Algebra Appl. , 436(10):3839–3863, 2012.[13] Susanne C. Brenner and L. Ridgway Scott.
The mathematical theory of finite element methods . Texts in appliedmathematics. Springer, New York, Berlin, Paris, 2002.[14] P.G. Ciarlet and P.-A. Raviart. Interpolation theory over curved elements, with applications to finite elementmethods.
Computer Methods in Applied Mechanics and Engineering , 1(2):217 – 249, 1972.[15] D. Colton and R. Kress.
Inverse Acoustic and Electromagnetic Scattering Theory . Springer-Verlag, Berlin, 1992.[16] C. Dettmann, G. V. Morozov, M. Sieber, and H. Waalkens. Internal and external resonances of dielectric disks.
Europhysics letters , 87(3), 2009.[17] C. P. Dettmann, G. V. Morozov, M. Sieber, and H. Waalkens. Internal and external resonances of dielectricdisks.
EPL (Europhysics Letters) , 87(3):34003, 2009.[18] Christian Engstr¨om and Luka Grubiˇsi´c. A subspace iteration algorithm for Fredholm valued functions.
Math.Probl. Eng. , pages Art. ID 459895, 14, 2015.[19] S. Fliss. A Dirichlet-to-Neumann approach for the exact computation of guided modes in photonic crystalwaveguides.
SIAM J. Sci. Comput. , 35(2):B438–B461, 2013.[20] I. C. Gohberg and M. G. Kre˘ın.
Introduction to the theory of linear nonselfadjoint operators , volume 18 of
Translations of mathematical monographs . American Mathematical Society, Providence, Rhode Island, 1969.[21] S. G¨uttel, R. Van Beeumen, K. Meerbergen, and W. Michiels. NLEIGS: a class of fully rational Krylov methodsfor nonlinear eigenvalue problems.
SIAM J. Sci. Comput. , 36(6):A2842–A2864, 2014.[22] Isaac Harari and Thomas J. R. Hughes. Analysis of continuous formulations underlying the computation oftime-harmonic acoustics in exterior domains.
Comput. Methods Appl. Mech. Eng. , 97(1):103–124, May 1992.[23] N. J. Higham.
Functions of Matrices. Theory and Computation . SIAM, 2008.
24] A D Humphrey and W L Barnes. Plasmonic surface lattice resonances in arrays of metallic nanoparticle dimers.
Journal of Optics , 18(3):035005, 2016.[25] T. Ikegami, T. Sakurai, and U. Nagashima. A filter diagonalization for generalized eigenvalue problems basedon the Sakurai-Sugiura projection method.
J. Comput. Appl. Math. , 233(8):1927–1936, 2010.[26] V. S. Ilchenko and A. B. Matsko. Optical resonators with whispering-gallery modes - Part II: Applications.
IEEE journal of selected topics in quantum electronics , 12(1), January 2006.[27] E. Jarlebring, G. Mele, and O. Runborg. The waveguide eigenvalue problem and the tensor infinite Arnoldimethod. Technical report, KTH Royal Institute of Technology, 2015. arxiv preprint.[28] E. Jarlebring, W. Michiels, and K. Meerbergen. A linear eigenvalue algorithm for the nonlinear eigenvalueproblem.
Numer. Math. , 122(1):169–195, 2012.[29] O. Karma. Approximation in eigenvalue problems for holomorphic Fredholm operator functions. I.
Numer.Funct. Anal. Optim. , 17(3–4):365–387, 1996.[30] O. Karma. Approximation in eigenvalue problems for holomorphic Fredholm operator functions. II.
Numer.Funct. Anal. Optim. , 17(3–4):389–408, 1996.[31] L. Kaufman. Eigenvalue problems in fiber optic design.
SIAM J. Matrix Anal. Appl. , 28(1):105–117, 2006.[32] D. Klindworth and K. Schmidt. Dirichlet-to-Neumann transparent boundary conditions for photonic crystalwaveguides.
IEEE Transactions on Magnetics , 50(2):7005204, 2014.[33] D. Kressner. A block Newton method for nonlinear eigenvalue problems.
Numer. Math. , 114(2):355–372, 2009.[34] M. Lenoir, M. Vullierme-Ledard, and C. Hazard. Variational formulations for the determination of resonantstates in scattering problems.
SIAM J. Math. Anal. , 23(3):579–608, 1992.[35] B.-S. Liao, Z. Bai, L.-Q. Lee, and K. Ko. Solving large scale nonlinear eigenvalue problems in next-generationaccelerator design. Technical Report SLAC-PUB-12137, Stanford University, 2006.[36] B.-S. Liao, Z. Bai, L.-Q. Lee, and K. Ko. Nonlinear Rayleigh-Ritz iterative method for solving large scalenonlinear eigenvalue problems.
Taiwanese Journal of Mathematics , 14(3):869–883, 2010.[37] Adi Ben-Israel Lily Yau. The newton and halley methods for complex roots.
The American MathematicalMonthly , 105(9):806–818, 1998.[38] Richard B. Melrose.
Geometric scattering theory . Stanford Lectures. Cambridge University Press, Cambridge,1995.[39] N. Asger Mortensen, Søren Raza, Martijn Wubs, Thomas Søndergaard, and Sergey I. Bozhevolnyi. A generalizednon-local optical response theory for plasmonic nanostructures.
Nature Communications , 5, 2014.[40] Hae Soo Oh and Ivo Babuˇska. The p -version of the finite element method for the elliptic boundary valueproblems with interfaces. Comput. Methods Appl. Mech. Engrg. , 97(2):211–231, 1992.[41] G. Rosolen and B. Maes. Asymmetric and connected graphene dimers for a tunable plasmonic response.
Phys.Rev. B , 92:205405, Nov 2015.[42] J. Sanchez Hubert and E. S´anchez-Palencia.
Vibration and coupling of continuous systems . Springer-Verlag,Berlin, 1989. Asymptotic methods.[43] S. Sauter. hp -finite elements for elliptic eigenvalue problems: Error estimates which are explicit with respect to λ , h , and p . SIAM J. Numer. Anal. , 48(1):95–108, 2010.[44] C. Schwab. p- and hp- Finite Element Methods: Theory and Applications in Solid and Fluid Mechanics . OxfordUniversity Press, 1998.[45] Pavel Solin, Karel Segeth, and Ivo Dolezel.
Higher-order finite element methods . Studies in advanced mathe-matics. Chapman &Hall/CRC, Boca Raton, London, 2004.[46] A. Spence and C. Poulton. Photonic band structure calculations using nonlinear eigenvalue techniques.
J.Comput. Phys. , 204(1):65–81, 2005.[47] D. B. Szyld and F. Xue. Preconditioned eigensolvers for large-scale nonlinear hermitian eigenproblems withvariational characterization. I. conjugate gradient methods. Technical report, Temple University, 2014.[48] J. Tausch. Computing Floquet-Bloch modes in biperiodic slabs with boundary elements.
J. Comput. Appl.Math. , 254:192–203, 2013.[49] R. Van Beeumen.
Rational Krylov Methods for Nonlinear Eigenvalue Problems . PhD thesis, Department ofComputer Science, KU Leuven, Leuven, Belgium, 2015.[50] H. Voss. An Arnoldi method for nonlinear eigenvalue problems.
BIT , 44:387 – 401, 2004.[51] H. Voss. Nonlinear eigenvalue problems. In L. Hogben, editor,
Handbook of Linear Algebra, Second Edition ,number 164 in Discrete Mathematics and Its Applications. Chapman and Hall/CRC, 2013.[52] G. N. Watson.
A treatise on the theory of Bessel functions . Cambridge Mathematical Library. CambridgeUniversity Press, Cambridge, 1995. Reprint of the second (1944) edition.[53] Maciej Zworski. Lectures on scattering resonances version 0 .
03. University Lecture, 2015.03. University Lecture, 2015.