Synopsis

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

Functions for obtaining both a single zero or root of the Bessel function, and placing multiple zeros into a container like std::vector by providing an output iterator.

The signature of the single value functions are:

template <class T>
T cyl_bessel_j_zero(
         T v,            // Floating-point value for Jv.
         int m);         // 1-based index of zero.

template <class T>
T cyl_neumann_zero(
         T v,            // Floating-point value for Jv.
         int m);         // 1-based index of zero.

and for multiple zeros:

template <class T, class OutputIterator>
OutputIterator cyl_bessel_j_zero(
                     T v,                       // Floating-point value for Jv.
                     int start_index,           // 1-based index of first zero.
                     unsigned number_of_zeros,  // How many zeros to generate.
                     OutputIterator out_it);    // Destination for zeros.

template <class T, class OutputIterator>
OutputIterator cyl_neumann_zero(
                     T v,                       // Floating-point value for Jv.
                     int start_index,           // 1-based index of zero.
                     unsigned number_of_zeros,  // How many zeros to generate
                     OutputIterator out_it);    // Destination for zeros.

There are also versions which allow control of the Policies for error handling and precision.

 template <class T>
 T cyl_bessel_j_zero(
          T v,            // Floating-point value for Jv.
          int m,          // 1-based index of zero.
          const Policy&); // Policy to use.

 template <class T>
 T cyl_neumann_zero(
          T v,            // Floating-point value for Jv.
          int m,          // 1-based index of zero.
          const Policy&); // Policy to use.


template <class T, class OutputIterator>
OutputIterator cyl_bessel_j_zero(
                     T v,                       // Floating-point value for Jv.
                     int start_index,           // 1-based index of first zero.
                     unsigned number_of_zeros,  // How many zeros to generate.
                     OutputIterator out_it,     // Destination for zeros.
                     const Policy& pol);        // Policy to use.

template <class T, class OutputIterator>
OutputIterator cyl_neumann_zero(
                     T v,                       // Floating-point value for Jv.
                     int start_index,           // 1-based index of zero.
                     unsigned number_of_zeros,  // How many zeros to generate.
                     OutputIterator out_it,     // Destination for zeros.
                     const Policy& pol);        // Policy to use.

Description

Every real order ν cylindrical Bessel and Neumann functions have an infinite number of zeros on the positive real axis. The real zeros on the positive real axis can be found by solving for the roots of

  J_ν(j_{ν, m}) = 0

  Y_ν(y_{ν, m}) = 0

Here, j_{ν, m} represents the m^th root of the cylindrical Bessel function of order ν, and y_{ν, m} represents the m^th root of the cylindrical Neumann function of order ν.

The zeros or roots (values of x where the function crosses the horizontal y = 0 axis) of the Bessel and Neumann functions are computed by two functions, cyl_bessel_j_zero and cyl_neumann_zero.

In each case the index or rank of the zero returned is 1-based, which is to say:

cyl_bessel_j_zero(v, 1);

returns the first zero of Bessel J.

Passing an start_index <= 0 results in a std::domain_error being raised.

For certain parameters, however, the zero'th root is defined and it has a value of zero. For example, the zero'th root of J[sub v](x) is defined and it has a value of zero for all values of v > 0 and for negative integer values of v = -n. Similar cases are described in the implementation details below.

The order v of J can be positive, negative and zero for the cyl_bessel_j and cyl_neumann functions, but not infinite nor NaN.

Examples of finding Bessel and Neumann zeros

This example demonstrates calculating zeros of the Bessel and Neumann functions. It also shows how Boost.Math and Boost.Multiprecision can be combined to provide a many decimal digit precision. For 50 decimal digit precision we need to include

#include <boost/multiprecision/cpp_dec_float.hpp>

and a typedef for float_type may be convenient (allowing a quick switch to re-compute at built-in double or other precision)

typedef boost::multiprecision::cpp_dec_float_50 float_type;

To use the functions for finding zeros of the functions we need

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

This file includes the forward declaration signatures for the zero-finding functions:

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

but more details are in the full documentation, for example at Boost.Math Bessel functions.

This example shows obtaining both a single zero of the Bessel function, and then placing multiple zeros into a container like std::vector by providing an iterator.

	Tip
	It is always wise to place code using Boost.Math inside try'n'catch blocks; this will ensure that helpful error messages are shown when exceptional conditions arise.

First, evaluate a single Bessel zero.

The precision is controlled by the float-point type of template parameter T of v so this example has double precision, at least 15 but up to 17 decimal digits (for the common 64-bit double).

//    double root = boost::math::cyl_bessel_j_zero(0.0, 1);
//    // Displaying with default precision of 6 decimal digits:
//    std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40483
//    // And with all the guaranteed (15) digits:
//    std::cout.precision(std::numeric_limits<double>::digits10);
//    std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40482555769577

