Boost C++ Libraries

...one of the most highly regarded and expertly designed C++ library projects in the world. Herb Sutter and Andrei Alexandrescu, C++ Coding Standards

PrevUpHomeNext

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_j(T1 v, T2 x);

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

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

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

The functions cyl_bessel_j and cyl_neumann return the result of the Bessel functions of the first and second kinds respectively:

cyl_bessel_j(v, x) = Jv(x)

cyl_neumann(v, x) = Yv(x) = Nv(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 cyclic nature of Jv:

The following graph shows the behaviour of Yv: this is also cyclic for large x, but tends to -∞ for small x:

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 comparisons to other libraries. Note that the cyclic nature of these functions means that they have an infinite number of irrational roots: in general these functions have arbitrarily large relative errors when the arguments are sufficiently close to a root. Of course the absolute error in such cases is always small. 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. Most of the gross errors exhibited by other libraries occur for very large arguments - you will need to drill down into the actual program output if you need more information on this.

Table 8.40. Error rates for cyl_bessel_j (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 J0: Mathworld Data (Integer Version)

Max = 6.55ε (Mean = 2.86ε)

(<cmath>: Max = 5.04ε (Mean = 1.78ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 1.12ε (Mean = 0.488ε))
(Rmath 3.2.3: Max = 0.629ε (Mean = 0.223ε) And other failures.)

Max = 6.55ε (Mean = 2.86ε)

Max = 2.52ε (Mean = 1.2ε)

(<math.h>: Max = 1.89ε (Mean = 0.988ε))

Bessel J0: Mathworld Data (Tricky cases) (Integer Version)

Max = 1.64e+08ε (Mean = 6.69e+07ε)

(<cmath>: Max = 4.79e+08ε (Mean = 1.96e+08ε))

Max = 8e+04ε (Mean = 3.27e+04ε)

(GSL 2.1: Max = 1e+07ε (Mean = 4.11e+06ε))
(Rmath 3.2.3: Max = 1.04e+07ε (Mean = 4.29e+06ε))

Max = 1.64e+08ε (Mean = 6.69e+07ε)

Max = 1e+07ε (Mean = 4.09e+06ε)

(<math.h>: Max = 2.54e+08ε (Mean = 1.04e+08ε))

Bessel J1: Mathworld Data (Integer Version)

Max = 3.59ε (Mean = 1.33ε)

(<cmath>: Max = 6.1ε (Mean = 2.95ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 1.89ε (Mean = 0.721ε))
(Rmath 3.2.3: Max = 0.946ε (Mean = 0.39ε) And other failures.)

Max = 1.44ε (Mean = 0.637ε)

Max = 1.73ε (Mean = 0.976ε)

(<math.h>: Max = 11.4ε (Mean = 4.15ε))

Bessel J1: Mathworld Data (tricky cases) (Integer Version)

Max = 2.18e+05ε (Mean = 9.76e+04ε)

(<cmath>: Max = 2.15e+06ε (Mean = 1.58e+06ε))

Max = 106ε (Mean = 47.5ε)

(GSL 2.1: Max = 1.26e+06ε (Mean = 6.28e+05ε))
(Rmath 3.2.3: Max = 2.93e+06ε (Mean = 1.7e+06ε))

Max = 2.18e+05ε (Mean = 9.76e+04ε)

Max = 3.23e+04ε (Mean = 1.45e+04ε)

(<math.h>: Max = 1.44e+07ε (Mean = 6.5e+06ε))

Bessel JN: Mathworld Data (Integer Version)

Max = 6.85ε (Mean = 3.35ε)

(<cmath>: Max = 2.13e+19ε (Mean = 5.16e+18ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 6.9e+05ε (Mean = 2.53e+05ε) And other failures.)
(Rmath 3.2.3: Max = +INFε (Mean = +INFε) And other failures.)

Max = 463ε (Mean = 112ε)

Max = 14.7ε (Mean = 5.4ε)

(<math.h>: Max = +INFε (Mean = +INFε) And other failures.)


Table 8.41. Error rates for cyl_bessel_j

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 J0: Mathworld Data

Max = 6.55ε (Mean = 2.86ε)

(<cmath>: Max = 5.04ε (Mean = 1.78ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 0.629ε (Mean = 0.223ε) And other failures.)
(Rmath 3.2.3: Max = 0.629ε (Mean = 0.223ε) And other failures.)

Max = 6.55ε (Mean = 2.86ε)

Max = 2.52ε (Mean = 1.2ε)

Bessel J0: Mathworld Data (Tricky cases)

Max = 1.64e+08ε (Mean = 6.69e+07ε)

(<cmath>: Max = 4.79e+08ε (Mean = 1.96e+08ε))

Max = 8e+04ε (Mean = 3.27e+04ε)

(GSL 2.1: Max = 6.5e+07ε (Mean = 2.66e+07ε))
(Rmath 3.2.3: Max = 1.04e+07ε (Mean = 4.29e+06ε))

Max = 1.64e+08ε (Mean = 6.69e+07ε)

Max = 1e+07ε (Mean = 4.09e+06ε)

Bessel J1: Mathworld Data

Max = 3.59ε (Mean = 1.33ε)

(<cmath>: Max = 6.1ε (Mean = 2.95ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 6.62ε (Mean = 2.35ε) And other failures.)
(Rmath 3.2.3: Max = 0.946ε (Mean = 0.39ε) And other failures.)

Max = 1.44ε (Mean = 0.637ε)

Max = 1.73ε (Mean = 0.976ε)

Bessel J1: Mathworld Data (tricky cases)

Max = 2.18e+05ε (Mean = 9.76e+04ε)

(<cmath>: Max = 2.15e+06ε (Mean = 1.58e+06ε))

Max = 106ε (Mean = 47.5ε)

(GSL 2.1: Max = 8.75e+05ε (Mean = 5.32e+05ε))
(Rmath 3.2.3: Max = 2.93e+06ε (Mean = 1.7e+06ε))

Max = 2.18e+05ε (Mean = 9.76e+04ε)

Max = 3.23e+04ε (Mean = 1.45e+04ε)

Bessel JN: Mathworld Data

Max = 6.85ε (Mean = 3.35ε)

(<cmath>: Max = 2.13e+19ε (Mean = 5.16e+18ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 6.9e+05ε (Mean = 2.15e+05ε) And other failures.)
(Rmath 3.2.3: Max = +INFε (Mean = +INFε) And other failures.)

Max = 463ε (Mean = 112ε)

Max = 14.7ε (Mean = 5.4ε)

Bessel J: Mathworld Data

Max = 14.7ε (Mean = 4.11ε)

(<cmath>: Max = 3.49e+05ε (Mean = 8.09e+04ε) And other failures.)

Max = 10ε (Mean = 2.24ε)

(GSL 2.1: Max = 2.39e+05ε (Mean = 5.37e+04ε) And other failures.)
(Rmath 3.2.3: Max = +INFε (Mean = +INFε) And other failures.)

Max = 14.7ε (Mean = 4.22ε)

Max = 14.9ε (Mean = 3.89ε)

Bessel J: Mathworld Data (large values)

Max = 607ε (Mean = 305ε)

(<cmath>: Max = 34.9ε (Mean = 17.4ε) And other failures.)

Max = 0.536ε (Mean = 0.268ε)

(GSL 2.1: Max = 4.91e+03ε (Mean = 2.46e+03ε) And other failures.)
(Rmath 3.2.3: Max = 5.9ε (Mean = 3.76ε))

Max = 607ε (Mean = 305ε)

Max = 9.31ε (Mean = 5.52ε)

Bessel JN: Random Data

Max = 50.8ε (Mean = 3.69ε)

(<cmath>: Max = 1.12e+03ε (Mean = 88.7ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 75.7ε (Mean = 5.36ε))
(Rmath 3.2.3: Max = 3.93ε (Mean = 1.22ε))

Max = 99.6ε (Mean = 22ε)

Max = 17.5ε (Mean = 1.46ε)

Bessel J: Random Data

Max = 11.4ε (Mean = 1.68ε)

(<cmath>: Max = 501ε (Mean = 52.3ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 15.5ε (Mean = 3.33ε) And other failures.)
(Rmath 3.2.3: Max = 6.74ε (Mean = 1.3ε))

Max = 260ε (Mean = 34ε)

Max = 9.24ε (Mean = 1.17ε)

Bessel J: Random Data (Tricky large values)

Max = 785ε (Mean = 94.2ε)

(<cmath>: Max = 5.01e+17ε (Mean = 6.23e+16ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 2.48e+05ε (Mean = 5.11e+04ε))
(Rmath 3.2.3: Max = 71.6ε (Mean = 11.7ε))

Max = 785ε (Mean = 97.4ε)

Max = 59.2ε (Mean = 8.67ε)


Table 8.42. Error rates for cyl_neumann (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

Y0: Mathworld Data (Integer Version)

Max = 5.53ε (Mean = 2.4ε)

(<cmath>: Max = 2.05e+05ε (Mean = 6.87e+04ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 6.46ε (Mean = 2.38ε))
(Rmath 3.2.3: Max = 167ε (Mean = 56.5ε))

Max = 5.53ε (Mean = 2.4ε)

Max = 4.61ε (Mean = 2.29ε)

(<math.h>: Max = 5.37e+03ε (Mean = 1.81e+03ε))

Y1: Mathworld Data (Integer Version)

Max = 6.33ε (Mean = 2.25ε)

(<cmath>: Max = 9.71e+03ε (Mean = 4.08e+03ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 1.51ε (Mean = 0.839ε))
(Rmath 3.2.3: Max = 193ε (Mean = 64.4ε))

Max = 6.33ε (Mean = 2.29ε)

Max = 4.75ε (Mean = 1.72ε)

(<math.h>: Max = 1.86e+04ε (Mean = 6.2e+03ε))

Yn: Mathworld Data (Integer Version)

Max = 55.2ε (Mean = 17.8ε)

(<cmath>: Max = 2.2e+20ε (Mean = 6.97e+19ε) And other failures.)

Max = 0.993ε (Mean = 0.314ε)

(GSL 2.1: Max = 2.41e+05ε (Mean = 7.62e+04ε))
(Rmath 3.2.3: Max = 1.24e+04ε (Mean = 4e+03ε))

Max = 55.2ε (Mean = 17.8ε)

Max = 35ε (Mean = 11.9ε)

(<math.h>: Max = 2.49e+05ε (Mean = 8.14e+04ε))


Table 8.43. Error rates for cyl_neumann

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

Y0: Mathworld Data

Max = 5.53ε (Mean = 2.4ε)

(<cmath>: Max = 2.05e+05ε (Mean = 6.87e+04ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 60.9ε (Mean = 20.4ε))
(Rmath 3.2.3: Max = 167ε (Mean = 56.5ε))

Max = 5.53ε (Mean = 2.4ε)

Max = 4.61ε (Mean = 2.29ε)

Y1: Mathworld Data

Max = 6.33ε (Mean = 2.25ε)

(<cmath>: Max = 9.71e+03ε (Mean = 4.08e+03ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 23.4ε (Mean = 8.1ε))
(Rmath 3.2.3: Max = 193ε (Mean = 64.4ε))

Max = 6.33ε (Mean = 2.29ε)

Max = 4.75ε (Mean = 1.72ε)

Yn: Mathworld Data

Max = 55.2ε (Mean = 17.8ε)

(<cmath>: Max = 2.2e+20ε (Mean = 6.97e+19ε) And other failures.)

Max = 0.993ε (Mean = 0.314ε)

(GSL 2.1: Max = 2.41e+05ε (Mean = 7.62e+04ε) And other failures.)
(Rmath 3.2.3: Max = 1.24e+04ε (Mean = 4e+03ε))

Max = 55.2ε (Mean = 17.8ε)

Max = 35ε (Mean = 11.9ε)

Yv: Mathworld Data

Max = 10.7ε (Mean = 4.93ε)

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

Max = 10ε (Mean = 3.02ε)

(GSL 2.1: Max = 1.07e+05ε (Mean = 3.22e+04ε) And other failures.)
(Rmath 3.2.3: Max = 243ε (Mean = 73.9ε))

Max = 10.7ε (Mean = 5.1ε)

Max = 7.89ε (Mean = 3.27ε)

Yv: Mathworld Data (large values)

Max = 1.7ε (Mean = 1.33ε)

(<cmath>: Max = 43.2ε (Mean = 16.3ε) And other failures.)

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 60.8ε (Mean = 23ε) And other failures.)
(Rmath 3.2.3: Max = 0.682ε (Mean = 0.335ε))

Max = 1.7ε (Mean = 1.33ε)

Max = 0.682ε (Mean = 0.423ε)

Y0 and Y1: Random Data

Max = 10.8ε (Mean = 3.04ε)

(<cmath>: Max = 2.59e+03ε (Mean = 500ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 34.4ε (Mean = 8.9ε))
(Rmath 3.2.3: Max = 83ε (Mean = 14.2ε))

Max = 10.8ε (Mean = 3.04ε)

Max = 4.17ε (Mean = 1.24ε)

Yn: Random Data

Max = 338ε (Mean = 27.5ε)

(<cmath>: Max = 4.01e+03ε (Mean = 348ε))

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 500ε (Mean = 47.8ε))
(Rmath 3.2.3: Max = 691ε (Mean = 67.9ε))

Max = 338ε (Mean = 27.5ε)

Max = 117ε (Mean = 10.2ε)

Yv: Random Data

Max = 2.08e+03ε (Mean = 149ε)

(<cmath>: Max = +INFε (Mean = +INFε) And other failures.)

Max = 1.53ε (Mean = 0.102ε)

(GSL 2.1: Max = 1.41e+06ε (Mean = 7.67e+04ε))
(Rmath 3.2.3: Max = 1.79e+05ε (Mean = 9.64e+03ε))

Max = 2.08e+03ε (Mean = 149ε)

Max = 1.23e+03ε (Mean = 69.9ε)


Note that for large x these functions are largely dependent on the accuracy of the std::sin and std::cos functions.

Comparison to GSL and Cephes is interesting: both Cephes and this library optimise the integer order case - leading to identical results - simply using the general case is for the most part slightly more accurate though, as noted by the better accuracy of GSL in the integer argument cases. This implementation tends to perform much better when the arguments become large, Cephes in particular produces some remarkably inaccurate results with some of the test data (no significant figures correct), and even GSL performs badly with some inputs to Jv. Note that by way of double-checking these results, the worst performing Cephes and GSL cases were recomputed using functions.wolfram.com, and the result checked against our test data: no errors in the test data were found.

The following error plot are based on an exhaustive search of the functions domain for J0 and Y0, MSVC-15.5 at double precision, other compilers and precisions are very similar - the plots simply illustrate the relatively large errors as you approach a zero, and the very low errors elsewhere.

Implementation

The implementation is mostly about filtering off various special cases:

When x is negative, then the order v must be an integer or the result is a domain error. If the order is an integer then the function is odd for odd orders and even for even orders, so we reflect to x > 0.

When the order v is negative then the reflection formulae can be used to move to v > 0:

Note that if the order is an integer, then these formulae reduce to:

J-n = (-1)nJn

Y-n = (-1)nYn

However, in general, a negative order implies that we will need to compute both J and Y.

When x is large compared to the order v then the asymptotic expansions for large x in M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions 9.2.19 are used (these were found to be more reliable than those in A&S 9.2.5).

When the order v is an integer the method first relates the result to J0, J1, Y0 and Y1 using either forwards or backwards recurrence (Miller's algorithm) depending upon which is stable. The values for J0, J1, Y0 and Y1 are calculated using the rational minimax approximations on root-bracketing intervals for small |x| and Hankel asymptotic expansion for large |x|. The coefficients are from:

W.J. Cody, ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of Special Function Routines and Test Drivers, ACM Transactions on Mathematical Software, vol 19, 22 (1993).

and

J.F. Hart et al, Computer Approximations, John Wiley & Sons, New York, 1968.

These approximations are accurate to around 19 decimal digits: therefore these methods are not used when type T has more than 64 binary digits.

When x is smaller than machine epsilon then the following approximations for Y0(x), Y1(x), Y2(x) and Yn(x) can be used (see: http://functions.wolfram.com/03.03.06.0037.01, http://functions.wolfram.com/03.03.06.0038.01, http://functions.wolfram.com/03.03.06.0039.01 and http://functions.wolfram.com/03.03.06.0040.01):

When x is small compared to v and v is not an integer, then the following series approximation can be used for Yv(x), this is also an area where other approximations are often too slow to converge to be used (see http://functions.wolfram.com/03.03.06.0034.01):

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

In the general case we compute Jv and Yv simultaneously.

To get the initial values, 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 Jμ(x), Jμ+1(x), Yμ(x), Yμ+1(x) and use them to obtain Jν(x), Yν(x).

The algorithm is called Steed's method, which needs two continued fractions as well as the Wronskian:

See: F.S. Acton, Numerical Methods that Work, The Mathematical Association of America, Washington, 1997.

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). Jμ, Jμ+1, Yμ, Yμ+1 can be calculated by

where

Jν and Yμ are then calculated using backward (Miller's algorithm) and forward recurrence respectively.

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

where

gk 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 ordinary Bessel function of the second kind, Journal of Computational Physics, vol 21, 343 (1976). Note Temme's series converge only for |μ| <= 1/2.

As the previous case, Yν is calculated from the forward recurrence, so is Yν+1. With these two values and fν, the Wronskian yields Jν(x) directly without backward recurrence.


PrevUpHomeNext