bessel/bessel.md

Code

The code is embedded in this document! VOC can run through this file and generate bessel.js.

I really dislike writing Math so I use:

var M = Math;

Horner Method

The methods use an approximating polynomial and evaluate using Horner's method:

function _horner(arr/*:Array<number>*/, v/*:number*/) { return arr.reduce(function(z,w){return v * z + w;},0); }

Recurrence

It can be shown that the four Bessel functions satisfy (on their support):

B_{n} (x) = \frac{2n}{x} B_{n-1}(x) - B_{n-2}(x)

So rather than go back and try to find solution for each order, we will build solutions for n=0 and n=1 and then apply the recurrence. The helper:

function _bessel_iter(x/*:number*/, n/*:number*/, f0/*:number*/, f1/*:number*/, sign/*:?number*/) {
  if(!sign) sign = -1;
  var tdx = 2 / x, f2;
  if(n === 0) return f0;
  if(n === 1) return f1;
  for(var o = 1; o != n; ++o) {
    f2 = f1 * o * tdx + sign * f0;
    f0 = f1; f1 = f2;
  }
  return f1;
}

We can directly generate the JS function given the basic solutions bessel0 and bessel1 by leveraging _bessel_iter from above. We have to add a few sanity checks since Y_n is undefined at 0 and K_n is real only when x>0

function _bessel_wrap(bessel0, bessel1, name, nonzero, sign) {
  return function bessel(x/*:number*/,n/*:number*/) {
    if(n === 0) return bessel0(x);
    if(n === 1) return bessel1(x);
    if(n < 0) throw name + ': Order (' + n + ') must be nonnegative';
    if(nonzero == 1 && x === 0) throw name + ': Undefined when x == 0';
    if(nonzero == 2 && x <= 0) throw name + ': Undefined when x <= 0';
    var b0 = bessel0(x), b1 = bessel1(x);
    return _bessel_iter(x, n, b0, b1, sign);
  };
}

Individual Solutions

To determine each individual solution, we first calculate a Chebyshev polynomial based on the regime (lower and higher-value approximations). This module uses the constants from the GNU Scientific Library, and from the venerable Numerical Recipes book but I have independently verified the constants using Mathematica.

var besselj = (function() {

The constants are named b[01]_[ab]([12][ab])? with the blocks corresponding to the function order (e.g. b0_ refers to order 0), variable name in the function and conditional.

  var b0_a1a = [57568490574.0,-13362590354.0,651619640.7,-11214424.18,77392.33017,-184.9052456].reverse();
  var b0_a2a = [57568490411.0,1029532985.0,9494680.718,59272.64853,267.8532712,1.0].reverse();
  var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse();
  var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse();

I noticed some strange oddities when leaning on Math.PI, so it is cached:

  var W = 0.636619772; // 2 / Math.PI

  function bessel0(x) {
    var a, a1, a2, y = x * x, xx = M.abs(x) - 0.785398164;

For small x, the direct Laurent approximation gives better results.

    if(M.abs(x) < 8) {
      a1 = _horner(b0_a1a, y);
      a2 = _horner(b0_a2a, y);
      a = a1/a2;
    }

For larger x, the Chebyshev approach is taken:

    else {
      y = 64 / y;
      a1 = _horner(b0_a1b, y);
      a2 = _horner(b0_a2b, y);
      a = M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x));
    }
    return a;
  }

A similar approach is taken for the first-order bessel function

  var b1_a1a = [72362614232.0,-7895059235.0,242396853.1,-2972611.439, 15704.48260, -30.16036606].reverse();
  var b1_a2a = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse();
  var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse();
  var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse();
  function bessel1(x) {
    var a, a1, a2, y = x*x, xx = M.abs(x) - 2.356194491;
    if(Math.abs(x)< 8) {
      a1 = x*_horner(b1_a1a, y);
      a2 = _horner(b1_a2a, y);
      a = a1 / a2;
    } else {
      y = 64 / y;
      a1=_horner(b1_a1b, y);
      a2=_horner(b1_a2b, y);
      a=M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x));
      if(x < 0) a = -a;
    }
    return a;
  }

For large values of x, the aforementioned iteration is fine, but for small values the expressions quickly blow up. Hence a more careful iteration is used:

  return function besselj(x/*:number*/, n/*:number*/) {
    n = Math.round(n);
    if(n === 0) return bessel0(M.abs(x));
    if(n === 1) return bessel1(M.abs(x));
    if(n < 0) throw 'BESSELJ: Order (' + n + ') must be nonnegative';
    if(M.abs(x) === 0) return 0;
  
    var ret, j, tox = 2 / M.abs(x), m, jsum, sum, bjp, bj, bjm;
    if(M.abs(x) > n) {
      ret = _bessel_iter(x, n, bessel0(M.abs(x)), bessel1(M.abs(x)),-1);
    } else {
      m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2);
      jsum=0;
      bjp=ret=sum=0.0;
      bj=1.0;
      for (j=m;j>0;j--) {
        bjm=j*tox*bj-bjp;
        bjp=bj;
        bj=bjm;
        if (M.abs(bj) > 1E10) {
          bj *= 1E-10;
          bjp *= 1E-10;
          ret *= 1E-10;
          sum *= 1E-10;
        }
        if (jsum) sum += bj;
        jsum=!jsum;
        if (j == n) ret=bjp;
      }
      sum=2.0*sum-bj;
      ret /= sum;
    }
    return x < 0 && (n%2) ? -ret : ret;
  };
})();

The second kind function Y is a bit less finicky:

var bessely = (function() {
  var b0_a1a = [-2957821389.0, 7062834065.0, -512359803.6, 10879881.29, -86327.92757, 228.4622733].reverse();
  var b0_a2a = [40076544269.0, 745249964.8, 7189466.438, 47447.26470, 226.1030244, 1.0].reverse();
  var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse();
  var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934945152e-7].reverse();
  
  var W = 0.636619772;
  function bessel0(x) {
    var a, a1, a2, y = x * x, xx = x - 0.785398164;
    if(x < 8) {
      a1 = _horner(b0_a1a, y);
      a2 = _horner(b0_a2a, y);
      a = a1/a2 + W * besselj(x,0) * M.log(x);
    } else {
      y = 64 / y;
      a1 = _horner(b0_a1b, y);
      a2 = _horner(b0_a2b, y);
      a = M.sqrt(W/x)*(M.sin(xx)*a1+M.cos(xx)*a2*8/x);
    }
    return a;
  }

  var b1_a1a = [-0.4900604943e13, 0.1275274390e13, -0.5153438139e11, 0.7349264551e9, -0.4237922726e7, 0.8511937935e4].reverse();
  var b1_a2a = [0.2499580570e14, 0.4244419664e12, 0.3733650367e10, 0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1].reverse();
  var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse();
  var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse();
  function bessel1(x) {
    var a, a1, a2, y = x*x, xx = x - 2.356194491;
    if(x < 8) {
      a1 = x*_horner(b1_a1a, y);
      a2 = _horner(b1_a2a, y);
      a = a1/a2 + W * (besselj(x,1) * M.log(x) - 1 / x);
    } else {
      y = 64 / y;
      a1=_horner(b1_a1b, y);
      a2=_horner(b1_a2b, y);
      a=M.sqrt(W/x)*(M.sin(xx)*a1+M.cos(xx)*a2*8/x);
    }
    return a;
  }

  return _bessel_wrap(bessel0, bessel1, 'BESSELY', 1, -1);
})();

And the modified Bessel functions are even easier:

var besseli = (function() {
  var b0_a = [1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732, 0.360768e-1, 0.45813e-2].reverse();
  var b0_b = [0.39894228, 0.1328592e-1, 0.225319e-2, -0.157565e-2, 0.916281e-2, -0.2057706e-1, 0.2635537e-1, -0.1647633e-1, 0.392377e-2].reverse();
  function bessel0(x) {
    if(x <= 3.75) return _horner(b0_a, x*x/(3.75*3.75));
    return M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b0_b, 3.75/M.abs(x));
  }

  var b1_a = [0.5, 0.87890594, 0.51498869, 0.15084934, 0.2658733e-1, 0.301532e-2, 0.32411e-3].reverse();
  var b1_b = [0.39894228, -0.3988024e-1, -0.362018e-2, 0.163801e-2, -0.1031555e-1, 0.2282967e-1, -0.2895312e-1, 0.1787654e-1, -0.420059e-2].reverse();
  function bessel1(x) {
    if(x < 3.75) return x * _horner(b1_a, x*x/(3.75*3.75));
    return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b1_b, 3.75/M.abs(x));
  }

  return function besseli(x/*:number*/, n/*:number*/) {
    n = Math.round(n);
    if(n === 0) return bessel0(x);
    if(n == 1) return bessel1(x);
    if(n < 0) throw 'BESSELI Order (' + n + ') must be nonnegative';
    if(M.abs(x) === 0) return 0;
  
    var ret, j, tox = 2 / M.abs(x), m, bip, bi, bim;
    m=2*M.round((n+M.round(M.sqrt(40*n)))/2);
    bip=ret=0.0;
    bi=1.0;
    for (j=m;j>0;j--) {
      bim=j*tox*bi + bip;
      bip=bi; bi=bim;
      if (M.abs(bi) > 1E10) {
        bi *= 1E-10;
        bip *= 1E-10;
        ret *= 1E-10;
      }
      if(j == n) ret = bip;
    }
    ret *= besseli(x, 0) / bi;
    return x < 0 && (n%2) ? -ret : ret;
  };
  
})();

var besselk = (function() {
  var b0_a = [-0.57721566, 0.42278420, 0.23069756, 0.3488590e-1, 0.262698e-2, 0.10750e-3, 0.74e-5].reverse();
  var b0_b = [1.25331414, -0.7832358e-1, 0.2189568e-1, -0.1062446e-1, 0.587872e-2, -0.251540e-2, 0.53208e-3].reverse();
  function bessel0(x) {
    if(x <= 2) return -M.log(x/2)*besseli(x,0) + _horner(b0_a, x*x/4);
    return M.exp(-x)/M.sqrt(x)*_horner(b0_b, 2/x);
  }

  var b1_a = [1.0, 0.15443144, -0.67278579, -0.18156897, -0.1919402e-1, -0.110404e-2, -0.4686e-4].reverse();
  var b1_b = [1.25331414, 0.23498619, -0.3655620e-1, 0.1504268e-1, -0.780353e-2, 0.325614e-2, -0.68245e-3].reverse();
  function bessel1(x) {
    if(x <= 2) return M.log(x/2)*besseli(x,1) + (1/x)*_horner(b1_a, x*x/4);
    return M.exp(-x)/M.sqrt(x)*_horner(b1_b, 2/x);
  }

  return _bessel_wrap(bessel0, bessel1, 'BESSELK', 2, 1);
})();

Export Magic

Now we need to export:

if(typeof exports !== "undefined") {
  exports.besselj = besselj;
  exports.bessely = bessely;
  exports.besseli = besseli;
  exports.besselk = besselk;
}

Testing

Fortunately, these functions are available in Excel, so you can test there (or use the same functions in Mathematica). Note that this uses Excel semantics, which are reverse from Mathematica:

BesselJ[n, x] (* Mathematica *)