Gaussian spectral rules for second order finite-difference schemes

Earlier the authors suggested an algorithm of grid optimization for a second order finite-difference approximation of a two-point problem. The algorithm yields exponential superconvergence of the Neumann-to-Dirichlet map (or the boundary impedance). Here we extend that approach to PDEs with piecewise-constant coefficients and rectangular homogeneous subdomains. Examples of the numerical solution of the 2-dimensional oscillatory Helmholtz equation exhibit exponential convergence at prescribed points, where the cost per grid node is close to that of the standard five-point finite-difference scheme. Our scheme demonstrates high accuracy with slightly more than two grid points per wavelength and reduces the grid size by more than three orders as compared with the standard scheme.