Boost C++ 库

...世界上最受推崇和专业设计的 C++ 库项目之一。 — Herb Sutter 和 Andrei Alexandrescu, C++ 编码标准

PrevUpHomeNext

第一类和第二类修正贝塞尔函数

概要

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

template <class T1, class T2>
BOOST_MATH_GPU_ENABLED calculated-result-type cyl_bessel_i(T1 v, T2 x);

template <class T1, class T2, class Policy>
BOOST_MATH_GPU_ENABLED calculated-result-type cyl_bessel_i(T1 v, T2 x, const Policy&);

template <class T1, class T2>
BOOST_MATH_GPU_ENABLED calculated-result-type cyl_bessel_k(T1 v, T2 x);

template <class T1, class T2, class Policy>
BOOST_MATH_GPU_ENABLED calculated-result-type cyl_bessel_k(T1 v, T2 x, const Policy&);
描述

函数 cyl_bessel_i 和 cyl_bessel_k 分别返回第一类和第二类修正贝塞尔函数的结果

cyl_bessel_i(v, x) = Iv(x)

cyl_bessel_k(v, x) = Kv(x)

其中

当 T1 和 T2 是不同类型时，这些函数的返回类型使用 结果类型计算规则 计算。这些函数也针对 T1 是整数的相对常见情况进行了优化。

最后的 策略 参数是可选的，可用于控制函数的行为：如何处理错误，使用什么级别的精度等。有关更多详细信息，请参阅 策略文档。

当结果未定义或为复数时，函数返回 domain_error。对于 cyl_bessel_j，当 x < 0 且 v 不是整数，或者当 x == 0 且 v != 0 时，会发生这种情况。对于 cyl_neumann，当 x <= 0 时，会发生这种情况。

下图说明了 Iv 的指数行为。

下图说明了 Kv 的指数衰减。

测试

有两组测试值：使用 functions.wolfram.com 计算的点值，以及使用此实现的简化版本（去除了所有特殊情况处理）计算的更大的一组测试。

精度

下表显示了这些函数在各种平台上的精度变化情况，以及与其他库的比较。请注意，仅给出系统上最宽浮点类型的结果，因为较窄的类型具有 实际上为零的误差。所有值都是以 epsilon 为单位的相对误差。请注意，我们的测试套件包含一些相当极端的输入，这导致其他库中出现大多数最坏的问题案例

表 8.44. cyl_bessel_i（整数阶）的错误率

GNU C++ 版本 7.1.0
linux
双精度

GNU C++ 版本 7.1.0
linux
长双精度

Sun 编译器版本 0x5150
Sun Solaris
长双精度

Microsoft Visual C++ 版本 14.1
Win32
双精度

Bessel I0：Mathworld 数据（整数版本）

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 0.79ε (平均值 = 0.482ε))
(Rmath 3.2.3: 最大值 = 1.52ε (平均值 = 0.622ε) 以及其他失败。)

最大值 = 1.95ε (平均值 = 0.738ε)

(<cmath>: 最大值 = 8.49ε (平均值 = 3.46ε) 以及其他失败。)

最大值 = 1.95ε (平均值 = 0.661ε)

最大值 = 0.762ε (平均值 = 0.329ε)

Bessel I1：Mathworld 数据（整数版本）

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 0.82ε (平均值 = 0.456ε))
(Rmath 3.2.3: 最大值 = 1.53ε (平均值 = 0.483ε) 以及其他失败。)

最大值 = 0.64ε (平均值 = 0.202ε)

(<cmath>: 最大值 = 5ε (平均值 = 2.15ε) 以及其他失败。)

最大值 = 0.64ε (平均值 = 0.202ε)

最大值 = 0.767ε (平均值 = 0.398ε)

Bessel In：Mathworld 数据（整数版本）

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 5.15ε (平均值 = 2.13ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 1.73ε (平均值 = 0.601ε) 以及其他失败。)

最大值 = 1.8ε (平均值 = 1.33ε)

(<cmath>: 最大值 = 430ε (平均值 = 163ε) 以及其他失败。)

最大值 = 463ε (平均值 = 140ε)

最大值 = 3.46ε (平均值 = 1.32ε)


表 8.45. cyl_bessel_i 的错误率

