An algorithm for the rapid numerical evaluation of Bessel functions of real orders and arguments (Q2631984): Difference between revisions

The reviewer would like to begin the review of this paper by first quoting the following beginning lines from a most recently published paper of \textit{P. Greengard} and \textit{V. Rokhlin} [Adv. Comput. Math. 45, No. 1, 23--49 (2019, Zbl 1428.33002)]: ``The evaluation of special functions is one of the most developed areas of numerical analysis. For some special functions, such as Bessel functions, the theory has been fairly complete for many decades; for others, such as Prolate Spheroidal Wave Functions, the theory is still an active area of research.'' These quotes highlight the importance of the subject dealt with in this paper where the author presents a very powerful algorithm which is capable of numerically evaluating `the Bessel functions of the first kind \(J_\nu\) and the second kind \(Y_\nu\) of nonnegative real orders [\dots] between 0 and 1,000,000,000 at essentially any positive real argument.' Some very well known schemes for the numerical evaluation of the Bessel functions are currently in vogue, like [\textit{D. E. Amos}, ACM Trans. Math. Softw. 12, 265--273 (1986; Zbl 0613.65013 )] and [\textit{G. Matviyenko}, Appl. Comput. Harmon. Anal. 1, No. 1, 116--135 (1993; Zbl 0802.33004)], yet the algorithm presented in this paper outshines them in certain respects. The most beautiful aspect of the algorithm of the author is that he presents here a scheme which represents `an incredibly straightforward approach' that has the potential of being applicable for the numerical evaluation of `a large class of special functions satisfying second-order ordinary differential equations', and he first of all ventures to apply this scheme only for the case of Bessel functions because, the only reason for that, according to him, is the availability of the previously well-known and widely used algorithms of Amos [loc. cit.] and Matviyenko [loc. cit.] with the results of which the outcomes of the proposed scheme could be compared to test its efficacy and appropriateness. The author very confidently and aptly remarks that `the results of applying this approach to other classes of special functions, such as the associated Legendre functions and prolate spheroidal wave functions, will be reported by the author at a later date.' Given that the scaled Bessel functions \({J_\nu }\left( t \right)\sqrt t \) and \({Y_\nu }\left( t \right)\sqrt t \) are solutions of the second-order linear ordinary differential equation \(y''\left( t \right) + \left( {1 - \frac{{{\nu ^2} - \frac{1}{4}}}{{{t^2}}}} \right)y\left( t \right) = 0\) for \(t \in \left( {0,\infty } \right),\) the author refers to the subsets \[ \mathcal{ O} = \left\{ {\left( {\nu ,t} \right):0 \le \nu \le \frac{1}{2}\;\text{{ and}}\;\text{{ }}t > 0} \right\} \cup \left\{ {\left( {\nu ,t} \right):\nu > \frac{1}{2}\;\text{{ and }}\;t \ge \sqrt {{\nu ^2} - \frac{1}{4}} } \right\} \] and \[\mathcal{ N} = \left\{ {\left( {\nu ,t} \right):\nu > \frac{1}{2}\;\text{{ and }}\;\text{{0 < }}t < \sqrt {{\nu ^2} - \frac{1}{4}} } \right\} \] of \(\mathbb{R} \times \mathbb{R}\) respectively as the oscillatory region and the nonoscillatory region. For large \(\nu\), the Bessel functions cannot be efficiently represented by polynomial expansions in \(\mathcal N\) as `they behave like rapidly increasing or decreasing exponentials' and the same is true for any large subset of \(\mathcal O\) where the Bessel polynomials oscillate. The long familiar fact of the Bessel differential equation admitting a nonoscillatory phase function (see, [\textit{G. N. Watson}, A treatise on the theory of Bessel functions. 2nd ed. New York: Cambridge University Press (1995; Zbl 0849.33001)]) guarantees the existence of a solution of Bessel's differential equation whose logarithm is nonoscillatory, therefore, this solution can be efficiently represented by polynomial expansions on large subsets of \(\mathcal O\), whereas, in the region \(\mathcal N\) `the logarithms of Bessel functions can be represented efficiently via polynomial expansions'. The key idea of the algorithm developed in this paper is to operate `by numerically precomputing the logarithms of certain (carefully chosen) solutions of Bessel's differential equation.' This approach, in fact, is a merit of the algorithm presented in this study because the recourse to this method avoids `many issues which arise from numerical overflow and underflow'. As `the nonoscilatory phase function is the imaginary part of the logarithm of a solution of Bessel's differential equation', these functions are represented `via piecewise bivariate Chebyshev expansions' and the resulting coefficients are arranged in a table which is `stored on the disk and loaded into memory when needed'. These precomputed expansions are supplemented by the known asymptotic and series expansions of Bessel expansions available in the literature for evaluation of the \(J_\nu\) and \(Y_\nu\) for nonegative values of \(\nu\) up to one billion and at positive real arguments. Wherever the asymptotic and series expansions of Bessel functions are not available in the literature, the author has made arrangements in the scheme of precomputed expansions to suitably cover the range of the parameter and argument so as to suffice to meet his objectives of writing this paper. One more noteworthy point is that for the very wide variety of numerical experiments reported in this paper, the size of the precomputed table is only about 1.3 megabytes. \par Talking about the structure of the paper, the necessary mathematical preliminaries and required numerical procedures are detailed in Section 2. A robust adaptive solver for nonlinear differential equations is developed in the third section which in turn is used by the algorithm described in Section 4. This algorithm proceeds through three stages to rapidly evaluate the numerical solutions of the Bessel's differential equation for a fixed value of \(\nu\) and is an essential precursor for the procedure described in Section 5 for the generation of the precomputed table, which `stores the coefficients in the piecewise compressed bivariate Chebyshev expansions of several functions'. The author's core algorithm is described in Section 6 and is available as a Fortran package for the interested reader from the websites \url{http://github.com/JamesCBremerJr/BesselEval} or, \url{http://www.math.ucdavis.edu/~bremer/code.html}. Extensive numerical experiments to describe the powerful properties of this novel algorithm are described in Section 7 where one very important fact, which the reviewer is most strongly inclined to point out here, is that the results of the Table 5 most clearly show that whereas the Amos algorithm [Amos, loc. cit.] may be slightly more accurate (with lesser maximum relative error) but at the expense of a slightly higher time of computation than compared with the author's algorithm for the values of \(\nu\) ranging from 0 to 100,000, but it fails altogether for numerically evaluating the Bessel functions for the values of \(\nu\) beyond 100,000 where it returns an abort message `and returns an error code in cases in which it is unable to evaluate the Bessel functions to at least 7-digit accuracy', the author's code works very well with an appreciable accuracy for values of \(\nu\) above 100,000 and up to 1,000,000,000! Another prominent merit of the algorithm of this study is that in the oscillatory region besides returning the values of the Bessel function it `also returns the values of a nonoscillatory phase function for Bessel's equation and its derivative', which is of vital importance not only for the computation of the zeros of special functions (see [\textit{J. Bremer}, SIAM J. Sci. Comput. 39, No. 1, A55--A82 (2017; Zbl 1355.65037)]) but also for applying special function transforms via the butterfly algorithm (see, for instance, [\textit{E. Candés} et al., Multiscale Model. Simul. 7, No. 4, 1727--1750 (2009; Zbl 1184.65125); \textit{Y. Li} et al., Multiscale Model. Simul. 13, No. 2, 714--732 (2015; Zbl 1317.44004); \textit{M. O'Neil} et al., Appl. Comput. Harmon. Anal. 28, No. 2, 203--226 (2010; Zbl 1191.65016)]). Commenting upon the ability of Amos' code [loc. cit.] regarding its power to numerically evaluate the Bessel functions of complex arguments, the author remarks that he would, `at a later date', announce an extension of this algorithm which would be capable of efficiently evaluating numerically the Bessel functions in the upper half of the complex plane, by exploiting the fact that `a phase function which is nonoscillatory on the upper half of the complex plane' is admitted by the Bessel's differential equation. He concludes the paper by saying that: ``The author will report on the use of the method of this paper to evaluate associated Legendre functions and prolate spheroidal wave functions at a later date, as well as on the rapid application of special function transforms using techniques related to those discussed here.'' \par In the opinion of the reviewer these pertinent remarks of the author which coincide so brilliantly with the opening quotes cited at the very beginning of this review should not come as a surprise to any reader because the author of this paper and the second author of the aforesaid paper (V. Rokhlin) have been recent research collaborators (see, for example, [\textit{Z. Heitman} et al., Appl. Comput. Harmon. Anal. 39, No. 2, 347--356 (2015; Zbl 1321.33007)]). The reviewer hereby most strongly recommends this excellent work and the very recent work of Greengard and Rokhlin [loc. cit.] to the researchers working in the field of numerical evaluation of special functions with his disclosure in the conclusion that this reviewer has also been strongly inspired and motivated by his study of these two remarkable works in his quest for the numerical zeroes of two special cases of the degenerate confluent hypergeometric function and the degenerate Gauss's hypergeometric function which are introduced by the present reviewer in his just concluded study.