Sheehan Olver

Tao Chen, Sheehan Olver

We show that the $QR$ factorization of a banded-plus-semiseparable (BPS) matrix is computable in optimal linear complexity with respect to the discretization size by showing that the intermediate stages of a $QR$ factorization as computed using Householder reflection maintain a specific structure which has optimal storage. This theoretical result enables the design of stable, linear-complexity algorithms for solving the associated linear systems. For symmetric BPS matrices, we further show that the $RQ$ product -- central to eigenvalue computations via the $QR$ algorithm -- also preserves the BPS structure, leading to a linear-complexity algorithm for each iteration. Numerical experiments validate the optimal linear complexity, confirm high numerical accuracy, and demonstrate substantial speedups compared with existing hierarchical approaches. The algorithms have been implemented in an open-source Julia package, providing an efficient and accessible platform for practical use.

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.

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, 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.

Jiajie Yao, Marco Fasondini, Sheehan Olver

We develop a sparse spectral method for solving partial differential equations on a class of two-dimensional geometries bounded by algebraic curves. The numerical method uses generalised bivariate Koornwinder polynomials which form a complete orthogonal basis, but one which is not graded in terms of polynomial degree. The polynomials are built from new families of univariate semiclassical orthogonal polynomials whose associated operator matrices (Jacobi matrices, raising matrices and differentiation matrices) are computed with optimal linear complexity in the number of basis functions. When used to discretise partial differential equations the resulting matrices are sparse enabling efficient numerical solution. Moreover, we develop fast transforms from values on a grid to expansion coefficients. The efficiency and accuracy of the resulting spectral method are illustrated through a series of numerical experiments on geometries whose boundaries are smooth and piecewise smooth including non-convex geometries. We observe algebraic convergence for geometries with corners, which accelerates to exponentially fast (spectral) convergence when the boundary is smooth.

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, 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})$.