Clean up besselj(nu,x) and add docs

Project.toml

@@ -3,11 +3,8 @@ uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 authors = ["Michael Helton <[email protected]> and contributors"]
 version = "0.1.0"
 
-[deps]
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-
 [compat]
-julia = "1.5"
+julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/Bessels.jl

@@ -1,7 +1,5 @@
 module Bessels
 
-using SpecialFunctions: loggamma
-
 export besselj0
 export besselj1
 export besselj

src/U_polynomials.jl

@@ -1,5 +1,19 @@
-besseljy_debye_cutoff(nu, x) = nu > 2.0 + 1.00035*x + Base.Math._approx_cbrt(Float64(302.681)*x) && x > 15
-# valid when x < v (uniform asymptotic expansions)
+#                 Debye asymptotic expansions
+#                `besseljy_debye`, `hankel_debye`
+#
+#  This file contains the debye asymptotic asymptotic expansions for large orders.
+#  These routines can be used to calculate `besselj`, `bessely`, `besselk`, `besseli`
+#  and the Hankel functions. These are uniform expansions for `besselk` and `besseli` for large orders.
+#  The forms used for `besselj` and `bessely` as well as the Hankel functions follow the notation by Matviyenko [1]
+#  but also see similar forms provided by NIST 10.19. All of these routines use a core routine for calculating
+#  the U-polynomials (NIST 10.41.E9) where the coefficients are dependent on each function. The forms for the modified
+#  Bessel functions `besselk` and `besseli` follow NIST 10.41 [3].
+#
+# [1] Matviyenko, Gregory. "On the evaluation of Bessel functions." 
+#     Applied and Computational Harmonic Analysis 1.1 (1993): 116-135.
+# [2] http://dlmf.nist.gov/10.41.E9
+# [3] https://dlmf.nist.gov/10.41
+
 """
     besseljy_debye(nu, x::T)
 
@@ -26,9 +40,9 @@ function besseljy_debye(v, x)
 
     return coef_Jn * Uk_Jn, coef_Yn * Uk_Yn
 end
-hankel_debye_cutoff(nu, x) = nu < 0.2 + x + Base.Math._approx_cbrt(-411*x)
 
-# valid when v < x (uniform asymptotic expansions)
+besseljy_debye_cutoff(nu, x) = nu > 2.0 + 1.00035*x + Base.Math._approx_cbrt(Float64(302.681)*x) && x > 15
+
 """
     hankel_debye(nu, x::T)
 
@@ -54,6 +68,8 @@ function hankel_debye(v, x::T) where T
     return coef_Yn * Uk_Yn
 end
 
+hankel_debye_cutoff(nu, x) = nu < 0.2 + x + Base.Math._approx_cbrt(-411*x)
+
 function Uk_poly_Jn(p, v, p2, x::T) where T <: Float64
     if v > 5.0 + 1.00033*x + Base.Math._approx_cbrt(1427.61*x)
         return Uk_poly10(p, v, p2)

src/asymptotics.jl

@@ -1,3 +1,13 @@
+#             Asymptotics expansions for x > nu
+#                `besseljy_large_argument`
+#
+#  This file contains the asymptotic asymptotic expansions for large arguments 
+#  which is accurate in the fresnel regime |x| > nu following the derivations given by Heitman [1].
+#  These routines can be used to calculate `besselj` and `bessely` using phase functions.
+#
+# [1] Heitman, Z., Bremer, J., Rokhlin, V., & Vioreanu, B. (2015). On the asymptotics of Bessel functions in the Fresnel regime. 
+#     Applied and Computational Harmonic Analysis, 39(2), 347-356.
+
 #besseljy_large_argument_min(::Type{Float32}) = 15.0f0
 besseljy_large_argument_min(::Type{Float64}) = 20.0
 besseljy_large_argument_min(::Type{T}) where T <: AbstractFloat = 40.0
