Problems with SpecialFunctions.jl

General Usage
,
1

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 :slight_smile:

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
2

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

2 Likes
3

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

4

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
5

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
6

Of course! I’m glad you found it useful :slight_smile:

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 :slight_smile:

3 Likes
7

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.

8

In that case, the taylor series will save you.

3 Likes
9

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