The power-series expansions given in §§10.2 and 10.8, together with the connection formulas of §10.4, can be used to compute the Bessel and Hankel functions when the argument or is sufficiently small in absolute value. In the case of the modified Bessel function see especially Temme (1975).
In other circumstances the power series are prone to slow convergence and heavy numerical cancellation.
If or is large compared with , then the asymptotic expansions of §§10.17(i)–10.17(iv) are available. Furthermore, the attainable accuracy can be increased substantially by use of the exponentially-improved expansions given in §10.17(v), even more so by application of the hyperasymptotic expansions to be found in the references in that subsection.
For large positive real values of the uniform asymptotic expansions of §§10.20(i) and 10.20(ii) can be used. Moreover, because of their double asymptotic properties (§10.41(v)) these expansions can also be used for large or , whether or not is large. It should be noted, however, that there is a difficulty in evaluating the coefficients , , , and , from the explicit expressions (10.20.10)–(10.20.13) when is close to owing to severe cancellation. Temme (1997) shows how to overcome this difficulty by use of the Maclaurin expansions for these coefficients or by use of auxiliary functions.
Similar observations apply to the computation of modified Bessel functions, spherical Bessel functions, and Kelvin functions. In the case of the spherical Bessel functions the explicit formulas given in §§10.49(i) and 10.49(ii) are terminating cases of the asymptotic expansions given in §§10.17(i) and 10.40(i) for the Bessel functions and modified Bessel functions. And since there are no error terms they could, in theory, be used for all values of ; however, there may be severe cancellation when is not large compared with .
A comprehensive and powerful approach is to integrate the differential equations (10.2.1) and (10.25.1) by direct numerical methods. As described in §3.7(ii), to insure stability the integration path must be chosen in such a way that as we proceed along it the wanted solution grows in magnitude at least as fast as all other solutions of the differential equation.
In the interval , needs to be integrated in the forward direction and in the backward direction, with initial values for the former obtained from the power-series expansion (10.2.2) and for the latter from asymptotic expansions (§§10.17(i) and 10.20(i)). In the interval either direction of integration can be used for both functions.
Similarly, to maintain stability in the interval the integration direction has to be forwards in the case of and backwards in the case of , with initial values obtained in an analogous manner to those for and .
For the function , for example, can always be computed in a stable manner in the sector by integrating along rays towards the origin.
Similar considerations apply to the spherical Bessel functions and Kelvin functions.
For further information, including parallel methods for solving the differential equations, see Lozier and Olver (1993).
For applications of generalized Gauss–Laguerre quadrature (§3.5(v)) to the evaluation of the modified Bessel functions for and see Gautschi (2002a). The integral representation used is based on (10.32.8).
If values of the Bessel functions , , or the other functions treated in this chapter, are needed for integer-spaced ranges of values of the order , then a simple and powerful procedure is provided by recurrence relations typified by the first of (10.6.1).
Suppose, for example, and . Then and can be generated by either forward or backward recurrence on when , but if then to maintain stability has to be generated by backward recurrence on , and has to be generated by forward recurrence on . In the case of , the need for initial values can be avoided by application of Olver’s algorithm (§3.6(v)) in conjunction with Equation (10.12.4) used as a normalizing condition, or in the case of noninteger orders, (10.23.15).
For applications of the continued-fraction expansions (10.10.1), (10.10.2), (10.33.1), and (10.33.2) to the computation of Bessel functions and modified Bessel functions see Gargantini and Henrici (1967), Amos (1974), Gautschi and Slavik (1978), Tretter and Walster (1980), Thompson and Barnett (1986), and Cuyt et al. (2008).
Newton’s rule (§3.8(i)) or Halley’s rule (§3.8(v)) can be used to compute to arbitrarily high accuracy the real or complex zeros of all the functions treated in this chapter. Necessary values of the first derivatives of the functions are obtained by the use of (10.6.2), for example. Newton’s rule is quadratically convergent and Halley’s rule is cubically convergent. See also Segura (1998, 2001).
Methods for obtaining initial approximations to the zeros include asymptotic expansions (§§10.21(vi)-10.21(ix)), graphical intersection of graphs in (e.g., §10.3(i)) with the -axis, or graphical intersection of complex-variable surfaces (e.g., §10.3(ii)) with the plane .
To ensure that no zeros are overlooked, standard tools are the phase principle and Rouché’s theorem; see §1.10(iv).
See Cornille (1972), Johansen and Sørensen (1979), Gabutti (1979), Gabutti and Minetti (1981), Candel (1981), Wong (1982), Lund (1985), Piessens and Branders (1985), Hansen (1985), Bezvoda et al. (1986), Puoskari (1988), Christensen (1990), Campos (1995), Lucas and Stone (1995), Barakat and Parshall (1996), Sidi (1997), Secada (1999).
The spherical Bessel transform is the Hankel transform (10.22.76) in the case when is half an odd positive integer.
See Ehrenmark (1995).