@@ -6,7 +16,12 @@ besseljy_large_argument_min(::Type{T}) where T <: AbstractFloat = 40.0
 besseljy_large_argument_cutoff(v, x::Float64) = (x > 1.65*v && x > besseljy_large_argument_min(Float64))
 besseljy_large_argument_cutoff(v, x::T) where T = (x > 4*v && x > besseljy_large_argument_min(T))
 
+"""
+    besseljy_large_argument(nu, x::T)
 
+Asymptotic expansions for large arguments valid when x > 1.6*nu and x > 20.0.
+Returns both (besselj(nu, x), bessely(nu, x)).
+"""
 function besseljy_large_argument(v, x::T) where T
     # gives both (besselj, bessely) for x > 1.6*v
     α, αp = _α_αp_asymptotic(v, x)
@@ -119,7 +134,6 @@ function _α_αp_poly_10(v, x::T) where T
     αp = evalpoly(xinv, (s0, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10))
     α = x * evalpoly(xinv, (s0, -s1, -s2/3, -s3/5, -s4/7, -s5/9, -s6/11, -s7/13, -s8/15, -s9/17, -s10/19))
     return α, αp
-    return α, αp
 end
 #=
 function _α_αp_poly_15(v, x::T) where T

src/besselj.jl

@@ -151,6 +151,47 @@ function besselj1(x::Float32)
     end
 end
 
+#                  Bessel functions of the first kind of order nu
+#                               besselj(nu, x)
+#
+#    A numerical routine to compute the Bessel function of the first kind J_{ν}(x) [1]
+#    for real orders and arguments of positive or negative value. The routine is based on several
+#    publications [2, 3, 4, 5] that calculate J_{ν}(x) for positive arguments and orders where
+#    reflection identities are used to compute negative arguments and orders.
+#
+#    In particular, the reflectance identities for negative integer orders J_{-n}(x) = (-1)^n * J_{n}(x) (Eq. 9.1.5; [6])
+#    and for negative noninteger orders J_{−ν}(x) = cos(πν) * J_{ν}(x) - sin(πν) * Y_{ν}(x) are used.
+#    For negative arguments of integer order, J_{n}(-x) = (-1)^n * J_{n}(x) is used and for
+#    noninteger orders, J_{ν}(−x) = exp(im*π*ν) * J_{ν}(x) is used. For negative orders and arguments the previous identities are combined.
+#
+#    The identities are computed by calling the `besselj_positive_args(nu, x)` function which computes J_{ν}(x)
+#    for positive arguments and orders. When x < ν and ν being reasonably large, the debye asymptotic expansion (Eq. 32; [3]) is used `besseljy_debye(nu, x)`.
+#    For large arguments x >> ν, the phase functions are used (Eq. 14 [4]) with `besseljy_large_argument(nu, x)`.
+#    When x > ν and x being reasonably large, the Hankel function is calculated from the debye expansion (Eq. 29; [3]) with `hankel_debye(nu, x)`
+#    and J_{n}(x) is calculated from the real part of the Hankel function. These expansions are not uniform so are not strictly used when the above inequalities are true, therefore, cutoffs
+#    were determined depending on the desired accuracy. For large arguments x >> ν, the phase functions are used (Eq. 15 [4]) with `besseljy_large_argument(nu, x)`.
+#    For small arguments, J_{ν}(x) is calculated from the power series (`bessely_power_series(nu, x`)
+#
+#    For values where the expansions for large arguments and orders are not valid, backward recurrence is employed after shifting the order up
+#    to where `besseljy_debye` is accurate then using downward recurrence. In general, the routine will be the slowest when ν ≈ x as all methods struggle at this point.
+#    
+# [1] http://dlmf.nist.gov/10.2.E2
+# [2] Bremer, James. "An algorithm for the rapid numerical evaluation of Bessel functions of real orders and arguments." 
+#     Advances in Computational Mathematics 45.1 (2019): 173-211.
+# [3] Matviyenko, Gregory. "On the evaluation of Bessel functions." 
+#     Applied and Computational Harmonic Analysis 1.1 (1993): 116-135.
+# [4] Heitman, Z., Bremer, J., Rokhlin, V., & Vioreanu, B. (2015). On the asymptotics of Bessel functions in the Fresnel regime. 
+#     Applied and Computational Harmonic Analysis, 39(2), 347-356.
+# [5] Ratis, Yu L., and P. Fernández de Córdoba. "A code to calculate (high order) Bessel functions based on the continued fractions method." 
+#     Computer physics communications 76.3 (1993): 381-388.
+# [6] Abramowitz, Milton, and Irene A. Stegun, eds. Handbook of mathematical functions with formulas, graphs, and mathematical tables. 
+#     Vol. 55. US Government printing office, 1964.
+#
+
+#####
+##### Generic routine for `besselj`
+#####
+
 function besselj(nu::Real, x::T) where T
     isinteger(nu) && return besselj(Int(nu), x)
     abs_nu = abs(nu)
