juliaintervals/intervalrootfinding.jl

Library for finding the roots of a function using interval arithmetic

intervalrootfinding.jl's Issues

BigFloats don't give enclosure

rts = roots(x->x^2 - 2, big(-∞..∞), 1e-70)

diam.(interval.(rts))

First diameter is 0, which is a bug, since it is not an enclosure of sqrt(2).

Clean up code

• Remove old stuff
• Split out contractors into new file
• Rename N (Newton operator) to 𝒩 for clarity

Krawczyk fails for functions with interval parameters

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)

Return semi-infinite intervals when cannot guarantee roots

e.g. for f(x) = (x-1)x(x+1) on (-Inf..Inf).

newton1d struggles with double roots

newton1d(x->(x^2-2)^2, -10..10)

seems to hang.

cc @eeshan9815

Use queue instead of stack

Make derivative another parameter of roots?

E.g.

roots(x->x^2 - 2, x->2x, X, tol)

This is how it used to be. This removes the need for keyword arguments and improves type inference.

Add centered form contractor

Done in #110.

Methods on IntervalBoxes should be moved to IntervalArithmetic

Cross-ref: JuliaIntervals/IntervalArithmetic.jl#158

Limit the maximum number of intervals

Eg finding roots of sin on -Inf..Inf

Use SVector for Newton

This should give a significant speed-up and less allocation.

Add hard test cases

There are some nice ones in

https://www.sciencedirect.com/science/article/pii/S0377042700007111

with results we can compare with.

Newton problem with Lagrange points

using StaticArrays

const G = 1.0
const ω = 1 
const μ1 = 0.3 # 1/2
const μ2 = one(μ1) - μ1
const xc1 = -μ2
const xc2 = μ1

fgx(x, y, xc, mu) = G * mu * (x - xc) / ((x - xc)^2+y^2)^(3/2)
fgy(x, y, xc, mu) = G * mu * y / ((x - xc)^2+y^2)^(3/2)

f1(X) = ( (x, y) = X; ω^2*x - fgx(x, y, xc1, μ1) - fgx(x, y, xc2, μ2) )
f2(X) = ( (x, y) = X; ω^2*y - fgy(x, y, xc1, μ1) - fgy(x, y, xc2, μ2) )

f(X) = SVector(f1(X), f2(X))

X = IntervalBox(-10..11, -10..11)

rts = roots(f, X, Bisection);

julia> roots(f, rts, Newton)
7-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([1.1232, 1.12321] × [0, 0], :unique)
 Root([-0.200001, -0.199999] × [0.866025, 0.866026], :unique)
 Root([0.299407, 0.300049] × [-0.000518799, 0.000122071], :unknown)
 Root([-0.28613, -0.286129] × [0, 0], :unique)
 Root([-0.700348, -0.699707] × [-0.000518799, 0.000122071], :unknown)
 Root([-0.200001, -0.199999] × [-0.866026, -0.866025], :unique)
 Root([-1.25674, -1.25673] × ∅, :unique)

The last root should have 0 (zero) instead of empty set...

Change to work with AbstractInterval

Bad behaviour with symmetric intervals

julia> f(xx) = ( (x,y) = xx; sin(x)*sin(y) )
f (generic function with 1 method)

julia> rts = roots(∇(f), IntervalBox(-5..5, 2), Newton)
1-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([0, 0] × [0, 0], :unique)

Make Newton work with infinite intervals

julia> roots(x->x^2 - 2, -Inf..Inf)
1-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([-∞, ∞], :unique)

This seems to be because N(X) = X

Repeating Newton gives :unknown

It should keep track of the :unique information.

julia> rts = roots(x->x^2 - 2, -10..10)
2-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root(Interval(1.414213562373095, 1.4142135623730951), :unique)
 Root(Interval(-1.4142135623730951, -1.414213562373095), :unique)

julia> roots(x->x^2 - 2, rts)
2-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root(Interval(1.414213562373095, 1.4142135623730951), :unknown)
 Root(Interval(-1.4142135623730951, -1.414213562373095), :unknown)

Newton fails for dummy example with continuous solutions

While exploring #100, I found that Newton fails in the dummy example presented there:

julia> f(xx) = SVector(xx[1] + xx[2], xx[1] + xx[2])
f (generic function with 1 method)