GNU C++ 版本 7.1.0
linux
双精度

GNU C++ 版本 7.1.0
linux
长双精度

Sun 编译器版本 0x5150
Sun Solaris
长双精度

Microsoft Visual C++ 版本 14.1
Win32
双精度

Bessel I0：Mathworld 数据

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 270ε (平均值 = 91.6ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 1.52ε (平均值 = 0.622ε) 以及其他失败。)

最大值 = 1.95ε (平均值 = 0.738ε)

(<cmath>: 最大值 = 8.49ε (平均值 = 3.46ε) 以及其他失败。)

最大值 = 1.95ε (平均值 = 0.661ε)

最大值 = 0.762ε (平均值 = 0.329ε)

Bessel I1：Mathworld 数据

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 128ε (平均值 = 41ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 1.53ε (平均值 = 0.483ε) 以及其他失败。)

最大值 = 0.64ε (平均值 = 0.202ε)

(<cmath>: 最大值 = 5ε (平均值 = 2.15ε) 以及其他失败。)

最大值 = 0.64ε (平均值 = 0.202ε)

最大值 = 0.767ε (平均值 = 0.398ε)

Bessel In：Mathworld 数据

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 2.31ε (平均值 = 0.838ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 1.73ε (平均值 = 0.601ε) 以及其他失败。)

最大值 = 1.8ε (平均值 = 1.33ε)

(<cmath>: 最大值 = 430ε (平均值 = 163ε) 以及其他失败。)

最大值 = 463ε (平均值 = 140ε)

最大值 = 3.46ε (平均值 = 1.32ε)

Bessel Iv：Mathworld 数据

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 5.95ε (平均值 = 2.08ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 3.53ε (平均值 = 1.39ε))

最大值 = 4.12ε (平均值 = 1.85ε)

(<cmath>: 最大值 = 616ε (平均值 = 221ε) 以及其他失败。)

最大值 = 4.12ε (平均值 = 1.95ε)

最大值 = 2.97ε (平均值 = 1.24ε)

Bessel In：随机数据

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 261ε (平均值 = 53.2ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 7.37ε (平均值 = 2.4ε))

最大值 = 4.62ε (平均值 = 1.06ε)

(<cmath>: 最大值 = 645ε (平均值 = 132ε))

最大值 = 176ε (平均值 = 39.1ε)

最大值 = 9.67ε (平均值 = 1.88ε)

Bessel Iv：随机数据

最大值 = 0.661ε (平均值 = 0.0441ε)

(GSL 2.1: 最大值 = 6.18e+03ε (平均值 = 1.55e+03ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 4.28e+08ε (平均值 = 2.85e+07ε))

最大值 = 8.35ε (平均值 = 1.62ε)

(<cmath>: 最大值 = 1.05e+03ε (平均值 = 224ε) 以及其他失败。)

最大值 = 283ε (平均值 = 88.4ε)

最大值 = 7.46ε (平均值 = 1.71ε)

Bessel Iv：Mathworld 数据（大值）

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 37ε (平均值 = 18ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 3.77e+168ε (平均值 = 2.39e+168ε) 以及其他失败。)

最大值 = 14.7ε (平均值 = 6.66ε)

(<cmath>: 最大值 = 118ε (平均值 = 57.2ε) 以及其他失败。)

最大值 = 14.7ε (平均值 = 6.59ε)

最大值 = 3.67ε (平均值 = 1.64ε)


表 8.46. cyl_bessel_k（整数阶）的错误率

GNU C++ 版本 7.1.0
linux
长双精度

GNU C++ 版本 7.1.0
linux
双精度

Sun 编译器版本 0x5150
Sun Solaris
长双精度

Microsoft Visual C++ 版本 14.1
Win32
双精度

Bessel K0：Mathworld 数据（整数版本）

最大值 = 0.833ε (平均值 = 0.436ε)

(<cmath>: 最大值 = 9.33ε (平均值 = 3.25ε))

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 1.2ε (平均值 = 0.733ε))
(Rmath 3.2.3: 最大值 = 0.833ε (平均值 = 0.601ε))

最大值 = 0.833ε (平均值 = 0.436ε)

最大值 = 0.833ε (平均值 = 0.552ε)

Bessel K1：Mathworld 数据（整数版本）

