Add bisection root-finder in arbitrary dimension

examples/complex_roots.jl

@@ -0,0 +1,49 @@
+# Example of finding roots of a complex function
+# f(z) = z^4 + z - 2  (book of Hansen-Walster, chap. 10)
+
+using IntervalRootFinding, IntervalArithmetic
+
+doc"""
+    complex_bisection(f, X)
+
+Find complex roots of $f: \mathbb{C} \to \mathbb{C}$.
+
+Inputs:
+
+- `f`: function $f: \mathbb{C} \to \mathbb{C}$, i.e. a function that accepts a complex number and returns another complex number.
+
+- `X`: An `IntervalBox` specifying the bounds on the real and imaginary parts
+of `z`.
+
+"""
+function complex_bisection(f, X::IntervalBox)
+
+    """
+    Make a 2D real version of the complex function `f` suitable for `bisection`,
+    i.e. that accepts an `IntervalBox` and returns an `IntervalBox`
+    """
+    function g(X::IntervalBox)
+        z = Complex(X...)
+        zz = f(z)
+
+        return IntervalBox(reim(zz))
+    end
+
+    roots = bisection(g, X)
+
+    return g, [Complex(root.interval...) for root in roots]
+end
+
+"""
+    complex_bisection(f, x, y)
+
+Version in which the bounds are specified as two separate `Interval`s.
+"""
+complex_bisection(f, x::Interval, y::Interval) = complex_bisection(f, x × y)
+
+f(z) = z^4 + z - 2
+
+L = 10
+g, roots = complex_bisection(f, -L..L, -L..L)
+
+println(roots)

src/IntervalRootFinding.jl

@@ -14,15 +14,17 @@ export
     Root, is_unique,
     find_roots,
     find_roots_midpoint,
+    bisection,
     bisect
 
 import Base: ⊆, show
 
 const derivative = ForwardDiff.derivative
 const D = derivative
 
-immutable Root{T<:Real}
-    interval::Interval{T}
+# Root object:
+immutable Root{T<:Union{Interval,IntervalBox}}
+    interval::T
     status::Symbol
 end
 
@@ -34,7 +36,26 @@ is_unique{T}(root::Root{T}) = root.status == :unique
 ⊆(a::Root, b::Root) = a.interval ⊆ b.interval
 
 
+# Common functionality:
+
+doc"""Returns the midpoint of the interval x, slightly shifted in case
+the midpoint is an exact root"""
+function guarded_mid{T}(f::Function, x::Interval{T})
+    m = mid(x)
+
+    if f(m) == zero(T)                      # 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
+
+    @assert m ∈ x
+
+    return m
+end
+
+
 include("bisect.jl")
+include("bisection.jl")
 include("newton.jl")
 include("krawczyk.jl")
 

src/bisection.jl

@@ -0,0 +1,37 @@
+
+doc"""
+    bisection(f, X; tolerance=1e-3)
+
+Find possible roots of the function `f` inside the `Interval` or `IntervalBox` `X`.
+"""
+function bisection{T<:Union{Interval,IntervalBox}}(f, X::T; tolerance=1e-3, debug=false)
+
+    image = f(X)
+
+    debug && @show X, image
+
+    if !(zero(X) ⊆ image)
+        return Root{typeof(X)}[]
+    end
+
+    if diam(X) < tolerance
+        return [Root(X, :unknown)]
+    end
+
+    X1, X2 = bisect(X)
+
+    if debug
+        @show X1, X2
+    end
+
+    return [bisection(f, X1, tolerance=tolerance);
+            bisection(f, X2, tolerance=tolerance)]
+
+end
+
+
+function bisection{T<:Union{Interval,IntervalBox}}(f, V::Vector{T}; tolerance=1e-3, debug=false)
+
+    return vcat([bisection(f, X, tolerance=tolerance, debug=debug) for X in V])
+
+end

src/krawczyk.jl

