Numerical routines for computing Bessel, Airy, and Hankel functions for real arguments. These routines are written in the Julia programming language and are self contained without any external dependencies.

The goal of the library is to provide high quality numerical implementations of Bessel functions with high accuracy without comprimising on computational time. In general, we try to match (and often exceed) the accuracy of other open source routines such as those provided by SpecialFunctions.jl. There are instances where we don't quite match that desired accuracy (within a digit or two) but in general will provide implementations that are 5-10x faster (see benchmarks).

The library currently supports Bessel functions, modified Bessel functions, Hankel functions, spherical Bessel functions, and Airy functions of the first and second kind for positive real arguments and integer and noninteger orders. Negative arguments are also supported only if the return value is real. Limited support is provided for complex arguments. An unexported gamma function is also provided.

# add the package
pkg> add https://github.com/JuliaMath/Bessels.jl.git

julia> using Bessels

julia> x = 12.3; nu = 1.3

julia> besselj(nu, x)
-0.2267581644816917

$$ J_{\nu}(x) = \sum_{m=0}^{\infty} \frac{(-1)^m}{m!\Gamma(m+\nu+1)}(\frac{x}{2})^{2m+\nu} $$

Bessel functions of the first kind, denoted as $J_{\nu}(x)$, can be called with besselj(nu, x) where nu is the order of the Bessel function with argument x. Routines are also available for orders 0 and 1 which can be called with besselj0(x) and besselj1(x).

julia> ν, x = 1.4, 12.3

# generic call for any order ν
julia> besselj(ν, x)
-0.22796228516266345

# ν = 0
julia> besselj0(x)
0.11079795030758544

# ν = 1
julia> besselj1(x)
-0.1942588480405914

$$ Y_{\nu}(x) = \frac{J_{\nu}(x) \cos(\nu \pi) - J_{-\nu}(x)}{\sin(\nu \pi)} $$

Bessel functions of the second kind, denoted as $Y_{\nu}(x)$, can be called with bessely(nu, x). Routines are also available for orders 0 and 1 which can be called with bessely0(x) and bessely1(x).

julia> ν, x = 1.4, 12.3

# generic call for any order ν
julia> bessely(ν, x)
0.00911009829832235

# ν = 0
julia> bessely0(x)
-0.19859309463502633

# ν = 1
julia> bessely1(x)
-0.11894840329926631

$$ I_{\nu}(x) = \sum_{m=0}^{\infty} \frac{1}{m!\Gamma(m+\nu+1)}(\frac{x}{2})^{2m+\nu} $$

Modified Bessel functions of the first kind, denoted as $I_{\nu}(x)$, can be called with besseli(nu, x) where nu is the order of the Bessel function with argument x. Routines are also available for orders 0 and 1 which can be called with besseli0(x) and besseli1(x). Exponentially scaled versions of these functions $I_{\nu}(x) \cdot e^{-x}$ are also provided which can be called with besseli0x(x), besseli1x(x), and besselix(nu, x).

julia> ν, x = 1.4, 12.3

# generic call for any order v
julia> besseli(ν, x)
23242.698263113296

# exponentially scaled version
julia> besselix(ν, x)
0.10579482312624018

# ν = 0
julia> besseli0(x)
25257.48759692308
julia> besseli0x(x)
0.11496562932068803

# ν = 1
julia> besseli1(x)
24207.933018435186
julia> besseli1x(x)
0.11018832507935208

$$ K_{\nu}(x) = \frac{\pi}{2} \frac{I_{-\nu}(x) - I_{\nu}(x)}{\sin(\nu \pi)} $$

Modified Bessel functions of the second kind, denoted as $K_{\nu}(x)$, can be called with besselk(nu, x). Routines are available for orders 0 and 1 which can be called with besselk0(x) and besselk1(x). Exponentially scaled versions of these functions $K_{\nu}(x) \cdot e^{x}$ are also provided which can be called with besselk0x(x), besselk1x(x), and besselkx(nu, x).

julia> ν, x = 1.4, 12.3

julia> besselk(ν, x)
1.739055243080153e-6

julia> besselk0(x)
1.6107849768886856e-6

julia> besselk1(x)
1.6750295538365835e-6

