search

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

Reciprocal of reciprocal misses root #131

Closed jwscook closed 5 years ago

jwscook commented 5 years ago

(Unnecessarily) taking the reciprocal of a reciprocal causes the package to miss a root:

julia> using IntervalArithmetic, IntervalRootFinding

julia> f(x) = (1-x)*(1+x)
f (generic function with 1 method)

julia> roots(x->f(x), interval(-2, 2))
2-element Array{Root{Interval{Float64}},1}:
 Root([0.999999, 1.00001], :unique)  
 Root([-1.00001, -0.999999], :unique)

julia> roots(x->1/(1/f(x)), interval(-2, 2))
1-element Array{Root{Interval{Float64}},1}:
 Root([-1.00001, -0.999999], :unknown)

Is this behaviour expected?

  [d1acc4aa] IntervalArithmetic v0.15.0
  [d2bf35a9] IntervalRootFinding v0.4.0
  [8197267c] IntervalSets v0.3.1
  [d8418881] Intervals v0.3.3
Version 1.0.4-pre.0 (2018-12-19)
release-1.0/5b7e8d9d4e (fork: 362 commits, 334 days)
dpsanders commented 5 years ago

1 / f(x) has poles (infinities) at plus and minus 1, so 1 / (1 / f(x)) is not really even defined there?

dpsanders commented 5 years ago

But maybe it's true that it should catch that "something funny happens" at each pole.

dpsanders commented 5 years ago

Actually on Julia 1.1 I get

julia> roots(x -> 1 / ( 1 / f(x)), -10..10)
3-element Array{Root{Interval{Float64}},1}:
 Root([0.999999, 1.00001], :unknown)
 Root([-1, -0.999999], :unknown)
 Root([-1.00001, -1], :unknown)
dpsanders commented 5 years ago

(With master (?) branches of the interval packages, if that makes a difference.)

Kolaru commented 5 years ago

The expected result is indeed to have the roots status to be :unkown.

With the stable versions of the packages (IntervalArithmetic 0.15.2 and IntervalRootFinding 0.4.0) and Julia 1.1, I get correct result too.

Anyway, it would not hurt to have this in the test suite.

jwscook commented 5 years ago

Ok, great!

dpsanders commented 5 years ago

Can you confirm that the behaviour is correct with Julia 1.1? Is there a reason you reopened the issue?

jwscook commented 5 years ago

One 1.1 on my mac I get

julia> using IntervalArithmetic, IntervalRootFinding

julia> f(x) = (1-x)*(1+x)
f (generic function with 1 method)

