Estimates on Neumann eigenfunctions at the boundary, and the "Method of Particular Solutions" for computing them

We consider the "Method of particular solutions" for numerically computing eigenvalues and eigenfunctions of the Laplacian $\Delta$ on a smooth, bounded domain Omega in RR^n with either Dirichlet or Neumann boundary conditions. This method constructs approximate eigenvalues E, and approximate eigenfunctions u that satisfy $\Delta u=Eu$ in Omega, but not the exact boundary condition. An inclusion bound is then an estimate on the distance of E from the actual spectrum of the Laplacian, in terms of (boundary data of) u. We prove operator norm estimates on certain operators on $L^2(\partial \Omega)$ constructed from the boundary values of the true eigenfunctions, and show that these estimates lead to sharp inclusion bounds in the sense that their scaling with $E$ is optimal. This is advantageous for the accurate computation of large eigenvalues. The Dirichlet case can be treated using elementary arguments and has appeared in SIAM J. Num. Anal. 49 (2011), 1046-1063, while the Neumann case seems to require much more sophisticated technology. We include preliminary numerical examples for the Neumann case.