BesselK Function in C#

Question

I am attempting to implement the BesselK method from Boost (a C++ library). The Boost method accepts two doubles and returns a double. (I have it implemented below as cyl_bessel_k .)

The equation I modeled this off of comes from Boosts documention: http://www.boost.org/doc/libs/1_45_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/mbessel.html

I have also been checking values against Wolfram: http://www.wolframalpha.com/input/?i=BesselK%283%2C1%29

I am able to match output from the Boost method when passing a positive non-integer value for "v". However, when an integer is passed, my output is severely off. So,there is an obvious discontinuity issue. From reading up on this, it seems that this issue arises from passing a negative integer to the gamma function. Somehow reflection comes into play here with the Bessel_I method, but I'm nearing the end of my math skillset.

1.) What needs to happen to the bessel_i method with reflection to make this work?

2.) I'm currently doing a partial sum approach. Boost uses a continuous fraction approach. How can I modify this to account for convergence?

Any input is appreciated! Thank you!

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

neonScarecrow · Answer 1 · 2015-07-22T20:46:43.107

After lots of refactoring and trying things that didn't work, this is what I came up with. It's mostly Boost logic that has been adapted and translated into C#.

It's not perfect though (likely due to rounding, precision,etc). Any improvements are welcome! Max error is 0.0000001926% between true Bessel_K value from Wolfram and my adapted method. This is occurs when parameter 'v' is an integer. For my purposes, this was close enough.

Link to fiddle: https://dotnetfiddle.net/QIYzK6

Hopefully it saves someone some headache.

BesselK Function in C#

1 Answers1