julia> roots(x->f(x), interval(-2, 2))
2-element Array{Root{Interval{Float64}},1}:
Error showing value of type Array{Root{Interval{Float64}},1}:
ERROR: UndefVarError: to_mpfr not defined
Stacktrace:
 [1] getproperty at ./sysimg.jl:13 [inlined]
 [2] round_string(::BigFloat, ::Int64, ::RoundingMode{:Down}) at /Users/james/.julia/packages/IntervalArithmetic/iO2gf/src/display.jl:114
 [3] round_string at /Users/james/.julia/packages/IntervalArithmetic/iO2gf/src/display.jl:125 [inlined]
 [4] basic_representation(::Interval{Float64}, ::Nothing) at /Users/james/.julia/packages/IntervalArithmetic/iO2gf/src/display.jl:143
 [5] representation at /Users/james/.julia/packages/IntervalArithmetic/iO2gf/src/display.jl:198 [inlined] (repeats 2 times)
 [6] show(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Interval{Float64}) at /Users/james/.julia/packages/IntervalArithmetic/iO2gf/src/display.jl:258
 [7] print(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Interval{Float64}) at ./strings/io.jl:31
 [8] print_to_string(::String, ::Vararg{Any,N} where N) at ./strings/io.jl:123
 [9] string at ./strings/io.jl:156 [inlined]
 [10] show(::IOContext{Base.GenericIOBuffer{Array{UInt8,1}}}, ::Root{Interval{Float64}}) at /Users/james/.julia/packages/IntervalRootFinding/H2vL9/src/root_object.jl:17
 [11] #sprint#340(::IOContext{REPL.Terminals.TTYTerminal}, ::Int64, ::Function, ::Function, ::Root{Interval{Float64}}) at ./strings/io.jl:99
 [12] #sprint at ./array.jl:0 [inlined]
 [13] alignment at ./show.jl:1761 [inlined]
 [14] alignment(::IOContext{REPL.Terminals.TTYTerminal}, ::Array{Root{Interval{Float64}},1}, ::UnitRange{Int64}, ::UnitRange{Int64}, ::Int64, ::Int64, ::Int64) at ./arrayshow.jl:68
 [15] print_matrix(::IOContext{REPL.Terminals.TTYTerminal}, ::Array{Root{Interval{Float64}},1}, ::String, ::String, ::String, ::String, ::String, ::String, ::Int64, ::Int64) at ./arrayshow.jl:186
 [16] print_matrix at ./arrayshow.jl:159 [inlined]
 [17] print_array at ./arrayshow.jl:308 [inlined]
 [18] show(::IOContext{REPL.Terminals.TTYTerminal}, ::MIME{Symbol("text/plain")}, ::Array{Root{Interval{Float64}},1}) at ./arrayshow.jl:345
 [19] display(::REPL.REPLDisplay, ::MIME{Symbol("text/plain")}, ::Any) at /Users/james/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:131
 [20] display(::REPL.REPLDisplay, ::Any) at /Users/james/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:135
 [21] display(::Any) at ./multimedia.jl:287
 [22] #invokelatest#1 at ./essentials.jl:742 [inlined]
 [23] invokelatest at ./essentials.jl:741 [inlined]
 [24] print_response(::IO, ::Any, ::Any, ::Bool, ::Bool, ::Any) at /Users/james/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:155
 [25] print_response(::REPL.AbstractREPL, ::Any, ::Any, ::Bool, ::Bool) at /Users/james/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:140
 [26] (::getfield(REPL, Symbol("#do_respond#38")){Bool,getfield(REPL, Symbol("##48#57")){REPL.LineEditREPL},REPL.LineEditREPL,REPL.LineEdit.Prompt})(::Any, ::Any, ::Any) at /Users/james/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:714
 [27] #invokelatest#1 at ./essentials.jl:742 [inlined]
 [28] invokelatest at ./essentials.jl:741 [inlined]
 [29] run_interface(::REPL.Terminals.TextTerminal, ::REPL.LineEdit.ModalInterface, ::REPL.LineEdit.MIState) at /Users/james/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/REPL/src/LineEdit.jl:2273
 [30] run_frontend(::REPL.LineEditREPL, ::REPL.REPLBackendRef) at /Users/james/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:1035
 [31] run_repl(::REPL.AbstractREPL, ::Any) at /Users/james/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/REPL/src/REPL.jl:192
 [32] (::getfield(Base, Symbol("##734#736")){Bool,Bool,Bool,Bool})(::Module) at ./client.jl:362
 [33] #invokelatest#1 at ./essentials.jl:742 [inlined]
 [34] invokelatest at ./essentials.jl:741 [inlined]
 [35] run_main_repl(::Bool, ::Bool, ::Bool, ::Bool, ::Bool) at ./client.jl:346
 [36] exec_options(::Base.JLOptions) at ./client.jl:284
 [37] _start() at ./client.jl:436
Version 1.1.2-pre.0 (2019-05-17)
release-1.1/2aee2debbb (fork: 117 commits, 214 days)
  [d1acc4aa] IntervalArithmetic v0.15.0
  [d2bf35a9] IntervalRootFinding v0.4.0
  [d8418881] Intervals v0.3.3

Which I suppose is a different issue

jwscook commented 5 years ago

I keep hitting the wrong button!

dpsanders commented 5 years ago

I think you need to update to the latest versions of the packages to solve that problem.

jwscook commented 5 years ago

I just did an update but to no avail

dpsanders commented 5 years ago

What is the result of ]st?

