1

我正在尝试从 Boost（一个 C++ 库）中实现 BesselK 方法。Boost 方法接受两个双精度并返回一个双精度。（我在下面将它实现为 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;
    }
4

1 回答 1

1

经过大量的重构和尝试不起作用的事情，这就是我想出的。它主要是 Boost 逻辑已被改编并翻译成 C#。

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

小提琴链接： https ://dotnetfiddle.net/QIYzK6

希望它可以避免一些人头疼。

于 2015-07-22T20:08:15.067 回答