JuliaIntervals / IntervalRootFinding.jl

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

Krawczyk fails with infinite interval #104

Closed dpsanders closed 5 years ago

dpsanders commented 5 years ago
julia> rts = roots(x->x^2 - 2, -∞..∞, Krawczyk)
0-element Array{Root{Interval{Float64}},1}

This is presumably related to the same issue about what happens when the derivative is zero at the midpoint. We should just explicitly check that and move the midpoint if necessary.

dpsanders commented 5 years ago

Why does it work for Newton?

ys1999 commented 5 years ago

Sir, Krawczyk is failed in this case because in the code contractors.jl contract function for Krawczyk and Newton for both is

function newtonlike_contract(op, C, X, tol)
    # use Bisection contractor for this:
    if !(contains_zero(C.f(X)))
        return :empty, X
    end

    # given that have the Jacobian, can also do mean value form

    NX = op(C.f, C.f′, X) ∩ X

    isempty(NX) && return :empty, X
    isinf(X) && return :unknown, NX  # force bisection

    if NX ⪽ X  # isinterior; know there's a unique root inside
        NX =  refine(X -> op(C.f, C.f′, X), NX, tol)
        return :unique, NX
    end

    return :unknown, NX
end

As we can see from code that if NX is empty than it will return empty status. 'op' function for Krawczyk ( single variable) is

function 𝒦(f, f′, X::Interval{T}) where {T}
    m = Interval(mid(X))
    Y = 1 / f′(m)

    m - Y*f(m) + (1 - Y*f′(X)) * (X - m)
end

so for the function f(x)=x^2-2 and interval X=-Inf..Inf will give m=[0,0] and if f′(m) become zero and due to this Y will become ∅ ( this below case ) and it will result in 'empty' status.

julia> X=0..0
[0, 0]

julia> p=1/X
∅

julia> isempty(p)
true
ys1999 commented 5 years ago

This problem can be solved by using this function for Krawczyk ( single variable) -

function 𝒦(f, f′, X::Interval{T}) where {T}
    m = Interval(mid(X , where_bisect))
    Y = 1 / f′(m)

    m - Y*f(m) + (1 - Y*f′(X)) * (X - m)
end

where_bisect is not exactly equal to 0.5 , so it will not give m=[0,0]. where_bisect =0.49609375(approx 0.5)

dpsanders commented 5 years ago

@ys1999 Yes, you are right. I believe this is fixed in #93.

Kolaru commented 5 years ago

Fixed by #93 using the proposition of @ys1999.