@@ -66,10 +66,10 @@ function krawczyk{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=
     # Maximum level of bisection
     level >= maxlevel && return [Root(x, :unknown)]
 
-    isempty(x) && return Root{T}[]  # [(x, :none)]
+    isempty(x) && return Root{typeof(x)}[]  # [(x, :none)]
     Kx = K(f, f_prime, x) ∩ x
 
-    isempty(Kx ∩ x) && return Root{T}[]  # [(x, :none)]
+    isempty(Kx ∩ x) && return Root{typeof(x)}[]  # [(x, :none)]
 
     if Kx ⪽ x  # isinterior
         debug && (print("Refining "); @show(x))

src/newton.jl

@@ -64,7 +64,7 @@ function newton{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
 
     debug && (print("Entering newton:"); @show(level); @show(x))
 
-    isempty(x) && return Root{T}[]  #[(x, :none)]
+    isempty(x) && return Root{typeof(x)}[]  #[(x, :none)]
 
     z = zero(x.lo)
     tolerance = abs(tolerance)
@@ -79,7 +79,7 @@ function newton{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
         Nx = N(f, x, deriv)
         debug && @show(Nx, Nx ⊆ x, Nx ∩ x)
 
-        isempty(Nx ∩ x) && return Root{T}[]
+        isempty(Nx ∩ x) && return Root{typeof(x)}[]
 
         if Nx ⊆ x
             debug && (print("Refining "); @show(x))

test/bisect.jl

@@ -2,7 +2,7 @@ using IntervalArithmetic, IntervalRootFinding
 using Base.Test
 
 
-@testset "Bisection tests" begin
+@testset "`bisect` function" begin
     X = 0..1
     @test bisect(X) == (0..0.5, 0.5..1)
     @test bisect(X, 0.25) == (0..0.25, 0.25..1)

test/bisection.jl

@@ -0,0 +1,34 @@
+using ValidatedNumerics, ValidatedNumerics.RootFinding
+using Base.Test
+
+@testset "Bisection tests" begin
+    @testset "Single variable" begin
+        f(x) = sin(x) + cos(x/2) + x/5
+
+        roots = bisection(f, -∞..∞, tolerance=1e-3)
+        roots2 = [root.interval for root in roots]
+
+        @test roots2 ==  [  Interval(-0.8408203124999999, -0.8398437499999999),
+                            Interval(3.6425781249999996, 3.6435546874999996),
+                            Interval(6.062499999999999, 6.063476562499999)
+                         ]
+
+    end
+
+    @testset "Two variables" begin
+        @intervalbox g(x, y) = (x^2 + y^2 - 1, y - x)
+
+        X = (-∞..∞) × (-∞..∞)
+        roots = bisection(g, X, tolerance=1e-3)
+        roots2 = [root.interval for root in roots]
+
+        @test roots2 == [IntervalBox(Interval(-0.7080078124999999, -0.7070312499999999), Interval(-0.7080078124999999, -0.7070312499999999)),
+ IntervalBox(Interval(-0.7080078124999999, -0.7070312499999999), Interval(-0.7070312499999999, -0.7060546874999999)),
+ IntervalBox(Interval(-0.7070312499999999, -0.7060546874999999), Interval(-0.7080078124999999, -0.7070312499999999)),
+ IntervalBox(Interval(0.7060546874999999, 0.7070312499999999), Interval(0.7070312499999999, 0.7080078124999999)),
+ IntervalBox(Interval(0.7070312499999999, 0.7080078124999999), Interval(0.7060546874999999, 0.7070312499999999)),
+ IntervalBox(Interval(0.7070312499999999, 0.7080078124999999), Interval(0.7070312499999999, 0.7080078124999999))]
+
+    end
+
+end

test/findroots.jl

@@ -86,7 +86,7 @@ function_list = [
                                             for i in 1:length(roots)
                                                 root = roots[i]
 
-                                                @test isa(root, Root{prec[1]})
+                                                @test isa(root, Root)
                                                 @test is_unique(root)
                                                 @test true_roots[i] ⊆ root.interval
                                             end