QuadGK evaluates function at upper bound about quadgk.jl HOT 4 OPEN

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)

from quadgk.jl.

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?

from quadgk.jl.

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

from quadgk.jl.

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

from quadgk.jl.