From 0bc27e35c51be463d3317e30d86651153647e62c Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Tue, 2 Aug 2022 10:32:48 -0400 Subject: [PATCH 1/4] update format for besselj --- src/besselj.jl | 96 +++++++++++++++++++++++++------------------- src/bessely.jl | 2 +- test/besselj_test.jl | 16 ++++---- 3 files changed, 63 insertions(+), 51 deletions(-) diff --git a/src/besselj.jl b/src/besselj.jl index 3e63cfc..222d575 100644 --- a/src/besselj.jl +++ b/src/besselj.jl @@ -186,59 +186,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 - - # 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)) + # use power series for small x and for when nu > x + besselj_series_cutoff(nu, x) && return besselj_power_series(nu, x) - debye_diff = debye_cutoff - nu - large_arg_diff = nu - x / 2.0 + # 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 - 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 + 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 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 +256,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 +277,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 +=# diff --git a/src/bessely.jl b/src/bessely.jl index 8b2d6dc..e676706 100644 --- a/src/bessely.jl +++ b/src/bessely.jl @@ -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. diff --git a/test/besselj_test.jl b/test/besselj_test.jl index a68cbfb..c7e89aa 100644 --- a/test/besselj_test.jl +++ b/test/besselj_test.jl @@ -109,8 +109,8 @@ 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(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) # need to test accuracy of negative orders and negative arguments and all combinations within @@ -118,14 +118,14 @@ end # 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) From b6c56856c20d16566f165d47870d8a18a7ba1b91 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Tue, 2 Aug 2022 10:51:15 -0400 Subject: [PATCH 2/4] add more docs --- src/besselj.jl | 41 +++++++++++++++++++++++++++++++++++++++++ src/bessely.jl | 2 +- 2 files changed, 42 insertions(+), 1 deletion(-) diff --git a/src/besselj.jl b/src/besselj.jl index 222d575..74e31ae 100644 --- a/src/besselj.jl +++ b/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) diff --git a/src/bessely.jl b/src/bessely.jl index e676706..2412d32 100644 --- a/src/bessely.jl +++ b/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. From d608e518d0211e970f4488e98e8795932df9bb0b Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Tue, 2 Aug 2022 11:11:48 -0400 Subject: [PATCH 3/4] add more docs and udpate compat --- Project.toml | 5 +---- src/Bessels.jl | 2 -- src/U_polynomials.jl | 24 ++++++++++++++++++++---- src/asymptotics.jl | 15 +++++++++++++++ 4 files changed, 36 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index f139f7c..3f3b7ee 100644 --- a/Project.toml +++ b/Project.toml @@ -3,11 +3,8 @@ uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38" authors = ["Michael Helton 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" diff --git a/src/Bessels.jl b/src/Bessels.jl index 280b75b..47c7310 100644 --- a/src/Bessels.jl +++ b/src/Bessels.jl @@ -1,7 +1,5 @@ module Bessels -using SpecialFunctions: loggamma - export besselj0 export besselj1 export besselj diff --git a/src/U_polynomials.jl b/src/U_polynomials.jl index 0096d2f..eec298d 100644 --- a/src/U_polynomials.jl +++ b/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) diff --git a/src/asymptotics.jl b/src/asymptotics.jl index 5f04b5f..c431fcb 100644 --- a/src/asymptotics.jl +++ b/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) From 4fd03bda2a7126587dac9e2809c1a7906774eb45 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Tue, 2 Aug 2022 11:39:13 -0400 Subject: [PATCH 4/4] fix coverage --- src/asymptotics.jl | 1 - test/besselj_test.jl | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/asymptotics.jl b/src/asymptotics.jl index c431fcb..e71117b 100644 --- a/src/asymptotics.jl +++ b/src/asymptotics.jl @@ -134,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 diff --git a/test/besselj_test.jl b/test/besselj_test.jl index c7e89aa..238df6e 100644 --- a/test/besselj_test.jl +++ b/test/besselj_test.jl @@ -109,8 +109,8 @@ 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