But note that because the parameter v controls the precision of the result, v must be a floating-point type. So if you provide an integer type, say 0, rather than 0.0, then it will fail to compile thus:

root = boost::math::cyl_bessel_j_zero(0, 1);

with this error message

error C2338: Order must be a floating-point type.

Optionally, we can use a policy to ignore errors, C-style, returning some value, perhaps infinity or NaN, or the best that can be done. (See user error handling).

To create a (possibly unwise!) policy ignore_all_policy that ignores all errors:

typedef boost::math::policies::policy<
  boost::math::policies::domain_error<boost::math::policies::ignore_error>,
  boost::math::policies::overflow_error<boost::math::policies::ignore_error>,
  boost::math::policies::underflow_error<boost::math::policies::ignore_error>,
  boost::math::policies::denorm_error<boost::math::policies::ignore_error>,
  boost::math::policies::pole_error<boost::math::policies::ignore_error>,
  boost::math::policies::evaluation_error<boost::math::policies::ignore_error>
            > ignore_all_policy;

Examples of use of this ignore_all_policy are

double inf = std::numeric_limits<double>::infinity();
double nan = std::numeric_limits<double>::quiet_NaN();

double dodgy_root = boost::math::cyl_bessel_j_zero(-1.0, 1, ignore_all_policy());
std::cout << "boost::math::cyl_bessel_j_zero(-1.0, 1) " << dodgy_root << std::endl; // 1.#QNAN
double inf_root = boost::math::cyl_bessel_j_zero(inf, 1, ignore_all_policy());
std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl; // 1.#QNAN
double nan_root = boost::math::cyl_bessel_j_zero(nan, 1, ignore_all_policy());
std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl; // 1.#QNAN

Another version of cyl_bessel_j_zero allows calculation of multiple zeros with one call, placing the results in a container, often std::vector. For example, generate and display the first five double roots of J_v for integral order 2, as column J₂(x) in table 1 of Wolfram Bessel Function Zeros.

unsigned int n_roots = 5U;
std::vector<double> roots;
boost::math::cyl_bessel_j_zero(2.0, 1, n_roots, std::back_inserter(roots));
std::copy(roots.begin(),
          roots.end(),
          std::ostream_iterator<double>(std::cout, "\n"));

Or we can use Boost.Multiprecision to generate 50 decimal digit roots of J_v for non-integral order v= 71/19 == 3.736842, expressed as an exact-integer fraction to generate the most accurate value possible for all floating-point types.

We set the precision of the output stream, and show trailing zeros to display a fixed 50 decimal digits.

std::cout.precision(std::numeric_limits<float_type>::digits10); // 50 decimal digits.
std::cout << std::showpoint << std::endl; // Show trailing zeros.

float_type x = float_type(71) / 19;
float_type r = boost::math::cyl_bessel_j_zero(x, 1); // 1st root.
std::cout << "x = " << x << ", r = " << r << std::endl;

r = boost::math::cyl_bessel_j_zero(x, 20U); // 20th root.
std::cout << "x = " << x << ", r = " << r << std::endl;

std::vector<float_type> zeros;
boost::math::cyl_bessel_j_zero(x, 1, 3, std::back_inserter(zeros));

std::cout << "cyl_bessel_j_zeros" << std::endl;
// Print the roots to the output stream.
std::copy(zeros.begin(), zeros.end(),
          std::ostream_iterator<float_type>(std::cout, "\n"));

Using Output Iterator to sum zeros of Bessel Functions

This example demonstrates summing zeros of the Bessel functions. To use the functions for finding zeros of the functions we need

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

We use the cyl_bessel_j_zero output iterator parameter out_it to create a sum of 1/zeros² by defining a custom output iterator:

template <class T>
struct output_summation_iterator
{
   output_summation_iterator(T* p) : p_sum(p)
   {}
   output_summation_iterator& operator*()
   { return *this; }
    output_summation_iterator& operator++()
   { return *this; }
   output_summation_iterator& operator++(int)
   { return *this; }
   output_summation_iterator& operator = (T const& val)
   {
     *p_sum += 1./ (val * val); // Summing 1/zero^2.
     return *this;
   }
private:
   T* p_sum;
};

The sum is calculated for many values, converging on the analytical exact value of 1/8.

using boost::math::cyl_bessel_j_zero;
double nu = 1.;
double sum = 0;
output_summation_iterator<double> it(&sum);  // sum of 1/zeros^2
cyl_bessel_j_zero(nu, 1, 10000, it);

double s = 1/(4 * (nu + 1)); // 0.125 = 1/8 is exact analytical solution.
std::cout << std::setprecision(6) << "nu = " << nu << ", sum = " << sum
  << ", exact = " << s << std::endl;
// nu = 1.00000, sum = 0.124990, exact = 0.125000

