Add slope interval newton method

Author: eeshan9815

src/IntervalRootFinding.jl

@@ -18,7 +18,8 @@ export
     derivative, jacobian,  # reexport derivative from ForwardDiff
     Root, is_unique,
     roots, find_roots,
-    bisect, newton1d, slope
+    bisect, newton1d, slope,
+    slope_newton1d
 
 export isunique, root_status
 

src/newton1d.jl

@@ -25,7 +25,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T};
             continue
         end
 
-        if 0 ∉ f′(X)
+        if 0 ∉ ForwardDiff.derivative(f, X)
 
             debug && println("0 ∉ f′(X)")
 
@@ -161,7 +161,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T};
 
 
             a = f(interval(expansion_pt))
-            b = f′(X)
+            b = ForwardDiff.derivative(f, X)
 
             if 0 < b.hi && 0 > b.lo && 0 ∉ a
                 if a.hi < 0
@@ -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 = newton1d(sin, -5..5)
+    rts2 = newton1d(f, -∞..∞)
+    rts3 = newton1d(f1, -10..10)
+    rts4 = newton1d(f2, -10..11)
+    rts5 = newton1d(f3, -10..10)
+    rts6 = newton1d(f4, -10..10)
+    rts7 = 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    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