# Problems with SpecialFunctions.jl

Hey guys, im having trouble using BigFloat variables on spherical bessel functions given by the SpecialFunctions,jl package. I’m using the following code and receiving the error below. Thanks beforehand

``````function partialwavexp(psiamp, axiconang, order, r, theta, phi)
k = 1000 * 2 * big(pi) / 1.54
kr = k * r
krho = k * sin(axiconang)
kz = k * cos(axiconang)
phi0 = z0 = rho0 = 0
nmax = Integer(ceil(kr + (big(405) / 100) * (kr^(1/3)) + 2))
psi = 0
for n in 0:nmax
for m in -n:n
BSC = bsccalc(n, m, axiconang, order, krho, kz, phi0, z0, rho0)
psi += BSC * sphericalbesselj(n, kr) * Plm(n, m, cos(theta)) * cis(m * phi)
end
end
return psi * psiamp
end
``````

And the error:

``````ERROR: LoadError: MethodError: no method matching besselj(::BigFloat, ::Complex{BigFloat})
Closest candidates are:
besselj(::T, ::Complex{T}) where T<:AbstractFloat at ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:590
besselj(::Real, ::Complex) at ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:584
besselj(::Real, ::AbstractFloat) at ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:482
...
Stacktrace:
[1] besselj(k::BigFloat, z::Complex{BigFloat})
@ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:590
[2] besselj(nu::BigFloat, x::BigFloat)
@ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:490
[3] sphericalbesselj(nu::BigInt, x::BigFloat)
@ SpecialFunctions ~/.julia/packages/SpecialFunctions/hefUc/src/bessel.jl:701
[4] partialwavexp(psiamp::Int64, axiconang::BigFloat, order::Int64, r::BigFloat, theta::Float64, phi::BigFloat)
@ Main ~/Documentos/BSCSim/src/naive/naive_implem.jl:30

``````
1 Like

if you need extended precision Bessel functions, they are provided by Arblib.jl.

2 Likes

Thanks for the advice , but i was not able to find spherical functions on the arb packages and documentation. Are there any other wrappers you could indicate please.

Just use the relation to the normal Bessel function http://dlmf.nist.gov/10.47.E3.

``````sphericalbesselj(nu, x::T) where T = sqrt(T(π)/2) * besselj(nu + one(T)/2, x) / sqrt(x)
``````

You could do much better for performance though by moving the `sphericalbesselj` outside of the `m` loop and using recurrence. You should probably also move the `sphericalbesselj` completely outside the loop and use recurrence.

This is only stable in the backward direction so you would need to compute `sphericalbesselj(nmax, x)` and `sphericalbesselj(nmax-1, x)` and sum from `for n in nmax:-1:0` then you would only need two bessel function evaluations… anyway that requires a bit more care in your implementation but the speed gains will be immense if the other calculations in your loop are cheap.

Edit: My previous advice if you are evaluating for complex numbers might need more care for stability in recurrence.

Edit2: I also usually use `ArbNumerics.jl` but the default precision is not always enough if `nu` or `x` is large so you might need to increase the precision to get fully accurate results.

2 Likes

Hello Helton, i cannot thank you enough for helping with this specific problem and giving valuable advice. bessel.jl developed by you and Oscar that also answered my question is a cornerstone on my new research project. I am still a julia newbie as you can see, but somehow in the future i still want to contribute to some JuliaMath repo due to how impressive they are. Again, thank you very much.

2 Likes

Of course! I’m glad you found it useful

Just to add another recommendation for you to consider. In general, you usually can avoid directly calculating Bessel (or Gamma) functions within the loop by using recurrence. Though, your problem sets up well for a different approach.

You could apply Miller’s recurrence algorithm. It will require some care for your desired precision and you need to determine how much higher your trial M should be over N. I’m not for sure if you need an arbitrary precision routine but if you only need a set amount (say `Float128` or `Double64`) this number can be determined before hand.

The main advantage of this is that you must calculate `sphericalbesselj(0, kr)` at the start (or end of the loop) which has a closed form expression (`j_0=sin(x)/x`) this is important because you can use that as your normalization condition at the end. So, now instead of calculating `sphericalbesselj` everytime within your loop you just have to calculate `sin(x)/x` at the end of your loop one time and normalize your sum. The advantage of this is you can find fast quadruple or double-double implementations of `sin(x)` in `DoubleFloats.jl` and `QuadMath.jl` if you don’t need an arbitrary precision routine. I’m away traveling currently so can’t code this up but would be happy to try and answer any questions if you want to give it a try! This approach essentially avoids computing the Bessel function in higher precision which will give you a significant speed gain. Of course if performance doesn’t matter then the above suggestion will be much easier

3 Likes

Heeeey, Miller’s algorithm seems great. I will implement it tomorrow, probably using float128, as a first approach i don’t think i’ll be needing more precision than that. My concern solely is when x approaches 0 in `j_0=sin(x)/x`, but i’ll try to work around that!
Thank you so much again and i hope you have a great traveling experience.

In that case, the taylor series will save you.

3 Likes

Nice, i’ll implement it using the taylor series expansion to some degree. I cannot stress how much you guys helped. Thanks a lot.

1 Like