@@ -186,59 +227,65 @@ function besselj(nu::Integer, x::T) where T
     end
 end
 
+"""
+    besselj_positive_args(nu, x::T) where T <: Float64
+
+Bessel function of the first kind of order nu, ``J_{nu}(x)``.
+nu and x must be real and nu and x must be positive.
+
+No checks on arguments are performed and should only be called if certain nu, x >= 0.
+"""
 function besselj_positive_args(nu::Real, x::T) where T
     nu == 0 && return besselj0(x)
     nu == 1 && return besselj1(x)
 
-    x < 4.0 && return besselj_small_arguments_orders(nu, x)
+    # x < ~nu branch see src/U_polynomials.jl
+    besseljy_debye_cutoff(nu, x) && return besseljy_debye(nu, x)[1]
 
-    large_arg_cutoff = 1.65*nu
-    (x > large_arg_cutoff && x > 20.0) && return besseljy_large_argument(nu, x)[1]
+    # large argument branch see src/asymptotics.jl
+    besseljy_large_argument_cutoff(nu, x) && return besseljy_large_argument(nu, x)[1]
 
+    # x > ~nu branch see src/U_polynomials.jl on computing Hankel function
+    hankel_debye_cutoff(nu, x) && return real(hankel_debye(nu, x))
 
-    debye_cutoff = 2.0 + 1.00035*x + Base.Math._approx_cbrt(302.681*Float64(x))
-    nu > debye_cutoff && return besseljy_debye(nu, x)[1]
-
-    if nu >= x
-        nu_shift = ceil(Int, debye_cutoff - nu)
-        v = nu + nu_shift
-        jnu = besseljy_debye(v, x)[1]
-        jnup1 = besseljy_debye(v+1, x)[1]
-        return besselj_down_recurrence(x, jnu, jnup1, v, nu)[1]
-    end
+    # use power series for small x and for when nu > x
+    besselj_series_cutoff(nu, x) && return besselj_power_series(nu, x)
 
