An algorithm for the rapid numerical evaluation of Bessel functions of real orders and arguments

Abstract

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. We supplement these precomputed expansions with two asymptotic expansions, one for large orders and extremely small arguments and the other for large orders and extremely large arguments, and with series expansions for small orders and small arguments. Our scheme is able to evaluate Bessel functions of orders between 0 and 1, 000, 000, 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 at least three 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.

Interpolative Decomposition Butterfly Factorization

“…It is demonstrated in [13] that (31) can be represented via a matvec u = Kg, where K satisfies the complementary low-rank property. An arbitrary entry of K can be calculated in O(1) operations [44] and hence IDBF is suitable for accelerating the matvec u = Kg.…”

A hierarchical butterfly LU preconditioner for two-dimensional electromagnetic scattering problems involving open surfaces

“…The analysis for other off-diagonal blocks is similar. By the theory of non-oscillatory phase functions for special functions [50,51,52], the Hankel function H…”

“…where α(x) is an analytic non-oscillatory phase function. For more details about the definition, computation, and asymptotic expansion of the non-oscillatory phase function, the reader is referred to Section 2.5 of [52], where the notation α(t) is used to represent the phase function as well. If we introduce a non-oscillatory matrix B and a purely oscillatory matrix C as follows,…”

“…since ρ is nearly isometric. Note that α(x) is non-oscillatory and analytic with the property (see Equations (57) in [52]) that lim x→∞ α (x) = 1.…”

A unified framework for oscillatory integral transforms: When to use NUFFT or butterfly factorization?

“…When the numerical rank of the phase function r is only larger than the dimension of the problem by one or two, NUFFT is usually faster than BF and hence it will be applied to compute Kf . Scenario 1 : There exists an algorithm for evaluating an arbitrary entry of the kernel matrix in O(1) operations [3,4,19,25].…”