最大值 = 0.786ε (平均值 = 0.329ε)

(<cmath>: 最大值 = 8.94ε (平均值 = 3.19ε))

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 0.626ε (平均值 = 0.333ε))
(Rmath 3.2.3: 最大值 = 0.894ε (平均值 = 0.516ε))

最大值 = 0.786ε (平均值 = 0.329ε)

最大值 = 0.786ε (平均值 = 0.39ε)

Bessel Kn：Mathworld 数据（整数版本）

最大值 = 2.6ε (平均值 = 1.21ε)

(<cmath>: 最大值 = 12.9ε (平均值 = 4.91ε) 以及其他失败。)

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 168ε (平均值 = 59.5ε))
(Rmath 3.2.3: 最大值 = 8.48ε (平均值 = 2.98ε))

最大值 = 2.6ε (平均值 = 1.21ε)

最大值 = 3.63ε (平均值 = 1.46ε)


表 8.47. cyl_bessel_k 的错误率

GNU C++ 版本 7.1.0
linux
长双精度

GNU C++ 版本 7.1.0
linux
双精度

Sun 编译器版本 0x5150
Sun Solaris
长双精度

Microsoft Visual C++ 版本 14.1
Win32
双精度

Bessel K0：Mathworld 数据

最大值 = 0.833ε (平均值 = 0.436ε)

(<cmath>: 最大值 = 9.33ε (平均值 = 3.25ε))

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 6.04ε (平均值 = 2.16ε))
(Rmath 3.2.3: 最大值 = 0.833ε (平均值 = 0.601ε))

最大值 = 0.833ε (平均值 = 0.436ε)

最大值 = 0.833ε (平均值 = 0.552ε)

Bessel K1：Mathworld 数据

最大值 = 0.786ε (平均值 = 0.329ε)

(<cmath>: 最大值 = 8.94ε (平均值 = 3.19ε))

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 6.26ε (平均值 = 2.21ε))
(Rmath 3.2.3: 最大值 = 0.894ε (平均值 = 0.516ε))

最大值 = 0.786ε (平均值 = 0.329ε)

最大值 = 0.786ε (平均值 = 0.39ε)

Bessel Kn：Mathworld 数据

最大值 = 2.6ε (平均值 = 1.21ε)

(<cmath>: 最大值 = 12.9ε (平均值 = 4.91ε) 以及其他失败。)

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 3.36ε (平均值 = 1.43ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 8.48ε (平均值 = 2.98ε))

最大值 = 2.6ε (平均值 = 1.21ε)

最大值 = 3.63ε (平均值 = 1.46ε)

Bessel Kv：Mathworld 数据

最大值 = 3.58ε (平均值 = 2.39ε)

(<cmath>: 最大值 = 13ε (平均值 = 4.81ε) 以及其他失败。)

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 5.47ε (平均值 = 2.04ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 3.15ε (平均值 = 1.35ε))

最大值 = 5.21ε (平均值 = 2.53ε)

最大值 = 4.78ε (平均值 = 2.19ε)

Bessel Kv：Mathworld 数据（大值）

最大值 = 42.3ε (平均值 = 21ε)

(<cmath>: 最大值 = 42.3ε (平均值 = 19.8ε) 以及其他失败。)

最大值 = 0ε (平均值 = 0ε)

