QuadGK evaluates function at upper bound

[urgomo]

This is probably related to #38, I am trying to integrate with quadgk a function from 1 to infinity which should be converging, and encounter this error (julia 1.8.5 and QuadGK v2.8.2).

using QuadGK
quadgk(x->1/x^1.1,1,Inf)
DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)

Stacktrace:
  [1] evalrule(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, a::Float64, b::Float64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(LinearAlgebra.norm))
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/evalrule.jl:37
  [2] adapt(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, segs::Vector{QuadGK.Segment{Float64, Float64, Float64}}, I::Float64, E::Float64, numevals::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, n::Int64, atol::Float64, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm))
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:48
  [3] do_quadgk(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, s::Tuple{Float64, Float64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:37
  [4] #46
    @ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:219 [inlined]
  [5] handle_infinities(workfunc::QuadGK.var"#46#47"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::var"#1#2", s::Tuple{Float64, Float64})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:110
  [6] quadgk(::Function, ::Float64, ::Vararg{Float64}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:218
  [7] quadgk(::Function, ::Float64, ::Vararg{Float64})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:216
  [8] quadgk(::Function, ::Int64, ::Vararg{Any}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:213
  [9] quadgk(::Function, ::Int64, ::Vararg{Any})
    @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:213
 [10] top-level scope
    @ In[2]:1

I think that the error is due to a call by quadgk of the integrand at the upper bound, which should not happen given the documentation of quadgk. https://docs.juliahub.com/QuadGK/hq5Ol/2.8.0/quadgk-examples/

This is confirmed by doing the change of variable done by quadgk manually (x=1+t/(1-t), dx = 1/(1-t)^2dt) as can be seen in the following code.

using QuadGK
function integrand(t)
    if t==1.0
        println(t)
    end

    return 1/(1+t/(1-t))^1.1 * 1/(1-t)^2
end
quadgk(integrand,0,1)
1.0
DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)

Stacktrace:
 [1] evalrule(f::typeof(integrand), a::Float64, b::Float64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(LinearAlgebra.norm))
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/evalrule.jl:37
 [2] adapt(f::typeof(integrand), segs::Vector{QuadGK.Segment{Float64, Float64, Float64}}, I::Float64, E::Float64, numevals::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, n::Int64, atol::Float64, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm))
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:48
 [3] do_quadgk(f::typeof(integrand), s::Tuple{Int64, Int64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:37
 [4] #46
   @ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:219 [inlined]
 [5] handle_infinities(workfunc::QuadGK.var"#46#47"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::typeof(integrand), s::Tuple{Int64, Int64})
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:118
 [6] quadgk(::Function, ::Int64, ::Vararg{Int64}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:218
 [7] quadgk(::Function, ::Int64, ::Vararg{Int64})
   @ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:216
 [8] top-level scope
   @ In[15]:1

[lxvm]

The change of variables only works reliably for 1/x^λ for λ>=2. If λ>1 the change of variables gives a singular integrand (with a converging integral) for which quadgk may or may not converge for a given λ depending on the error tolerances. Here is a figure showing what happens to the transformed problem

julia> function integrand(t, λ)
          return 1/(1+t/(1-t))^λ * 1/(1-t)^2
       end
integrand (generic function with 2 methods)

julia> using Plots

julia> plt = plot(; xguide="t", yguide="1/(1+t/(1-t))^λ * 1/(1-t)^2", xlims=(0,1), ylims=(0,2));

julia> for λ in [1.1, 1.5, 1.9, 2.0, 2.1, 2.5]; plot!(plt, range(0, 1, length=1000), t -> integrand(t, λ), label="λ=$λ"); end

julia> plt

So the fix is to apply a change of variables of the form x = y^a to your integrand that removes the endpoint singularity so that dx/x^λ = a*y^(a-1)*dy/y^(λ*a) = a*dy/y^((λ-1)*a+1) with (λ-1)*a+1 > 2 or a > 1/(λ-1). Here is your solution

julia> quadgk(y->10*y^9/(y^10)^1.1,1,Inf)
(9.999999999999993, 1.7763568394002505e-15)

[stevengj]

Basically, if you have a badly behaved integrand and request too low an error tolerance (e.g. the default is rtol ≈ 1e-8), then it may subdivide the integration integral closer and closer to the endpoint until roundoff errors cause it to evaluate at the endpoint, or close enough to the endpoint to cause an error..

It might be worth including a check for this (quadrature point == endpoint) in the rule evaluation, and if so stop refining (since at this point the floating-point precision will make it impossible to go further, and also we should never evaluate integrands exactly at endpoints).

However, in this case, it is evaluating at 0.9999999999999964, not exactly at the endpoint?

[lxvm]

I think 0.9999999999999964 is just the midpoint of the most refined interval, (0.9999999999999929, 1.0), and since the GK nodes cluster at the boundaries of the interval there is at least one node that rounds off to the endpoint

julia> using QuadGK

julia> x_end = -QuadGK.cachedrule(Float64, 7)[1][1] # outermost point in GK rule on [-1,1]
0.9914553711208126

julia> t_a, t_b = (0.9999999999999929, 1.0) # problem interval in transformed domain
(0.9999999999999929, 1.0)

julia> t_end = t_a + ((t_b - t_a)/2)*(1+x_end) # roundoff of outermost GK node to endpoint
1.0

[lxvm]

With #91 the result of the OP would become

julia> quadgk(x->1/x^1.1,1,Inf)
(9.772220920299873, 0.011179284983380486)

which has underestimated error