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

Krawczyk fails for functions with interval parameters #99

Closed Kolaru closed 5 years ago

Kolaru commented 5 years ago
julia> a = 2..2.1
[2, 2.10001]

julia> roots(x -> x^2 - a, -4..4, Krawczyk)
0-element Array{Root{Interval{Float64}},1}

Newton seems to be fine

julia> roots(x -> x^2 - a, -4..4, Newton)
2-element Array{Root{Interval{Float64}},1}:
 Root([1.4141, 1.44946], :unique)
 Root([-1.44946, -1.4141], :unique)
dpsanders commented 5 years ago

Thanks for the report.

Krawczyk hasn't received much attention. Presumably it's supposed to work for cases like this?

Bisection also does badly, since it ends up bisecting the interval that is the solution :/

Kolaru commented 5 years ago

The problem is not the interval parameters, but the fact that our favorite example currently fails with Krawczyk, because the derivative is 0 at the middle point:

julia> roots(x -> x^2 - 2, -4..4, Krawczyk)
0-element Array{Root{Interval{Float64}},1}

The reason is given in the new issue #102.

If the interval is not symmetric, it does just fine:

julia> roots(x -> x^2 - (2..2.1), -4..4.0001, Krawczyk)
2-element Array{Root{Interval{Float64}},1}:
 Root([1.4141, 1.44947], :unique)
 Root([-1.44947, -1.4141], :unique)

Theoretically it is supposed to work fine in this case because for an interval A, the interval extension of x^2 - A is also a (bad) interval extension for x^2 - a for any a in A. Then, all methods may bisect the solution if they are unable to conclude that it exists for all a in A. In that case it is the parameter interval A that should be bisected.

On the up side, I need to do this, so I will be implementing an algorithm doing exactly this tonight :D

dpsanders commented 5 years ago

This seems to be fixed now:

julia> a = 2..2.1
[2, 2.10001]

julia> roots(x -> x^2 - a, -4..4, Krawczyk)
2-element Array{Root{Interval{Float64}},1}:
 Root([1.41409, 1.44947], :unique)
 Root([-1.44946, -1.4141], :unique)
dpsanders commented 5 years ago

Needs a test.

Kolaru commented 5 years ago

I've added it to the (ever growing) list of missing tests in #115.