Julia: How to change the type of an argument in the Bessel function?

user_14426118 0
Viewed 67

Here is my code:

using Plots
using SpecialFunctions
using QuadGK

kappa = 1
B = 1
xi = (kappa/B)^4

function correlation(x)
    quadgk(q -> q * SpecialFunctions.besselj0(x*q)/(q^4 + xi), 0, 1e6)[1]/kappa 
end

r = range(-20, 20, length = 1001)

plot(r, correlation(r))

I get an error on the Bessel function. I get that the argument is the problem and that it should be of the format ::BigFloat, ::Float16, or ::Float32, but I don't know how to do it. I tried to write x .* q instead, but the problem remains the same, I get the error:

ERROR: MethodError: no method matching besselj0(::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}})

Also, I'm searching for a way to write +infinity instead of 1e6.

julia
1 Answers

Just replace correlation(r) with correlation.(r) in your code to use broadcasting, as is explained here.

The core of your problem that in Julia functions are not broadcastable by default - you usually have to opt-in to have it (especially when you work with numerical code). Here is a basic example:

julia> sin(1)
0.8414709848078965

julia> sin([1])
ERROR: MethodError: no method matching sin(::Array{Int64,1})
Closest candidates are:
  sin(::BigFloat) at mpfr.jl:727
  sin(::Missing) at math.jl:1197
  sin(::Complex{Float16}) at math.jl:1145
  ...
Stacktrace:
 [1] top-level scope at REPL[2]:1

julia> sin.([1])
1-element Array{Float64,1}:
 0.8414709848078965

However, in your case the correlation function is quite expensive. In such a case I usually use ProgressMeter.jl to monitor the progress of the computations (it shows how long you can expect for the computations to finish). So you can write:

using ProgressMeter
result = @showprogress map(correlation, r)

and use the map function to apply correlation function to all elements of r (in this case the result will be the same as for broadcasting).

Finally, your computations will be much faster if you do not use global variables in quadgk. It is better to pass kappa and xi as arguments to the function like this:

function correlation(x, kappa, xi)
    quadgk(q -> q * SpecialFunctions.besselj0(x*q)/(q^4 + xi), 0, 1e6)[1]/kappa 
end

result = @showprogress map(x -> correlation(x, kappa, xi), r)
user_12695670
Related