Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

Modified Bessel Functions of the First and Second Kinds

Synopsis

#include <boost/math/special_functions/bessel.hpp>

template <class T1, class T2>
calculated-result-type cyl_bessel_i(T1 v, T2 x);

template <class T1, class T2, class Policy>
calculated-result-type cyl_bessel_i(T1 v, T2 x, const Policy&);

template <class T1, class T2>
calculated-result-type cyl_bessel_k(T1 v, T2 x);

template <class T1, class T2, class Policy>
calculated-result-type cyl_bessel_k(T1 v, T2 x, const Policy&);
Description

The functions cyl_bessel_i and cyl_bessel_k return the result of the modified Bessel functions of the first and second kind respectively:

cyl_bessel_i(v, x) = Iv(x)

cyl_bessel_k(v, x) = Kv(x)

where:

The return type of these functions is computed using the result type calculation rules when T1 and T2 are different types. The functions are also optimised for the relatively common case that T1 is an integer.

The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.

The functions return the result of domain_error whenever the result is undefined or complex. For cyl_bessel_j this occurs when x < 0 and v is not an integer, or when x == 0 and v != 0. For cyl_neumann this occurs when x <= 0.

The following graph illustrates the exponential behaviour of Iv.

The following graph illustrates the exponential decay of Kv.

Testing

There are two sets of test values: spot values calculated using functions.wolfram.com, and a much larger set of tests computed using a simplified version of this implementation (with all the special case handling removed).

Accuracy

The following tables show how the accuracy of these functions varies on various platforms, along with comparison to other libraries. Note that only results for the widest floating-point type on the system are given, as narrower types have effectively zero error. All values are relative errors in units of epsilon. Note that our test suite includes some fairly extreme inputs which results in most of the worst problem cases in other libraries:

Table 8.44. Error rates for cyl_bessel_i (integer orders)

GNU C++ version 7.1.0
linux
double

GNU C++ version 7.1.0
linux
long double

Sun compiler version 0x5150
Sun Solaris
long double

Microsoft Visual C++ version 14.1
Win32
double

Bessel I0: Mathworld Data (Integer Version)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 0.79ε (Mean = 0.482ε))
(Rmath 3.2.3: Max = 1.52ε (Mean = 0.622ε) And other failures.)

Max = 1.95ε (Mean = 0.738ε)

(<cmath>: Max = 8.49ε (Mean = 3.46ε) And other failures.)

Max = 1.95ε (Mean = 0.661ε)

Max = 0.762ε (Mean = 0.329ε)

Bessel I1: Mathworld Data (Integer Version)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 0.82ε (Mean = 0.456ε))
(Rmath 3.2.3: Max = 1.53ε (Mean = 0.483ε) And other failures.)

Max = 0.64ε (Mean = 0.202ε)

(<cmath>: Max = 5ε (Mean = 2.15ε) And other failures.)

Max = 0.64ε (Mean = 0.202ε)

Max = 0.767ε (Mean = 0.398ε)

Bessel In: Mathworld Data (Integer Version)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 5.15ε (Mean = 2.13ε) And other failures.)
(Rmath 3.2.3: Max = 1.73ε (Mean = 0.601ε) And other failures.)

Max = 1.8ε (Mean = 1.33ε)

(<cmath>: Max = 430ε (Mean = 163ε) And other failures.)

Max = 463ε (Mean = 140ε)

Max = 3.46ε (Mean = 1.32ε)


Table 8.45. Error rates for cyl_bessel_i

GNU C++ version 7.1.0
linux
double

GNU C++ version 7.1.0
linux
long double

Sun compiler version 0x5150
Sun Solaris
long double

Microsoft Visual C++ version 14.1
Win32
double

Bessel I0: Mathworld Data

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 270ε (Mean = 91.6ε) And other failures.)
(Rmath 3.2.3: Max = 1.52ε (Mean = 0.622ε) And other failures.)

Max = 1.95ε (Mean = 0.738ε)

(<cmath>: Max = 8.49ε (Mean = 3.46ε) And other failures.)

Max = 1.95ε (Mean = 0.661ε)

Max = 0.762ε (Mean = 0.329ε)

Bessel I1: Mathworld Data

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 128ε (Mean = 41ε) And other failures.)
(Rmath 3.2.3: Max = 1.53ε (Mean = 0.483ε) And other failures.)

Max = 0.64ε (Mean = 0.202ε)

(<cmath>: Max = 5ε (Mean = 2.15ε) And other failures.)

Max = 0.64ε (Mean = 0.202ε)

Max = 0.767ε (Mean = 0.398ε)

Bessel In: Mathworld Data

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 2.31ε (Mean = 0.838ε) And other failures.)
(Rmath 3.2.3: Max = 1.73ε (Mean = 0.601ε) And other failures.)

Max = 1.8ε (Mean = 1.33ε)

(<cmath>: Max = 430ε (Mean = 163ε) And other failures.)

Max = 463ε (Mean = 140ε)

