Are these functions correct?


(David Rutten) #1

I’m working on a visual programming platform which uses and re-exposes a bunch of MathNet functionality as part of the standard toolkit. I state this merely to indicate I’m not solving specific problems at this point, but just trying to figure out what parts of MathNet would be useful for my own users to have access to.

As part of my test project I’ve started generating graphs of various functions and some of them look nothing like the other examples I see online. Specifically the various kinds of Bessel functions don’t have that wavy character like the graphics on Wikipedia show.

For example here’s Second kind, order=0, exponential:

The DiGamma inverse also doesn’t look at all like the inverse of the DiGamma results.

I’ll be the first to admit I don’t know what one would use these functions for, and maybe I’m just not looking at them right. Here’s a collection of the graphs my test project is outputting, in case someone wants to have a look.


(Peter Vanderwaart) #2

I can’t help you but (1) I was able to reproduce your result above for BesselK0e, (2) the graph of BesselK0 does not agree with Wikipedia, and (3) it’s hard to reconcile the results of BesselK0e with those of BesselK0.


(Peter Vanderwaart) #3

It looks like BesselK0 is the Modified Bessel Function. Compare to Figure 4.4 in the pdf found here. The documentation does say “modified Bessel function.”

The Modified Bessel Functions are very different of the Bessel Functions, not just an alternate representation. It’s explained fairly clearly in the material at the link.


(David Rutten) #4

Indeed, seems like it was right all along (except maybe for the exponential versions, I’m still not sure whether that discontinuity at x=2 should be there). Thanks for looking into it.


(Peter Vanderwaart) #5

I saw a reference to using one approximation for (0,x) and second model for (x,infinity) for some x. I think that’s what is causing the discontinuity: two models that don’t agree exactly at the cutoff point.


(David Rutten) #6

Here’s a Bessel First Kind order=n implementation that seems to work:

private static double BesselFirstKind(int alpha, double x)
{
  if (alpha < 0)
    throw new ArgumentException(nameof(alpha) + " must be strictly positive.");

  double sum = 0.0;
  for (int m = 0; m < 1000; m++)
  {
    BigInteger gamma = Factorial(m + alpha);
    BigInteger factorial = Factorial(m);

    double denominator = (double)(gamma * factorial);
    if (double.IsInfinity(denominator))
      break;

    double firstPart = m.IsEven() ?
      +1.0 / denominator :
      -1.0 / denominator;

    double secondPart = Math.Pow(0.5 * x, 2 * m + alpha);
    double newSum = sum + (firstPart * secondPart);
    ulong proximity = sum.NumbersBetween(newSum);
    double accuracy = Math.Abs(sum - newSum);
    sum = newSum;

    if (proximity < 3 || accuracy < 1e-64)
      break;
  }

  return sum;
}


(David Rutten) #7

Yeah never mind it seeming to work. Only within a very limited domain…

Apologies for keeping on amending these posts. I’ve now settled on the following, which is an approximate solution but does not suffer from the same numeric instability:

public static double BesselFirstKindIntegral(int order, double x)
{
  if (order < 0)
    throw new ArgumentException(nameof(order) + " may not be negative.");

  double samplingAccuracy = _besselSamples[1] - _besselSamples[0];

  double integral = 0.0;
  for (int i = 0; i < _besselSamples.Length - 1; i++)
  {
    double t = _besselSamples[i];
    double y = Math.Cos((order * t) - x * Math.Sin(t));
    integral += y * samplingAccuracy;
  }

  return OnePiInverse * integral;
}
private static readonly double[] _besselSamples = (0.0).ArrayTo(OnePi, 1000);

(Peter Vanderwaart) #8

FWIW, you can use the functions in Numerics.Integration to get a more sophisticated method of integration without much trouble. For example:

public double BesselOne(int order, double x)
        {
            var f0 = new Func<int, double, double, double>((o, v, t) => Math.Cos(o * t - v * Math.Sin(t)));  // integrand
            double sum = NewtonCotesTrapeziumRule.IntegrateComposite(t => f0(order, x, t), 0.0, Math.PI, 250);

            return sum / Math.PI;
        }

I think that no matter what you do, it’s going to go wonky at for some large enough value of order or x, so test it to be sure it’s working in the desired range.