Add slope interval newton method

src/IntervalRootFinding.jl

@@ -22,7 +22,7 @@ export
     gauss_seidel_interval, gauss_seidel_interval!,
     gauss_seidel_contractor, gauss_seidel_contractor!,
     gauss_elimination_interval, gauss_elimination_interval!,
-    slope
+    slope, slope_newton1d
 
 export isunique, root_status
 

src/newton1d.jl

@@ -204,3 +204,8 @@ and a `debug` boolean argument that prints out diagnostic information."""
 
 newton1d{T}(f::Function, x::Interval{T};  args...) =
     newton1d(f, x->D(f,x), x; args...)
+
+function slope_newton1d(f::Function, x::Interval{T};
+                        reltol=eps(T), abstol=eps(T), debug=false, debugroot=false) where {T}
+    newton1d(f, x->slope(f, x), x, reltol=reltol, abstol=abstol, debug=debug, debugroot=debugroot)
+end

src/slopes.jl

@@ -3,7 +3,7 @@ using IntervalArithmetic, ForwardDiff
 import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan
 import IntervalArithmetic: mid, interval
 
-function slope(f::Function, x::Interval, c::Number)
+function slope(f::Function, x::Interval, c::Number = mid(x))
     f(slope_var(x, c)).s
 end
 

test/newton1d.jl

@@ -84,3 +84,43 @@ three_halves_pi = 3*big_pi/2
 
     end
 end
+
+@testset "Testing slope newton1d" begin
+
+    f(x) = exp(x^2) - cos(x)    #double root
+    f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24  #four unique roots
+    f2(x) = 4567x^2 - 9134x + 4567  #double root
+    f3(x) = (x^2 - 2)^2 #two double roots
+    f4(x) = sin(x) - x  #triple root at 0
+    f5(x) = (x^2 - 1)^4 * (x^2 - 2)^4 #two quadruple roots
+
+    rts1 = slope_newton1d(sin, -5..5)
+    rts2 = slope_newton1d(f, -∞..∞)
+    rts3 = slope_newton1d(f1, -10..10)
+    rts4 = slope_newton1d(f3, -10..10)
+    rts5 = slope_newton1d(f4, -10..10)
+
+    @test length(rts1) == 3
+    L = [ -pi_interval(Float64), 0..0, pi_interval(Float64)]
+    for i = 1:length(rts1)
+        @test L[i] in rts1[i].interval && :unique == rts1[i].status
+    end
+
+    @test length(rts2) == 1
+    @test (0..0) == rts2[1].interval && :unknown == rts2[1].status
+
+    @test length(rts3) == 4
+    L = [1, 2, 3, 4]
+    for i = 1:length(rts3)
+        @test L[i] in rts3[i].interval && :unique == rts3[i].status
+    end
+
+    L1 = [-sqrt(2), sqrt(2)]
+    for i = 1:length(rts4)
+        @test L1[i] in rts4[i].interval && :unknown == rts4[i].status
+    end
+
+    @test length(rts5) == 1
+    @test 0 in rts5[1].interval && :unknown == rts5[1].status
+
+end