Skip to content

Commit 00cd0ec

Browse files
committed
add docs
1 parent cd97b85 commit 00cd0ec

File tree

2 files changed

+63
-2
lines changed

2 files changed

+63
-2
lines changed

src/besseli.jl

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -130,6 +130,35 @@ function besseli1x(x::T) where T <: Union{Float32, Float64}
130130
return z
131131
end
132132

133+
# Modified Bessel functions of the first kind of order nu
134+
# besseli(nu, x)
135+
#
136+
# A numerical routine to compute the modified Bessel function of the first kind I_{ν}(x) [1]
137+
# for real orders and arguments of positive or negative value. The routine is based on several
138+
# publications [2-6] that calculate I_{ν}(x) for positive arguments and orders where
139+
# reflection identities are used to compute negative arguments and orders.
140+
#
141+
# In particular, the reflectance identities for negative noninteger orders I_{−ν}(x) = I_{ν}(x) + 2 / πsin(πν)*Kν(x)
142+
# and for negative integer orders I_{−n}(x) = I_n(x) are used.
143+
# For negative arguments of integer order, In(−x) = (−1)^n In(x) is used and for
144+
# noninteger orders, Iν(−x) = exp(iπν) Iν(x) is used. For negative orders and arguments the previous identities are combined.
145+
#
146+
# The identities are computed by calling the `besseli_positive_args(nu, x)` function which computes I_{ν}(x)
147+
# for positive arguments and orders. For large orders, Debye's uniform asymptotic expansions are used where large arguments (x>>nu)
148+
# a large argument expansion is used. The rest of the values are computed using the power series.
149+
150+
# [1] https://dlmf.nist.gov/10.40.E1
151+
# [2] Amos, Donald E. "Computation of modified Bessel functions and their ratios." Mathematics of computation 28.125 (1974): 239-251.
152+
# [3] Gatto, M. A., and J. B. Seery. "Numerical evaluation of the modified Bessel functions I and K."
153+
# Computers & Mathematics with Applications 7.3 (1981): 203-209.
154+
# [4] Temme, Nico M. "On the numerical evaluation of the modified Bessel function of the third kind."
155+
# Journal of Computational Physics 19.3 (1975): 324-337.
156+
# [5] Amos, DEv. "Algorithm 644: A portable package for Bessel functions of a complex argument and nonnegative order."
157+
# ACM Transactions on Mathematical Software (TOMS) 12.3 (1986): 265-273.
158+
# [6] Segura, Javier, P. Fernández de Córdoba, and Yu L. Ratis. "A code to evaluate modified bessel functions based on thecontinued fraction method."
159+
# Computer physics communications 105.2-3 (1997): 263-272.
160+
#
161+
133162
"""
134163
besseli(nu, x::T) where T <: Union{Float32, Float64}
135164

src/besselk.jl

Lines changed: 34 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -146,6 +146,40 @@ function besselk1x(x::T) where T <: Union{Float32, Float64}
146146
end
147147
end
148148

149+
# Modified Bessel functions of the second kind of order nu
150+
# besselk(nu, x)
151+
#
152+
# A numerical routine to compute the modified Bessel function of the second kind K_{ν}(x) [1]
153+
# for real orders and arguments of positive or negative value. The routine is based on several
154+
# publications [2-8] that calculate K_{ν}(x) for positive arguments and orders where
155+
# reflection identities are used to compute negative arguments and orders.
156+
#
157+
# In particular, the reflectance identities for negative orders I_{−ν}(x) = I_{ν}(x).
158+
# For negative arguments of integer order, Kn(−x) = (−1)^n Kn(x) − im * π In(x) is used and for
159+
# noninteger orders, Kν(−x) = exp(−iπν)*Kν(x) − im π Iν(x) is used. For negative orders and arguments the previous identities are combined.
160+
#
161+
# The identities are computed by calling the `besseli_positive_args(nu, x)` function which computes K_{ν}(x)
162+
# for positive arguments and orders. For large orders, Debye's uniform asymptotic expansions are used.
163+
# For small value and when nu > ~x the power series is used. The rest of the values are computed using slightly different methods.
164+
# The power series for besseli is modified to give both I_{v} and I_{v-1} where the ratio K_{v+1} / K_{v} is computed using continued fractions [8].
165+
# The wronskian connection formula is then used to compute K_v.
166+
167+
# [1] http://dlmf.nist.gov/10.27.E4
168+
# [2] Amos, Donald E. "Computation of modified Bessel functions and their ratios." Mathematics of computation 28.125 (1974): 239-251.
169+
# [3] Gatto, M. A., and J. B. Seery. "Numerical evaluation of the modified Bessel functions I and K."
170+
# Computers & Mathematics with Applications 7.3 (1981): 203-209.
171+
# [4] Temme, Nico M. "On the numerical evaluation of the modified Bessel function of the third kind."
172+
# Journal of Computational Physics 19.3 (1975): 324-337.
173+
# [5] Amos, DEv. "Algorithm 644: A portable package for Bessel functions of a complex argument and nonnegative order."
174+
# ACM Transactions on Mathematical Software (TOMS) 12.3 (1986): 265-273.
175+
# [6] Segura, Javier, P. Fernández de Córdoba, and Yu L. Ratis. "A code to evaluate modified bessel functions based on thecontinued fraction method."
176+
# Computer physics communications 105.2-3 (1997): 263-272.
177+
# [7] Geoga, Christopher J., et al. "Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation."
178+
# arXiv preprint arXiv:2201.00090 (2022).
179+
# [8] Cuyt, A. A., Petersen, V., Verdonk, B., Waadeland, H., & Jones, W. B. (2008).
180+
# Handbook of continued fractions for special functions. Springer Science & Business Media.
181+
#
182+
149183
"""
150184
besselk(x::T) where T <: Union{Float32, Float64}
151185
@@ -305,8 +339,6 @@ end
305339
# A modified form appears below. See more discussion at https://github.com/heltonmc/Bessels.jl/pull/29
306340
# This is only accurate for noninteger orders (nu) and no checks are performed.
307341
#
308-
# [7] Geoga, Christopher J., et al. "Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation."
309-
# arXiv preprint arXiv:2201.00090 (2022).
310342
"""
311343
besselk_power_series(nu, x::T) where T <: Float64
312344

0 commit comments

Comments
 (0)