  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.