(GSL 2.1: 最大值 = 308ε (平均值 = 142ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 84.6ε (平均值 = 37.8ε))

最大值 = 42.3ε (平均值 = 21ε)

最大值 = 59.8ε (平均值 = 26.9ε)

Bessel Kn：随机数据

最大值 = 4.55ε (平均值 = 1.12ε)

(<cmath>: 最大值 = 13.9ε (平均值 = 2.91ε))

最大值 = 0.764ε (平均值 = 0.0348ε)

(GSL 2.1: 最大值 = 8.71ε (平均值 = 1.76ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 7.47ε (平均值 = 1.34ε))

最大值 = 4.55ε (平均值 = 1.12ε)

最大值 = 9.34ε (平均值 = 1.7ε)

Bessel Kv：随机数据

最大值 = 7.88ε (平均值 = 1.48ε)

(<cmath>: 最大值 = 13.6ε (平均值 = 2.68ε) 以及其他失败。)

最大值 = 0.507ε (平均值 = 0.0313ε)

(GSL 2.1: 最大值 = 9.71ε (平均值 = 1.47ε) 以及其他失败。)
(Rmath 3.2.3: 最大值 = 7.37ε (平均值 = 1.49ε))

最大值 = 7.88ε (平均值 = 1.47ε)

最大值 = 8.33ε (平均值 = 1.62ε)


以下误差图基于对 I0、I1、K0 和 K1 函数域的详尽搜索，MSVC-15.5 在 double 精度下，以及 GCC-7.1/Ubuntu 在 long double 和 __float128 精度下。

实现

以下情况首先作为特殊情况处理

当计算 x < 0 的 Iv 时，则 ν 必须是整数，否则会发生域错误。如果 ν 是整数，则如果 ν 是奇数，则函数为奇函数，如果 ν 是偶数，则函数为偶函数，我们可以反射到 x > 0。

对于 v 等于 0、1 或 0.5 的 Iv，作为特殊情况处理。

0 和 1 的情况在有限和无限区间上使用多项式近似。近似形式基于 Pavel Holoborodko 的 “第一类修正贝塞尔函数的有理逼近 - 用于双精度计算的 I0(x)”，我们将其扩展到处理高达 128 位精度（每种目标精度使用不同的近似值）。

类似地，我们有

0.5 的情况是一个简单的三角函数

I0.5(x) = sqrt(2 / πx) * sinh(x)

对于 v 为整数的 Kv，结果使用递推关系计算

从 K0 和 K1 开始，它们使用上面的有理近似计算。这些有理近似的精度约为 19 位数字，因此仅在 T 的精度不超过 64 个二进制数字时使用。

当 x 相对于 v 较小时，最好直接从级数计算 Ivx

在一般情况下，我们首先借助反射公式将 ν 归一化为 [0, [inf])

令 μ = ν - floor(ν + 1/2)，则 μ 是 ν 的小数部分，使得 |μ| <= 1/2（我们稍后需要它来收敛）。我们的想法是计算 Kμ(x) 和 Kμ+1(x)，并使用它们来获得 Iν(x) 和 Kν(x)。

该算法由 Temme 在以下文献中提出

N.M. Temme, 关于第三类修正贝塞尔函数的数值评估, Journal of Computational Physics, vol 19, 324 (1975),

这需要两个连分数以及 Wronskian 行列式

连分数使用改进的 Lentz 方法计算

(W.J. Lentz, 使用连分数生成米散射计算中的贝塞尔函数, Applied Optics, vol 15, 668 (1976)).

它们的收敛速度取决于 x，因此对于大的 x 和小的 x，我们需要不同的策略。

x > v, CF1 需要 O(x) 次迭代才能收敛，CF2 收敛迅速。

x <= v, CF1 收敛迅速，当 x -> 0 时，CF2 无法收敛。

当 x 较大时（x > 2），两个连分数都收敛（对于非常大的 x，CF1 可能很慢）。Kμ 和 Kμ+1 可以通过以下方式计算

其中

S 也是一个与 CF2 一起求和的级数，请参阅

I.J. Thompson 和 A.R. Barnett, 实阶和复宗量的修正贝塞尔函数 I_v 和 K_v，达到选定的精度, Computer Physics Communications, vol 47, 245 (1987).

当 x 较小时（x <= 2），CF2 收敛可能失败（但 CF1 工作得非常好）。这里的解决方案是 Temme 的级数

其中

fk 和 hk 也通过递归（涉及 gamma 函数）计算，但公式有点复杂，读者可以参考

N.M. Temme, 关于第三类修正贝塞尔函数的数值评估, Journal of Computational Physics, vol 19, 324 (1975).

注意：Temme 的级数仅对 |μ| <= 1/2 收敛。

然后从正向递推计算 Kν(x)，Kν+1(x) 也是如此。有了这两个值和 fν，Wronskian 行列式直接产生 Iν(x)。

版权所有 © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle Walker 和 Xiaogang Zhang

根据 Boost 软件许可协议 1.0 版分发。（请参阅随附文件 LICENSE_1_0.txt 或在 https://boost.ac.cn/LICENSE_1_0.txt 复制）

PrevUpHomeNext