...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, // 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.
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 mth root of the cylindrical Bessel function of order ν, and yν, m represents the mth 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.
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 Jv for integral order 2, as column J2(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
Jv 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"));
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/zeros2 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 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
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 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 0th 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 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(-ζ) + f1(-ζ/ν)
where -ζ = ν-2/3am and am is the absolute value of the mth 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,
f1(-ζ) = ½x(-ζ) {h(-ζ)}² ⋅ b0(-ζ)
where
h(-ζ) = {4(-ζ) / (x² - 1)}4
and
b0(-ζ) = -5/(48ζ²) + 1/(-ζ)½ ⋅ { 5/(24(x2-1)3/2) + 1/(8(x2-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 mth root of the cylindrical Bessel function j-ν,m is bracketed by the mth root and the (m+1)th root of the Bessel function of corresponding positive integer order. In other words,
jnν,m < j-ν,m < jnν,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.
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.