Sheehan Olver

Sheehan Olver, Alex Townsend

A spectral method is developed for the direct solution of linear ordinary differential equations with variable coefficients. The method leads to matrices which are almost banded, and a numerical solver is presented that takes O(m^2n) operations, where m is the number of Chebyshev points needed to resolve the coefficients of the differential operator and n is the number of Chebyshev coefficients needed to resolve the solution to the differential equation. We prove stability of the method by relating it to a diagonally preconditioned system which has a bounded condition number, in a suitable norm. For Dirichlet boundary conditions, this implies stability in the standard 2-norm. An adaptive QR factorization is developed to efficiently solve the resulting linear system and automatically choose the optimal number of Chebyshev coefficients needed to represent the solution. The resulting algorithm can efficiently and reliably solve for solutions that require as many as a million unknowns.

Percy Deift, Govind Menon, Sheehan Olver et al.

The authors present evidence for universality in numerical computations with random data. Given a (possibly stochastic) numerical algorithm with random input data, the time (or number of iterations) to convergence (within a given tolerance) is a random variable, called the halting time. Two-component universality is observed for the fluctuations of the halting time, i.e., the histogram for the halting times, centered by the sample average and scaled by the sample variance, collapses to a universal curve, independent of the input data distribution, as the dimension increases. Thus, up to two components, the sample average and the sample variance, the statistics for the halting time are universally prescribed. The case studies include six standard numerical algorithms, as well as a model of neural computation and decision making. A link to relevant software is provided in for the reader who would like to do computations of his'r own.

Sheehan Olver, Alex Townsend

We describe a framework for solving a broad class of infinite-dimensional linear equations, consisting of almost banded operators, which can be used to resepresent linear ordinary differential equations with general boundary conditions. The framework contains a data structure on which row operations can be performed, allowing for the solution of linear equations by the adaptive QR approach. The algorithm achieves $O(n^{\rm opt})$ complexity, where $n^{\rm opt}$ is the number of degrees of freedom required to achieve a desired accuracy, which is determined adaptively. In addition, special tensor product equations, such as partial differential equations on rectangles, can be solved by truncating the operator in the $y$-direction with $n_y$ degrees of freedom and using a generalized Schur decomposition to upper triangularize, before applying the adaptive QR approach to the $x$-direction, requiring $O(n_y^2 n_x^{\rm opt})$ operations. The framework is implemented in the ApproxFun package written in the Julia programming language, which achieves highly competitive computational costs by exploiting unique features of Julia.

Alex Townsend, Marcus Webb, Sheehan Olver

Many standard conversion matrices between coefficients in classical orthogonal polynomial expansions can be decomposed using diagonally-scaled Hadamard products involving Toeplitz and Hankel matrices. This allows us to derive $\smash{\mathcal{O}(N(\log N)^2)}$ algorithms, based on the fast Fourier transform, for converting coefficients of a degree $N$ polynomial in one polynomial basis to coefficients in another. Numerical results show that this approach is competitive with state-of-the-art techniques, requires no precomputational cost, can be implemented in a handful of lines of code, and is easily adapted to extended precision arithmetic.

Richard Mikael Slevinsky, Sheehan Olver

We develop a spectral method for solving univariate singular integral equations over unions of intervals by utilizing Chebyshev and ultraspherical polynomials to reformulate the equations as almost-banded infinite-dimensional systems. This is accomplished by utilizing low rank approximations for sparse representations of the bivariate kernels. The resulting system can be solved in ${\cal O}(m^2n)$ operations using an adaptive QR factorization, where $m$ is the bandwidth and $n$ is the optimal number of unknowns needed to resolve the true solution. The complexity is reduced to ${\cal O}(m n)$ operations by pre-caching the QR factorization when the same operator is used for multiple right-hand sides. Stability is proved by showing that the resulting linear operator can be diagonally preconditioned to be a compact perturbation of the identity. Applications considered include the Faraday cage, and acoustic scattering for the Helmholtz and gravity Helmholtz equations, including spectrally accurate numerical evaluation of the far- and near-field solution. The Julia software package SingularIntegralEquations.jl implements our method with a convenient, user-friendly interface.

Sheehan Olver, Thomas Trogdon