We also provide support for besselj(nu::M, x::T), bessely(nu::M, x::T), besseli(nu::M, x::T), besselk(nu::M, x::T), besseli(nu::M, x::T), besselh(nu::M, k, x::T) when M is some AbstractRange and T is some float.

julia> besselj(0:10, 1.0)
11-element Vector{Float64}:
 0.7651976865579666
 0.44005058574493355
 0.11490348493190049
 0.019563353982668407
 0.0024766389641099553
 0.00024975773021123444
 2.0938338002389273e-5
 1.5023258174368085e-6
 9.422344172604502e-8
 5.249250179911876e-9
 2.630615123687453e-10

In general, this provides a fast way to generate a sequence of Bessel functions for many orders.

julia> @btime besselj(0:100, 50.0)
  398.095 ns (1 allocation: 896 bytes)

This function will allocate so it is recommended that you calculate the Bessel functions at the top level of your function outside any hot loop. You can also call the mutating function on your preallocated vector Bessels.besselj!(out, nu, x)

a = zeros(10)
out = Bessels.besselj!(a, 1:10, 1.0)

Support for complex numbers is only provided for the Airy functions (airyai, airyaiprime, airybi, airybiprime) and the Bessel functions of the first kind with orders 0 and 1 (besselj0, besselj1, besseli0, besseli1).

Support is provided for negative arguments and orders only if the return value is real. A domain error will be thrown if the return value is complex. See #30 for more details.

julia> ν, x = 13.0, -1.0
julia> besseli(ν, x)
-1.9956316782072005e-14

julia> ν, x = -14.0, -9.9
julia> besseli(ν, x)
0.2892290867115618

julia> ν, x = 12.6, -3.0
julia> besseli(ν, x)
ERROR: DomainError with -3.0:
Complex result returned for real arguments. Complex arguments are currently not supported
Stacktrace:
 [1] _besseli(nu::Float64, x::Float64)
   @ Bessels ~/.julia/packages/Bessels/OBoYU/src/besseli.jl:181
 [2] besseli(nu::Float64, x::Float64)
   @ Bessels ~/.julia/packages/Bessels/OBoYU/src/besseli.jl:167
 [3] top-level scope
   @ REPL[62]:1

We also provide an unexported gamma function for real arguments that can be called with Bessels.gamma(x).

We report the relative errors (abs(1 - Bessels.f(x)/ArbNumerics.f(ArbFloat(x)))) compared to ArbNumerics.jl when computed in a higher precision. The working precision was set to setworkingprecision(ArbFloat, 500); setextrabits(128) for the calculations in arbitrary precision. We generate a thousand random points for $x \in (0, 100)$ and compute the mean and maximum absolute relative errors.

In general the largest relative errors are observed near the zeros of Bessel functions for besselj and bessely. Accuracy might also be slightly worse for very large arguments when using Float64 precision.

We give brief performance comparisons to the implementations provided by SpecialFunctions.jl. In general, special functions are computed with separate algorithms in different domains leading to computational time being dependent on argument. For these comparisons we show the relative speed increase for computing random values between 0 and 100 for x and order nu. In some ranges, performance may be significantly better while others will be more similar.

Benchmarks were run using Julia Version 1.7.2 on an Apple M1 using Rosetta.

• besselj0(x)
• besselj1(x)
• besselj(nu, x)
• bessely0(x)
• bessely1(x)
• bessely(nu, x)
• besseli0(x)
• besseli1(x)
• besseli(nu, x)
• besselk0(x)
• besselk1(x)
• besselk(nu, x)
• besselh(nu, k, x)
• hankelh1(nu, x)
• hankelh2(nu, x)
• sphericalbesselj(nu, x)
• sphericalbessely(nu, x)
• Bessels.sphericalbesseli(nu, x)
• Bessels.sphericalbesselk(nu, x)
• airyai(x)
• airyaiprime(x)
• airybi(x)
• airybiprime(x)
• Bessels.gamma(x)

• Support for higher precision Double64, Float128
• Support for complex arguments (x and nu)
• Support for derivatives with respect to argument and order

Contributions are very welcome, as are feature requests, suggestions or general discussions. Please open an issue for discussion on newer implementations, share papers, new features, or if you encounter any problems. Our goal is to provide high quality implementations of Bessel functions that match or exceed the accuracy of the implementations provided by SpecialFunctions.jl. Please let us know if you encounter any accuracy or performance differences.