newton1d - debug mode, robustness, harder test cases (#50)

* Add debug mode and bug fixes

* Add tests

* Major robustness improvements and hard test cases added

* Documentation and minor fixes

Author: eeshan9815

src/IntervalRootFinding.jl

@@ -23,7 +23,7 @@ export
 export isunique, root_status
 
 
-import IntervalArithmetic.interval
+import IntervalArithmetic: interval, wideinterval
 
 
 

src/newton1d.jl

@@ -5,70 +5,151 @@ Optional keyword arguments give the tolerances `reltol` and `abstol`.
 and a `debug` boolean argument that prints out diagnostic information."""
 
 function newton1d{T}(f::Function, f′::Function, x::Interval{T};
-                    reltol=eps(T), abstol=eps(T), debug=false)
+                    reltol=eps(T), abstol=eps(T), debug=false, debugroot=false)
 
-    L = Interval{T}[]
+    L = Interval{T}[] # Array to hold the intervals still to be processed
 
-    R = Root{Interval{T}}[]
+    R = Root{Interval{T}}[] # Array to hold the `root` objects obtained
+    reps = reps1 = 0
 
-    push!(L, x)
+    push!(L, x) # Initialize
+    initial_width = ∞
+    X = emptyinterval(T) # Initialize
+    while !isempty(L)   # Until all intervals have been processed
+        X = pop!(L)     # Process next interval
+
+        debug && (print("Current interval popped: "); @show X)
 
-    while !isempty(L)
-        X = pop!(L)
         m = mid(X)
         if (isempty(X))
             continue
         end
 
-        if 0 ∉ f′(X)
+        if 0 ∉ f′(X)    # if 0 ∉ f′(X), Newton steps can be made normally until either X becomes empty, or is known to contain a unique root
+
+            debug && println("0 ∉ f′(X)")
+
             while true
+
                 m = mid(X)
-                N = m - (f(Interval(m)) / f′(X))
+                N = m - (f(interval(m)) / f′(X))    # Newton step
+
+                debug && (print("Newton step1: "); @show (X, X ∩ N))
+                if X == X ∩ N   # Checking if Newton step was redundant
+                    reps1 += 1
+                    if reps1 > 20
+                        reps1 = 0
+                        break
+                    end
+                end
                 X = X ∩ N
 
-                if isempty(X)
+                if (isempty(X))         # No root in X
                     break
 
-                elseif 0 ∈ f(Interval(prevfloat(m), nextfloat(m)))
-                    push!(R, Root(X, :unique))
+                elseif 0 ∈ f(wideinterval(mid(X)))  # Root guaranteed to be in X
+                    n = fa = fb = 0
+                    root_exist = true
+                    while (n < 4 && (fa == 0 || fb == 0))   # Narrowing the interval further
+                        if fa == 0
+                            if 0 ∈ f(wideinterval(X.lo))
+                                fa = 1
+                            else
+                                N = X.lo - (f(interval(X.lo)) / f′(X))
+                                X = X ∩ N
+                                if (isempty(X))
+                                    root_exist = false
+                                    break
+                                end
+                            end
+                        end
+                        if fb == 0
+                            if 0 ∈ f(wideinterval(X.hi))
+                                fb = 1
+                            else
+                                if 0 ∈ f(wideinterval(mid(X)))
+                                    N = X.hi - (f(interval(X.hi)) / f′(X))
+                                else
+                                    N = mid(X) - (f(interval(mid(X))) / f′(X))
+                                end
+                                X = X ∩ N
+                                if (isempty(X))
+                                    root_exist = false
+                                    break
+                                end
+                            end
+                        end
+                        N = mid(X) - (f(interval(mid(X))) / f′(X))
+                        X = X ∩ N
+                        if (isempty(X))
+                            root_exist = false
+                            break
+                        end
+                        n += 1
+                    end
+                    if root_exist
+                        push!(R, Root(X, :unique))
+                        debugroot && @show "Root found", X  # Storing determined unique root
+                    end
+
                     break
                 end
             end
 
-        else
-            # 0 ∈ f'(X)
+        else                                # 0 ∈ f′(X)
+            if diam(X) == initial_width     # if no improvement occuring for a number of iterations
+                reps += 1
+                if reps > 10
+                    push!(R, Root(X, :unknown))
+                    debugroot && @show "Repeated root found", X
+                    reps = 0
+                    continue
+                end
+            end
+            initial_width = diam(X)
+            debug && println("0 ∈ f'(X)")
+
             expansion_pt = Inf
             # expansion point for the newton step might be m, X.lo or X.hi according to some conditions
 
-            if 0 ∈ f(Interval(mid(X)))
+            if 0 ∈ f(wideinterval(mid(X)))  # if root in X, narrow interval further
                 # 0 ∈ fⁱ(x)
-                # Step 7
 
-                if 0 ∉ f(Interval(X.lo))
+                debug && println("0 ∈ fⁱ(x)")
+
+                if 0 ∉ f(wideinterval(X.lo))
                     expansion_pt = X.lo
 
-                elseif 0 ∉ f(Interval(X.hi))
+                elseif 0 ∉ f(wideinterval(X.hi))
                     expansion_pt = X.hi
 
                 else
-                    x1 = mid(Interval(X.lo, mid(X)))
-                    x2 = mid(Interval(mid(X), X.hi))
-                    if 0 ∉ f(Interval(x1)) || 0 ∉ f(Interval(x2))
-                        push!(L, Interval(X.lo, m))
-                        push!(L, Interval(m, X.hi))
+                    x1 = mid(interval(X.lo, mid(X)))
+                    x2 = mid(interval(mid(X), X.hi))
+                    if 0 ∉ f(wideinterval(x1)) || 0 ∉ f(wideinterval(x2))
+                        push!(L, interval(X.lo, m))
+                        push!(L, interval(m, X.hi))
                         continue
 
                     else
-                        push!(R, Root(X, :unique))
+                        push!(R, Root(X, :unknown))
+
+                        debugroot && @show "Multiple root found", X
+
                         continue
                     end
                 end
 
             else
                 # 0 ∉ fⁱ(x)
 
-                if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol
+                debug && println("0 ∉ fⁱ(x)")
+
+                if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol     # checking if X is still within tolerances
                     push!(R, Root(X, :unknown))
+
+                    debugroot && @show "Tolerance root found", X
+
                     continue
                 end
             end
@@ -78,26 +159,28 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T};
                 expansion_pt = mid(X)
             end
 
-            initial_width = diam(X)
 
-            a = f(Interval(expansion_pt))
+            a = f(interval(expansion_pt))
             b = f′(X)
-
+            # Newton steps with extended division creating two intervals
             if 0 < b.hi && 0 > b.lo && 0 ∉ a
                 if a.hi < 0
-                    push!(L, X ∩ (expansion_pt - Interval(-Inf, a.hi / b.hi)))
-                    push!(L, X ∩ (expansion_pt - Interval(a.hi / b.lo, Inf)))
+                    push!(L, X ∩ (expansion_pt - interval(-Inf, a.hi / b.hi)))
+                    push!(L, X ∩ (expansion_pt - interval(a.hi / b.lo, Inf)))
 
                 elseif a.lo > 0
-                    push!(L, X ∩ (expansion_pt - Interval(-Inf, a.lo / b.lo)))
-                    push!(L, X ∩ (expansion_pt - Interval(a.lo / b.hi, Inf)))
+                    push!(L, X ∩ (expansion_pt - interval(-Inf, a.lo / b.lo)))
+                    push!(L, X ∩ (expansion_pt - interval(a.lo / b.hi, Inf)))
 
                 end
 
                 continue
 
             else
-                N = expansion_pt - (f(Interval(expansion_pt))/f′(X))
+                N = expansion_pt - (f(interval(expansion_pt))/f′(X))
+
+                debug && (print("Newton step2: "); @show (X, X ∩ N))
+
                 X = X ∩ N
                 m = mid(X)
 
@@ -106,12 +189,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T};
                 end
             end
 
-            if diam(X) > initial_width/2
-                push!(L, Interval(m, X.hi))
-                X = Interval(X.lo, m)
-            end
-
-            push!(L, X)
+            push!(L, X) # Pushing X into L to be processed again
         end
     end
 

test/newton1d.jl

@@ -18,36 +18,69 @@ three_halves_pi = 3*big_pi/2
 
 @testset "Testing newton1d" begin
 
-    f(x) = exp(x^2) - cos(x)
+    f(x) = exp(x^2) - cos(x)    #double root
     f′(x) = 2*x*exp(x^2) + sin(x)
-    f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24
+    f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24  #four unique roots
     f1′(x) = 4x^3 - 30x^2 + 70x - 50
-
+    f2(x) = 4567x^2 - 9134x + 4567  #double root
+    f2′(x) = 9134x - 9134
+    f3(x) = (x^2 - 2)^2 #two double roots
+    f3′(x) = 4x * (x^2 - 2)
+    f4(x) = sin(x) - x  #triple root at 0
+    f4′(x) = cos(x) - 1
+    f5(x) = (x^2 - 1)^4 * (x^2 - 2)^4 #two quadruple roots
+    f5′(x) = 8x * (-3 + 2x^2) * (2 - 3x^2 + x^4)^3
     for autodiff in (false, true)
         if autodiff
             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)
 
         else
             rts1 = newton1d(sin, cos, -5..5)
             rts2 = newton1d(f, f′, -∞..∞)
             rts3 = newton1d(f1, f1′, -10..10)
+            rts4 = newton1d(f2, f2′, -10..11)
+            rts5 = newton1d(f3, f3′, -10..10)
+            rts6 = newton1d(f4, f4′, -10..10)
+            rts7 = newton1d(f5, f5′, -10..10, reltol=0)
         end
 
         @test length(rts1) == 3
         L = [ -pi_interval(Float64), 0..0, pi_interval(Float64)]
         for i = 1:length(rts1)
-            @test L[i] == rts1[i].interval && :unique == rts1[i].status
+            @test L[i] in rts1[i].interval && :unique == rts1[i].status
         end
 
         @test length(rts2) == 1
-        @test (0..0) == rts2[1].interval && :unique == rts2[1].status
+        @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
 end