The effective and efficient numerical solution of Riemann-Hilbert problems has been demonstrated in recent work. With the aid of ideas from the method of nonlinear steepest descent for Riemann-Hilbert problems, the resulting numerical methods have been shown numerically to retain accuracy as values of certain parameters become arbitrarily large. Remarkably, this numerical approach does not require knowledge of local parametrices; rather, the deformed contour is scaled near stationary points at a specific rate. The primary aim of this paper is to prove that this observed asymptotic accuracy is indeed achieved. To do so, we first construct a general theoretical framework for the numerical solution of Riemann-Hilbert problems. Second, we demonstrate the precise link between nonlinear steepest descent and the success of numerics in asymptotic regimes. In particular, we prove sufficient conditions for numerical methods to retain accuracy. Finally, we compute solutions to the homogeneous Painlevé II equation and the modified Korteweg-de Vries equations to explicitly demonstrate the practical validity of the theory.

Sheehan Olver, Yuan Xu

Orthogonal polynomials with respect to a weight function defined on a wedge in the plane are studied. A basis of orthogonal polynomials is explicitly constructed for two large class of weight functions and the convergence of Fourier orthogonal expansions is studied. These are used to establish analogous results for orthogonal polynomials on the boundary of the square. As an application, we study the statistics of the associated determinantal point process and use the basis to calculate Stieltjes transforms.

Sheehan Olver, Alex Townsend, Geoff Vasil

In this paper, we demonstrate that many of the computational tools for univariate orthogonal polynomials have analogues for a family of bivariate orthogonal polynomials on the triangle, including Clenshaw's algorithm and sparse differentiation operators. This allows us to derive a practical spectral method for solving linear partial differential equations on triangles with sparse discretizations. We can thereby rapidly solve partial differential equations using polynomials with degrees in the thousands, resulting in sparse discretizations with as many as several million degrees of freedom.

Sheehan Olver, Raj Rao Nadakuditi

We develop a numerical approach for computing the additive, multiplicative and compressive convolution operations from free probability theory. We utilize the regularity properties of free convolution to identify (pairs of) `admissible' measures whose convolution results in a so-called `invertible measure' which is either a smoothly-decaying measure supported on the entire real line (such as the Gaussian) or square-root decaying measure supported on a compact interval (such as the semi-circle). This class of measures is important because these measures along with their Cauchy transforms can be accurately represented via a Fourier or Chebyshev series expansion, respectively. Thus, knowledge of the functional inverse of their Cauchy transform suffices for numerically recovering the invertible measure via a non-standard yet well-behaved Vandermonde system of equations. We describe explicit algorithms for computing the inverse Cauchy transform alluded to and recovering the associated measure with spectral accuracy. Convergence is guaranteed under broad assumptions on the input measures.

Sheehan Olver, Thomas Trogdon

In recent developments, a general approach for solving Riemann--Hilbert problems numerically has been developed. We review this numerical framework, and apply it to the calculation of orthogonal polynomials on the real line. Combining this numerical algorithm with an approach to compute Fredholm determinants, we are able to calculate level densities and gap statistics for general finite-dimensional unitary ensembles. We also include a description of how to compute the Hastings--McLeod solution of the homogeneous Painlevé II equation.

Nicholas Hale, Sheehan Olver

A fast algorithm (linear in the degrees of freedom) for the solution of linear variable-coefficient rational-order fractional integral and differential equations is described. The approach is related to the ultraspherical method for ordinary differential equations [Olver & Townsend 2013], and involves constructing two different bases, one for the domain of the operator and one for the range of the operator. The bases are constructed from direct sums of suitably weighted ultraspherical or Jacobi polynomial expansions, for which explicit representations of fractional integrals and derivatives are known, and are carefully chosen so that the resulting operators are banded or almost-banded. Geometric convergence is demonstrated for numerous model problems when the variable coefficients and right-hand side are sufficiently smooth.

Tom Claeys, Sheehan Olver

We study a family of distributions that arise in critical unitary random matrix ensembles. They are expressed as Fredholm determinants and describe the limiting distribution of the largest eigenvalue when the dimension of the random matrices tends to infinity. The family contains the Tracy-Widom distribution and higher order analogues of it. We compute the distributions numerically by solving a Riemann-Hilbert problem numerically, plot the distributions, and discuss several properties that they appear to exhibit.

