I am writing a program to perform a spherical Bessel transformation which requires some one dimensional integrals but I am having some trouble with the Integrals.jl package.
Naturally my first attempt to do this was to use QuadGK.jl but I also wanted to try using the Integrals.jl package since I might want to make use of it’s in place or batching features in the future. I am aware that the QuadGK.jl package also has in place features with quadgk! but even then I’d prefer to understand what I am doing wrong.
using QuadGK
using Integrals
using SpecialFunctions
using Plots
When I do this with the QuadGK package my code looks as follows:
function spherical_bessel_transform(f::Function, root::Float64, nu::Int64)::Float64
norm = 2 / sphericalbesselj(nu + 1, root)^2;
func = (x) -> x^2 * f(x) * sphericalbesselj(nu, root * x);
return norm * quadgk(func, 0, 1; atol=1e-14, rtol=1e-14)[1];
end
When I do this with the Integrals package my code looks as follows:
function spherical_bessel_transform(f::Function, root::Float64, nu::Int64)::Float64
norm = 2 / sphericalbesselj(nu + 1, root)^2;
func = (x, p) -> x^2 * f(x) * sphericalbesselj(nu, root * x);
prob = IntegralProblem(func, 0, 1; abstol = 1e-14, reltol = 1e-14);
sol = solve(prob, QuadGKJL());
return norm * sol.u;
end
I plot both of these methods for a test function as follows:
test_func = (x) -> exp(-x^2 * 10000);
test_roots = collect(1.0:500.0) * pi
plot(spherical_bessel_transform.(test_func, test_roots, 0))
They should both give the same results but the result using “Integrals.jl” does not look right.
I do not understand why these would give different results since they are supposedly using the same machinery under the hood.