close #454; some guard on runaway roots (#435)

* close #454; some guard on runaway roots

* fix doc

Author: jverzani

README.md

@@ -206,7 +206,7 @@ f (generic function with 1 method)
 julia> x0 = 0.1147
 0.1147
 
-julia> find_zero(f, x0, Roots.Order1()) ≈ 5.075844588445206 # stopped as |f(xₙ)| ≤ |xₙ|ϵ
+julia> find_zero(f, x0, Roots.Order5()) ≈ 5.936596662527689 # stopped as |f(xₙ)| ≤ |xₙ|ϵ
 true
 
 julia> find_zero(f, x0, Roots.Order1(), atol=0.0, rtol=0.0) # error as no check on `|f(xn)|`

src/convergence.jl

@@ -199,6 +199,18 @@ function iszero_Δx(
     isapprox(a, b, atol=δₐ, rtol=δᵣ)
 end
 
+# test when fconverged to ensure not runawa
+function is_small_Δx(M::AbstractUnivariateZeroMethod,
+                     state::AbstractUnivariateZeroState,
+                     options)
+    δ = abs(state.xn1 - state.xn0)
+    δₐ, δᵣ = options.xabstol, options.xreltol
+    Δₓ = max(_unitless(δₐ), _unitless(abs(state.xn1)) * δᵣ)
+    Δₓ = sqrt(sqrt(sqrt((abs(_unitless(Δₓ)))))) # faster than x^(1/8)
+    return δ ≤ Δₓ
+end
+
+
 isnan_f(M::AbstractBracketingMethod, state) = isnan(state.fxn1) || isnan(state.fxn0)
 isnan_f(M::AbstractNonBracketingMethod, state) = isnan(state.fxn1)
 
@@ -279,6 +291,7 @@ function decide_convergence(
 ) where {T,S}
     xn0, xn1 = state.xn0, state.xn1
     fxn1 = state.fxn1
+
     val ∈ (:f_converged, :exact_zero, :converged) && return xn1
 
     ## XXX this could be problematic
@@ -299,7 +312,8 @@ function decide_convergence(
         elseif val == :not_converged
             # this is the case where runaway can happen
             ## XXX Need a good heuristic to catch that
-            is_approx_zero_f(M, state, options, :relaxed) && return xn1
+            is_approx_zero_f(M, state, options, :relaxed) &&
+                is_small_Δx(M, state, options) && return xn1
         end
     end
 

src/find_zeros.jl

@@ -107,8 +107,8 @@ end
 
 # adjust what we mean by x1 ~ x2 for purposes of adding a new zero
 function approx_close(z1, z2, xatol, xrtol)
-    tol = max(sqrt(xatol), max(abs(z1), abs(z2)) * sqrt(xrtol))
-    abs(z1 - z2) < tol
+    z₁,z₂,δ,ϵ = _unitless.((z1, z2, xatol, xrtol))
+    return isapprox(z₁, z₂; atol=sqrt(δ), rtol=sqrt(ϵ))
 end
 
 # is proposed not near xs? (false and we add proposed)

test/test_bracketing.jl

@@ -292,20 +292,23 @@ end
         Roots.Bisection(),
     ]
     results = [run_tests((f, b) -> find_zero(f, b, M), name="$M") for M in Ms]
-    maxfailures = maximum([length(result.failures) for result in results])
-    maxresidual = maximum([result.maxresidual for result in results])
+    maxfailures = maximum(length(result.failures) for result in results)
+    maxresidual = maximum(result.maxresidual for result in results)
     cnts = [result.evalcount for result in results]
+
     @test maxfailures == 0
     @test maxresidual <= 5e-13
     @test avg(cnts) <= 4700
 
-    ## False position has larger residuals (and failures until maxsteps is increased)
+    ## False position has larger residuals
+    ## Fn #13 fails on numbers 2 and 4 until maxsteps is increased; 100 works
     Ms = [Roots.FalsePosition(i) for i in 1:12]
     results = [run_tests((f, b) -> find_zero(f, b, M), name="$M") for M in Ms]
-    maxfailures = maximum([length(result.failures) for result in results])
-    maxresidual = maximum([result.maxresidual for result in results])
+    maxfailures = maximum(length(result.failures) for result in results)
+    maxresidual = maximum(result.maxresidual for result in results)
     cnts = [result.evalcount for result in results]
-    @test maxfailures <= 0
+
+    @test maxfailures <= 1
     @test maxresidual <= 1e-5
     @test avg(cnts) <= 3000
 

test/test_composable.jl

@@ -40,6 +40,10 @@ using ForwardDiff
         xrts = find_zeros(y, 0s, 10s)
         @test length(xrts) == 1
         @test xrts[1] ≈ 1.886053370668014s
+
+        # issue #434
+        xzs1=find_zeros(x -> cos(x / 1u"m"), -1.6u"m",2u"m")
+        @test length(xzs1) == 2 && maximum(xzs1) ≈ 1.5707963267948966 * u"m"
     end
 
     # Polynomials