juliaintervals/intervalrootfinding.jl

Library for finding the roots of a function using interval arithmetic

IntervalRootFinding

IntervalRootFinding.jl is a Julia package for finding the roots of functions, i.e. solutions to the equation $f(x) = 0$. To do so, it uses interval arithmetic from the IntervalArithmetic library.

Documentation

The official documentation is available online: https://juliaintervals.github.io/IntervalRootFinding.jl/stable.

Installation

The IntervalRootFinding.jl package requires to install Julia.

Then, start Julia and execute the following command in the REPL:

using Pkg; Pkg.add("IntervalRootFinding")

Basic usage examples

The basic function is roots. Given a standard Julia function and an interval, the roots function returns a list of intervals containing all roots of the function located in the prescribed interval.

julia> using IntervalArithmetic, IntervalRootFinding

julia> f(x) = sin(x) - 0.1*x^2 + 1
f (generic function with 1 method)

julia> roots(f, -10..10)
4-element Array{Root{Interval{Float64}},1}:
 Root([3.14959, 3.1496], :unique)
 Root([-4.42654, -4.42653], :unique)
 Root([-1.08205, -1.08204], :unique)
 Root([-3.10682, -3.10681], :unique)

The :unique status indicates that each listed interval contains exactly one root. The other possible status is :unknown, which corresponds to intervals that may contain zero, one, or more roots (no guarantee is provided for these intervals).

These results are represented in the following plot, the region containing roots being in green.

intervalrootfinding.jl's Issues

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

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...

Limit the maximum number of intervals

Eg finding roots of sin on -Inf..Inf

Need absolute *and* relative tolerance

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.

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

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.

newton1d struggles with double roots

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

seems to hang.

cc @eeshan9815

roots hangs with double roots

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

seems to hang

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

Run femtocleaner

roots cannot handle functions that return intervals

Also seems to be true of newton1d.

cc @eeshan9815

Remove `Dual` / `Interval` promotion tests

These don't seem to actually ever be used.

Document new roots syntax and iterator interface.

• Syntax for roots
• Iterator interface
• Remove old docs

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

Use SVector for Newton

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

Break up `branch_and_prune.jl`

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

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

Change to work with AbstractInterval

Methods on IntervalBoxes should be moved to IntervalArithmetic

Cross-ref: JuliaIntervals/IntervalArithmetic.jl#158

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

Use queue instead of stack

Move lexcmp on Intervals to IntervalArithmetic.jl

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 old code

Remove old code:

• bisection.jl

• move out definitions of N used in branch_and_prune

• krawczyk.jl

• newton.jl

• find_roots stuff

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

Newton of BigFloat should iterate to completion

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.

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)

Clean up code

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

Add test with gradient

Done in #88

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

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.

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

Allow passing explicit derivative to Newton

Add debug mode

That gives information about what is going on.

Add centered form contractor

Done in #110.

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 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.

Implement Krawczyk in N dimensions

Tolerance in Newton is not respected

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

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.

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)

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)

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

Remove NewtonContractor

I believe we can remove NewtonContractor by just using

roots(f, f_prime, Newton)

Am I incorrect, @Kolaru?

Add hard test cases

There are some nice ones in

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

with results we can compare with.

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)

Return semi-infinite intervals when cannot guarantee roots

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

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)