JuliaMath / QuadGK.jl

adaptive 1d numerical Gauss–Kronrod integration in Julia
MIT License
268 stars 37 forks source link

quadgk with identical lower and upper bound is broken #97

Closed bclyons12 closed 8 months ago

bclyons12 commented 8 months ago

I used to be able to do

julia> g(y) = quadgk(x -> x, y, 1.0)
g (generic function with 1 method)
julia> g(0.0)
(0.5, 0.0)
julia> g(1.0)
(0.0, 0.0)

Doing an integral where the upper and lower bound were identical gave the answer 0.0, as expected. After upgrading to v2.9.2, I get the error:

julia> g(1.0)
ERROR: DomainError with 1.0:
roundoff error detected near endpoint of the initial interval (1.0, 1.0)
Stacktrace:
  [1] throw_endpoint_error(x::Float64, a::Float64, b::Float64)
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:186
  [2] #check_endpoint_roundoff#48
    @ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:181 [inlined]
  [3] check_endpoint_roundoff
    @ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:179 [inlined]
  [4] #6
    @ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:15 [inlined]
  [5] ntuple
    @ ./ntuple.jl:48 [inlined]
  [6] do_quadgk(f::var"#5#6", s::Tuple{…}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:13
  [7] #51
    @ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:257 [inlined]
  [8] handle_infinities(workfunc::QuadGK.var"#51#52"{…}, f::var"#5#6", s::Tuple{…})
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:146
  [9] #quadgk#50
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:256 [inlined]
 [10] quadgk
    @ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:254 [inlined]
 [11] g(y::Float64)
    @ Main ./REPL[2]:1
 [12] top-level scope
    @ REPL[3]:1

which seems to be caused by https://github.com/JuliaMath/QuadGK.jl/pull/91 .

Is there anyway to restore the functionality to accept identical lower and upper bounds? The above function is just an MWE, but it's not uncommon to have an integral inside a function where the lower bound is an argument and can vary to the upper bound.

stevengj commented 8 months ago

Yes, this is a bug.

stevengj commented 8 months ago

I should have a fix merged and a new version tagged shortly, thanks!

bclyons12 commented 8 months ago

There's still an issue if the bounds are approximately machine precision apart. Perhaps you need isapprox instead of ==?

DomainError with 0.9999999999999996:
  roundoff error detected near endpoint of the initial interval (0.9999999999999996, 1.0)
bclyons12 commented 7 months ago

@stevengj I can't reopen this issue is seems. Let me know if you'd like me to make a new one

stevengj commented 7 months ago

Can you open a new one?

stevengj commented 7 months ago

Nevermind, I just filed a PR.