We describe a method for the rapid numerical evaluation of the Bessel functions of the first and second kinds of nonnegative real orders and positive arguments. Our algorithm makes use of the well-known observation that although the Bessel functions themselves are expensive to represent via piecewise polynomial expansions, the logarithms of certain solutions of Bessel's equation are not. We exploit this observation by numerically precomputing the logarithms of carefully chosen Bessel functions and representing them with piecewise bivariate Chebyshev expansions. Our scheme is able to evaluate Bessel functions of orders between 0 and $1\sep,000\sep,000\sep,000$ at essentially any positive real argument. In that regime, it is competitive with existing methods for the rapid evaluation of Bessel functions and has several advantages over them. First, our approach is quite general and can be readily applied to many other special functions which satisfy second order ordinary differential equations. Second, by calculating the logarithms of the Bessel functions rather than the Bessel functions themselves, we avoid many issues which arise from numerical overflow and underflow. Third, in the oscillatory regime, our algorithm calculates the values of a nonoscillatory phase function for Bessel's differential equation and its derivative. These quantities are useful for computing the zeros of Bessel functions, as well as for rapidly applying the Fourier-Bessel transform. The results of extensive numerical experiments demonstrating the efficacy of our algorithm are presented. A Fortran package which includes our code for evaluating the Bessel functions as well as our code for all of the numerical experiments described here is publically available.

We describe a method for calculating the roots of special functions satisfying second order linear ordinary differential equations. It exploits the recent observation that the solutions of a large class of such equations can be represented via nonoscillatory phase functions, even in the high-frequency regime. Our algorithm achieves near machine precision accuracy and the time required to compute one root of a solution is independent of the frequency of oscillations of that solution. Moreover, despite its great generality, our approach is competitive with specialized, state-of-the-art methods for the construction of Gaussian quadrature rules of large orders when it used in such a capacity. The performance of the scheme is illustrated with several numerical experiments and a Fortran implementation of our algorithm is available at the author's website.

We describe an algorithm for the numerical solution of second order linear differential equations in the highly-oscillatory regime. It is founded on the recent observation that the solutions of equations of this type can be accurately represented using nonoscillatory phase functions. Unlike standard solvers for ordinary differential equations, the running time of our algorithm is independent of the frequency of oscillation of the solutions. We illustrate the performance of the method with several numerical experiments.

We express a certain complex-valued solution of Legendre's differential equation as the product of an oscillatory exponential function and an integral involving only nonoscillatory elementary functions. By calculating the logarithmic derivative of this solution, we show that Legendre's differential equation admits a nonoscillatory phase function. Moreover, we derive from our expression an asymptotic expansion useful for evaluating Legendre functions of the first and second kinds of large orders, as well as the derivative of the nonoscillatory phase function. Numerical experiments demonstrating the properties of our asymptotic expansion are presented.

Recently, it was observed that solutions of a large class of highly oscillatory second order linear ordinary differential equations can be approximated using nonoscillatory phase functions. In particular, under mild assumptions on the coefficients and wavenumber $\lambda$ of the equation, there exists a function whose Fourier transform decays as $\exp(-\mu |\xi|)$ and which represents solutions of the differential equation with accuracy on the order of $\lambda^{-1} \exp(-\mu \lambda)$. In this article, we establish an improved existence theorem for nonoscillatory phase functions. Among other things, we show that solutions of second order linear ordinary differential equations can be represented with accuracy on the order of $\lambda^{-1} \exp(-\mu \lambda)$ using functions in the space of rapidly decaying Schwartz functions whose Fourier transforms are both exponentially decaying and compactly supported. These new observations play an important role in the analysis of a method for the numerical solution of second order ordinary differential equations whose running time is independent of the parameter $\lambda$. This algorithm will be reported at a later date.

We observe that solutions of a large class of highly oscillatory second order linear ordinary differential equations can be approximated using nonoscillatory phase functions. In addition, we describe numerical experiments which illustrate important implications of this fact. For example, that many special functions of great interest --- such as the Bessel functions $J_\nu$ and $Y_\nu$ --- can be evaluated accurately using a number of operations which is $O(1)$ in the order $\nu$. The present paper is devoted to the development of an analytical apparatus. Numerical aspects of this work will be reported at a later date.

We introduce a version of the asymptotic expansions for Bessel functions $J_\nu(z)$, $Y_\nu(z)$ that is valid whenever $|z| > \nu$ (which is deep in the Fresnel regime), as opposed to the standard expansions that are applicable only in the Fraunhofer regime (i.e. when $|z| > \nu^2$). As expected, in the Fraunhofer regime our asymptotics reduce to the classical ones. The approach is based on the observation that Bessel's equation admits a non-oscillatory phase function, and uses classical formulas to obtain an asymptotic expansion for this function; this in turn leads to both an analytical tool and a numerical scheme for the efficient evaluation of $J_\nu(z)$, $Y_\nu(z)$, as well as various related quantities. The effectiveness of the technique is demonstrated via several numerical examples. We also observe that the procedure admits far-reaching generalizations to wide classes of second order differential equations, to be reported at a later date.