julia> roots(f, IntervalBox(-1..1, 2), Newton)
0-element Array{Root{IntervalBox{2,Float64}},1}

There are obvious solution in that interval box though, (0, 0) for example.

Remove old code

Remove old code:

• bisection.jl

• move out definitions of N used in branch_and_prune

• krawczyk.jl

• newton.jl

• find_roots stuff

Problems when using current master of IntervalArithmetics

Working on #178 with current master of IntervalArithmetics.jl tests are broken. This is is related with JuliaIntervals/IntervalArithmetic.jl#162, but after that, it is still yielding errors that look like this:

2D roots: Error During Test
  Got an exception of type MethodError outside of a @test
  MethodError: no method matching ForwardDiff.Dual{ForwardDiff.Tag{#f#6,IntervalArithmetic.IntervalBox{2,Float64}},V,N} where N where V<:Real(::IntervalArithmetic.IntervalBox{2,Float64}, ::ForwardDiff.Chunk{1}, ::Val{1})
  Stacktrace:
   [1] macro expansion at /Users/benet/.julia/v0.6/ForwardDiff/src/apiutils.jl:26 [inlined]
   [2] dualize at /Users/benet/.julia/v0.6/ForwardDiff/src/apiutils.jl:22 [inlined]
   [3] static_dual_eval at /Users/benet/.julia/v0.6/ForwardDiff/src/apiutils.jl:30 [inlined]
   [4] vector_mode_jacobian at /Users/benet/.julia/v0.6/ForwardDiff/src/jacobian.jl:172 [inlined]
   [5] jacobian at /Users/benet/.julia/v0.6/ForwardDiff/src/jacobian.jl:81 [inlined]
   [6] #40 at /Users/benet/.julia/v0.6/IntervalRootFinding/src/roots.jl:146 [inlined]
   [7] 𝒩(::#f#6, ::IntervalRootFinding.##40#41{#f#6}, ::IntervalArithmetic.IntervalBox{2,Float64}) at /Users/benet/.julia/v0.6/IntervalRootFinding/src/contractors.jl:105
   [8] newtonlike_contract(::IntervalRootFinding.#𝒩, ::IntervalRootFinding.Newton{#f#6,IntervalRootFinding.##40#41{#f#6}}, ::IntervalArithmetic.IntervalBox{2,Float64}, ::Float64) at /Users/benet/.julia/v0.6/IntervalRootFinding/src/contractors.jl:47
...

Any idea where is the problem?

Add debug mode

That gives information about what is going on.

Integration with Symbolic packages

Currently calling lambdify on a symbolic expression does not work as nicely as explicitly defining a function

I believe this is why I got the following error:

DomainError:
log will only return a complex result if called with a complex argument. Try log(complex(x)).

From,

• https://discourse.julialang.org/t/way-to-find-positive-roots-by-treating-complex-numbers-as-negatives-symengine-roots/11688/10?u=djsegal

Incorrect result for newton1d(x->sin(x)-x)

This has a triple (?) root at 0.

julia> newton1d(x -> sin(x) - x, -10..10)
4-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([-2.55949e-08, 9.21914e-10], :unique)
 Root([9.21913e-10, 1.0455e-08], :unique)
 Root([1.04549e-08, 1.99881e-08], :unique)
 Root([2.74386e-08, 3.41862e-08], :unique)

This result is incorrect: there are not unique roots in those intervals.

cc @eeshan9815

Krawczyk fails when the Jacobian is singular at the center of starting region

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

while

julia> roots(x -> x^2 - 2, -4..4.0001, Krawczyk)
2-element Array{Root{Interval{Float64}},1}:
 Root([1.41389, 1.41454], :unique)
 Root([-1.41445, -1.41398], :unique)

It can be hotfixed by shifting the point where the Krawczyk operator take the Jacobian (similar to what is done when bisecting an interval), but it can also be more seriously fixed by implementing more advanced Krawczyk methods described in the litterature and that claim to avoid the problem altogether (and also faster convergence).

Document new roots syntax and iterator interface.

• Syntax for roots
• Iterator interface
• Remove old docs

Add test with gradient

Done in #88

Need absolute *and* relative tolerance

Cannot use complex functions

julia> roots(sin, Complex(-5..5, -5..5))
ERROR: MethodError: no method matching sin(::IntervalRootFinding.Compl{IntervalArithmetic.Interval{Float64}})
Closest candidates are:
  sin(::Float64) at math.jl:419
  sin(::Float32) at math.jl:420
  sin(::Float16) at math.jl:950

There was an issue that meant that standard Julia complex numbers could not be used, hence the reimplementation in complex.jl as Compl that is cited here.

roots cannot handle functions that return intervals

Also seems to be true of newton1d.

cc @eeshan9815

find_roots_midpoint should gracefully handle cases with no roots

julia> using IntervalRootFinding

julia> f(x)=x^0
f (generic function with 1 method)

julia> IntervalRootFinding.find_roots_midpoint(f, -5, 5)
ERROR: BoundsError: attempt to access 0-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1} at index [1]
Stacktrace:
 [1] getindex(::Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}, ::Int64) at ./array.jl:520
 [2] #find_roots_midpoint#61(::Float64, ::Bool, ::Int64, ::Int64, ::Function, ::Function, ::Int64, ::Int64, ::IntervalRootFinding.#newton) at /home/pkm/.julia/v0.6/IntervalRootFinding/src/IntervalRootFinding.jl:107
 [3] find_roots_midpoint(::Function, ::Int64, ::Int64) at /home/pkm/.julia/v0.6/IntervalRootFinding/src/IntervalRootFinding.jl:105

julia> IntervalRootFinding.find_roots(f, -5, 5)
0-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}

It shouldn't throw an error, but return an empty result. Or error with a better error :)

Remove NewtonContractor

I believe we can remove NewtonContractor by just using

roots(f, f_prime, Newton)

Am I incorrect, @Kolaru?

Add new root status when the presence (but not uniqueness) of a solution is known

For Krawczyk and Newton method, the fact that a contracted interval is subset of the original interval is a sufficient condition for the presence of a solution.

The proposition here is to add a new root status :exist that correspoond to this information.

This can be usefull when solving for equations with a continuous set of solution, for which this situation is common, it allows to prove the presence of a solution in some region.

See a dummy example below.

julia> import ForwardDiff: jacobian

julia> using StaticArrays

julia> f(xx) = SVector(xx[1] + xx[2], xx[1] + xx[2])
f (generic function with 1 method)

julia> rts = roots(f, IntervalBox(-1..1, 2), Krawczyk)
4559-element Array{Root{IntervalBox{2,Float64}},1}:
 Root([0.00721069, 0.00813498] × [-0.0078125, -0.00690255], :unknown)
 Root([-0.0078125, -0.00690255] × [0.00721069, 0.00813498], :unknown)
 Root([-0.000359588, 0.000564693] × [-0.000359588, 0.000564693], :unknown)
 Root([0.00339598, 0.00433482] × [-0.00411516, -0.00319087], :unknown)
 Root([0.00528843, 0.00624206] × [-0.00597828, -0.00505399], :unknown)
 Root([0.00624205, 0.0072107] × [-0.00690256, -0.00597827], :unknown)
 Root([0.00624205, 0.0072107] × [-0.0078125, -0.00690255], :unknown)
 Root([0.00528843, 0.00624206] × [-0.00690256, -0.00597827], :unknown)
 Root([0.00433481, 0.00528844] × [-0.005054, -0.00411515], :unknown)
 Root([0.00433481, 0.00528844] × [-0.00597828, -0.00505399], :unknown)
 ⋮
 Root([-0.324746, -0.323745] × [0.32355, 0.324535], :unknown)
 Root([-0.323243, -0.32273] × [0.322566, 0.323551], :unknown)
 Root([-0.323746, -0.323242] × [0.322566, 0.323551], :unknown)
 Root([-0.325745, -0.324745] × [0.324534, 0.325535], :unknown)
 Root([-0.332648, -0.331663] × [0.332452, 0.332957], :unknown)
 Root([-0.331664, -0.330664] × [0.331453, 0.332453], :unknown)
 Root([-0.332648, -0.331663] × [0.331453, 0.332453], :unknown)
 Root([-0.333632, -0.332647] × [0.332956, 0.333468], :unknown)
 Root([-0.333632, -0.332647] × [0.332452, 0.332957], :unknown)
 Root([-0.331664, -0.330664] × [0.330453, 0.331454], :unknown)

julia> X = interval(rts[1])
[0.00433481, 0.00528844] × [-0.00597828, -0.00505399]

julia> C = Krawczyk(f, x -> jacobian(f, x))
Krawczyk{typeof(f),getfield(Main, Symbol("##9#10"))}(f, getfield(Main, Symbol("##9#10"))())

julia> status, CX = C(X, 1e-15)

julia> issubset(CX, X)
true

PS: I'm actually volunteer to implement it when I found a bit of time, but I prefer to first raise it as an issue as I suspect I may have overlooked something preventing this.

Move lexcmp on Intervals to IntervalArithmetic.jl

Example 5.5 from Smiley and Chun (2001) solved inconsistently

The example in test/smiley_examples.jl has 41 roots, however, only 20 are found on some setups. On other setups all 41 are identified correctly. This inconsistency was first discovered and discussed in #70.

Implement Krawczyk in N dimensions

Remove `Dual` / `Interval` promotion tests

These don't seem to actually ever be used.

cannot run basic example on readme

hi there,

I'm getting

               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: https://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.6.0 (2017-06-19 13:05 UTC)
 _/ |\__'_|_|_|\__'_|  |  Official http://julialang.org/ release
|__/                   |  x86_64-apple-darwin13.4.0

julia> using IntervalArithmetic, IntervalRootFinding

julia> rts = roots(x->x^2 - 2, -10..10, Bisection)
ERROR: UndefVarError: roots not defined

julia> rts = IntervalRootFinding.roots(x->x^2 - 2, -10..10, Bisection)
ERROR: UndefVarError: roots not defined

julia> 

cheers

Newton of BigFloat should iterate to completion

Newton fails with abs

julia> f(p) = abs(1 / ( (1+p)^30 ) * 10_000 - 100)
f (generic function with 1 method)

julia> roots(f, -1000..1000, Newton)
2-element Array{Root{Interval{Float64}},1}:
 Root(∅, :unique)
 Root(∅, :unique)

roots should not return "unique" for multiple roots

Seems to be a special case when the root is exact:

julia> f5(x) = (x^2 - 1)^4 * (x^2 - 2)^4             #two quadruple roots
f5 (generic function with 1 method)

julia>

julia> roots(f5, -10..10)
4-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([1.41421, 1.41422], :unknown)
 Root([1, 1], :unique)
 Root([-1, -1], :unique)
 Root([-1.41422, -1.41421], :unknown)

Break up `branch_and_prune.jl`

There is far too much going on in this file. Split it up into pieces.

Allow passing explicit derivative to Newton

Roots disappear

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)!

