Fast Computation of High‐Frequency Dirichlet Eigenmodes via Spectral Flow of the Interior Neumann‐to‐Dirichlet Map

We present a new algorithm for numerical computation of large eigenvalues and associated eigenfunctions of the Dirichlet Laplacian in a smooth, star‐shaped domain in ?^d, d ≥ 2. Conventional boundary‐based methods require a root search in eigenfrequency k, hence take O(N³) effort per eigenpair found, where N = O(k^d?1) is the number of unknowns required to discretize the boundary. Our method is O(N) faster, achieved by linearizing with respect to k the spectrum of a weighted interior Neumann‐to‐Dirichlet (NtD) operator for the Helmholtz equation. Approximations to the square roots k_j of all O(N) eigenvalues lying in [k ? ?, k], where ? = O(1), are found with O(N³) effort. We prove an error estimate with C independent of k. We present a higher‐order variant with eigenvalue error scaling empirically as O(?⁵) and eigenfunction error as O(?³), the former improving upon the “scaling method” of Vergini and Saraceno. For planar domains (d = 2), with an assumption of absence of spectral concentration, we also prove rigorous error bounds that are close to those numerically observed. For d = 2 we compute robustly the spectrum of the NtD operator via potential theory, Nyström discretization, and the Cayley transform. At high frequencies (400 wavelengths across), with eigenfrequency relative error 10^?10, we show that the method is 10³ times faster than standard ones based upon a root search. © 2014 Wiley Periodicals, Inc.