About the Project
10 Bessel FunctionsComputation

§10.74 Methods of Computation

  1. §10.74(i) Series Expansions
  2. §10.74(ii) Differential Equations
  3. §10.74(iii) Integral Representations
  4. §10.74(iv) Recurrence Relations
  5. §10.74(v) Continued Fractions
  6. §10.74(vi) Zeros and Associated Values
  7. §10.74(vii) Integrals
  8. §10.74(viii) Functions of Imaginary Order

§10.74(i) Series Expansions

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ν(z) 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 |ν|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 ν 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 ν is large. It should be noted, however, that there is a difficulty in evaluating the coefficients Ak(ζ), Bk(ζ), Ck(ζ), and Dk(ζ), 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 n2.

§10.74(ii) Differential Equations

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 0<x<ν, Jν(x) needs to be integrated in the forward direction and Yν(x) 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 ν<x< either direction of integration can be used for both functions.

Similarly, to maintain stability in the interval 0<x< the integration direction has to be forwards in the case of Iν(x) and backwards in the case of Kν(x), with initial values obtained in an analogous manner to those for Jν(x) and Yν(x).

For z the function Hν(1)(z), for example, can always be computed in a stable manner in the sector 0phzπ 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).

§10.74(iii) Integral Representations

For evaluation of the Hankel functions Hν(1)(z) and Hν(2)(z) for complex values of ν 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ν(z) for 0<ν<1 and 0<x< see Gautschi (2002a). The integral representation used is based on (10.32.8).

For evaluation of Kν(z) from (10.32.14) with ν=n and z complex, see Mechel (1966).

§10.74(iv) Recurrence Relations

If values of the Bessel functions Jν(z), Yν(z), 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, ν=n0,1,2,, and x(0,). Then Jn(x) and Yn(x) can be generated by either forward or backward recurrence on n when n<x, but if n>x then to maintain stability Jn(x) has to be generated by backward recurrence on n, and Yn(x) has to be generated by forward recurrence on n. In the case of Jn(x), 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).

§10.74(v) Continued Fractions

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).

§10.74(vi) Zeros and Associated Values

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).

Real Zeros

See Olver (1960, pp. xvi–xxix), Grad and Zakrajšek (1973), Temme (1979a), Ikebe et al. (1991), Zafiropoulos et al. (1996), Vrahatis et al. (1997a), Ball (2000), and Gil and Segura (2003).

Complex Zeros

See Leung and Ghaderpanah (1979), Kerimov and Skorokhodov (1984b, c, 1985a, 1985b), Skorokhodov (1985), Modenov and Filonov (1986), Vrahatis et al. (1997b), and Segura (2013).

Multiple Zeros

See Kerimov and Skorokhodov (1985c, 1986, 1987, 1988).

§10.74(vii) Integrals

Hankel Transform

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).

Fourier–Bessel Expansion

For the computation of the integral (10.23.19) see Piessens and Branders (1983, 1985), Lewanowicz (1991), and Zhileĭkin and Kukarkin (1995).

Spherical Bessel Transform

The spherical Bessel transform is the Hankel transform (10.22.76) in the case when ν is half an odd positive integer.

See Lehman et al. (1981), Puoskari (1988), and Sharafeddin et al. (1992).

Kontorovich–Lebedev Transform

See Ehrenmark (1995).


For infinite integrals involving products of Bessel functions of the first kind, see Linz and Kropp (1973), Gabutti (1980), Ikonomou et al. (1995), Lucas (1995), and Van Deun and Cools (2008). For infinite integrals involving products of Bessel functions of the first and second kinds, see Ratnanather et al. (2014).

§10.74(viii) Functions of Imaginary Order

For the computation of the functions I~ν(x) and K~ν(x) defined by (10.45.2) see Temme (1994c) and Gil et al. (2002d, 2003a, 2004b).