Add slope interval newton method

Author: eeshan9815

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{T}(f::Function, x::Interval{T};
+                    reltol=eps(T), abstol=eps(T), debug=false, debugroot=false)
+    newton1d(f, x->slope(f, x), x, reltol=reltol, abstol=abstol, debug=debug, debugroot=debugroot)
+end

test/newton1d.jl

@@ -84,3 +84,55 @@ 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(f2, -10..11)
+    rts5 = slope_newton1d(f3, -10..10)
+    rts6 = slope_newton1d(f4, -10..10)
+    rts7 = slope_newton1d(f5, -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
+
+    @test length(rts4) == 1
+    @test 1 in rts4[1].interval && :unknown == rts4[1].status
+
+    L1 = [-sqrt(2), sqrt(2)]
+    for i = 1:length(rts5)
+        @test L1[i] in rts5[i].interval && :unknown == rts5[i].status
+    end
+
+
+    @test length(rts6) == 1
+    @test 0 in rts6[1].interval && :unknown == rts6[1].status
+
+    @test length(rts7) == 4
+    L = [-sqrt(2), -1, 1, sqrt(2)]
+    for i = 1:length(rts7)
+        @test L[i] in rts7[i].interval && :unknown == rts7[i].status
+    end
+
+end