Sheehan Olver

This paper constructs polynomial bases that capture the structure of the de Rham complex with boundary conditions in disks and cylinders (both periodic and finite) in a way that respects rotational symmetry. The starting point is explicit constructions of vector and matrix orthogonal polynomials on the unit disk that are analogous to the (scalar) generalised Zernike polynomials. We use these to build new orthogonal polynomials with respect to a matrix weight that forces vector polynomials to be normal on the boundary of the disk. The resulting weighted vector orthogonal polynomials have a simple connection to the gradient of weighted generalised Zernike polynomials, and their curl (i.e. vorticity or rot) is a constant multiple of the standard Zernike polynomials which are orthogonal with respect to $L^2$ on the disk. This construction naturally leads to bases in cylinders with simple recurrences relating their gradient, curl and divergence. These bases decouple the de Rham complex into small exact sub-complexes.

Geoffrey M. Vasil, Keaton J. Burns, Daniel Lecoanet et al.

Spectral methods are an efficient way to solve partial differential equations on domains possessing certain symmetries. The utility of a method depends strongly on the choice of spectral basis. In this paper we describe a set of bases built out of Jacobi polynomials, and associated operators for solving scalar, vector, and tensor partial differential equations in polar coordinates on a unit disk. By construction, the bases satisfy regularity conditions at r=0 for any tensorial field. The coordinate singularity in a disk is a prototypical case for many coordinate singularities. The work presented here extends to other geometries. The operators represent covariant derivatives, multiplication by azimuthally symmetric functions, and the tensorial relationship between fields. These arise naturally from relations between classical orthogonal polynomials, and form a Heisenberg algebra. Other past work uses more specific polynomial bases for solving equations in polar coordinates. The main innovation in this paper is to use a larger set of possible bases to achieve maximum bandedness of linear operations. We provide a series of applications of the methods, illustrating their ease-of-use and accuracy.

John W. Pearson, Sheehan Olver, Mason A. Porter

The two most commonly used hypergeometric functions are the confluent hypergeometric function and the Gauss hypergeometric function. We review the available techniques for accurate, fast, and reliable computation of these two hypergeometric functions in different parameter and variable regimes. The methods that we investigate include Taylor and asymptotic series computations, Gauss-Jacobi quadrature, numerical solution of differential equations, recurrence relations, and others. We discuss the results of numerical experiments used to determine the best methods, in practice, for each parameter and variable regime considered. We provide 'roadmaps' with our recommendation for which methods should be used in each situation.

Folkmar Bornemann, Alexander Its, Sheehan Olver et al.

As a basic example in nonlinear theories of discrete complex analysis, we explore various numerical methods for the accurate evaluation of the discrete map $Z^a$ introduced by Agafonov and Bobenko. The methods are based either on a discrete Painlevé equation or on the Riemann-Hilbert method. In the latter case, the underlying structure of a triangular Riemann-Hilbert problem with a non-triangular solution requires special care in the numerical approach. Complexity and numerical stability are discussed, the results are illustrated by numerical examples

Alex Townsend, Sheehan Olver

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.

Alex Townsend, Thomas Trogdon, Sheehan Olver

A fast and accurate algorithm for the computation of Gauss-Hermite and generalized Gauss-Hermite quadrature nodes and weights is presented. The algorithm is based on Newton's method with carefully selected initial guesses for the nodes and a fast evaluation scheme for the associated orthogonal polynomial. In the Gauss-Hermite case the initial guesses and evaluation scheme rely on explicit asymptotic formulas. For generalized Gauss-Hermite, the initial guesses are furnished by sampling a certain equilibrium measure and the associated polynomial evaluated via a Riemann-Hilbert reformulation. In both cases the $n$-point quadrature rule is computed in $\mathcal{O}(n)$ operations to an accuracy that is close to machine precision. For sufficiently large $n$, some of the quadrature weights have a value less than the smallest positive normalized floating-point number in double precision and we exploit this fact to achieve a complexity as low as $\mathcal{O}(\sqrt{n})$.