Source author record

Sheehan Olver

Sheehan Olver appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

21works
14topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

21 published item(s)

preprint2020arXiv

Non-homogeneous wave equation on a cone

The wave equation $\left(\partial_{tt} - c^2 Δ_x\right) u(x,t) = e^{-t} f(x,t)$ is shown to have a unique solution if $u$ and its partial derivatives in $x$ are in $L^2(e^{-t})$ on the cone, and the solution can be explicit given in the Fourier series of orthogonal polynomials on the cone. This provides a particular solution for the boundary value problems of the non-homogeneous wave equation on the cone, which can be combined with a solution to the homogeneous wave equation in the cone to obtain the full solution.

preprint2020arXiv

Orthogonal structure on a quadratic curve

Orthogonal polynomials on quadratic curves in the plane are studied. These include orthogonal polynomials on ellipses, parabolas, hyperbolas, and two lines. For an integral with respect to an appropriate weight function defined on any quadratic curve, an explicit basis of orthogonal polynomials is constructed in terms of two families of orthogonal polynomials in one variable. Convergence of the Fourier orthogonal expansions is also studied in each case. As an application, we see that the resulting bases can be used to interpolate functions on the real line with singularities of the form $|x|$, $\sqrt{x^2+ ε^2}$, or $1/x$, with exponential convergence.

preprint2020arXiv

Sparse spectral and p-finite element methods for partial differential equations on disk slices and trapeziums

Sparse spectral methods for solving partial differential equations have been derived in recent years using hierarchies of classical orthogonal polynomials on intervals, disks, and triangles. In this work we extend this methodology to a hierarchy of non-classical orthogonal polynomials on disk slices (e.g. a half-disk) and trapeziums. This builds on the observation that sparsity is guaranteed due to the boundary being defined by an algebraic curve, and that the entries of partial differential operators can be determined using formulae in terms of (non-classical) univariate orthogonal polynomials. We apply the framework to solving the Poisson, variable coefficient Helmholtz, and Biharmonic equations.

preprint2016arXiv

Fast polynomial transforms based on Toeplitz and Hankel matrices

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.

preprint2016arXiv

Tensor calculus in polar coordinates using Jacobi polynomials

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.

preprint2015arXiv

Numerical Methods for the Computation of the Confluent and Gauss Hypergeometric Functions

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.

preprint2015arXiv

Numerical Methods for the Discrete Map $Z^a$

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

preprint2015arXiv

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 $(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.

preprint2014arXiv

A practical framework for infinite-dimensional linear algebra

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.

preprint2014arXiv

Fast computation of Gauss quadrature nodes and weights on the whole real line

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

preprint2014arXiv

On the use of conformal maps for the acceleration of convergence of the trapezoidal rule and Sinc numerical methods

We investigate the use of conformal maps for the acceleration of convergence of the trapezoidal rule and Sinc numerical methods. The conformal map is a polynomial adjustment to the $\sinh$ map, and allows the treatment of a finite number of singularities in the complex plane. In the case where locations are unknown, the so-called Sinc-Padé approximants are used to provide approximate results. This adaptive method is shown to have almost the same convergence properties. We use the conformal maps to generate high accuracy solutions to several challenging integrals, nonlinear waves, and multidimensional integrals.

preprint2014arXiv

Sampling unitary invariant ensembles

We develop an algorithm for sampling from the unitary invariant random matrix ensembles. The algorithm is based on the representation of their eigenvalues as a determinantal point process whose kernel is given in terms of orthogonal polynomials. Using this algorithm, statistics beyond those known through analysis are calculable through Monte Carlo simulation. Unexpected phenomena are observed in the simulations.

preprint2014arXiv

Universality in Numerical Computations with Random Data. Case Studies

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.

preprint2013arXiv

A Riemann--Hilbert approach to Jacobi operators and Gaussian quadrature

The computation of the entries of Jacobi operators associated with orthogonal polynomials has important applications in numerical analysis. From truncating the operator to form a Jacobi matrix, one can apply the Golub--Welsh algorithm to compute the Gaussian quadrature weights and nodes. Furthermore, the entries of the Jacobi operator are the coefficients in the three-term recurrence relationship for the polynomials. This provides an efficient method for evaluating the orthogonal polynomials. Here, we present an $\mathcal O(N)$ method to compute the first $N$ rows of Jacobi operators from the associated weight. The method exploits the Riemann--Hilbert representation of the polynomials by solving a deformed Riemann--Hilbert problem numerically. We further adapt this computational approach to certain entire weights that are beyond the reach of current asymptotic Riemann--Hilbert techniques.

preprint2013arXiv

Change of variable formulæ for regularizing slowly decaying and oscillatory Cauchy and Hilbert transforms

Formulæ are derived for expressing Cauchy and Hilbert transforms of a function $f$ in terms of Cauchy and Hilbert transforms of $f(x^r)$. When $r$ is an integer, this corresponds to evaluating the Cauchy transform of $f(x^r)$ at all choices of $z^{1/r}$. Related formuæ for rational $r$ result in a reduction to a generalized Cauchy transform living on a Riemann surface, which in turn is reducible to the standard Cauchy transform. These formulæ are used to regularize the behaviour of functions that are slowly decaying or oscillatory, in order to facilitate numerical computation and extend asymptotic results.

preprint2013arXiv

Fast inverse transform sampling in one and two dimensions

We develop a computationally efficient and robust algorithm for generating pseudo-random samples from a broad class of smooth probability distributions in one and two dimensions. The algorithm is based on inverse transform sampling with a polynomial approximation scheme using Chebyshev polynomials, Chebyshev grids, and low rank function approximation. Numerical experiments demonstrate that our algorithm outperforms existing approaches.

preprint2013arXiv

Numerical computation of convolutions in free probability theory

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.

preprint2012arXiv

A fast and well-conditioned spectral method

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.

preprint2012arXiv

Nonlinear steepest descent and the numerical solution of Riemann-Hilbert problems

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.

preprint2012arXiv

Numerical solution of Riemann--Hilbert problems: random matrix theory and orthogonal polynomials

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.

preprint2011arXiv

Numerical study of higher order analogues of the Tracy-Widom distribution

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.