CoulombFunctions.jl
This library provide efficient and accurate implementations of the regular and irregular spherical Bessel functions [$j_n(z)$ and $y_n(z)$] and Coulomb functions [$F_\lambda(\eta,z)$ and $G_\lambda(\eta,z)$], and their derivatives. The former are solutions to the radial part of the Helmholtz equation in spherical coordinates:
\[\tag{DLMF10.47.1} z^2w'' + 2zw' + [z^2 - n(n+1)] = 0,\]
whereas the latter obey
\[\tag{DLMF33.2.1} w'' + \left[ 1 - \frac{2\eta}{z} - \frac{\lambda(\lambda+1)}{z^2} \right]w = 0.\]
The spherical Bessel functions are related to the Coulomb functions as
\[\tag{DLMF33.5.3} \begin{aligned} j_n(z) &= \frac{F_n(0,z)}{z}, & y_n(z) &= -\frac{G_n(0,z)}{z}. \end{aligned}\]
Usage
x = range(-12, stop=12, length=1000)
nℓ = 10
j, j′, y, y′ = bessels(x, nℓ)
x = range(0, stop=12, length=1000)
nℓ = 10
η = -1.0
F, F′, G, G′ = coulombs(x, η, nℓ)
Accuracy
To check the accuracy, we compare with SpecialFunctions.jl, which however does not provide the spherical Bessel functions but the ordinary (cylindrical) ones. They are however related as
\[\tag{DLMF10.47.\{3,4\}} \begin{aligned} j_n(z) &\equiv \sqrt{\frac{\pi}{2z}} J_{n+\frac{1}{2}}(z),\\ y_n(z) &\equiv \sqrt{\frac{\pi}{2z}} Y_{n+\frac{1}{2}}(z). \end{aligned}\]
Furthermore, SpecialFunctions.jl does not provide the derivatives out-of-the-box, but they are easily found using the recurrence relation.:
\[f'_n(z) = f_{n-1}(z) - \frac{n+1}{z} f_n(z), \qquad f_n = j_n,y_n.\]
This time, we investigate a larger domain of parameters, but avoid smaller values of $x$ than $0.1$ since that does not seem to work in SpecialFunctions.jl (CoulombFunctions.jl works at $x\leq0$ as well):
nx = 1001
x = 10 .^ range(-1, stop=4, length=nx)
nℓ = 105
j, j′, y, y′ = bessels(x, nℓ)
We note that CoulombFunctions.jl seems to be ~3 times faster than SpecialFunctions.jl when evaluating all Bessel functions for a fixed value of $x$, most likely due to extra processing taking place when SpecialFunctions.jl does not compute all values simultaneously, but one order at a time:
CoulombFunctions.jl:
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 76.331 μs (0.00% GC)
median time: 77.057 μs (0.00% GC)
mean time: 80.337 μs (0.00% GC)
maximum time: 197.715 μs (0.00% GC)
--------------
samples: 10000
evals/sample: 1
SpecialFunctions.jl:
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 170.295 μs (0.00% GC)
median time: 185.132 μs (0.00% GC)
mean time: 193.312 μs (0.00% GC)
maximum time: 442.129 μs (0.00% GC)
--------------
samples: 10000
evals/sample: 1
Finally, the agreement for the (irregular) Neumann functions for $x < 100$ is terrible, unclear why. But they are irregular (diverging) as $x$ tends to zero.