Merge pull request #39 from JuliaMath/airy

support for Airy functions

Author: heltonmc

NEWS.md

@@ -14,6 +14,7 @@ For bug fixes, performance enhancements, or fixes to unexported functions we wil
  - 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)
  - add exported methods for spherical bessel function `sphericalbesselj(nu, x)`, `sphericalbesselj(nu, x)`, (#38)
+ - add exported methods for airy functions `airyai(x)`, `airyaiprime(x)`, `airybi(x)`, `airybiprime(x)`, (#39)
 
 ### Fixed
  - fix cutoff in `bessely` to not return error for integer orders and small arguments (#33)

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 and Hankel functions for real arguments. These routines are written in the Julia programming language and are self contained without any external dependencies.
+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](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 supports Bessel functions, modified Bessel functions, Hankel functions and spherical 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, 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. We plan to support complex arguments in the future. An unexported gamma function is also provided.
 
 # Quick start
 
@@ -170,7 +170,6 @@ We report the relative errors (`abs(1 - Bessels.f(x)/ArbNumerics.f(ArbFloat(x)))
 | besseli(92.12, x)  | 9e-15   | 7e-14  |
 | Bessels.gamma(x)   | 1.3e-16  | 5e-16
 
-
 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.
 
 # Benchmarks
@@ -189,7 +188,6 @@ We give brief performance comparisons to the implementations provided by [Specia
 | besselk(nu, x)   | 4x  |
 | Bessels.gamma(x)   | 5x  |
 
-
 Benchmarks were run using Julia Version 1.7.2 on an Apple M1 using Rosetta. 
 
 # API
@@ -211,13 +209,16 @@ Benchmarks were run using Julia Version 1.7.2 on an Apple M1 using Rosetta.
 - `hankelh2(nu, x)`
 - `sphericalbesselj(nu, x)`
 - `sphericalbessely(nu, x)`
+- `airyai(x)`
+- `airyaiprime(x)`
+- `airybi(x)`
+- `airybiprime(x)`
 - `Bessels.gamma(x)`
 
 # Current Development Plans
 
 - Support for higher precision `Double64`, `Float128`
 - Support for complex arguments (`x` and `nu`)
-- Airy functions
 - Support for derivatives with respect to argument and order
 
 # Contributing 

src/Bessels.jl

@@ -29,11 +29,17 @@ export besselh
 export hankelh1
 export hankelh2
 
+export airyai
+export airyaiprime
+export airybi
+export airybiprime
+
 include("besseli.jl")
 include("besselj.jl")
 include("besselk.jl")
 include("bessely.jl")
 include("hankel.jl")
+include("airy.jl")
 include("sphericalbessel.jl")
 
 include("Float128/besseli.jl")

src/airy.jl

@@ -0,0 +1,107 @@
+#                           Airy functions
+#
+#                       airyai(x), airybi(nu, x)
+#                       airyaiprime(x), airybiprime(nu, x)
+#
+#    A numerical routine to compute the airy functions and their derivatives.
+#    These routines use their relations to other special functions using https://dlmf.nist.gov/9.6.
+#    Specifically see (NIST 9.6.E1 - 9.6.E9) for computation from the defined bessel functions.
+#    For negative arguments these definitions are prone to some cancellation leading to higher errors.
+#    In the future, these could be replaced with more custom routines as they depend on a single variable.
+#
+    
+"""
+    airyai(x)
+Airy function of the first kind ``\\operatorname{Ai}(x)``.
+"""
+airyai(x::Real) = _airyai(float(x))
+
+_airyai(x::Float16) = Float16(_airyai(Float32(x)))
+
+function _airyai(x::T) where T <: Union{Float32, Float64}
+    isnan(x) && return x
+    x_abs = abs(x)
+    z = 2 * x_abs^(T(3)/2) / 3
+
+    if x > zero(T)
+        return isinf(z) ? zero(T) : sqrt(x / 3) * besselk(one(T)/3, z) / π
+    elseif x < zero(T)
+        Jv, Yv = besseljy_positive_args(one(T)/3, z)
+        Jmv = (Jv - sqrt(T(3)) * Yv) / 2
+        return isinf(z) ? throw(DomainError(x, "airyai(x) is not defined.")) : sqrt(x_abs) * (Jmv + Jv) / 3
+    elseif iszero(x)
+        return T(0.3550280538878172)
+    end
+end
+
+"""
+    airyaiprime(x)
+Derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x)``.
+"""
+airyaiprime(x::Real) = _airyaiprime(float(x))
+
+_airyaiprime(x::Float16) = Float16(_airyaiprime(Float32(x)))
+
+function _airyaiprime(x::T) where T <: Union{Float32, Float64}
+    isnan(x) && return x
+    x_abs = abs(x)
+    z = 2 * x_abs^(T(3)/2) / 3
+
+    if x > zero(T)
+        return isinf(z) ? zero(T) : -x * besselk(T(2)/3, z) / (sqrt(T(3)) * π)
+    elseif x < zero(T)
+        Jv, Yv = besseljy_positive_args(T(2)/3, z)
+        Jmv = -(Jv + sqrt(T(3))*Yv) / 2
+        return isinf(z) ? throw(DomainError(x, "airyaiprime(x) is not defined.")) : x_abs * (Jv - Jmv) / 3
+    elseif iszero(x)
+        return T(-0.2588194037928068)
+    end
+end
+
+"""
+    airybi(x)
+Airy function of the second kind ``\\operatorname{Bi}(x)``.
+"""
+airybi(x::Real) = _airybi(float(x))
+
+_airybi(x::Float16) = Float16(_airybi(Float32(x)))
+
+function _airybi(x::T) where T <: Union{Float32, Float64}
+    isnan(x) && return x
+    x_abs = abs(x)
+    z = 2 * x_abs^(T(3)/2) / 3
+
+    if x > zero(T)
+        return isinf(z) ? z : sqrt(x / 3) * (besseli(-one(T)/3, z) + besseli(one(T)/3, z))
+    elseif x < zero(T)
+        Jv, Yv = besseljy_positive_args(one(T)/3, z)
+        Jmv = (Jv - sqrt(T(3)) * Yv) / 2
+        return isinf(z) ? throw(DomainError(x, "airybi(x) is not defined.")) : sqrt(x_abs/3) * (Jmv - Jv)
+    elseif iszero(x)
+        return T(0.6149266274460007)
+    end
+end
+
+"""
+    airybiprime(x)
+Derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x)``.
+"""
+airybiprime(x::Real) = _airybiprime(float(x))
+
+_airybiprime(x::Float16) = Float16(_airybiprime(Float32(x)))
+
+function _airybiprime(x::T) where T <: Union{Float32, Float64}
+    isnan(x) && return x
+    x_abs = abs(x)
+    z = 2 * x_abs^(T(3)/2) / 3
+
+    if x > zero(T)
+        return isinf(z) ? z : x * (besseli(T(2)/3, z) + besseli(-T(2)/3, z)) / sqrt(T(3))
+    elseif x < zero(T)
+        Jv, Yv = besseljy_positive_args(T(2)/3, z)
+        Jmv = -(Jv + sqrt(T(3))*Yv) / 2
+        return isinf(z) ? throw(DomainError(x, "airybiprime(x) is not defined.")) : x_abs * (Jv + Jmv) / sqrt(T(3))
+    elseif iszero(x)
+        return T(0.4482883573538264)
+    end
+end

src/besseli.jl

@@ -207,6 +207,7 @@ Modified Bessel function of the first kind of order nu, ``I_{nu}(x)`` for positi
 function besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64}
     iszero(nu) && return besseli0(x)
     isone(nu) && return besseli1(x)
+    isinf(x) && return T(Inf)
 
     # use large argument expansion if x >> nu
     besseli_large_argument_cutoff(nu, x) && return besseli_large_argument(nu, x)

src/besselk.jl

@@ -217,6 +217,7 @@ Modified Bessel function of the second kind of order nu, ``K_{nu}(x)`` valid for
 """
 function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64}
     iszero(x) && return T(Inf)
+    isinf(x) && return zero(T)
 
     # dispatch to avoid uniform expansion when nu = 0 
     iszero(nu) && return besselk0(x)

test/airy_test.jl

@@ -0,0 +1,29 @@
+# test the airy functions
+# these are prone to some cancellation as they are indirectly calculated from relations to bessel functions
+# this is amplified for negative arguments
+# the implementations provided by SpecialFunctions.jl suffer from same inaccuracies
+for x in [0.0; 1e-17:0.1:100.0]
+    @test isapprox(airyai(x), SpecialFunctions.airyai(x), rtol=1e-12)
+    @test isapprox(airyai(-x), SpecialFunctions.airyai(-x), rtol=1e-9)
+
+    @test isapprox(airyaiprime(x), SpecialFunctions.airyaiprime(x), rtol=1e-12)
+    @test isapprox(airyaiprime(-x), SpecialFunctions.airyaiprime(-x), rtol=1e-9)
+
+    @test isapprox(airybi(x), SpecialFunctions.airybi(x), rtol=1e-12)
+    @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=5e-8)
+
+    @test isapprox(airybiprime(x), SpecialFunctions.airybiprime(x), rtol=1e-12)
+    @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=1e-9)
+end
+
+# test Inf
+@test iszero(airyai(Inf))
+@test iszero(airyaiprime(Inf))
+@test isinf(airybi(Inf))
+@test isinf(airybiprime(Inf))
+
+# test Float16 types
+@test airyai(Float16(1.2)) isa Float16
+@test airyaiprime(Float16(1.9)) isa Float16
+@test airybi(Float16(1.2)) isa Float16
+@test airybiprime(Float16(1.9)) isa Float16

test/besseli_test.jl

@@ -102,6 +102,9 @@ for v in nu, xx in x
     @test isapprox(besseli(v, xx), Float64(sf), rtol=2e-13)
 end
 
+# test Inf
+@test isinf(besseli(2, Inf))
+
  ### tests for negative arguments
 
 (v, x) = 12.0, 3.2

test/besselk_test.jl

@@ -120,6 +120,9 @@ for v in nu, xx in x
     @test isapprox(besselk(v, xx), Float64(sf), rtol=2e-13)
 end
 
+# test Inf
+@test iszero(besselk(2, Inf))
+
 ### tests for negative arguments
 
 (v, x) = 12.0, 3.2

test/runtests.jl

@@ -8,4 +8,5 @@ import SpecialFunctions
 @time @testset "bessely" begin include("bessely_test.jl") end
 @time @testset "hankel" begin include("hankel_test.jl") end
 @time @testset "gamma" begin include("gamma_test.jl") end
+@time @testset "airy" begin include("airy_test.jl") end
 @time @testset "sphericalbessel" begin include("sphericalbessel_test.jl") end