Bessel function of second with complex order

General Usage
1

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.

2

Your best bet might be Nemo.jl.

3

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.

4

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)
   end
end

then

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

besselj(order, argument)
# 0.6415000256096106 - 0.2518459064184279im

bessely(order, argument)
# 0.12578491928293956 + 0.17534238733104368im
5

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)
			s+=r

			if ((abs(real(r/s))<precision)&&(abs(imag(r/s))<precision))
				break
			end
		end
	else
		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
			s+=r

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

	return s

end

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

Cheers!

3 Likes