A Fast Algorithm for the Calculation of the Roots of Special Functions

Abstract

We describe a procedure for the determination of the roots of functions satisfying second-order ordinary differential equations, including the classical special functions. The scheme is based on a combination of the Prüfer transform with the classical Taylor series method for the solution of ordinary differential equations and requires O(1) operations for the determination of each root. When the functions in question are classical orthogonal polynomials (Legendre, Hermite, etc.), the techniques presented here also provide tools for the evaluation of the weights for the associated Gaussian quadratures. The performance of the scheme for several classical special functions (prolate spheroidal wave functions, Bessel functions, and Legendre, Hermite, and Laguerre polynomials) is illustrated with numerical examples. Introduction.Zeros of classical special functions play an important role in the mathematical sciences, being related to Gaussian quadratures (both classical and generalized), resonances in mechanical and electrical systems, quantum mechanical calculations, etc. As a result, algorithms for the determination of roots of functions of this type are quite well developed; most of them are based on the existence of three-term recurrence relations for the classical special functions and the associated connection to the eigenvalues of certain tridiagonal matrices (see, for example, [6,9,2]). The latter are dealt with via well-understood and reliable techniques, such as the QR iteration.While such methods are quite efficient for small-scale problems, they require order n 2 operations, where n is the number of roots to be computed. Relatively recently, algorithms for the diagonalization of tridiagonal matrices were introduced with CPU time requirements proportional to n · log(n) (see [7]); a detailed description of the resulting fast root-finding scheme is given in [13]. Another relatively recent development involves straightforward root-finding algorithms applied to the orthogonal polynomials (see [15,11]); this leads to more efficient schemes producing nodes and weights to higher precision. However, the cost of the resulting procedures is still of the order n 2 , albeit with a smaller constant. The scheme of this paper is a modification of this latter approach using the solution of the underlying ordinary differential equation (ODE) as a tool for the evaluation of the special function whose roots are to be found.The subject of this paper is the observation that the roots of classical special functions (and of many others) can be found for a cost proportional to n, using the

Fast and Accurate Computation of Gauss--Legendre and Gauss--Jacobi Quadrature Nodes and Weights

“…The current state of the art is the Glasier-Liu-Rohklin (GLR) algorithm [21], which computes all the nodes and weights of the n-point quadrature rule in a total of O(n) operations. This also employs Newton's method, but here the function and derivative evaluations are computed with n-independent complexity by sequentially hopping from one node to the next using a local 30-term Taylor series approximation generated from the second-order differential equation satisfied by the orthogonal polynomial.…”

Linear barycentric rational quadrature

“…This may seem expensive at first sight. However the Chebfun [3,28] command legpts, an implementation of the method introduced in [16], provides a very fast algorithm. Moreover, we had to compute the Legendre points and weights only once for all the examples we investigated.…”

Computation of connection coefficients and measure modifications for orthogonal polynomials

“…Notable exceptions are special Chebyshev-like families when explicit, simple expressions are known for the quadrature nodes and weights. However, asymptotic O(N) algorithms exist for any ω, at least for the determination of the nodes [13].…”