第一类和第二类贝塞尔函数的零点查找

概要

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

用于获取贝塞尔函数的单个零点或根，以及通过提供输出迭代器将多个零点放入像 std::vector 这样的容器中的函数。

单值函数的签名是

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.

以及多个零点的签名是

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.

还有一些版本允许控制用于错误处理和精度的策略。

 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.

描述

每个实数阶 ν 柱贝塞尔函数和诺伊曼函数在正实轴上都有无数个零点。正实轴上的实数零点可以通过求解以下方程的根来找到

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

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

这里，j_{ν, m} 表示阶数为 ν 的柱贝塞尔函数的第 m^th 个根，而 y_{ν, m} 表示阶数为 ν 的柱诺伊曼函数的第 m^th 个根。

贝塞尔函数和诺伊曼函数的零点或根（函数穿过水平 y = 0 轴的 x 值）由两个函数 cyl_bessel_j_zero 和 cyl_neumann_zero 计算。

在每种情况下，返回的零点的索引或秩都是从 1 开始的，也就是说

cyl_bessel_j_zero(v, 1);

返回贝塞尔 J 的第一个零点。

传递 start_index <= 0 会导致引发 std::domain_error。

但是，对于某些参数，定义了零阶根，其值为零。例如，对于所有 v > 0 的值以及 v = -n 的负整数值，J[sub v](x) 的零阶根被定义，并且其值为零。下面在实现细节中描述了类似的情况。

对于 cyl_bessel_j 和 cyl_neumann 函数，J 的阶数 v 可以是正数、负数和零，但不能是无穷大或 NaN。

查找贝塞尔和诺伊曼零点的示例

此示例演示了计算贝塞尔函数和诺伊曼函数的零点。它还展示了如何将 Boost.Math 和 Boost.Multiprecision 结合起来以提供多位十进制数字精度。对于 50 位十进制数字精度，我们需要包含

#include <boost/multiprecision/cpp_dec_float.hpp>

以及 typedef 用于 float_type 可能很方便（允许快速切换以在内置 double 或其他精度下重新计算）

typedef boost::multiprecision::cpp_dec_float_50 float_type;

要使用查找函数零点的函数，我们需要

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

此文件包含查找零点函数的前向声明签名

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

但更多详细信息请参见完整文档，例如 Boost.Math 贝塞尔函数。

此示例显示了如何获取贝塞尔函数的单个零点，然后通过提供迭代器将多个零点放入像 std::vector 这样的容器中。

	提示
	始终明智的做法是将使用 Boost.Math 的代码放在 try'n'catch 块中；这将确保在出现异常情况时显示有用的错误消息。

首先，评估单个贝塞尔零点。

精度由 v 的模板参数 T 的浮点类型控制，因此此示例具有 double 精度，至少 15 位但最多 17 位十进制数字（对于常见的 64 位 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

但请注意，由于参数 v 控制结果的精度，v 必须是浮点类型。因此，如果您提供整数类型，例如 0，而不是 0.0，则编译将失败，如下所示

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

带有此错误消息

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

或者，我们可以使用策略来忽略错误，C 风格，返回一些值，可能是无穷大或 NaN，或尽力而为。（参见用户错误处理）。

要创建忽略所有错误的（可能不明智的！）策略 ignore_all_policy

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;

使用此 ignore_all_policy 的示例是

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

cyl_bessel_j_zero 的另一个版本允许一次调用计算多个零点，并将结果放入容器中，通常是 std::vector。例如，生成并显示积分阶 2 的 J_v 的前五个 double 根，如 Wolfram 贝塞尔函数零点表 1 中的列 J₂(x)。

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"));

或者我们可以使用 Boost.Multiprecision 生成非积分阶 v= 71/19 == 3.736842 的 J_v 的 50 位十进制数字根，表示为精确整数分数，以便为所有浮点类型生成最精确的值。

我们设置输出流的精度，并显示尾随零以显示固定的 50 位十进制数字。

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"));

使用输出迭代器求贝塞尔函数零点之和

此示例演示了求贝塞尔函数零点之和。要使用查找函数零点的函数，我们需要

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

我们使用 cyl_bessel_j_zero 输出迭代器参数 out_it 通过定义自定义输出迭代器来创建 1/zeros² 的总和

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;
};

总和是为许多值计算的，收敛于 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

计算诺伊曼函数的零点。

此示例还展示了如何将 Boost.Math 和 Boost.Multiprecision 结合起来以提供多位十进制数字精度。对于 50 位十进制数字精度，我们需要包含

#include <boost/multiprecision/cpp_dec_float.hpp>

以及 typedef 用于 float_type 可能很方便（允许快速切换以在内置 double 或其他精度下重新计算）

typedef boost::multiprecision::cpp_dec_float_50 float_type;

要使用查找 cyl_neumann 函数零点的函数，我们需要

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

诺伊曼（贝塞尔 Y）函数零点的评估方式非常相似

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

来自“错误”输入的错误消息

另一个示例演示了计算贝塞尔函数的零点，显示了通过抛出异常处理来自“错误”输入的错误消息。

要使用查找函数零点的函数，我们需要

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

	提示
	始终明智的做法是将所有使用 Boost.Math 的代码放在 try'n'catch 块中；这将确保在出现异常情况时可以显示有用的错误消息。

下面的示例显示了来自几个“错误”参数的消息，这些参数抛出了 domain_error 异常。

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;
}

	注意
	错误消息中显示的类型是提升后的类型，使用来自 `float` 到 `double` 的精度策略和内部提升策略，在本例中。

