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:
bessel(order, x) = sum(n, 0, infinity, ((-1) ** n / (factorial(n) * gamma(n + order + 1))) * (x / 2) ** (2 * n + order))
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.
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))
(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))))
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)
extracted and translated from SheetJs/bessel
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