The purpose of this article is to provide some insight on the efficiency and accuracy of a few available algorithms for the computation of integer-order Bessel functions First, the computation of integer-order Bessel hnctions of the first kind (J,,(z)), using the Fast Fourier Transform (FFT) algorithm as opposed to recurrence techniques, is investigated It is shown that recurrence techniques are superior to the FFT technique, both in accuracy and speed efficiency Second, an erroneous algorithm, suggested in the literature and used by commercially available software, specifically MA TLAR 3.5 and MATHFNA AKA 1.2, for computing integer-order Bessel functions of the second kind ( < , ( z ) ) , is revealed by comparing these routines with an algorithm developed by the author Catastrophic errors result from the use of the erroneous algorithm, for the computation of large orders with non-real arguments
IntroductionPhysical problems, formulated in a cylindrical-or sphericalcoordinate system, have prompted the development of various methods for the computation of the Bessel hnctions J , , ( z ) and T,(z), for integer order, I?, and a complex or real argument, z. To be able to use them successfully, it is important to know the numerical limitations of these algorithms. In terms of flexibility and efficiency, the methods may be divided roughly into two groups. The first group describes methods where, due to numerical problems, severe restrictions are placed on the magnitude of the argument and the maximum order [ I , 2, 3; 4, Appendix 11. The second group describes methods which allow for the computation of a range of orders, using either recurrence techniques [5, p. 385; 6, 7, 8, 9, 10, 1 1 , 121 or alternative methods, of which the use of the FFT algorithm [13, p. 2901, in particular, is an interesting example [14, IS].In Section 3, the accuracy and efficiency of the FFT method is compared to recurrence techniques. The problems of determining the starting point of backward recurrence, and the minimum number of sample points required by the FFT method, are also discussed.Finally, the problem of computing Y,,(z) for complex arguments is addressed in Section 4, and the accuracy of a few commercial mathematical routines, namely MATLAB 3.5 and MA THE-MATICA versions 1.2 and 2.0, is investigated. Forward recurrence is usually recommended [5,6,7, 81, to compute all orders of E, from U, and E;. Recurrence is stable in the direction of increasing function values only, hence this procedure is valid for real arguments. For complex arguments, however, forward recurrence becomes unstable, since I~, ( z ) l decreases initially with increasing ii, before a turning point is reached (see Figure 4), and it increases only for larger orders beyond the turning point. Apparently, this erroneous algorithm is employed by MATZAB 3.5 and MATHE-MATICA I . 2 . Since the answers are presumably accepted in good faith by the scientific community, it is evidently quite a serious problem.
Comparison between the FFT method and recurrence technique...