Open dpsanders opened 5 years ago
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)
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
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.
I think this will require directly specifying how the derivative of
abs
on anInterval
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 Interval
s may include x_0
somewhere in the middle.
I hope this may help.
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.
This should now be fixable after https://github.com/JuliaDiff/DiffRules.jl/pull/33.