Max = 3.46ε (Mean = 1.32ε)

Bessel Iv: Mathworld Data

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 5.95ε (Mean = 2.08ε) And other failures.)
(Rmath 3.2.3: Max = 3.53ε (Mean = 1.39ε))

Max = 4.12ε (Mean = 1.85ε)

(<cmath>: Max = 616ε (Mean = 221ε) And other failures.)

Max = 4.12ε (Mean = 1.95ε)

Max = 2.97ε (Mean = 1.24ε)

Bessel In: Random Data

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 261ε (Mean = 53.2ε) And other failures.)
(Rmath 3.2.3: Max = 7.37ε (Mean = 2.4ε))

Max = 4.62ε (Mean = 1.06ε)

(<cmath>: Max = 645ε (Mean = 132ε))

Max = 176ε (Mean = 39.1ε)

Max = 9.67ε (Mean = 1.88ε)

Bessel Iv: Random Data

Max = 0.661ε (Mean = 0.0441ε)

(GSL 2.1: Max = 6.18e+03ε (Mean = 1.55e+03ε) And other failures.)
(Rmath 3.2.3: Max = 4.28e+08ε (Mean = 2.85e+07ε))

Max = 8.35ε (Mean = 1.62ε)

(<cmath>: Max = 1.05e+03ε (Mean = 224ε) And other failures.)

Max = 283ε (Mean = 88.4ε)

Max = 7.46ε (Mean = 1.71ε)

Bessel Iv: Mathworld Data (large values)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 37ε (Mean = 18ε) And other failures.)
(Rmath 3.2.3: Max = 3.77e+168ε (Mean = 2.39e+168ε) And other failures.)

Max = 14.7ε (Mean = 6.66ε)

(<cmath>: Max = 118ε (Mean = 57.2ε) And other failures.)

Max = 14.7ε (Mean = 6.59ε)

Max = 3.67ε (Mean = 1.64ε)


Table 8.46. Error rates for cyl_bessel_k (integer orders)

GNU C++ version 7.1.0
linux
long double

GNU C++ version 7.1.0
linux
double

Sun compiler version 0x5150
Sun Solaris
long double

Microsoft Visual C++ version 14.1
Win32
double

Bessel K0: Mathworld Data (Integer Version)

Max = 0.833ε (Mean = 0.436ε)

