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

Reenable and fix example 5.5 from Smiley et al. #153

Closed gwater closed 4 years ago

gwater commented 4 years ago

fix for #86

gwater commented 4 years ago

unfortunately, the example takes to long to solve and some CI tests time out. So we might have to create a faster test for CI

gwater commented 4 years ago

implements the idea described in https://github.com/JuliaArrays/StaticArrays.jl/issues/450

gwater commented 4 years ago

Ok, given that there apparently is no way to test whether LU decomposition will work, I dont see a solution. The fundamental problem remains that IntervalArithmetic tolerates silent failures in external algorithms. I spent a few hours generalizing (and cleaning up) IntervalRootFinding to support NumberIntervals and then the solution is easy:

function 𝒩(f::Function, jacobian::Function, X::AbstractVector{T}, α) where T # multidimensional Newton operator
    m = T.(mid.(X, α))
    J = jacobian(X)
    try
        return convert(typeof(X), m .- (J \ f(m)))
    catch
        return X .± Inf
    end
end

I've tested it: The exception is triggered (disturbingly often) and eventually all roots in Smiley 5.5 are discovered.

gwater commented 4 years ago

ok, this version uses the internal algorithm and finally finds all the roots for Smiley 5.5

gwater commented 4 years ago

not sure I trust the internal solver – for the matrix discussed above, it basically just returns this hard-coded value from https://github.com/JuliaIntervals/IntervalRootFinding.jl/blob/80328b460deae07f447f83c8a6dbefad21837b92/src/linear_eq.jl#L90

gwater commented 4 years ago

not sure I trust the internal solver – for the matrix discussed above, it basically just returns this hard-coded value from

https://github.com/JuliaIntervals/IntervalRootFinding.jl/blob/80328b460deae07f447f83c8a6dbefad21837b92/src/linear_eq.jl#L90

It's also breaking the Julia naming convention – this function is marked as mutating but it doesn't https://github.com/JuliaIntervals/IntervalRootFinding.jl/blob/80328b460deae07f447f83c8a6dbefad21837b92/src/linear_eq.jl#L102

gwater commented 4 years ago

cleaned up the commit log a bit. this PR contains a few more changes than strictly necessary, lmk if I should separate them into other PRs

dpsanders commented 4 years ago

Thanks a lot!