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. however, 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)!))
where sum relates to the capital-sigma notation and has the parameters variable, from, to, and body.
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)
this version also works well for larger values. it was extracted and translated from SheetJs/bessel.
bessel_j = do -> M = Math W = 2 / M.PI # Coefficients for small and large x approximations for J0(x) and J1(x) coeff_j0_small_x_numerator = [57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456].reverse() coeff_j0_small_x_denominator = [57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0].reverse() coeff_j0_large_x_cosine = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse() coeff_j0_large_x_sine = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse() coeff_j1_small_x_numerator = [72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606].reverse() coeff_j1_small_x_denominator = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse() coeff_j1_large_x_cosine = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse() coeff_j1_large_x_sine = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse() # Horner's method for polynomial evaluation horner = (coefficients, variable) -> result = 0 for coefficient in coefficients result = variable * result + coefficient result # Bessel function recurrence relation bessel_recurrence = (x, n, j0, j1, sign) -> return j0 if n == 0 return j1 if n == 1 two_over_x = 2 / x current_j = j1 for order in [2..n] current_j = j1 * (order - 1) * two_over_x + sign * j0 j0 = j1 j1 = current_j current_j # Bessel function J0(x) computation bessel_j0 = (x) -> y = x * x if x < 8 numerator = horner(coeff_j0_small_x_numerator, y) denominator = horner(coeff_j0_small_x_denominator, y) numerator / denominator else xx = x - 0.785398164 # x - pi/4 y = 64 / y cosine_term = horner(coeff_j0_large_x_cosine, y) sine_term = horner(coeff_j0_large_x_sine, y) M.sqrt(W / x) * (M.cos(xx) * cosine_term - (M.sin(xx) * sine_term * 8 / x)) # Bessel function J1(x) computation bessel_j1 = (x) -> y = x * x if M.abs(x) < 8 numerator = x * horner(coeff_j1_small_x_numerator, y) denominator = horner(coeff_j1_small_x_denominator, y) numerator / denominator else xx = M.abs(x) - 2.356194491 # abs(x) - 3*pi/4 y = 64 / y cosine_term = horner(coeff_j1_large_x_cosine, y) sine_term = horner(coeff_j1_large_x_sine, y) result = M.sqrt(W / M.abs(x)) * (M.cos(xx) * cosine_term - (M.sin(xx) * sine_term * 8 / M.abs(x))) result = -result if x < 0 result # Main Bessel function J_n(x) (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) * bessel_j(x, -n) if x < 0 return (if n % 2 then -1 else 1) * bessel_j(-x, n) return bessel_j0(x) if n == 0 return bessel_j1(x) if n == 1 return 0 if x == 0 if x > n bessel_recurrence(x, n, bessel_j0(x), bessel_j1(x), -1) else m = 2 * M.floor((n + M.floor(M.sqrt(40 * n))) / 2) sum_sign = false previous_j = 0.0 sum_j = 0.0 current_j = 1.0 previous_j_minus_1 = 0.0 two_over_x = 2 / x for j in [m..1] previous_j_minus_1 = j * two_over_x * current_j - previous_j previous_j = current_j current_j = previous_j_minus_1 if M.abs(current_j) > 1e10 current_j *= 1e-10 previous_j *= 1e-10 sum_j *= 1e-10 if sum_sign sum_j += current_j sum_sign = !sum_sign if j == n previous_j = previous_j_minus_1 sum_j = 2.0 * sum_j - current_j previous_j / sum_j