Discrepancies with `QuadGK`

[arthurmloureiro]

Hi @xzackli !

I continue to push a bit your package, I hope this is ok!

I am doing a "classic" weak lensing kernel integration where I have an interpolated kernel: which I want to integrate as $\Psi{\ell}(k) = \int \rm{d}\chi j{\ell}(k\chi)\psi^i(\chi)$

When I compare the results using OscillatoryIntegralsODE.jl with what I get simply using QuadGK I get weird results or high-values of $\ell$ I get several warnings and weird results:

┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/gTrkJ/src/integrator_interface.jl:536

Here are the functions I am using to generate this results:

QuadGK

function lensing_pot_integrand(χ::Float64, k::Float64, ℓ::Int64)
    return sphericalbesselj(ℓ, k * χ) * kernel_spline(χ)
end

function Ψℓk(k_::Float64, ℓ_::Int64, χ₁::Float64, χ₂::Float64)
    return quadgk(χ -> lensing_pot_integrand(χ, k_, ℓ_), χ₁, χ₂)[1]
end

function Ψℓk(k_::Vector{Float64}, ℓ_::Int64, χ₁::Float64, χ₂::Float64)
    results = similar(k_)
    Threads.@threads for i in 1:length(k_)
        results[i] = Ψℓk(k_[i], ℓ_, χ₁, χ₂)
    end
    return results
end

# ell=2 is in the first plot:
test_quad = Ψℓk(klog, 2, first(chi_of_z), last(chi_of_z))
# ell = 20 is in the second:
test_quad = Ψℓk(klog, 20, first(chi_of_z), last(chi_of_z))

Levi Integrate

function Ψℓk_levi(k_vec::Vector{Float64}, dχ::Vector{Float64}, ℓ::Int64)
    results = similar(k_vec)
    @inbounds Threads.@threads for i in 1:length(k_vec)
        kk = k_vec[i]
        jl = SphericalBesselJIntegral{Int64, Float64}(ℓ, kk)
        results[i] = levintegrate(jl, kernel_spline, first(dχ), last(dχ))
    end
    return results
end

# ell=2 is in the first plot:
test_levi = Ψℓk_levi(klog, chi_of_z, 2)
# ell = 20 is in the second:
test_levi = Ψℓk_levi(klog, chi_of_z, 20)

Am I doing something obviously wrong?

Cheers!

Edit: The low-ell problem was a bug on my side!

[arthurmloureiro]

Updating on some tests following a suggestion by @ntessore.

klog = np.logspace(np.log10(1e-3), np.log10(1e0), 150);
χ_test = reverse(2*π ./ klog)

# quadGK
f(x) = x*exp(-x)

function integrand(k::Float64, x::Float64, l::Int64)
    return sphericalbesselj(l, k * x) * f(x)
end

function integral(k::Float64, x::Vector{Float64}, l::Int64)
    return quadgk(x -> integrand(k, x, l), first(x), last(x))[1]
end

function integral(k::Vector{Float64}, x::Vector{Float64}, l::Int64)
    results = similar(k)
    Threads.@threads for i in 1:length(k)
        results[i] = integral(k[i], x, l)
    end
    return results
end

#Levi     
function integral_levi(k::Vector{Float64}, x::Vector{Float64}, l::Int64)
    results = similar(k)
    l_float = convert(Float64, l)
    Threads.@threads for i in 1:length(k)
        kk = k[i]
        jl = SphericalBesselJIntegral{Float64}(l_float, kk)
        results[i] = levintegrate(jl, f, first(x), last(x))
    end
    return results
end

int_quad = integral(klog,  χ_test, 0);
int_levi = integral_levi(klog, χ_test, 0);

# exact
F(k, a, b) = (exp(-a)*(k*cos(a*k)+sin(a*k)) - exp(-b)*(k*cos(b*k)+sin(b*k)))/(k+k^3)
res_exact = F.(klog, first(χ_test), last(χ_test))

plot(klog, int_quad, label="QuadGK")
plot!(klog, int_levi, label="Levi Method")
plot!(klog, res_exact, label="Exact")
plot!(xaxis=:log10, legend=:bottomleft)

But with higher ell:

int_quad = integral(klog,  χ_test, 10);
int_levi = integral_levi(klog, χ_test, 10);

plot(klog, int_quad, label="QuadGK")
plot!(klog, int_levi, label="Levi Method")
plot!(xaxis=:log10, legend=:bottomleft)