The automatic solution of partial differential equations using a global spectral method

A spectral method for solving linear partial differential equations (PDEs) with variable coefficients and general boundary conditions defined on rectangular domains is described, based on separable representations of partial differential operators and the one-dimensional ultraspherical spectral method. If a partial differential operator is of splitting rank 2, such as the operator associated with Poisson or Helmholtz, the corresponding PDE is solved via a generalized Sylvester matrix equation, and a bivariate polynomial approximation of the solution of degree (nx,ny) is computed in O((nxny)3/2) operations. Partial differential operators of splitting rank ≥3 are solved via a linear system involving a block-banded matrix in O(min⁡(nx3ny,nxny3)) operations. Numerical examples demonstrate the applicability of our 2D spectral method to a broad class of PDEs, which includes elliptic and dispersive time-evolution equations. The resulting PDE solver is written in Matlab and is publicly available as part of Chebfun. It can resolve solutions requiring over a million degrees of freedom in under 60 seconds. An experimental implementation in the Julia language can currently perform the same solve in 10 seconds.