Interval Newton Method (#39)

* newton1d for 0 not in f_prime(X)

* add other cases

* Add review changes

* add extended interval division

* add tests

* minor changes

* fix AD bug

Author: eeshan9815

src/IntervalRootFinding.jl

@@ -15,7 +15,7 @@ export
     derivative, jacobian,  # reexport derivative from ForwardDiff
     Root, is_unique,
     roots, find_roots,
-    bisect
+    bisect, newton1d
 
 import Base: ⊆, show
 
@@ -62,6 +62,8 @@ include("newton.jl")
 include("krawczyk.jl")
 include("complex.jl")
 include("branch_and_prune.jl")
+include("newton1d.jl")
+
 
 
 gradient(f) = X -> ForwardDiff.gradient(f, SVector(X))

src/newton1d.jl

@@ -0,0 +1,128 @@
+"""`newton1d` performs the interval Newton method on the given function `f`
+with its derivative `f′` and initial interval `x`.
+Optional keyword arguments give the tolerances `reltol` and `abstol`.
+`reltol` is the tolerance on the relative error whereas `abstol` is the tolerance on |f(X)|,
+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)
+
+    L = Interval{T}[]
+
+    R = Root{Interval{T}}[]
+
+    push!(L, x)
+
+    while !isempty(L)
+        X = pop!(L)
+        m = mid(X)
+        if (isempty(X))
+            continue
+        end
+
+        if 0 ∉ f′(X)
+            while true
+                m = mid(X)
+                N = m - (f(Interval(m)) / f′(X))
+                X = X ∩ N
+
+                if isempty(X)
+                    break
+
+                elseif 0 ∈ f(Interval(prevfloat(m), nextfloat(m)))
+                    push!(R, Root(X, :unique))
+                    break
+                end
+            end
+
+        else
+            # 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)))
+                # 0 ∈ fⁱ(x)
+                # Step 7
+
+                if 0 ∉ f(Interval(X.lo))
+                    expansion_pt = X.lo
+
+                elseif 0 ∉ f(Interval(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))
+                        continue
+
+                    else
+                        push!(R, Root(X, :unique))
+                        continue
+                    end
+                end
+
+            else
+                # 0 ∉ fⁱ(x)
+
+                if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol
+                    push!(R, Root(X, :unknown))
+                    continue
+                end
+            end
+            # Step 8
+
+            if isinf(expansion_pt)
+                expansion_pt = mid(X)
+            end
+
+            initial_width = diam(X)
+
+            a = f(Interval(expansion_pt))
+            b = f′(X)
+
+            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)))
+
+                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)))
+
+                end
+
+                continue
+
+            else
+                N = expansion_pt - (f(Interval(expansion_pt))/f′(X))
+                X = X ∩ N
+                m = mid(X)
+
+                if isempty(X)
+                    continue
+                end
+            end
+
+            if diam(X) > initial_width/2
+                push!(L, Interval(m, X.hi))
+                X = Interval(X.lo, m)
+            end
+
+            push!(L, X)
+        end
+    end
+
+    return R
+end
+
+
+"""`newton1d` performs the interval Newton method on the given function `f` and initial interval `x`.
+Optional keyword arguments give the tolerances `reltol` and `abstol`.
+`reltol` is the tolerance on the relative error whereas `abstol` is the tolerance on |f(X)|,
+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...)

test/newton1d.jl

@@ -0,0 +1,53 @@
+
+using IntervalArithmetic, IntervalRootFinding
+using ForwardDiff
+using Base.Test
+
+const D = IntervalRootFinding.derivative
+
+setprecision(Interval, Float64)
+float_pi = @interval(pi)
+
+setprecision(Interval, 10000)
+big_pi = @interval(pi)
+# Using precision "only" 256 leads to overestimation of the true roots for `cos`
+# i.e the Newton method gives more accurate results!
+
+half_pi = big_pi / 2
+three_halves_pi = 3*big_pi/2
+
+@testset "Testing newton1d" begin
+
+    f(x) = exp(x^2) - cos(x)
+    f′(x) = 2*x*exp(x^2) + sin(x)
+    f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24
+    f1′(x) = 4x^3 - 30x^2 + 70x - 50
+
+    for autodiff in (false, true)
+        if autodiff
+            rts1 = newton1d(sin, -5..5)
+            rts2 = newton1d(f, -∞..∞)
+            rts3 = newton1d(f1, -10..10)
+
+        else
+            rts1 = newton1d(sin, cos, -5..5)
+            rts2 = newton1d(f, f′, -∞..∞)
+            rts3 = newton1d(f1, f1′, -10..10)
+        end
+
+        @test length(rts1) == 3
+        L = [(-π.. -π), (0..0), (π..π)]
+        for i = 1:length(rts1)
+            @test L[i] == rts1[i].interval && :unique == rts1[i].status
+        end
+
+        @test length(rts2) == 1
+        @test (0..0) == rts2[1].interval && :unique == 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
+    end
+end

test/runtests.jl

@@ -4,3 +4,4 @@ using Base.Test
 include("bisect.jl")
 include("findroots.jl")
 include("roots.jl")
+include("newton1d.jl")