Extra roots bug for large intervals

[dpsanders]

https://discourse.julialang.org/t/question-about-interval-root-finding/57124

f(E) = tan(√(2E)) + √(E / (3π^2 - E))

roots(f, -1e100..1e100)

is buggy.

[lucaferranti]

Also, now that I look at the function better E=0 should also be a root, but that goes undetected with Newton 🤔. Krawczyk and Bisection can find it at the cost of a couple of extra roots.

julia> roots(f, 0..30, Bisection)
5-element Vector{Root{Interval{Float64}}}:
 Root(Interval(1.2337005023828356, 1.233700557778276), :unknown)
 Root(Interval(15.064561967948055, 15.064562025101973), :unknown)
 Root(Interval(3.844642092262015, 3.8446421485298243), :unknown)
 Root(Interval(11.103304945220835, 11.103305008980778), :unknown)
 Root(Interval(0.0, 8.972354260755779e-8), :unknown)

julia> roots(f, 0..30, Krawczyk)
5-element Vector{Root{Interval{Float64}}}:
 Root(Interval(15.064562002776027, 15.064562002781207), :unique)
 Root(Interval(11.103304945220835, 11.103305008980778), :unknown)
 Root(Interval(1.2337005023828356, 1.233700557778276), :unknown)
 Root(Interval(3.8446421166049274, 3.844642116613573), :unique)
 Root(Interval(0.0, 8.972354260755779e-8), :unknown)

Also, those extra intervals contain or are close to vertical asymptotes of the function (pi^2/8 and 9/8*pi^2), which I guess explains why the algorithm cannot through them away.

[dpsanders]

The root at 0 cannot be proved because there is not an open interval around it, since sqrt is undefined for negative numbers.

[dpsanders]

But the root at 0 should still be reported as :unknown. That seems like a serious bug! I am wondering if that is what is giving the empty interval.

[lucaferranti]

yeah I also wondered whether there was a relation between the returned empty set and the root at zero. It seems that starting with "small" intervals in the form [0, a] the root 0 is correctly found as :unknown , for bigger intervals the root at zero remains undetected. Furthermore, if 0 is in the interior of the starting interval, then sometimes an empty set is returned.

julia> f(E) = tan(sqrt(2*E)) + sqrt(E/(3*pi^2-E))
f (generic function with 1 method)

julia> roots(f, 0..1)
1-element Vector{Root{Interval{Float64}}}:
 Root([0, 9.95329e-08], :unknown)

julia> roots(f, 0..2)
2-element Vector{Root{Interval{Float64}}}:
 Root([0, 9.87553e-08], :unknown)
 Root([1.2337, 1.23371], :unknown)

julia> roots(f, 0..5)
3-element Vector{Root{Interval{Float64}}}:
 Root([1.2337, 1.23371], :unknown)
 Root([0, 6.07614e-08], :unknown)
 Root([3.84464, 3.84465], :unique)

julia> roots(f, 0..10)
1-element Vector{Root{Interval{Float64}}}:
 Root([3.84464, 3.84465], :unique)

julia> roots(f, 0..30)
2-element Vector{Root{Interval{Float64}}}:
 Root([15.0645, 15.0646], :unique)
 Root([3.84464, 3.84465], :unique)

julia> roots(f, 0..100)
2-element Vector{Root{Interval{Float64}}}:
 Root([3.84464, 3.84465], :unique)
 Root([15.0645, 15.0646], :unique)

julia> roots(f, -1..1)
1-element Vector{Root{Interval{Float64}}}:
 Root(∅, :unique)

julia> roots(f, -2..2)
2-element Vector{Root{Interval{Float64}}}:
 Root([1.2337, 1.23371], :unknown)
 Root(∅, :unique)

julia> roots(f, -3..3)
2-element Vector{Root{Interval{Float64}}}:
 Root([1.2337, 1.23371], :unknown)
 Root(∅, :unique)

julia> roots(f, -5..5)
3-element Vector{Root{Interval{Float64}}}:
 Root([3.84464, 3.84465], :unique)
 Root([1.2337, 1.23371], :unknown)
 Root(∅, :unique)

julia> roots(f, -10..10)
1-element Vector{Root{Interval{Float64}}}:
 Root([3.84464, 3.84465], :unique)

julia> roots(f, -100..100)
2-element Vector{Root{Interval{Float64}}}:
 Root([3.84464, 3.84465], :unique)
 Root([15.0645, 15.0646], :unique)

julia> roots(f, -1e100..1e100)
2-element Vector{Root{Interval{Float64}}}:
 Root([15.0645, 15.0646], :unique)
 Root([3.84464, 3.84465], :unique)

julia> roots(f, -1e200..1e200)
4-element Vector{Root{Interval{Float64}}}:
 Root([15.0645, 15.0646], :unique)
 Root([11.1033, 11.1034], :unknown)
 Root(∅, :unique)
 Root([3.84464, 3.84465], :unique)

[Kolaru]

I suspect the empty interval appears because it is a subset of any interval, so some inclusion test somewhere does not work as intended (in many places empty intervals are explicitly special cased, otherwise this example would be way more broken ^^').

This together with JuliaIntervals/IntervalArithmetic.jl#441 makes me wonder if we should use decorated itnerval as default, or at least convert everything to decorated interval in IntervalRootFinding to catch that kind of problems.

Moreover, looking at the graph of the function, it even may be that Newton should not be able to exclude the singularities, since the image of any interval containing one of the singularity will contain zero. The function may violate some of the hypothesis of the Newton method we take for granted here. Decorations could help for that as well.