2023-01-05

bessel function of the first kind

wikipedia: bessel function

there are multiple kinds of bessel functions that are grouped together just because they are similar. they are called first kind, second kind, and so on.

x values for which the bessel function of the first kind crosses zero:

  • order 0: 2.404825557695773, 5.520078110286311, 8.653727912911013
  • order 1: 3.8317059702075125, 7.015586669815619, 10.173468135062722
  • order 2: 5.135622301840682556301, 8.417244140399864857784, 11.61984117214905942709
  • order 3: 6.380161895923983506237, 9.761023129981669678545, 13.01520072169843441983

computation

formula

bessel(order, x) = sum(n, 0, infinity,
  ((-1) ** n / (factorial(n) * gamma(n + order + 1)))
  * (x / 2) ** (2 * n + order))
  • for integers, the gamma function can be replaced with the factorial
  • the most direct translation to programming code quickly reaches the limits of numerical representation because of the exponential (2 * n + order)

version for small values

this version is simple and compact, but in software it only gives accurate values for about x < 36 and then starts to rapidly move towards infinity.

formula

bessel(order, x) = (x / 2) ** order
  * sum(n, 0, infinity, (-1 ** n * (x ** 2 / 4) ** n) / (n! * (n + order)!))

the ratio between consecutively summed terms is

-1 * (x ** 2 / 4) / (n * (order + n))

source

scheme

(define (factorial n) (if (> 1 n) 1 (* n (factorial (- n 1)))))

(define (bessel order x term-count)
  (* (expt (/ x 2) order)
    (let loop ((n 1) (term (/ 1 (factorial order))))
      (if (< n term-count)
        (+ term (loop (+ 1 n) (* term -1 (/ (/ (* x x) 4) (* n (+ order n))))))
        term))))

coffeescript

factorial = (n) -> if 1 > n then 1 else n * factorial(n - 1)

bessel = (order, x, term_count) ->
  term = 1 / factorial order
  sum = term
  for n in [1...term_count]
    term = term * -1 * (x * x / 4) / (n * (order + n))
    sum += term
  sum * Math.pow(x / 2, order)

general version

extracted and translated from SheetJs/bessel

coffeescript

besselj = do ->
  M = Math
  W = 2 / M.PI

  b0_a1a = [57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456].reverse()
  b0_a2a = [57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0].reverse()
  b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse()
  b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse()
  b1_a1a = [72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606].reverse()
  b1_a2a = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse()
  b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse()
  b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse()

  horner = (arr, v) ->
    i = 0
    z = 0
    while i < arr.length
      z = v * z + arr[i]
      i += 1
    z

  bessel_iter = (x, n, f0, f1, sign) ->
    return f0 if n == 0
    return f1 if n == 1
    tdx = 2 / x
    f2 = f1
    o = 1
    while o < n
      f2 = f1 * o * tdx + sign * f0
      f0 = f1
      f1 = f2
      o += 1
    f2

  bessel0 = (x) ->
    a = 0
    a1 = 0
    a2 = 0
    y = x * x
    if x < 8
      a1 = horner(b0_a1a, y)
      a2 = horner(b0_a2a, y)
      a = a1 / a2
    else
      xx = x - 0.785398164
      y = 64 / y
      a1 = horner(b0_a1b, y)
      a2 = horner(b0_a2b, y)
      a = M.sqrt(W / x) * (M.cos(xx) * a1 - (M.sin(xx) * a2 * 8 / x))
    a

  bessel1 = (x) ->
    a = 0
    a1 = 0
    a2 = 0
    y = x * x
    xx = M.abs(x) - 2.356194491
    if M.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
    a

  (n, x) ->
    n = M.round(n)
    if !isFinite(x)
      return if isNaN(x) then x else 0
    if n < 0
      return (if n % 2 then -1 else 1) * besselj(x, -n)
    if x < 0
      return (if n % 2 then -1 else 1) * besselj(-x, n)
    return bessel0(x) if n == 0
    return bessel1(x) if n == 1
    return 0 if x == 0
    ret = 0.0
    if x > n
      ret = bessel_iter(x, n, bessel0(x), bessel1(x), -1)
    else
      m = 2 * M.floor((n + M.floor(M.sqrt(40 * n))) / 2)
      jsum = false
      bjp = 0.0
      sum = 0.0
      bj = 1.0
      bjm = 0.0
      tox = 2 / x
      j = m
      while j > 0
        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
        j -= 1
      sum = 2.0 * sum - bj
      ret /= sum
    ret