From b05047720ac54fc231ad186e93171a3c743554f1 Mon Sep 17 00:00:00 2001
From: "David P. Sanders" <dpsanders@gmail.com>
Date: Thu, 16 Mar 2017 09:02:02 -0600
Subject: [PATCH] Fix 1D Newton by checking correctly whether midpoint is
 exactly a root (#254)

* Fix 1D Newton by checking correctly whether midpoint is exactly a root

* Fix Krawczyk with new definition of guarded_mid

* Fix Krawczyk and logistic tests

* Fix Krawczyk
---
 src/root_finding/krawczyk.jl         |  9 +++++----
 src/root_finding/newton.jl           | 11 ++++++-----
 test/root_finding_tests/findroots.jl | 18 ++++++++++++++++++
 3 files changed, 29 insertions(+), 9 deletions(-)

diff --git a/src/root_finding/krawczyk.jl b/src/root_finding/krawczyk.jl
index 2a0b189..99fbf0e 100644
--- a/src/root_finding/krawczyk.jl
+++ b/src/root_finding/krawczyk.jl
@@ -10,13 +10,14 @@ does not contain zero, and the second is the inverse of its derivative"""
 function guarded_derivative_midpoint{T}(f::Function, f_prime::Function, x::Interval{T})
 
     α = convert(T, 0.46875)   # close to 0.5, but exactly representable as a floating point
-    m = Interval( guarded_mid(x) )
+    # m = Interval( guarded_mid(f, x) )
+    m = Interval(mid(x))
 
     C = inv(f_prime(m))
 
     # Check that 0 is not in C; if so, consider another point rather than m
     i = 0
-    while zero(T) ∈ C
+    while zero(T) ∈ C || isempty(C)
         m = Interval( α*x.lo + (one(T)-α)*x.hi )
         C = inv(f_prime(m))
 
@@ -70,12 +71,12 @@ function krawczyk{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=
 
     isempty(Kx ∩ x) && return Root{T}[]  # [(x, :none)]
 
-    if Kx ⪽ x
+    if Kx ⪽ x  # isinterior
         debug && (print("Refining "); @show(x))
         return krawczyk_refine(f, f_prime, Kx, tolerance=tolerance, debug=debug)
     end
 
-    m = guarded_mid(x)
+    m = guarded_mid(f, x)
 
     debug && @show(x,m)
 
diff --git a/src/root_finding/newton.jl b/src/root_finding/newton.jl
index c412681..cdce49b 100644
--- a/src/root_finding/newton.jl
+++ b/src/root_finding/newton.jl
@@ -4,10 +4,11 @@
 
 # What is this guarded_mid for? Shouldn't it be checking if f(m)==0?
 doc"""Returns the midpoint of the interval x, slightly shifted in case
-it is zero"""
-function guarded_mid{T}(x::Interval{T})
+the midpoint is an exact root"""
+function guarded_mid{T}(f, x::Interval{T})
     m = mid(x)
-    if m == zero(T)                     # midpoint exactly 0
+
+    if f(m) == 0                      # midpoint exactly a root
         α = convert(T, 0.46875)      # close to 0.5, but exactly representable as a floating point
         m = α*x.lo + (one(T)-α)*x.hi   # displace to another point in the interval
     end
@@ -19,7 +20,7 @@ end
 
 
 function N{T}(f::Function, x::Interval{T}, deriv::Interval{T})
-    m = guarded_mid(x)
+    m = guarded_mid(f, x)
     m = Interval(m)
 
     m - (f(m) / deriv)
@@ -85,7 +86,7 @@ function newton{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
             return newton_refine(f, f_prime, Nx, tolerance=tolerance, debug=debug)
         end
 
-        m = guarded_mid(x)
+        m = guarded_mid(f, x)
 
         debug && @show(x,m)
 
diff --git a/test/root_finding_tests/findroots.jl b/test/root_finding_tests/findroots.jl
index 72b94bf..4243f31 100644
--- a/test/root_finding_tests/findroots.jl
+++ b/test/root_finding_tests/findroots.jl
@@ -147,6 +147,24 @@ end
     end
 end
 
+@testset "Iterated logistic function fixed points" begin
+    ∘(f, g) = x -> f(g(x))
+    iterate(f, n) = n == 1 ? f : f ∘ iterate(f, n-1)
+
+    f(x) = 4x*(1-x)  # logistic map
+
+    for n in 1:10
+        g = x -> iterate(f, n)(x) - x  # look for fixed points of f^n
+
+        roots = newton(g, 0..1)
+        @test length(roots) == 2^n
+
+        roots = krawczyk(g, 0..1)
+        @test length(roots) == 2^n
+
+    end
+end
+
 
 
 # Example of a function with a double root at 0 from Burden & Faires, 9th ed, p.84