-    # at this point x > nu and  x < nu * 1.65
-    # in this region forward recurrence is stable
-    # we must decide if we should do backward recurrence if we are closer to debye accuracy
-    # or if we should do forward recurrence if we are closer to large argument expansion
-    debye_cutoff = 5.0 + 1.00033*x + Base.Math._approx_cbrt(1427.61*Float64(x))
-
-    debye_diff = debye_cutoff - nu
-    large_arg_diff = nu - x / 2.0
-
-    if (debye_diff > large_arg_diff && x > 20.0)
-        nu_shift = ceil(Int, large_arg_diff)
-        v2 = nu - nu_shift
-        jnu = besseljy_large_argument(v2, x)[1]
-        jnum1 = besseljy_large_argument(v2 - 1, x)[1]
-        return besselj_up_recurrence(x, jnu, jnum1, v2, nu)[1]
-    else
-        nu_shift = ceil(Int, debye_diff)
-        v = nu + nu_shift
-        jnu = besseljy_debye(v, x)[1]
-        jnup1 = besseljy_debye(v+1, x)[1]
-        return besselj_down_recurrence(x, jnu, jnup1, v, nu)[1]
-    end
+    # At this point we must fill the region when x ≈ v with recurrence
+    # Backward recurrence is always stable and forward recurrence is stable when x > nu
+    # However, we only use backward recurrence by shifting the order up and using `besseljy_debye` to generate start values
+    # Both `besseljy_debye` and `hankel_debye` get more accurate for large orders,
+    # however `besseljy_debye` is slightly more efficient (no complex variables) and we need no branches if only consider one direction.
+    # On the other hand, shifting the order down avoids any concern about underflow for large orders
+    # Shifting the order too high while keeping x fixed could result in numerical underflow
+    # Therefore we need to shift up only until the `besseljy_debye` is accurate and need to test that no underflow occurs
+    # Shifting the order up decreases the value substantially for high orders and results in a stable forward recurrence
+    # as the values rapidly increase
+
+    debye_cutoff = 2.0 + 1.00035*x + Base.Math._approx_cbrt(302.681*Float64(x))
+    nu_shift = ceil(Int, debye_cutoff - nu)
+    v = nu + nu_shift
+    jnu = besseljy_debye(v, x)[1]
+    jnup1 = besseljy_debye(v+1, x)[1]
+    return besselj_down_recurrence(x, jnu, jnup1, v, nu)[1]
 end
 
-# generally can only use for x < 4.0
-# this needs a better way to sum these as it produces large errors
+#####
+##### Power series for J_{nu}(x)
+#####
+
+# accurate for x < 7.0 or nu > 2+ 0.109x + 0.062x^2 for Float64
+# accurate for x < 20.0 or nu > 14.4 - 0.455x + 0.027x^2 for Float32 (when using F64 precision)
 # only valid in non-oscillatory regime (v>1/2, 0<t<sqrt(v^2 - 0.25))
-# power series has premature underflow for large orders
-function besselj_small_arguments_orders(v, x::T) where T
-    v > 60 && return log_besselj_small_arguments_orders(v, x)
+# power series has premature underflow for large orders though use besseljy_debye for large orders
+"""
+    besselj_power_series(nu, x::T) where T <: Float64
 
-    MaxIter = 2000
+Computes ``J_{nu}(x)`` using the power series.
+In general, this is most accurate for small arguments and when nu > x.
+"""
+function besselj_power_series(v, x::T) where T
+    MaxIter = 3000
     out = zero(T)
     a = (x/2)^v / gamma(v + one(T))
     t2 = (x/2)^2
@@ -250,11 +297,16 @@ function besselj_small_arguments_orders(v, x::T) where T
     return out
 end
 
+besselj_series_cutoff(v, x::Float64) = (x < 7.0) || v > (2 + x*(0.109 + 0.062x))
+besselj_series_cutoff(v, x::Float32) = (x < 20.0) || v > (14.4 + x*(-0.455 + 0.027x))
+
+#=
 # this needs a better way to sum these as it produces large errors
 # use when v is large and x is small
-# need for bessely 
+# though when v is large we should use the debye expansion instead
+# also do not have a julia implementation of loggamma so will not use for now
 function log_besselj_small_arguments_orders(v, x::T) where T
-    MaxIter = 2000
+    MaxIter = 3000
     out = zero(T)
     a = one(T)
     x2 = (x/2)^2
@@ -266,3 +318,4 @@ function log_besselj_small_arguments_orders(v, x::T) where T
     logout = -loggamma(v + 1) + fma(v, log(x/2), log(out))
     return exp(logout)
 end
+=#

src/bessely.jl

@@ -190,7 +190,7 @@ end
 #    reflection identities are used to compute negative arguments and orders.
 #
 #    In particular, the reflectance identities for negative integer orders Y_{-n}(x) = (-1)^n * Y_{n}(x) (Eq. 9.1.5; [6])
