...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#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, // Floatingpoint value for Jv. int m); // 1based index of zero. template <class T> T cyl_neumann_zero( T v, // Floatingpoint value for Jv. int m); // 1based index of zero.
and for multiple zeros:
template <class T, class OutputIterator> OutputIterator cyl_bessel_j_zero( T v, // Floatingpoint value for Jv. int start_index, // 1based 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, // Floatingpoint value for Jv. int start_index, // 1based 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, // Floatingpoint value for Jv. int m, // 1based index of zero. const Policy&); // Policy to use. template <class T> T cyl_neumann_zero( T v, // Floatingpoint value for Jv. int m, // 1based index of zero. const Policy&); // Policy to use. template <class T, class OutputIterator> OutputIterator cyl_bessel_j_zero( T v, // Floatingpoint value for Jv. int start_index, // 1based 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, // Floatingpoint value for Jv. int start_index, // 1based index of zero. unsigned number_of_zeros, // How many zeros to generate. OutputIterator out_it, // Destination for zeros. const Policy& pol); // Policy to use.
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 1based, 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.
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 recompute at builtin 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 zerofinding 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 floatpoint 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 64bit 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 floatingpoint 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 floatingpoint type.
Optionally, we can use a policy to ignore errors, Cstyle, 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_{2}(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 nonintegral order v= 71/19 == 3.736842
,
expressed as an exactinteger fraction to generate the most accurate value
possible for all floatingpoint 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"));
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^{2} 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
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 recompute at builtin 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
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 
In this example the promotion goes:
float
and
int
.
int
"as if"
it were a double
, so arguments
are float
and double
.
double

so that's the precision we want (and the type that will be returned).
double
for full float
precision.
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.
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 NewtonRaphson iteration from Boost.Math's rootfinding 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 halfinteger 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β)^{3}
32(μ1)(83μ²  982μ + 3779) / 15(8β)^{5}
64(μ1)(6949μ^{3}  153855μ² + 1585743μ 6277237) / 105(8a)^{7}
 … (5)
where μ = 4ν^{2} 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 easytounderstand 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_{1}(ζ/ν)
where ζ = ν^{2/3}a_{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_{1}(ζ) = ½x(ζ) {h(ζ)}² ⋅ b_{0}(ζ)
where
h(ζ) = {4(ζ) / (x²  1)}^{4}
and
b_{0}(ζ) = 5/(48ζ²) + 1/(ζ)^{½} ⋅ { 5/(24(x^{2}1)^{3/2}) + 1/(8(x^{2}1)^{½)}}
When solving for x(ζ) in Eq. 7 above, the righthandside 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 NewtonRaphson 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 halfinteger 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 rootfinding via bisection are used to localize j_{ν,1} and y_{ν,1}.
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.