JuliaIntervals / IntervalRootFinding.jl

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

LinearAlgebra.SingularException(4) #140

Closed newptcai closed 4 years ago

newptcai commented 4 years ago

The following code triggers an LinearAlgebra.SingularException(4)

using IntervalArithmetic, IntervalRootFinding, ForwardDiff
g((a, b, c, l))=-3*3^(1/3) + (1 + 2*a*b)^(1/3) + (1 + 2*a*c)^(1/3) + (1 + 2*b*c)^(1/3)+l*(a + b + c + a*b*c - 4)
∇g = ∇(g)
box = IntervalBox(0..4,3)×(-1..0)
rts = roots(∇g, box, Newton, 1e-5)
println(rts)

The error messages are

ERROR: LinearAlgebra.SingularException(4)
Stacktrace:
 [1] macro expansion at /home/myusername/.julia/packages/StaticArrays/DBECI/src/triangular.jl:404 [inlined]
 [2] _A_ldiv_B at /home/myusername/.julia/packages/StaticArrays/DBECI/src/triangular.jl:379 [inlined]
 [3] \ at /home/myusername/.julia/packages/StaticArrays/DBECI/src/triangular.jl:25 [inlined]
 [4] macro expansion at /home/myusername/.julia/packages/StaticArrays/DBECI/src/solve.jl:56 [inlined]
 [5] _solve at /home/myusername/.julia/packages/StaticArrays/DBECI/src/solve.jl:46 [inlined]
 [6] \ at /home/myusername/.julia/packages/StaticArrays/DBECI/src/solve.jl:1 [inlined]
 [7] 𝒩(::IntervalRootFinding.var"#69#70"{typeof(g)}, ::IntervalRootFinding.var"#43#44"{IntervalRootFinding.var"#69#70"{typeof(g)}}, ::IntervalBox{4,Float64}, ::Float64) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/contractors.jl:26
 [8] #35 at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/contractors.jl:115 [inlined]
 [9] determine_region_status(::IntervalRootFinding.var"#35#36"{Newton{IntervalRootFinding.var"#69#70"{typeof(g)},IntervalRootFinding.var"#43#44"{IntervalRootFinding.var"#69#70"{typeof(g)}}},Float64}, ::IntervalRootFinding.var"#69#70"{typeof(g)}, ::Root{IntervalBox{4,Float64}}) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/contractors.jl:153
 [10] (::Newton{IntervalRootFinding.var"#69#70"{typeof(g)},IntervalRootFinding.var"#43#44"{IntervalRootFinding.var"#69#70"{typeof(g)}}})(::Root{IntervalBox{4,Float64}}, ::Float64, ::Float64) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/contractors.jl:116
 [11] Newton at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/contractors.jl:115 [inlined]
 [12] process(::DepthFirstSearch{IntervalBox{4,Float64},Newton{IntervalRootFinding.var"#69#70"{typeof(g)},IntervalRootFinding.var"#43#44"{IntervalRootFinding.var"#69#70"{typeof(g)}}},Float64}, ::Root{IntervalBox{4,Float64}}) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/roots.jl:56
 [13] iterate(::DepthFirstSearch{IntervalBox{4,Float64},Newton{IntervalRootFinding.var"#69#70"{typeof(g)},IntervalRootFinding.var"#43#44"{IntervalRootFinding.var"#69#70"{typeof(g)}}},Float64}, ::IntervalRootFinding.BBTree{Root{IntervalBox{4,Float64}}}) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/branch_and_bound.jl:307
 [14] iterate(::DepthFirstSearch{IntervalBox{4,Float64},Newton{IntervalRootFinding.var"#69#70"{typeof(g)},IntervalRootFinding.var"#43#44"{IntervalRootFinding.var"#69#70"{typeof(g)}}},Float64}) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/branch_and_bound.jl:303
 [15] branch_and_prune(::Root{IntervalBox{4,Float64}}, ::Newton{IntervalRootFinding.var"#69#70"{typeof(g)},IntervalRootFinding.var"#43#44"{IntervalRootFinding.var"#69#70"{typeof(g)}}}, ::Type, ::Float64) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/roots.jl:86
 [16] _roots(::Function, ::Function, ::Root{IntervalBox{4,Float64}}, ::Type{Newton}, ::Type{DepthFirstSearch}, ::Float64) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/roots.jl:184
 [17] _roots at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/roots.jl:178 [inlined]
 [18] roots(::Function, ::IntervalBox{4,Float64}, ::Type{Newton}, ::Float64) at /home/myusername/.julia/packages/IntervalRootFinding/rDvpg/src/roots.jl:136
 [19] top-level scope at REPL[7]:1
newptcai commented 4 years ago

If seems that if I slightly shrink the box to avoid including $0$ solves the problem.

dpsanders commented 4 years ago

Maybe it works better with Krawczyk?

newptcai commented 4 years ago

Maybe it works better with Krawczyk?

Amazing! Problem solved.