JuliaIntervals / IntervalRootFinding.jl

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

Empty interval results on discontinuous functions #184

Open Timeroot opened 1 year ago

Timeroot commented 1 year ago

When given discontinuous functions, the results mess up, badly. For instance, this test function gives an interval labelled :unique that definitely does not contain a root:

function f( x )
  return x + sign(x)
end

X = -5 .. 5

rt = roots(f, X)
println(rt)
println("f returns: ",f(rt[1].interval))

gives

Root{Interval{Float64}}[Root([0.480468, 0.480469], :unique)]
f returns: [1.48046, 1.48047]

so evaluating on the function itself makes clear that there is no root here. Things get even sillier in 2D, where I get an empty set that is claimed to contain a root! 😄

function f( (x, y) )
  return SVector(
    x + 0.1 * sign(sin(10*y)),
    x - 0.1 * sign(cos(10*y)) + y - 0.5,
  )
end

X = -5 .. 5

rt = roots(f, X × X)

returns

6-element Vector{Root{IntervalBox{2, Float64}}}:
 Root(∅², :unique)
 Root([-0.0283185, -0.0283184] × [0.628318, 0.628319], :unknown)
 Root([0.1, 0.100001] × [0.5, 0.5], :unique)
 Root([-0.0283186, -0.0283184] × [0.628318, 0.628319], :unknown)
 Root([-0.100001, -0.1] × [0.699999, 0.700001], :unique)
 Root([0.1, 0.100001] × [0.471238, 0.471239], :unknown)

Now I'm not expecting great results exactly, as discontinuous functions are of course very hard to say anything meaningful about -- the fixed-point theorem doesn't apply, of course. But maybe some kind of sanity-checking (that contains_zero(f(X)), at least?) would be helpful. And if that sanity check fails, to give a warning to the user that their functions appear nondifferentiable and problematic.

This is similar in spirit to #146 , but that stems from max providing incorrect derivative information. I don't know if this "issue" is firmly just out of scope. :)

dpsanders commented 1 year ago

Decorated intervals are designed to solve this problem: they will tell you that "something strange happened during the computation".