Calculating zeros of the Neumann function.

This example also shows how Boost.Math and Boost.Multiprecision can be combined to provide a many decimal digit precision. For 50 decimal digit precision we need to include

#include <boost/multiprecision/cpp_dec_float.hpp>

and a typedef for float_type may be convenient (allowing a quick switch to re-compute at built-in double or other precision)

typedef boost::multiprecision::cpp_dec_float_50 float_type;

To use the functions for finding zeros of the cyl_neumann function we need:

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

The Neumann (Bessel Y) function zeros are evaluated very similarly:

using boost::math::cyl_neumann_zero;
double zn = cyl_neumann_zero(2., 1);
std::cout << "cyl_neumann_zero(2., 1) = " << zn << std::endl;

std::vector<float> nzeros(3); // Space for 3 zeros.
cyl_neumann_zero<float>(2.F, 1, nzeros.size(), nzeros.begin());

std::cout << "cyl_neumann_zero<float>(2.F, 1, ";
// Print the zeros to the output stream.
std::copy(nzeros.begin(), nzeros.end(),
          std::ostream_iterator<float>(std::cout, ", "));

std::cout << "\n""cyl_neumann_zero(static_cast<float_type>(220)/100, 1) = "
  << cyl_neumann_zero(static_cast<float_type>(220)/100, 1) << std::endl;
// 3.6154383428745996706772556069431792744372398748422

Error messages from 'bad' input

Another example demonstrates calculating zeros of the Bessel functions showing the error messages from 'bad' input is handled by throwing exceptions.

To use the functions for finding zeros of the functions we need:

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

	Tip
	It is always wise to place all code using Boost.Math inside try'n'catch blocks; this will ensure that helpful error messages can be shown when exceptional conditions arise.

Examples below show messages from several 'bad' arguments that throw a domain_error exception.

try
{ // Try a zero order v.
  float dodgy_root = boost::math::cyl_bessel_j_zero(0.F, 0);
  std::cout << "boost::math::cyl_bessel_j_zero(0.F, 0) " << dodgy_root << std::endl;
  // Thrown exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int):
  // Requested the 0'th zero of J0, but the rank must be > 0 !
}
catch (std::exception& ex)
{
  std::cout << "Thrown exception " << ex.what() << std::endl;
}

	Note
	The type shown in the error message is the type after promotion, using precision policy and internal promotion policy, from float to double in this case.

In this example the promotion goes:

See full code for other examples that promote from double to long double.

Other examples of 'bad' inputs like infinity and NaN are below. Some compiler warnings indicate that 'bad' values are detected at compile time.

try
{ // order v = inf
   std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << std::endl;
   double inf = std::numeric_limits<double>::infinity();
   double inf_root = boost::math::cyl_bessel_j_zero(inf, 1);
   std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl;
   // Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):
   // Order argument is 1.#INF, but must be finite >= 0 !
}
catch (std::exception& ex)
{
  std::cout << "Thrown exception " << ex.what() << std::endl;
}

try
{ // order v = NaN, rank m = 1
   std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << std::endl;
   double nan = std::numeric_limits<double>::quiet_NaN();
   double nan_root = boost::math::cyl_bessel_j_zero(nan, 1);
   std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl;
   // Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):
   // Order argument is 1.#QNAN, but must be finite >= 0 !
}
catch (std::exception& ex)
{
  std::cout << "Thrown exception " << ex.what() << std::endl;
}

The output from other examples are shown appended to the full code listing.

The full code (and output) for these examples is at Bessel zeros, Bessel zeros iterator, Neumann zeros, Bessel error messages.

Implementation

Various methods are used to compute initial estimates for j_{ν, m} and y_{ν, m} ; these are described in detail below.

After finding the initial estimate of a given root, its precision is subsequently refined to the desired level using Newton-Raphson iteration from Boost.Math's root-finding with derivatives utilities combined with the functions cyl_bessel_j and cyl_neumann.

Newton iteration requires both J_ν(x) or Y_ν(x) as well as its derivative. The derivatives of J_ν(x) and Y_ν(x) with respect to x are given by M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, NBS (1964). In particular,

  d/_dx J_ν(x) = J_ν-1(x) - ν J_ν(x) / x

  d/_dx Y_ν(x) = Y_ν-1(x) - ν Y_ν(x) / x

Enumeration of the rank of a root (in other words the index of a root) begins with one and counts up, in other words m,=1,2,3,… The value of the first root is always greater than zero.

For certain special parameters, cylindrical Bessel functions and cylindrical Neumann functions have a root at the origin. For example, J_ν(x) has a root at the origin for every positive order ν > 0, and for every negative integer order ν = -n with n ∈ ℕ ⁺ and n ≠ 0.

