Newton fails with abs

[dpsanders]

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)

[Kolaru]

It looks like the fact an interval should not be bisected on its center was lost on recent update

julia> roots(f, -1000..1000.001, Newton)
2-element Array{Root{Interval{Float64}},1}:
 Root([0.165881, 0.165943], :unique)
 Root([-2.16602, -2.16504], :unique)

[Kolaru]

Still not fixed despite the corrections in #93. Seems that the problem is indeed the absolute value:

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)

julia> f3(p) = (1 / ( (1+p)^30 ) * 10_000 - 100)^2
f3 (generic function with 1 method)

julia> roots(f3, -1000..1000, Newton)
2-element Array{Root{Interval{Float64}},1}:
 Root([0.16551, 0.166413], :unknown)
 Root([-2.16653, -2.16564], :unknown)

(squaring instead of taking the absolute value should yield a very similar problem)

The abs function together with Interval confuses ForwardDiff:

julia> fp(x) = derivative(f, x)
fp (generic function with 1 method)

julia> fp((-2.167..(-2.165)))
[2378.6, 2770.9]

julia> fp(-2.167)
-2499.911977412291

julia> fp(-2.165)
2636.4380353856686

[dpsanders]

I think this will require directly specifying how the derivative of abs on an Interval behaves; but I don't know how to do so.

[lbenet]

I think this will require directly specifying how the derivative of abs on an Interval behaves; but I don't know how to do so.

In TaylorSeries.jl, we define |f(x)| as the Taylor expansion of f(x) if f(x_0)>0, otherwise as as the Taylor expansion of -f(x) if f(x_0)<0, where x_0 is the expansion point. If f(x_0)=0 there is no first derivative of |f(x)| at x_0, so we throw an error.

Note that the above answer is valid for small enough δx = x-x_0, such that f(x) does not change sign. The problem here is that the Intervals may include x_0 somewhere in the middle.

I hope this may help.

[dpsanders]

What I mean is that we should define the derivative of abs on an interval as something like

function abs_derivative(x::Interval{T}) where {T}
    0 ∈ x && return Interval{T}(-1, 1)
    x.hi < 0 && return Interval{T}(-1)
    return Interval{T}(1)
end

But I don't know how to plug this into ForwardDiff's structure.

[dpsanders]

This should now be fixable after JuliaDiff/DiffRules.jl#33.