boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-28 19:32:08 +00:00
Code Issues Projects Releases Wiki Activity

Corrected cyl_neumann_zero() for negative order.

Browse Source 
Improved the clarity of source level comments in the Bessel zero codes.

[SVN r83109]
This commit is contained in:
Christopher Kormanyos
2013-02-23 20:54:18 +00:00
parent 712cbabdb1
commit f2be55beef
2 changed files with 51 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
include/boost/math/special_functions/bessel.hpp
View File

@@ -363,12 +363,13 @@ inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)
const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
// Handle the non-finite order.
// Handle non-finite order.
if (!(boost::math::isfinite)(v) )
{
return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
}
// Handle negative rank.
if(m < 0)
{
// Zeros of Jv(x) with negative rank are not defined and requesting one raises a domain error.
@@ -415,12 +416,14 @@ inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)
// Account for the radix of number representations having non-two radix!
const int my_digits2 = policies::digits<T, Policy>();
const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));
// Perform the root-finding using Newton-Raphson iteration from Boost.Math.
const T jvm =
boost::math::tools::newton_raphson_iterate(
boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv_and_jv_prime<T, Policy>((order_is_integer ? vv : v), order_is_zero, pol),
guess_root,
(std::max)(T(guess_root - 0.2F), T(0)),
T(guess_root - delta_lo),
T(guess_root + 0.2F),
my_digits2,
number_of_iterations);
@@ -447,7 +450,7 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
}
// Handle if the zero'th root or a negative root is requested.
// Handle negative rank.
if(m < 0)
{
return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", m, pol);
@@ -459,10 +462,19 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
const bool order_is_negative = (v < 0);
const T vv((!order_is_negative) ? v : -v);
// Check if the order is very close to a negative half-integer.
const bool order_is_negative_half_integer = (order_is_negative && (((vv * 2U) - floor(vv * 2U)) < boost::math::tools::epsilon<T>()));
const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);
// The zero'th zero of Jv(x) for v < 0 is not defined
// For negative integers, use reflection to positive integer order.
if(order_is_negative && order_is_integer)
return boost::math::detail::cyl_neumann_zero_imp(vv, m, pol);
// Check if the order is very close to a negative half-integer.
const T delta_half_integer(vv - (floor(vv) + 0.5F));
const bool order_is_negative_half_integer =
(order_is_negative && ((delta_half_integer > -half_epsilon) && (delta_half_integer < +half_epsilon)));
// The zero'th zero of Yv(x) for v < 0 is not defined
// unless the order is a negative integer.
if((m == 0) && (!order_is_negative_half_integer))
{
@@ -470,7 +482,7 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Yv for negative, non-half-integer order, but the rank must be > 0 !", m, pol);
}
// For nrgative half-integers, use the corresponding
// For negative half-integers, use the corresponding
// spherical Bessel function of positive half-integer order.
if(order_is_negative_half_integer)
return boost::math::detail::cyl_bessel_j_zero_imp(vv, m, pol);
@@ -487,12 +499,14 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
// Account for the radix of number representations having non-two radix!
const int my_digits2 = policies::digits<T, Policy>();
const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));
// Perform the root-finding using Newton-Raphson iteration from Boost.Math.
const T yvm =
boost::math::tools::newton_raphson_iterate(
boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv_and_yv_prime<T, Policy>(v, pol),
guess_root,
(std::max)(T(guess_root - 0.2F), T(0)),
T(guess_root - delta_lo),
T(guess_root + 0.2F),
my_digits2,
number_of_iterations);

35
include/boost/math/special_functions/detail/bessel_jy_zero.hpp
View File

@@ -9,9 +9,11 @@
//
// This header contains implementation details for estimating the zeros
// of cylindrical Bessel and Neumann functions on the positive real axis.
// Support is included for both positive as well as negative order.
// Various methods are used to estimate the roots. These include
// empirical curve fitting and McMahon's asymptotic approximation
// for small order, and uniform asymptotic expansion for large order.
// for small order, uniform asymptotic expansion for large order,
// and iteration and root interlacing for negative order.
//
#ifndef _BESSEL_JY_ZERO_2013_01_18_HPP_
#define _BESSEL_JY_ZERO_2013_01_18_HPP_
@@ -283,7 +285,14 @@
root_hi = root_lo;
// Decrease the lower end of the bracket using an adaptive algorithm.
(root_lo > 0.5F) ? root_lo -= 0.5F : root_lo *= 0.85F;
if(root_lo > 0.5F)
{
root_lo -= 0.5F;
}
else
{
root_lo *= 0.85F;
}
}
}
else
@@ -437,7 +446,7 @@
template<class T, class Policy>
T initial_guess(const T& v, const int m, const Policy& pol)
{
BOOST_MATH_STD_USING // ADL of std names, needed for ceil.
BOOST_MATH_STD_USING // ADL of std names, needed for floor.
// Compute an estimate of the m'th root of cyl_neumann.
@@ -448,12 +457,17 @@
{
// Create the positive order and extract its positive floor and ceiling integer parts.
const T vv(-v);
const T vv_ceil (ceil (vv));
const T vv_floor(floor(vv));
// The to-be-found root is bracketed by the roots of the
// Bessel function whose reflected, positive integer order
// is greater than, but nearest to vv.
// is less than, but nearest to vv.
// The special case of negative, half-integer order uses
// the relation between Yv and spherical Bessel functions
// in order to obtain the bracket for the root.
// In these special cases, cyl_neumann(-n/2, x) = sph_bessel_j(+n/2, x)
// for v = -n/2.
T root_hi;
T root_lo;
@@ -462,6 +476,8 @@
{
// The estimate of the first root for negative order is found using
// an adaptive range-searching algorithm.
// Take special precautions for the discontinuity at negative,
// half-integer orders and use different brackets above and below these.
if(T(vv - vv_floor) < 0.5F)
{
root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
@@ -487,7 +503,14 @@
root_hi = root_lo;
// Decrease the lower end of the bracket using an adaptive algorithm.
(root_lo > 0.5F) ? root_lo -= 0.5F : root_lo *= 0.85F;
if(root_lo > 0.5F)
{
root_lo -= 0.5F;
}
else
{
root_lo *= 0.85F;
}
}
}
else