JuliaIntervals / IntervalRootFinding.jl

Library for finding the roots of a function using interval arithmetic
https://juliaintervals.github.io/IntervalRootFinding.jl/
Other
127 stars 26 forks source link

Roots disappear #90

Closed dpsanders closed 5 years ago

dpsanders commented 6 years ago
f(xx) = ( (x, y) = xx; SVector( log(y/x) + 3x, y -2x) )

X = IntervalBox(-100..100, 2)
roots(f, X)

Gives a single root for X = IntervalBox(-10..10, 2) but no root for (-100..100)!

Kolaru commented 6 years ago

The culprit is the log function that can return an empty interval for intervals contained outside its domain. For example we have


julia> log(Interval(-0.2))
∅

Here is a small example showing the problem in the case you showed

julia> f(xx) = ( (x, y) = xx; SVector( log(y/x) + 3x, y -2x) )
f (generic function with 1 method)

julia> Y = IntervalBox(-100..100, 2)
[-100, 100] × [-100, 100]

julia> Y1, Y2 = bisect(Y)
([-100, -0.78125] × [-100, 100], [-0.78125, 100] × [-100, 100])

# Mid of the IntervalBox as found in the contractors
julia> m = IntervalBox(Interval.(mid(Y2, where_bisect)))
[49.2156, 49.2157] × [-0.78125, -0.78125]

julia> f(m)
2-element StaticArrays.SArray{Tuple{2},IntervalArithmetic.Interval{Float64},1,2}:
    ∅
 [-99.2127, -99.2126]

For X = IntervalBox(-10..10, 2), the empty intervals created just appear to never actually come from intervals containing the solution.

However, I have no idea how to solve the poblem rigorously.

dpsanders commented 6 years ago

Ouch. So with Bisection it should work I guess.

Now we are aware of this possibility, I guess we need to explicitly check for the possibility of an empty f(m) and bisect in that case.

dpsanders commented 6 years ago

This should also be helped by incorporating constraint propagation from IntervalConstraintProgramming.jl

Kolaru commented 6 years ago

Switching contractor would be delicate as the call f(m) appears in the deepest possible in the algorithm: this would imply carrying the information all the way back to nearly the higher level.

After thinking about it, I think we can safely throw an DomainError here, if we consider that this is the user responsibility to only search root where the function is defined.

dpsanders commented 6 years ago

The user can't know where the function will be defined; I think that's the package's job! For the moment simply checking if the result is empty and bisecting if so should be enough?

Kolaru commented 6 years ago

You are right I misunderstood you: we can do simple checks and bisect the interval (for some reason I thought you wanted to change the method used to the bisection method). I'll make a PR to solve this.