Merge pull request #33 from JuliaMath/hankel

Support Hankel functions and efficient methods to get Jv and Yv

Author: heltonmc

NEWS.md

@@ -0,0 +1,27 @@
+# Changelog
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+For pre-1.0.0 releases we will increment the minor version when we export any new functions or alter exported API.
+For bug fixes, performance enhancements, or fixes to unexported functions we will increment the patch version.
+**Note**: The exported API can be considered very stable and likely will not change without serious consideration.
+
+# Unreleased
+
+### Added
+ - add an unexport method (`Bessels.besseljy(nu, x)`) for faster computation of `besselj` and `bessely` (#33)
+ - add exported methods for Hankel functions `besselh(nu, k, x)`, `hankelh1(nu, x)`, `hankelh2(nu, x)` (#33)
+
+### Fixed
+ - fix cutoff in `bessely` to not return error for integer orders and small arguments (#33)
+
+# Version 0.1.0
+
+Initial release of Bessels.jl
+
+### Added
+
+ - support bessel functions (`besselj`, `bessely`, `besseli`, `besselk`) for real arguments
+ - provide a gamma function (`Bessels.gamma`) for real arguments

README.md

@@ -2,11 +2,11 @@
 [![Build Status](https://github.com/heltonmc/Bessels.jl/actions/workflows/CI.yml/badge.svg?branch=master)](https://github.com/heltonmc/Bessels.jl/actions/workflows/CI.yml?query=branch%3Amaster)
 [![Coverage](https://codecov.io/gh/heltonmc/Bessels.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/heltonmc/Bessels.jl)
 
-Numerical routines for computing Bessel functions and modified Bessel functions of the first and second kind. These routines are written in the Julia programming language and are self contained without any external dependencies.
+Numerical routines for computing Bessel 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](https://github.com/JuliaMath/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](https://github.com/heltonmc/Bessels.jl/edit/update_readme/README.md#benchmarks)).
 
-The library currently only supports Bessel functions and modified Bessel 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. We plan to support complex arguments in the future. An unexported gamma function is also provided.
+The library currently supports Bessel functions, modified Bessel functions and Hankel 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. We plan to support complex arguments in the future. An unexported gamma function is also provided.
 
 # Quick start
 
@@ -206,12 +206,14 @@ Benchmarks were run using Julia Version 1.7.2 on an Apple M1 using Rosetta.
 - `besselk0(x)`
 - `besselk1(x)`
 - `besselk(nu, x)`
+- `besselh(nu, k, x)`
+- `hankelh1(nu, x)`
+- `hankelh2(nu, x)`
 - `Bessels.gamma(x)`
 
 # Current Development Plans
 
 - Support for higher precision `Double64`, `Float128`
-- Hankel functions
 - Support for complex arguments (`x` and `nu`)
 - Airy functions
 - Support for derivatives with respect to argument and order

src/Bessels.jl

@@ -22,18 +22,23 @@ export besselk0x
 export besselk1
 export besselk1x
 
+export besselh
+export hankelh1
+export hankelh2
+
 include("besseli.jl")
 include("besselj.jl")
 include("besselk.jl")
 include("bessely.jl")
-include("constants.jl")
+include("hankel.jl")
 
 include("Float128/besseli.jl")
 include("Float128/besselj.jl")
 include("Float128/besselk.jl")
 include("Float128/bessely.jl")
 include("Float128/constants.jl")
 
+include("constants.jl")
 include("math_constants.jl")
 include("U_polynomials.jl")
 include("recurrence.jl")
@@ -42,6 +47,4 @@ include("Polynomials/besselj_polys.jl")
 include("asymptotics.jl")
 include("gamma.jl")
 
-#include("hankel.jl")
-
 end

src/U_polynomials.jl

@@ -41,7 +41,7 @@ function besseljy_debye(v, x)
     return coef_Jn * Uk_Jn, coef_Yn * Uk_Yn
 end
 
-besseljy_debye_cutoff(nu, x) = nu > 2.0 + 1.00035*x + Base.Math._approx_cbrt(Float64(302.681)*x) && x > 15
+besseljy_debye_cutoff(nu, x) = nu > 2.0 + 1.00035*x + Base.Math._approx_cbrt(Float64(302.681)*x) && nu > 15
 
 """
     hankel_debye(nu, x::T)

src/bessely.jl

@@ -228,13 +228,13 @@ end
 #####
 ##### Generic routine for `bessely`
 #####
+
 """
     bessely(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 where nu and x can be positive or negative.
 """
-
 function bessely(nu::Real, x::T) where T
     isinteger(nu) && return bessely(Int(nu), x)
     abs_nu = abs(nu)
@@ -310,10 +310,10 @@ function bessely_positive_args(nu, x::T) where T
     hankel_debye_cutoff(nu, x) && return imag(hankel_debye(nu, x))
 
     # use power series for small x and for when nu > x
-    bessely_series_cutoff(nu, x) && return bessely_power_series(nu, x)
+    bessely_series_cutoff(nu, x) && return bessely_power_series(nu, x)[1]
 
     # for x ∈ (6, 19) we use Chebyshev approximation and forward recurrence
-    besseljy_chebyshev_cutoff(nu, x) && return bessely_chebyshev(nu, x)
+    besseljy_chebyshev_cutoff(nu, x) && return bessely_chebyshev(nu, x)[1]
 
     # at this point x > 19.0 (for Float64) and fairly close to nu
     # shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
@@ -340,6 +340,7 @@ end
 
 Computes ``Y_{nu}(x)`` using the power series when nu is not an integer.
 In general, this is most accurate for small arguments and when nu > x.
+Outpus both (Y_{nu}(x), J_{nu}(x)).
 """
 function bessely_power_series(v, x::T) where T
     MaxIter = 3000
@@ -358,7 +359,7 @@ function bessely_power_series(v, x::T) where T
         b *= -inv((-v + i + one(T)) * (i + one(T))) * t2
     end
     s, c = sincospi(v)
-    return (out*c - out2) / s
+    return (out*c - out2) / s, out
 end
 bessely_series_cutoff(v, x) = (x < 7.0) || v > 1.35*x - 4.5
 
@@ -375,7 +376,7 @@ Forward recurrence is used to fill orders starting at low orders ν ∈ (0, 2).
 function bessely_chebyshev(v, x)
     v_floor, _ = modf(v)
     Y0, Y1 = bessely_chebyshev_low_orders(v_floor, x)
-    return besselj_up_recurrence(x, Y1, Y0, v_floor + 1, v)[1]
+    return besselj_up_recurrence(x, Y1, Y0, v_floor + 1, v)
 end
 
 # only implemented for Float64 so far

src/hankel.jl

@@ -1,14 +1,207 @@
-#= Aren't tested yet
-function besselh(nu::Float64, k::Integer, x::AbstractFloat)
+#                           Hankel functions
+#
+#                           besseljy(nu, x)
+#            besselh(nu, x), hankelh1(nu, x), hankelh2(nu, x)
+#
+#    A numerical routine to compute both Bessel functions of the first J_{ν}(x) and second kind Y_{ν}(x)
+#    for real orders and arguments of positive or negative value. Please see notes in src/besselj.jl and src/bessely.jl
+#    for details on implementation as most of the routine is similar. A key difference is when the methods to compute bessely
+#    don't also give besselj. We then rely on a continued fraction approach to compute J_{ν}(x)/J_{ν-1}(x) to get J_{ν}(x)
+#    from  Y_{ν}(x) and Y_{ν-1}(x)[1]. The continued fraction approach is more quickly converging when nu and x are small in magnitude.
+#    When x and nu are large and nu = x + ϵ, we fall back to computing J_{ν}(x) and Y_{ν}(x) separately as this was found to be more efficient.
+#    
+# [1] 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.
+#
+
+#####
+##### Generic routine for `besseljy`
+#####
+
+"""
+    besseljy(nu, x::T) where T <: Float64
+
+Return both Bessel functions of the first ``J_{nu}(x)`` and
+second ``Y_{nu}(x)`` kind. This method will be faster than calling (besselj(nu, x), bessely(nu, x)) separately,
+unless nu is slightly larger than x.
+
+Results may be slightly different than individual functions in some domains due to using different algorithms.
+"""
+function besseljy(nu::Real, x::T) where T
+    isinteger(nu) && return besseljy(Int(nu), x)
+    abs_nu = abs(nu)
+    abs_x = abs(x)
+
+    Jnu, Ynu =  besseljy_positive_args(abs_nu, abs_x)
+   
+    if nu >= zero(T)
+        if x >= zero(T)
+            return Jnu, Ynu
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return Ynu * cispi(-nu) + 2im * besselj_positive_args(abs_nu, abs_x) * cospi(abs_nu)
+        end
+    else
+        spi, cpi = sincospi(abs_nu)
+        out = Jnu * cpi - Ynu * spi
+        if x >= zero(T)
+            return out, Ynu * cpi + Jnu * spi
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return cpi * (Ynu * cispi(nu) + 2im * Jnu * cpi) + Jnu * spi * cispi(abs_nu)
+        end
+    end
+end
+
+function besseljy(nu::Integer, x::T) where T
+    abs_nu = abs(nu)
+    abs_x = abs(x)
+    sg = iseven(abs_nu) ? 1 : -1
+
+    Jnu, Ynu =  besseljy_positive_args(abs_nu, abs_x)
+
+    if nu >= zero(T)
+        if x >= zero(T)
+            return Jnu, Ynu
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #return Jnu * sg, Ynu * sg + 2im * sg * besselj_positive_args(abs_nu, abs_x)
+        end
+    else
+        if x >= zero(T)
+            return Jnu * sg, Ynu * sg
+        else
+            return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
+            #spi, cpi = sincospi(abs_nu)
+            #return (cpi*Jnu - spi*Ynu) * sg, Ynu + 2im * besselj_positive_args(abs_nu, abs_x)
+        end
+    end
+end
+
+#####
+#####  `besseljy` for positive arguments and orders
+#####
+
+function besseljy_positive_args(nu::Real, x::T) where T
+    nu == 0 && return (besselj0(x), bessely0(x))
+    nu == 1 && return (besselj1(x), bessely1(x))
+    iszero(x) && return (besselj_positive_args(nu, x), bessely_positive_args(nu, x))
+
+    # x < ~nu branch see src/U_polynomials.jl
+    besseljy_debye_cutoff(nu, x) && return besseljy_debye(nu, x)
+
+    # large argument branch see src/asymptotics.jl
+    besseljy_large_argument_cutoff(nu, x) && return besseljy_large_argument(nu, x)
+
+    # x > ~nu branch see src/U_polynomials.jl on computing Hankel function
+    if hankel_debye_cutoff(nu, x)
+        H = hankel_debye(nu, x)
+        return real(H), imag(H)
+    end
+    # use forward recurrence if nu is an integer up until the continued fraction becomes inefficient
+    if isinteger(nu) && nu < 150
+        Y0 = bessely0(x)
+        Y1 = bessely1(x)
+        Ynm1, Yn = besselj_up_recurrence(x, Y1, Y0, 1, nu-1)
+
+        ratio_Jv_Jvm1 = besselj_ratio_jnu_jnum1(nu, x)
+        Jn = 2 / (π*x * (Ynm1 - Yn / ratio_Jv_Jvm1))
+        return Jn, Yn
+    end
+
+    # use power series for small x and for when nu > x
+    if bessely_series_cutoff(nu, x) && besselj_series_cutoff(nu, x)
+        Yn, Jn = bessely_power_series(nu, x)
+        return Jn, Yn
+    end
+
+    # for x ∈ (6, 19) we use Chebyshev approximation and forward recurrence
+    if besseljy_chebyshev_cutoff(nu, x)
+        Yn, Ynp1 = bessely_chebyshev(nu, x)
+        ratio_Jvp1_Jv = besselj_ratio_jnu_jnum1(nu+1, x)
+        Jnp1 = 2 / (π*x * (Yn - Ynp1 / ratio_Jvp1_Jv))
+        return Jnp1 / ratio_Jvp1_Jv, Yn
+    end
+
+    # at this point x > 19.0 (for Float64) and fairly close to nu
+    # shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
+    nu_shift = floor(nu) - ceil(Int, -1.5 + x + Base.Math._approx_cbrt(-411*x))
+    v2 = maximum((nu - nu_shift, modf(nu)[1] + 1))
+
+    Hnu = hankel_debye(v2, x)
+    Hnum1 = hankel_debye(v2 - 1, x)
+
+    # forward recurrence is stable for Hankel when x >= nu
+    if x >= nu
+        H = besselj_up_recurrence(x, Hnu, Hnum1, v2, nu)[1]
+        return real(H), imag(H)
+    else
+        # At this point besselj can not be calculated with forward recurrence
+        # We could calculate it from bessely using the continued fraction approach like the following
+        #        Yn, Ynp1 = besselj_up_recurrence(x, imag(Hnu), imag(Hnum1), v2, nu)
+        #        ratio_Jvp1_Jv = besselj_ratio_jnu_jnum1(nu+1, x)
+        #        Jnp1 = 2 / (π*x * (Yn - Ynp1 / ratio_Jvp1_Jv))
+        #        return Jnp1 / ratio_Jvp1_Jv, Yn
+        # However, the continued fraction approach is slowly converging for large arguments
+        # We will fall back to computing besselj separately instead
+        return besselj_positive_args(nu, x), besselj_up_recurrence(x, imag(Hnu), imag(Hnum1), v2, nu)[1]
+    end
+end
+
+#####
+#####  Continued fraction for J_{ν}(x)/J_{ν-1}(x)
+#####
+
+# implements continued fraction to compute ratio of J_{ν}(x)/J_{ν-1}(x)
+# using equation 22 and 24 of [1]
+# in general faster converging for small magnitudes of x and nu and nu >> x
+function besselj_ratio_jnu_jnum1(n, x::T) where T
+    MaxIter = 5000
+    xinv = inv(x)
+    xinv2 = 2 * xinv
+    d = x / (n + n)
+    a = d
+    h = a
+    b = muladd(2, n, 2) * xinv
+    for _ in 1:MaxIter
+        d = inv(b - d)
+        a *= muladd(b, d, -1)
+        h = h + a
+        b = b + xinv2
+        abs(a / h) <= eps(T) && break
+    end
+    return h
+end
+
+#####
+#####  Hankel functions
+#####
+
+"""
+    besselh(nu, [k=1,] x)
+
+Bessel function of the third kind of order `nu` (the Hankel function). `k` is either 1 or 2,
+selecting [`hankelh1`](@ref) or [`hankelh2`](@ref), respectively.
+"""
+function besselh(nu::Real, k::Integer, x)
+    Jn, Yn = besseljy(nu, x)
     if k == 1
-        return complex(besselj(nu, x), bessely(nu, x))
+        return complex(Jn, Yn)
     elseif k == 2
-        return complex(besselj(nu, x), -bessely(nu, x))
+        return complex(Jn, -Yn)
     else
         throw(ArgumentError("k must be 1 or 2"))
     end
 end
 
-hankelh1(nu, z) = besselh(nu, 1, z)
-hankelh2(nu, z) = besselh(nu, 2, z)
-=#
+"""
+    hankelh1(nu, x)
+Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``.
+"""
+hankelh1(nu, x) = besselh(nu, 1, x)
+
+"""
+    hankelh2(nu, x)
+Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x)``.
+"""
+hankelh2(nu, x) = besselh(nu, 2, x)