Fast Poisson solvers for spectral methods

Abstract

Poisson's equation is the canonical elliptic partial differential equation. While there exist fast Poisson solvers for finite difference and finite element methods, fast Poisson solvers for spectral methods have remained elusive. Here, we derive spectral methods for solving Poisson's equation on a square, cylinder, solid sphere, and cube that have an optimal complexity (up to polylogarithmic terms) in terms of the degrees of freedom required to represent the solution. Whereas FFT-based fast Poisson solvers exploit structured eigenvectors of finite difference matrices, our solver exploits a separated spectra property that holds for our spectral discretizations. Without parallelization, we can solve Poisson's equation on a square with 100 million degrees of freedom in under two minutes on a standard laptop. Throughout this paper, "optimal complexity" means a computational complexity that is optimal up to polylogarithmic factors. 1

On the singular values of matrices with high displacement rank

“…If E and G are disks in the complex plane, the required single shift parameter (α, β) is given by Theorem 2, a rotation mapping, and the formula in [21]. When E and G are closed real intervals, we refer the reader to the formulas in [17], as well as the MATLAB code in [12,Appendix A]. For most other choices of E and G, heuristic shift selection strategies must be employed (see Section 5.1).…”

“…In [12], spectral discretizations are developed so that the ADI method can be used to solve Poisson's equation on a variety of domains in optimal computational complexity (up to polylogarithmic factors). Combining these ideas with FI-ADI leads to highly efficient Poisson solvers that construct low rank approximations to solutions.…”

“…A standard numerical approach for finding u involves discretizing (38) using second-order finite differences. 11 To achieve spectral accuracy, we instead apply the method in [12], where (38) is discretized in a way that leads to the matrix equation AX −XB =DFD T .…”

“…Here,X andF contain scaled expansion coefficients for expressing u and f , respectively, in a particular ultraspherical polynomial basis, andD is diagonal. We refer the reader to [12,Sec. 3] for further details.…”

“…In this paper, we relax the assumption that rank(F ) is small and assume instead that the singular values of F decay rapidly, so that rank (F ) is small. Such scenarios occur in numerical computing [12,13,14], where (1) arises from the discretization of certain PDEs and F is associated with a smooth 2D function (see Section 6). Current methods only bound singular values of X with indices that are multiples of ρ (see Section 2), and therefore provide little to no information in this setting.…”