(<cmath>: Max = 9.33ε (Mean = 3.25ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 1.2ε (Mean = 0.733ε))
(Rmath 3.2.3: Max = 0.833ε (Mean = 0.601ε))

Max = 0.833ε (Mean = 0.436ε)

Max = 0.833ε (Mean = 0.552ε)

Bessel K1: Mathworld Data (Integer Version)

Max = 0.786ε (Mean = 0.329ε)

(<cmath>: Max = 8.94ε (Mean = 3.19ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 0.626ε (Mean = 0.333ε))
(Rmath 3.2.3: Max = 0.894ε (Mean = 0.516ε))

Max = 0.786ε (Mean = 0.329ε)

Max = 0.786ε (Mean = 0.39ε)

Bessel Kn: Mathworld Data (Integer Version)

Max = 2.6ε (Mean = 1.21ε)

(<cmath>: Max = 12.9ε (Mean = 4.91ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 168ε (Mean = 59.5ε))
(Rmath 3.2.3: Max = 8.48ε (Mean = 2.98ε))

Max = 2.6ε (Mean = 1.21ε)

Max = 3.63ε (Mean = 1.46ε)


Table 8.47. Error rates for cyl_bessel_k

GNU C++ version 7.1.0
linux
long double

GNU C++ version 7.1.0
linux
double

Sun compiler version 0x5150
Sun Solaris
long double

Microsoft Visual C++ version 14.1
Win32
double

Bessel K0: Mathworld Data

Max = 0.833ε (Mean = 0.436ε)

(<cmath>: Max = 9.33ε (Mean = 3.25ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 6.04ε (Mean = 2.16ε))
(Rmath 3.2.3: Max = 0.833ε (Mean = 0.601ε))

Max = 0.833ε (Mean = 0.436ε)

Max = 0.833ε (Mean = 0.552ε)

Bessel K1: Mathworld Data

Max = 0.786ε (Mean = 0.329ε)

(<cmath>: Max = 8.94ε (Mean = 3.19ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 6.26ε (Mean = 2.21ε))
(Rmath 3.2.3: Max = 0.894ε (Mean = 0.516ε))

Max = 0.786ε (Mean = 0.329ε)

Max = 0.786ε (Mean = 0.39ε)

Bessel Kn: Mathworld Data

Max = 2.6ε (Mean = 1.21ε)

(<cmath>: Max = 12.9ε (Mean = 4.91ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 3.36ε (Mean = 1.43ε) And other failures.)
(Rmath 3.2.3: Max = 8.48ε (Mean = 2.98ε))

Max = 2.6ε (Mean = 1.21ε)

Max = 3.63ε (Mean = 1.46ε)

Bessel Kv: Mathworld Data

Max = 3.58ε (Mean = 2.39ε)

(<cmath>: Max = 13ε (Mean = 4.81ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 5.47ε (Mean = 2.04ε) And other failures.)
(Rmath 3.2.3: Max = 3.15ε (Mean = 1.35ε))

Max = 5.21ε (Mean = 2.53ε)

Max = 4.78ε (Mean = 2.19ε)

Bessel Kv: Mathworld Data (large values)

Max = 42.3ε (Mean = 21ε)

(<cmath>: Max = 42.3ε (Mean = 19.8ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 308ε (Mean = 142ε) And other failures.)
(Rmath 3.2.3: Max = 84.6ε (Mean = 37.8ε))

Max = 42.3ε (Mean = 21ε)

Max = 59.8ε (Mean = 26.9ε)

Bessel Kn: Random Data

Max = 4.55ε (Mean = 1.12ε)

(<cmath>: Max = 13.9ε (Mean = 2.91ε))

Max = 0.764ε (Mean = 0.0348ε)

(GSL 2.1: Max = 8.71ε (Mean = 1.76ε) And other failures.)
(Rmath 3.2.3: Max = 7.47ε (Mean = 1.34ε))

Max = 4.55ε (Mean = 1.12ε)

Max = 9.34ε (Mean = 1.7ε)

Bessel Kv: Random Data

Max = 7.88ε (Mean = 1.48ε)

(<cmath>: Max = 13.6ε (Mean = 2.68ε) And other failures.)

Max = 0.507ε (Mean = 0.0313ε)

(GSL 2.1: Max = 9.71ε (Mean = 1.47ε) And other failures.)
(Rmath 3.2.3: Max = 7.37ε (Mean = 1.49ε))

Max = 7.88ε (Mean = 1.47ε)

Max = 8.33ε (Mean = 1.62ε)


The following error plot are based on an exhaustive search of the functions domain for I0, I1, K0, and K1, MSVC-15.5 at double precision, and GCC-7.1/Ubuntu for long double and __float128.

Implementation

The following are handled as special cases first:

When computing Iv for x < 0, then ν must be an integer or a domain error occurs. If ν is an integer, then the function is odd if ν is odd and even if ν is even, and we can reflect to x > 0.

For Iv with v equal to 0, 1 or 0.5 are handled as special cases.

The 0 and 1 cases use polynomial approximations on finite and infinite intervals. The approximating forms are based on "Rational Approximations for the Modified Bessel Function of the First Kind - I0(x) for Computations with Double Precision" by Pavel Holoborodko, extended by us to deal with up to 128-bit precision (with different approximations for each target precision).

Similarly we have:

The 0.5 case is a simple trigonometric function:

I0.5(x) = sqrt(2 / πx) * sinh(x)

For Kv with v an integer, the result is calculated using the recurrence relation:

starting from K0 and K1 which are calculated using rational the approximations above. These rational approximations are accurate to around 19 digits, and are therefore only used when T has no more than 64 binary digits of precision.

When x is small compared to v, Ivx is best computed directly from the series:

In the general case, we first normalize ν to [0, [inf]) with the help of the reflection formulae:

Let μ = ν - floor(ν + 1/2), then μ is the fractional part of ν such that |μ| <= 1/2 (we need this for convergence later). The idea is to calculate Kμ(x) and Kμ+1(x), and use them to obtain Iν(x) and Kν(x).

The algorithm is proposed by Temme in

N.M. Temme, On the numerical evaluation of the modified bessel function of the third kind, Journal of Computational Physics, vol 19, 324 (1975),

which needs two continued fractions as well as the Wronskian:

The continued fractions are computed using the modified Lentz's method

(W.J. Lentz, Generating Bessel functions in Mie scattering calculations using continued fractions, Applied Optics, vol 15, 668 (1976)).

Their convergence rates depend on x, therefore we need different strategies for large x and small x.

x > v, CF1 needs O(x) iterations to converge, CF2 converges rapidly.

x <= v, CF1 converges rapidly, CF2 fails to converge when x -> 0.

When x is large (x > 2), both continued fractions converge (CF1 may be slow for really large x). Kμ and Kμ+1 can be calculated by

where

S is also a series that is summed along with CF2, see

I.J. Thompson and A.R. Barnett, Modified Bessel functions I_v and K_v of real order and complex argument to selected accuracy, Computer Physics Communications, vol 47, 245 (1987).

When x is small (x <= 2), CF2 convergence may fail (but CF1 works very well). The solution here is Temme's series:

where

fk and hk are also computed by recursions (involving gamma functions), but the formulas are a little complicated, readers are referred to

N.M. Temme, On the numerical evaluation of the modified Bessel function of the third kind, Journal of Computational Physics, vol 19, 324 (1975).

Note: Temme's series converge only for |μ| <= 1/2.

Kν(x) is then calculated from the forward recurrence, as is Kν+1(x). With these two values and fν, the Wronskian yields Iν(x) directly.


PrevUpHomeNext