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

TL;DR

This paper introduces a fast, spectral method-based solver for linear PDEs on rectangular domains, leveraging separable representations and matrix equations, capable of handling large-scale problems efficiently in MATLAB and Julia.

Contribution

The paper presents a novel spectral method for solving PDEs with variable coefficients, utilizing separable operators and Sylvester equations for efficient computation, and provides a publicly available MATLAB implementation.

Findings

Solves PDEs of splitting rank 2 using Sylvester matrix equations.

Achieves solutions with over a million degrees of freedom in under 60 seconds in MATLAB.

Julia implementation can perform the same solve in about 10 seconds.

Abstract

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 $(n_x,n_y)$ is computed in $\mathcal{O}((n_x n_y)^{3/2})$ operations. Partial differential operators of splitting rank $\geq 3$ are solved via a linear system involving a block-banded matrix in $\mathcal{O}(\min(n_x^{3} n_y,n_x n_y^{3}))$ 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.