c# - C# 中的 BesselK 函数

Question

我正在尝试从 Boost(一个 C++ 库)实现 BesselK 方法。 Boost 方法接受两个 double 值并返回一个 double 值。 (我在下面将其实现为 cyl_bessel_k 。)

我对此建模的方程式来自 Boosts 文档: http://www.boost.org/doc/libs/1_45_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/mbessel.html

我也一直在检查 Wolfram 的值: http://www.wolframalpha.com/input/?i=BesselK%283%2C1%29

当为“v”传递一个正的非整数值时，我能够匹配 Boost 方法的输出。但是，当传递一个整数时，我的输出严重关闭。因此，存在明显的不连续性问题。通过阅读这篇文章，这个问题似乎是由将负整数传递给 Gamma 函数引起的。不知何故，反射在这里与 Bessel_I 方法一起发挥作用，但我的数学技能已接近尾声。

1.) 带反射的 bessel_i 方法需要做什么才能使这项工作正常进行？

2.) 我目前正在使用部分求和法。 Boost 使用连续分数法。我如何修改它以考虑收敛？

欢迎任何意见!谢谢!

    static double cyl_bessel_k(double v, double x)
    {
        if (v > 0)
        {
            double iNegativeV = cyl_bessel_i(-v, x);
            double iPositiveV = cyl_bessel_i(v, x);
            double besselSecondKind = (Math.PI / 2) * ((iNegativeV - iPositiveV ) / (Math.Sin(Math.PI * v)));
            return besselSecondKind;
        }
        else
        {
           //error handling
        }
    }

    static double cyl_bessel_i(double v, double x)
    {
        if (x == 0)
        {
            return 0;
        }
        double summed = 0;
        double a = Math.Pow((0.5d * x), v);
        for (double k = 0; k < 10; k++) //how to account for convergence? 10 is arbitrary
        {
            double b = Math.Pow(0.25d * Math.Pow(x, 2), k);
            double kFactorial = SpecialFunctions.Factorial((int)k); //comes from MathNet.Numerics (Nuget)
            double gamma = SpecialFunctions.Gamma(v + k + 1); //comes from MathNet.Numerics
            summed += b / (kFactorial * gamma);
        }
        return a * summed;
    }

Answer 1

经过大量重构和尝试无效的方法后，这就是我想出的。它主要是已被改编并翻译成 C# 的 Boost 逻辑。

虽然它并不完美(可能是由于四舍五入、精度等原因)。欢迎任何改进! Wolfram 的真实 Bessel_K 值与我采用的方法之间的最大误差为 0.0000001926%。当参数“v”是整数时会发生这种情况。就我的目的而言，这已经足够接近了。

fiddle 链接: https://dotnetfiddle.net/QIYzK6

希望它能减轻一些人的头痛。