(Complex?) bisection doesn't work with BigFloat

Seems to be working fine in the new version.

julia> x = -5..6
[-5, 6]

julia> Xc = Complex(big(x), big(x))
[-5, 6]₂₅₆ + [-5, 6]₂₅₆*im

julia> rts = roots(f, Xc, Newton)
3-element Array{IntervalRootFinding.Root{Complex{IntervalArithmetic.Interval{BigFloat}}},1}:
 Root([0.999999, 1.00001]₂₅₆ + [-1.10032e-19, 1.10551e-19]₂₅₆*im, :unique)
 Root([-0.500001, -0.499999]₂₅₆ + [0.866025, 0.866026]₂₅₆*im, :unique)
 Root([-0.500001, -0.499999]₂₅₆ - [0.866025, 0.866026]₂₅₆*im, :unique)

Tolerance in Newton is not respected

The tolerance setting doesn't seem to be respected by Newton.

Always use autodiff for derivatives / gradients in Newton

This will simplify the interface, at a very small cost in efficiency.
(And it's surely pretty rare that people will really want to give the derivative by hand -- and they shouldn't be doing so.)

Run femtocleaner

Bisection of atomic interval

One solution for the bisection of an atomic interval [l, h] would be to return three intervals:

• the point l
• the interval [l, h]
• the point h

Probably the easiest solution is not to bisect it.
Note that this will cause problems since currently it is expected that bisect returns two intervals.
Could return the same interval twice, together with an :atomic flag.

roots hangs with double roots

 rts = roots(x->exp(x^2) - cos(x), -10..10)

seems to hang

Atomic intervals should be treated specially

If an Interval is atomic (unable to be bisected) it should be marked as unknown and not processed further.

This will be important for infinite intervals if the function is such that we cannot exclude the interval (prevfloat(Inf), Inf).