An adaptive fast solver for the modified Helmholtz equation in two dimensions

In this paper, we present a fast multipole-accelerated integral equation method for solving the modified Helmholtz equation Δu(x→)-β2u(x→)=f(x→) in two dimensions. The method is direct, and unlike classical FFT based fast solvers, it allows for adaptive mesh refinement but with comparable amount of work per grid point. When the computational domain is rectangular, Dirichlet, Neumann, periodic, and free-space boundary conditions can be imposed analytically without the need to solve a system of linear equations. Several important features of the algorithm are discussed, including the use of precomputed tables, diagonal translation operators, and lattice sums to impose periodic boundary conditions. Numerical experiments show that, for a wide range of the parameter β, the algorithm is stable and high-order accurate.