In addition, Y_ν(x) has a root at the origin for every negative half-integer order ν = -n/2, with n ∈ ℕ ⁺ and and n ≠ 0.

For these special parameter values, the origin with a value of x = 0 is provided as the 0^th root generated by cyl_bessel_j_zero() and cyl_neumann_zero().

When calculating initial estimates for the roots of Bessel functions, a distinction is made between positive order and negative order, and different methods are used for these. In addition, different algorithms are used for the first root m = 1 and for subsequent roots with higher rank m ≥ 2. Furthermore, estimates of the roots for Bessel functions with order above and below a cutoff at ν = 2.2 are calculated with different methods.

Calculations of the estimates of j_ν,1 and y_ν,1 with 0 ≤ ν < 2.2 use empirically tabulated values. The coefficients for these have been generated by a computer algebra system.

Calculations of the estimates of j_ν,1 and y_ν,1 with ν≥ 2.2 use Eqs.9.5.14 and 9.5.15 in M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, NBS (1964).

In particular,

  j_ν,1 ≅ ν + 1.85575 ν^⅓ + 1.033150 ν^-⅓ - 0.00397 ν^-1 - 0.0908 ν^-5/3 + 0.043 ν^-7/3 + …

and

  y_ν,1 ≅ ν + 0.93157 ν^⅓ + 0.26035 ν^-⅓ + 0.01198 ν^-1 - 0.0060 ν^-5/3 - 0.001 ν^-7/3 + …

Calculations of the estimates of j_{ν, m} and y_{ν, m} with rank m > 2 and 0 ≤ ν < 2.2 use McMahon's approximation, as described in M. Abramowitz and I. A. Stegan, Section 9.5 and 9.5.12. In particular,

  j_ν,m, y_ν,m ≅ β - (μ-1) / 8β

      - 4(μ-1)(7μ - 31) / 3(8β)³

      -32(μ-1)(83μ² - 982μ + 3779) / 15(8β)⁵

      -64(μ-1)(6949μ³ - 153855μ² + 1585743μ- 6277237) / 105(8a)⁷

      - …     (5)

where μ = 4ν² and β = (m + ½ν - ¼)π for j_ν,m and β = (m + ½ν -¾)π for y_ν,m.

Calculations of the estimates of j_{ν, m} and y_{ν, m} with ν ≥ 2.2 use one term in the asymptotic expansion given in Eq.9.5.22 and top line of Eq.9.5.26 combined with Eq. 9.3.39, all in M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, NBS (1964) explicit and easy-to-understand treatment for asymptotic expansion of zeros. The latter two equations are expressed for argument (x) greater than one. (Olver also gives the series form of the equations in §10.21(vi) McMahon's Asymptotic Expansions for Large Zeros - using slightly different variable names).

In summary,

  j_{ν, m} ∼ νx(-ζ) + f₁(-ζ/ν)

where -ζ = ν^-2/3a_m and a_m is the absolute value of the m^th root of Ai(x) on the negative real axis.

Here x = x(-ζ) is the inverse of the function

  ⅔(-ζ)^3/2 = √(x² - 1) - cos⁻¹(1/x)     (7)

Furthermore,

  f₁(-ζ) = ½x(-ζ) {h(-ζ)}² ⋅ b₀(-ζ)

where

  h(-ζ) = {4(-ζ) / (x² - 1)}⁴

and

  b₀(-ζ) = -5/(48ζ²) + 1/(-ζ)^½ ⋅ { 5/(24(x²-1)^3/2) + 1/(8(x²-1)^½)}

When solving for x(-ζ) in Eq. 7 above, the right-hand-side is expanded to order 2 in a Taylor series for large x. This results in

  ⅔(-ζ)^3/2 ≈ x + 1/2x - π/2

The positive root of the resulting quadratic equation is used to find an initial estimate x(-ζ). This initial estimate is subsequently refined with several steps of Newton-Raphson iteration in Eq. 7.

Estimates of the roots of cylindrical Bessel functions of negative order on the positive real axis are found using interlacing relations. For example, the m^th root of the cylindrical Bessel function j_-ν,m is bracketed by the m^th root and the (m+1)^th root of the Bessel function of corresponding positive integer order. In other words,

  j_nν,m < j_-ν,m < j_nν,m+1

where m > 1 and n_ν represents the integral floor of the absolute value of |-ν|.

Similar bracketing relations are used to find estimates of the roots of Neumann functions of negative order, whereby a discontinuity at every negative half-integer order needs to be handled.

Bracketing relations do not hold for the first root of cylindrical Bessel functions and cylindrical Neumann functions with negative order. Therefore, iterative algorithms combined with root-finding via bisection are used to localize j_-ν,1 and y_-ν,1.

Testing

The precision of evaluation of zeros was tested at 50 decimal digits using cpp_dec_float_50 and found identical with spot values computed by Wolfram Alpha.