Блог

Главная - Вопросы по программированию - Джулия: Как изменить тип аргумента в функции Бесселя?

Джулия: Как изменить тип аргумента в функции Бесселя?

#julia

#Юлия

Вопрос:

Вот мой код:

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

Я получаю сообщение об ошибке в функции Бесселя. Я понимаю, что аргумент является проблемой и что он должен иметь формат ::BigFloat , ::Float16 , или ::Float32 , но я не знаю, как это сделать. Я попытался написать x .* q вместо этого, но проблема остается той же, я получаю ошибку:

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

Кроме того, я ищу способ написать infinity вместо 1e6 .

Комментарии:

1. Ответ на второй вопрос, вероятно Inf , <strike> , но я не знаю, принимает ли QuadGK это как границу интеграции.</strike> Похоже, это так .

Ответ №1:

Просто замените correlation(r) на correlation.(r) в вашем коде, чтобы использовать широковещательную передачу, как объясняется здесь .

Суть вашей проблемы в том, что функции Julia по умолчанию недоступны для трансляции — обычно для этого нужно зарегистрироваться (особенно когда вы работаете с числовым кодом). Вот базовый пример:

 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

Однако в вашем случае correlation функция довольно дорогая. В таком случае я обычно использую ProgressMeter.jl для отслеживания хода вычислений (он показывает, сколько времени вы можете ожидать завершения вычислений). Так что вы можете написать:

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

и используйте map функцию, чтобы применить correlation функцию ко всем элементам r (в этом случае результат будет таким же, как для трансляции).

Наконец, ваши вычисления будут намного быстрее, если вы не будете использовать глобальные переменные в quadgk . Лучше передавать kappa and xi в качестве аргументов функции вот так:

 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)