Fast computation of high frequency Dirichlet eigenmodes via the spectral flow of the interior Neumann-to-Dirichlet map

TL;DR

This paper introduces a fast, linearized algorithm for computing high-frequency Dirichlet eigenmodes in smooth, star-shaped domains, significantly reducing computational effort compared to traditional boundary-based methods.

Contribution

The authors develop a novel spectral flow method that linearizes the interior Neumann-to-Dirichlet operator spectrum, enabling O(N) faster eigenvalue computation for the Dirichlet Laplacian.

Findings

Achieves O(N) speedup over traditional methods.

Provides rigorous error bounds for eigenvalues and eigenfunctions.

Demonstrates high-frequency eigenmode computation with 10^-10 accuracy, 1000 times faster.

Abstract

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 $\mathbb{R}^d$, $d\ge 2$. Conventional boundary-based methods require a root-search in eigenfrequency $k$, hence take $O(N^3)$ effort per eigenpair found, using dense linear algebra, 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 $\hat{k}_j$ to the square-roots $k_j$ of all O(N) eigenvalues lying in $[k - \epsilon, k]$, where $\epsilon=O(1)$, are found with $O(N^3)$ effort. We prove an error estimate $$ |\hat k_j - k_j| \leq C \Big(\frac{\epsilon^2}{k} + \epsilon^3 \Big), $$ with $C$ independent of $k$. We present a higher-order variant with eigenvalue error scaling empirically as $O(\epsilon^5)$ and eigenfunction error as $O(\epsilon^3)$, the former improving upon the 'scaling method' of Vergini--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\"{o}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^3$ times faster than standard ones based upon a root-search.