-#    and for negative noninteger orders Y_{−ν}(x) = cos(πν) * Y_{ν}(x) + sin(πν) * J_{ν}(x) (Eq. 9.1.2; [6]) are used.
+#    and for negative noninteger orders Y_{−ν}(x) = cos(πν) * Y_{ν}(x) + sin(πν) * J_{ν}(x) are used.
 #    For negative arguments of integer order, Y_{n}(-x) = (-1)^n * Y_{n}(x) + (-1)^n * 2im * J_{n}(x) is used and for
 #    noninteger orders, Y_{ν}(−x) = exp(−im*π*ν) * Y_{ν}(x) + 2im * cos(πν) * J_{ν}(x) is used.
 #    For negative orders and arguments the previous identities are combined.
@@ -277,7 +277,7 @@ end
 """
     bessely_positive_args(nu, x::T) where T <: Float64
 
-Bessel function of the first kind of order nu, ``Y_{nu}(x)``.
+Bessel function of the second kind of order nu, ``Y_{nu}(x)``.
 nu and x must be real and nu and x must be positive.
 
 No checks on arguments are performed and should only be called if certain nu, x >= 0.

test/besselj_test.jl

@@ -109,23 +109,23 @@ end
 @test isapprox(besselj(10.0, 150.0), SpecialFunctions.besselj(10.0, 150.0), rtol=1e-12)
 
 # test BigFloat for single point
-@test isapprox(besselj(big"2000", big"1500.0"), SpecialFunctions.besselj(big"2000", big"1500"), rtol=5e-20)
-@test isapprox(besselj(big"20", big"1500.0"), SpecialFunctions.besselj(big"20", big"1500"), rtol=5e-20)
+@test isapprox(Bessels.besseljy_debye(big"2000", big"1500.0")[1], SpecialFunctions.besselj(big"2000", big"1500"), rtol=5e-20)
+@test isapprox(Bessels.besseljy_large_argument(big"20", big"1500.0")[1], SpecialFunctions.besselj(big"20", big"1500"), rtol=5e-20)
 
 
 # need to test accuracy of negative orders and negative arguments and all combinations within
 # SpecialFunctions.jl doesn't provide these so will hand check against hard values
 # values taken from https://keisan.casio.com/exec/system/1180573474 which match mathematica
 # need to also account for different branches when nu isa integer
 nu = -9.102; x = -12.48
-@test isapprox(besselj(nu, x), 0.09842356047575545808128 -0.03266486217437818486161im, rtol=1e-14)
+@test isapprox(besselj(nu, x), 0.09842356047575545808128 -0.03266486217437818486161im, rtol=8e-14)
 nu = -5.0; x = -5.1
-@test isapprox(besselj(nu, x), 0.2740038554704588327387, rtol=1e-14)
+@test isapprox(besselj(nu, x), 0.2740038554704588327387, rtol=8e-14)
 nu = -7.3; x = 19.1
-@test isapprox(besselj(nu, x), 0.1848055978553359009813, rtol=1e-14)
+@test isapprox(besselj(nu, x), 0.1848055978553359009813, rtol=8e-14)
 nu = -14.0; x = 21.3
-@test isapprox(besselj(nu, x), -0.1962844898264965120021, rtol=1e-14)
+@test isapprox(besselj(nu, x), -0.1962844898264965120021, rtol=8e-14)
 nu = 13.0; x = -8.5
-@test isapprox(besselj(nu, x), -0.006128034621313167000171, rtol=1e-14)
+@test isapprox(besselj(nu, x), -0.006128034621313167000171, rtol=8e-14)
 nu = 17.45; x = -16.23
-@test isapprox(besselj(nu, x), -0.01607335977752705869797 -0.1014831996412783806255im, rtol=1e-14)
+@test isapprox(besselj(nu, x), -0.01607335977752705869797 -0.1014831996412783806255im, rtol=8e-14)