Open urgomo opened 1 year ago
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)
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?
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
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
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).
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.