JuliaIntervals / IntervalRootFinding.jl

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

Shouldn't this return a unique root? #141

Closed newptcai closed 4 years ago

newptcai commented 4 years ago

Try this.

julia> roots(x -> x^2-(1-cos(x)), -1/2..1/2)
7-element Array{Root{Interval{Float64}},1}:
 Root([1.4021e-07, 1.41988e-07], :unknown)   
 Root([-3.11849e-08, 2.56487e-08], :unknown) 
 Root([-8.71372e-08, -3.11848e-08], :unknown)
 Root([-3.76699e-05, -3.76244e-05], :unknown)
 Root([1.03555e-06, 1.09957e-06], :unknown)  
 Root([2.56486e-08, 8.24821e-08], :unknown)  
 Root([5.84386e-07, 5.8855e-07], :unknown)  
Kolaru commented 4 years ago

Derivative at the root is zero, therefore the method can't determine the status of the interval containing the root (that's an absolute limitation derived from the theory), and it may happen that some intervals close by also end up with unknown status (that's a numeric limitation coming from finite precision arithmetic).

Both of these are expected.

Best we could do here is to provide a hull function for arbitrary number of roots/intervals that would allow to easily fuse these intervals back together (it can't be done automatically though, since the returned roots are not contiguous).

newptcai commented 4 years ago

Cannot we check (automatically) that the 2nd derivative is positive (2 in this case), so there can be only one unique root near zero?

Kolaru commented 4 years ago

Checking the second derivative won't help, as far as I can tell.

Consider the following example:

julia> roots(x -> (x^2 - 1e-16) - (1 - cos(x)), -0.5 .. 0.5)
7-element Array{Root{Interval{Float64}},1}:
 Root([1.4021e-07, 1.42267e-07], :unknown)   
 Root([-3.11849e-08, 2.56487e-08], :unknown) 
 Root([-8.71372e-08, -3.11848e-08], :unknown)
 Root([-3.76699e-05, -3.76244e-05], :unknown)
 Root([1.03555e-06, 1.09961e-06], :unknown)  
 Root([2.56486e-08, 8.24821e-08], :unknown)  
 Root([5.84386e-07, 5.88617e-07], :unknown)  

The second derivative is still 2, but the number of solutions is not that clear (my guess is that there are two). Modifying it slightly further, we have

julia> roots(x -> (x^2 - 1e-16) + (1 - cos(x)), -0.5 .. 0.5)
1-element Array{Root{Interval{Float64}},1}:
 Root([-3.11849e-08, 2.56487e-08], :unknown)

where there are now definitively two distinct solutions, one between -1e-8 and 0, and the other between 0 and 1e-8 (the other roots disappears though, since the function is less flat at 0 than the original one).

dpsanders commented 4 years ago

That would tell you that there is either no root, one double root, or two roots. There is no way using interval methods to prove that there is a double root. You can prove that there are no roots or 2 unique distinct roots. Basically there is no way to check exact equality using interval methods (or indeed any other purely numerical method). (Unless the zero occurs at an exactly representable floating point number.)

dpsanders commented 4 years ago

You should be able to use BigFloat to solve the perturbed case

Kolaru commented 4 years ago

True

julia> roots(x -> (x^2 - 1e-16) - (1 - cos(x)), big(-0.5 .. 0.5), Newton, 1e-12)

2-element Array{Root{Interval{BigFloat}},1}:
 Root([1.41421e-08, 1.41422e-08]₂₅₆, :unique) 
 Root([-1.41423e-08, -1.4142e-08]₂₅₆, :unique)