jwscook commented 5 years ago 
(v1.1) pkg> st
    Status `~/.julia/environments/v1.1/Project.toml`
  [7d9fca2a] Arpack v0.3.0
  [6e4b80f9] BenchmarkTools v0.4.1
  [861a8166] Combinatorics v0.7.0
  [2a437133] ConcreteAbstractions v0.0.0 #master (https://github.com/tbreloff/ConcreteAbstractions.jl.git)
  [d38c429a] Contour v0.5.1
  [8a292aeb] Cuba v1.0.0
  [667455a9] Cubature v1.2.2+ #master (https://github.com/stevengj/Cubature.jl.git)
  [e30172f5] Documenter v0.20.0
  [7a1cc6ca] FFTW v0.2.4
  [442a2c76] FastGaussQuadrature v0.3.2
  [f6369f11] ForwardDiff v0.10.1
  [4d00f742] GeometryTypes v0.7.1
  [19dc6840] HCubature v1.1.2
  [d1acc4aa] IntervalArithmetic v0.15.0
  [d2bf35a9] IntervalRootFinding v0.4.0
  [8197267c] IntervalSets v0.3.1
  [d8418881] Intervals v0.3.3
  [76087f3c] NLopt v0.5.1
  [2774e3e8] NLsolve v3.0.1
  [5fb14364] OhMyREPL v0.4.0
  [429524aa] Optim v0.17.2
  [91a5bcdd] Plots v0.21.0
  [c46f51b8] ProfileView v0.4.0
  [438e738f] PyCall v1.18.5
  [d330b81b] PyPlot v2.6.3
  [1fd47b50] QuadGK v2.0.2
  [f2b01f46] Roots v0.7.4
  [22ff4e98] SeriesAccelerators v0.2.1 #master (https://github.com/jwscook/SeriesAccelerators.jl.git)
  [699a6c99] SimpleTraits v0.8.0
  [e094c991] SingularIntegralEquations v0.4.1
  [b85f4697] SoftGlobalScope v1.0.7
  [276daf66] SpecialFunctions v0.7.2
  [7957ca88] StatProfilerHTML v0.5.0 #master (https://github.com/tkluck/StatProfilerHTML.jl)
  [90137ffa] StaticArrays v0.10.0
  [fd094767] Suppressor v0.1.1
  [24249f21] SymPy v0.8.3
  [37e2e46d] LinearAlgebra 
  [9abbd945] Profile 
  [9a3f8284] Random 
  [8dfed614] Test
dpsanders commented 5 years ago

That's not the latest published version of IntervalArithmetic.

dpsanders commented 5 years ago

What happens when you do ]up?

dpsanders commented 5 years ago

Or try removing and re-adding the packages maybe?

jwscook commented 5 years ago 
(v1.1) pkg> up
  Updating registry at `~/.julia/registries/General`
┌ Warning: Some registries failed to update:
│     — `~/.julia/registries/General` — registry dirty
└ @ Pkg.Types ~/Documents/code/builds/julia1.1/usr/share/julia/stdlib/v1.1/Pkg/src/Types.jl:1269
  Updating git-repo `https://github.com/stevengj/Cubature.jl.git`
  Updating git-repo `https://github.com/tkluck/StatProfilerHTML.jl`
  Updating git-repo `https://github.com/jwscook/SeriesAccelerators.jl.git`
  Updating git-repo `https://github.com/tbreloff/ConcreteAbstractions.jl.git`
 Resolving package versions...
  Updating `~/.julia/environments/v1.1/Project.toml`
 [no changes]
  Updating `~/.julia/environments/v1.1/Manifest.toml`
 [no changes]
dpsanders commented 5 years ago

It looks like there's some problem with your installation. You should probably try asking on slack. (Or just delete .Julia and reinstall...)

jwscook commented 5 years ago

I'm not in a position to clobber my .Julia. I cloned the latest IntervalArithmetic.jl (v0.16.0 #master) and it works!

lbenet commented 5 years ago

Happy that tagging solved this issue!

jwscook commented 5 years ago

thanks to everyone for looking into it!