Return sequences of Bessel functions for many orders

[heltonmc]

See similar issue at JuliaMath/SpecialFunctions.jl#26.

This is a feature request from discourse. Essentially it would be nice to take arguments like besselj(nu:Vector, x::T) that returns a sequence of Bessel functions with step size equal to 1. Essentially a fast way to generate besselj(n, x) for n = 0, 1, 2, 3, ....

I'd probably pose we have two options a non allocating version where the user can preallocate the memory like besselj!(out, nu, x) and then one that allocates for besselj(nu, x) if nu is some AbstractVector.

Of course these are all filled by recurrence relations. Downward recurrence for besselj and besseli and upward recurrence using bessely and besselk.

There's a couple things to make sure we get right...

• Downward recurrence could have start values (nu[end]) be zero if nu >> x. Not sure we want to just compute until there is a finite value as we move down the vector as that is inefficient...
• Upward recurrence could return infinity which then could return NaN as we get a subtraction of infinities in the recurrence
• We also need to be careful about start and end nu values on either side of zero. Something like nu=-10:10 shouldn't be allowed as it will not be numerically stable crossing zero.
• There is also the other concern of during some routines we are already generating values from recurrence so we already have two starting values but not always. It would be faster to modify the whole routine but I'm not for sure that is worth it because as length(nu) >> 1 the computation time of a two bessel calls will be a fraction of total runtime....

Just some quick thoughts on things to be careful about. The code is already here though just need to figure out finer details.

[krcools]

I also seem to remember that one direction of recursion is numerical unstable for Bessels and the other for Neumanns/Hankels. You probably have much deeper insight in this than me but I thought I'd mention it for completeness...

[heltonmc]

Yes thank you! It should be noted that since the regular Bessel functions are oscillatory there are some regions where they are stable in both directions. For instance besselj can actually be computed with forward recurrence if x < nu but x~nu is a turning point in that function and becomes unstable afterwards. bessely can also be computed with downward recurrence in that region as well but that is never useful.

We use this a lot already for particularly difficult regions (

Bessels.jl/src/hankel.jl

Lines 126 to 138 in 372120a

	# at this point x > 19.0 (for Float64) and fairly close to nu
	# shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
	nu_shift = ceil(nu) - floor(Int, -1.5 + x + Base.Math._approx_cbrt(-411*x)) + 2
	v2 = maximum((nu - nu_shift, modf(nu)[1] + 1))
	Hnu = hankel_debye(v2, x)
	Hnum1 = hankel_debye(v2 - 1, x)
	# forward recurrence is stable for Hankel when x >= nu
	if x >= nu
	H = besselj_up_recurrence(x, Hnu, Hnum1, v2, nu)[1]
	return real(H), imag(H)
	else

). The Hankel function though because it is defined as complex(besselj, bessely) needs some additional care. It's a hard function to just arbitrarily use recurrence. The continued fraction approach will be useful here.