Bessel function of second with complex order


Is there a math library in Julia which can compute a modified Bessel function of the second kind with complex order and argument? I am only aware of SpecialFunctions.jl which contains a method to compute a modified Bessel function of the second kind with a complex argument of real order.


Your best bet might be Nemo.jl.


I had a student working on a related problem at one point — my recollection is that we couldn’t find any free code for this, and very little literature. He developed some prototypes, but as I recall he ran into tricky accuracy issues with the second-kind Bessel functions for complex arguments that were close to integers. We also wanted the derivative with respect to the order, however, so maybe the problem was easier if you don’t need that. I’ll bug him to see whether he has any notes he can post.


If you want to use a package that is in development,
add SpecialFunctions and
add ArbNumerics (must use version 0.3.5+, just merged).

using ArbNumerics
import SpecialFunctions: besselj, bessely

for F in (:besselj, :bessely)
   @eval function $F(order::Complex{T}, argument::Complex{T}) where {T}
       result = $F(ArbComplex(order), ArbComplex(argument))
       return Complex{T}(result)


order = 0.0 + 0.5im; argument = 1.0 + 1.0im;

besselj(order, argument)
# 0.6415000256096106 - 0.2518459064184279im

bessely(order, argument)
# 0.12578491928293956 + 0.17534238733104368im


Thank you all for your response. I finally ended up constructing my own function using SpecialFunctions.jl. The lines of codes below seem to be performing very well.

using SpecialFunctions

# Define the modified Bessel Function of the second kind with complex order & argument
function CBesselK(nu, z)

	Knu = exp(lgamma(nu))  * (z/2)^(-nu)
	Kmu = exp(lgamma(-nu)) * (z/2)^nu

	bound = 20

	precision = 1e-10

	s = 0

	if (abs(z) < 10)
		for i = 0:bound
			var = i
			facto = factorial(i)

			r = 0.5 * Knu * ((z/2)^(2*var) * exp(-lgamma(1+var-nu)+lgamma(1-nu))/facto) + 0.5 * Kmu * ((z/2)^(2*var) * exp(-lgamma(1+var+nu)+lgamma(1+nu))/facto)

			if ((abs(real(r/s))<precision)&&(abs(imag(r/s))<precision))
		for i = 0:bound
			var = i
			facto = factorial(i)

			r = sqrt(pi/(2*z)) * exp(-z) * exp(lgamma(0.5+var+nu)+lgamma(0.5+var-nu)-lgamma(0.5+nu)-lgamma(0.5-nu)/facto) * (-1/2*z)^var

			if ((abs(real(r/s))<precision)&&(abs(imag(r/s))<precision))

	return s


Comparing the performance with Mathematica:

# CBesselK
order   = exp(1-3im)
arg     = (4+6im)*exp(2*im*pi)
CBesselK(order, arg)
0.011090297922744 - 0.005887726316203im

# Mathematica 
N[BesselK[Exp[1 - 3 I], (4 + 6 I)*Exp[2*I*Pi]], 15]
0.011090297922749 - 0.005887726316251 I