在此示例中，提升过程如下

参数是 float 和 int。
将 int “视为”它是一个 double，因此参数是 float 和 double。
公共类型是 double - 所以这就是我们想要的精度（以及将返回的类型）。
在内部评估为 double 以获得完整的 float 精度。

有关从 double 提升到 long double 的其他示例，请参见完整代码。

下面是其他“错误”输入（如无穷大和 NaN）的示例。一些编译器警告表明在编译时检测到“错误”值。

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;
}

其他示例的输出显示在完整代码列表的末尾。

这些示例的完整代码（和输出）位于贝塞尔零点、贝塞尔零点迭代器、诺伊曼零点、贝塞尔错误消息。

实现

使用各种方法来计算 j_{ν, m} 和 y_{ν, m} 的初始估计值；这些在下面详细描述。

在找到给定根的初始估计值后，随后使用 Boost.Math 的使用导数的求根实用程序结合函数 cyl_bessel_j 和 cyl_neumann，将精度细化到所需的水平，使用牛顿-拉夫森迭代法。

牛顿迭代法需要 J_ν(x) 或 Y_ν(x) 及其导数。J_ν(x) 和 Y_ν(x) 关于 x 的导数由 M. Abramowitz 和 I. A. Stegun，《数学函数手册》，NBS (1964) 给出。特别是，

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

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

根的秩的枚举（换句话说，根的索引）从 1 开始并向上计数，换句话说 m,=1,2,3,… 第一个根的值始终大于零。

对于某些特殊参数，柱贝塞尔函数和柱诺伊曼函数在原点处有一个根。例如，对于每个正阶 ν > 0，以及对于每个负整数阶 ν = -n 且 n ∈ ℕ ⁺ 和 n ≠ 0，J_ν(x) 在原点处有一个根。

此外，对于每个负半整数阶 ν = -n/2，且 n ∈ ℕ ⁺ 和 n ≠ 0，Y_ν(x) 在原点处有一个根。

对于这些特殊的参数值，原点的值为 x = 0，作为 cyl_bessel_j_zero() 和 cyl_neumann_zero() 生成的 0^th 根提供。

在计算贝塞尔函数根的初始估计值时，正阶和负阶之间存在区别，并且对它们使用不同的方法。此外，对于第一个根 m = 1 和随后的更高秩的根 m ≥ 2，使用不同的算法。此外，对于阶数高于和低于 ν = 2.2 的贝塞尔函数根的估计值，采用不同的方法计算。

使用经验列表值计算 0 ≤ ν < 2.2 的 j_ν,1 和 y_ν,1 的估计值。这些系数是由计算机代数系统生成的。

使用 M. Abramowitz 和 I. A. Stegun，《数学函数手册》，NBS (1964) 中的公式 9.5.14 和 9.5.15 计算 ν≥ 2.2 的 j_ν,1 和 y_ν,1 的估计值。

特别是，

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

和

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

使用 McMahon 近似值计算秩 m > 2 和 0 ≤ ν < 2.2 的 j_{ν, m} 和 y_{ν, m} 的估计值，如 M. Abramowitz 和 I. A. Stegan，第 9.5 节和 9.5.12 节中所述。特别是，

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)

其中 μ = 4ν² 和 β = (m + ½ν - ¼)π 对于 j_ν,m，β = (m + ½ν -¾)π 对于 y_ν,m。

使用 M. Abramowitz 和 I. A. Stegun，《数学函数手册》，NBS (1964) 中的公式 9.5.22 和公式 9.5.26 的顶行以及公式 9.3.39 中的渐近展开式中的一项，计算 ν ≥ 2.2 的 j_{ν, m} 和 y_{ν, m} 的估计值，其中明确且易于理解地处理零点的渐近展开。后两个方程针对参数 (x) 大于 1 的情况表示。（Olver 还在 §10.21(vi) McMahon 大型零点的渐近展开式中给出了方程的级数形式 - 使用略有不同的变量名）。

总之，

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

其中 -ζ = ν^-2/3a_m，a_m 是负实轴上 Ai(x) 的第 m^th 个根的绝对值。

这里 x = x(-ζ) 是函数的反函数

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

(7)

此外，

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

其中

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

和

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

在求解公式 7 中的 x(-ζ) 时，对于大的 x，右侧展开为 2 阶泰勒级数。这导致

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

所得二次方程的正根用于查找初始估计值 x(-ζ)。随后在公式 7 中使用牛顿-拉夫森迭代的几个步骤来细化此初始估计值。

使用交错关系找到正实轴上负阶柱贝塞尔函数根的估计值。例如，柱贝塞尔函数 j_-ν,m 的第 m^th 个根由相应正整数阶的贝塞尔函数的第 m^th 个根和第 (m+1)^th 个根括起来。换句话说，

j_{nν, m} < j_{-ν, m} < j_{nν, m+1}

其中 m > 1，n_ν 表示 |-ν| 的绝对值的整数底。

类似的括号关系用于查找负阶诺伊曼函数根的估计值，由此需要处理每个负半整数阶的不连续性。

括号关系不适用于负阶柱贝塞尔函数和柱诺伊曼函数的第一个根。因此，迭代算法与二分法求根相结合用于定位 j_-ν,1 和 y_-ν,1。

测试

使用 cpp_dec_float_50 在 50 位十进制数字的精度下测试了零点评估的精度，发现与 Wolfram Alpha 计算的点值相同。