A direct solver for elliptic PDEs in three dimensions based on hierarchical merging of Poincaré–Steklov operators

A numerical method for variable coefficient elliptic PDEs on three dimensional domains is described. The method is designed for problems with smooth solutions, and is based on a multidomain spectral collocation discretization scheme. The resulting system of linear equations can very efficiently be solved using a nested dissection style direct (as opposed to iterative) solver. This makes the scheme particularly well suited to solving problems for which iterative solvers struggle; in particular for problems with oscillatory solutions. A principal feature of the scheme is that once the solution operator has been constructed, the actual solve is extremely fast. An upper bound on the asymptotic cost of the build stage of O(N4/3)O(N4/3) is proved (for the case where the PDE is held fixed as NN increases). The solve stage has close to linear complexity. The scheme requires a relatively large amount of storage per degree of freedom, but since it is a high order scheme, a small number of degrees of freedom is sufficient to achieve high accuracy. The method is presented for the case where there is no body load present, but it can with little difficulty be generalized to the non-homogeneous case. Numerical experiments demonstrate that the scheme is capable of solving Helmholtz-type equations on a domain of size 20×20×2020×20×20 wavelengths to three correct digits on a modest personal workstation, with N≈2⋅106N≈2⋅106.