§10.74 Methods of Computation

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 $x$ or $z$ is sufficiently small in absolute value. In the case of the modified Bessel function $K_{\nu }\left(z\right)$ see especially Temme (1975).

In other circumstances the power series are prone to slow convergence and heavy numerical cancellation.

If $x$ or $|z|$ is large compared with ${|\nu |}^2$, 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 $\nu$ 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 $x$ or $|z|$, whether or not $\nu$ is large. It should be noted, however, that there is a difficulty in evaluating the coefficients $A_k(\zeta )$, $B_k(\zeta )$, $C_k(\zeta )$, and $D_k(\zeta )$, from the explicit expressions (10.20.10)–(10.20.13) when $z$ is close to $1$ 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 $z$; however, there may be severe cancellation when $|z|$ is not large compared with $n^2$.

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 , $J_{\nu }\left(x\right)$ needs to be integrated in the forward direction and $Y_{\nu }\left(x\right)$ 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 $I_{\nu }\left(x\right)$ and backwards in the case of $K_{\nu }\left(x\right)$, with initial values obtained in an analogous manner to those for $J_{\nu }\left(x\right)$ and $Y_{\nu }\left(x\right)$.

For $z\in ℂ$ the function $H_{\nu }^{(1)}\left(z\right)$, for example, can always be computed in a stable manner in the sector $0\le \mathrm{ph}z\le \pi$ 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 evaluation of the Hankel functions $H_{\nu }^{(1)}\left(z\right)$ and $H_{\nu }^{(2)}\left(z\right)$ for complex values of $\nu$ and $z$ based on the integral representations (10.9.18) see Remenets (1973).

For applications of generalized Gauss–Laguerre quadrature (§3.5(v)) to the evaluation of the modified Bessel functions $K_{\nu }\left(z\right)$ for and see Gautschi (2002a). The integral representation used is based on (10.32.8).

For evaluation of $K_{\nu }\left(z\right)$ from (10.32.14) with $\nu =n$ and $z$ complex, see Mechel (1966).

If values of the Bessel functions $J_{\nu }\left(z\right)$, $Y_{\nu }\left(z\right)$, or the other functions treated in this chapter, are needed for integer-spaced ranges of values of the order $\nu$, then a simple and powerful procedure is provided by recurrence relations typified by the first of (10.6.1).

Suppose, for example, $\nu =n\in 0,1,2,\mathrm{\dots }$, and $x\in (0,\mathrm{\infty })$. Then $J_n\left(x\right)$ and $Y_n\left(x\right)$ can be generated by either forward or backward recurrence on $n$ when , but if $n>x$ then to maintain stability $J_n\left(x\right)$ has to be generated by backward recurrence on $n$, and $Y_n\left(x\right)$ has to be generated by forward recurrence on $n$. In the case of $J_n\left(x\right)$, 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 further information see Gautschi (1967), Olver and Sookne (1972), Temme (1975), Campbell (1980), and Kerimov and Skorokhodov (1984a).

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 $2D$ graphs in $ℝ$ (e.g., §10.3(i)) with the $x$-axis, or graphical intersection of $3D$ complex-variable surfaces (e.g., §10.3(ii)) with the plane $z=0$.

To ensure that no zeros are overlooked, standard tools are the phase principle and Rouché’s theorem; see §1.10(iv).

For the computation of the functions ${\tilde{I}}_{\nu }\left(x\right)$ and ${\tilde{K}}_{\nu }\left(x\right)$ defined by (10.45.2) see Temme (1994c) and Gil et al. (2002d, 2003a, 2004b).