Add solvers for quadratic interval equations (#68)

* Add solvers for quadratic interval equations

* Add tests

* interval arithmetic to bound rounding errors and refine tests

* Documentation

* optimization

Author: eeshan9815

src/IntervalRootFinding.jl

@@ -18,7 +18,7 @@ export
     derivative, jacobian,  # reexport derivative from ForwardDiff
     Root, is_unique,
     roots, find_roots,
-    bisect, newton1d
+    bisect, newton1d, quadratic_roots
 
 export isunique, root_status
 
@@ -58,7 +58,7 @@ include("complex.jl")
 include("contractors.jl")
 include("roots.jl")
 include("newton1d.jl")
-
+include("quadratic.jl")
 
 
 gradient(f) = X -> ForwardDiff.gradient(f, SVector(X))

src/quadratic.jl

@@ -0,0 +1,72 @@
+"""
+Helper function for quadratic_interval that computes roots of a
+real quadratic using interval arithmetic to bound rounding errors.
+"""
+function quadratic_helper!{T}(a::Interval{T}, b::Interval{T}, c::Interval{T}, L::Array{Interval{T}})
+
+    Δ = b^2 - 4*a*c
+
+    Δ.hi < 0 && return
+
+    Δ = sqrt(Δ)
+
+    if (b.lo >= 0)
+        x0 = -0.5 * (b + Δ)
+
+    else
+        x0 = -0.5 * (b - Δ)
+    end
+
+    if (0 ∈ x0)
+        push!(L, x0)
+
+    else
+        x1 = c / x0
+        x0 = x0 / a
+        push!(L, x0, x1)
+    end
+
+end
+
+
+"""
+Function to solve a quadratic equation where the coefficients are intervals.
+Returns an array of intervals of the roots.
+Arguments `a`, `b` and `c` are interval coefficients of `x²`, `x` and `1` respectively.
+The interval case differs from the non-interval case in that
+there might be three disjoint interval roots. In the third
+case, one interval root extends to −∞ and another extends to +∞.
+This algorithm finds the set of points where `F.lo(x) ≥ 0` and the set
+of points where `F.hi(x) ≤ 0` and takes the intersection of these two sets.
+Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 8
+"""
+function quadratic_roots{T}(a::Interval{T}, b::Interval{T}, c::Interval{T})
+
+    L = Interval{T}[]
+    R = Interval{T}[]
+
+    quadratic_helper!(Interval(a.lo), Interval(b.lo), Interval(c.lo), L)
+    quadratic_helper!(Interval(a.hi), Interval(b.hi), Interval(c.hi), L)
+    quadratic_helper!(Interval(a.lo), Interval(b.hi), Interval(c.lo), L)
+    quadratic_helper!(Interval(a.hi), Interval(b.lo), Interval(c.hi), L)
+
+    if (length(L) == 8)
+        resize!(L, 4)
+    end
+
+    if (a.lo < 0 || (a.lo == 0 && b.hi == 0) || (a.lo == 0 && b.hi == 0 && c.lo ≤ 0))
+        push!(L, Interval(-∞))
+    end
+
+    if (a.lo < 0 || (a.lo == 0 && b.lo == 0) || (a.lo == 0 && b.lo == 0 && c.lo ≤ 0))
+        push!(L, Interval(∞))
+    end
+
+    sort!(L, by = x -> x.lo)
+
+    for i in 1:2:length(L)
+        push!(R, Interval(L[i].lo, L[i+1].hi))
+    end
+
+    R
+end

test/quadratic.jl

@@ -0,0 +1,21 @@
+using IntervalArithmetic, IntervalRootFinding
+
+struct Quad{T}
+    a::Interval{T}
+    b::Interval{T}
+    c::Interval{T}
+    x::Vector{Interval{T}}
+end
+
+
+@testset "Quadratic Interval Equations" begin
+    rts = Quad{Float64}[]
+    push!(rts, Quad(interval(1, 5), interval(2, 4), interval(0, 2), [interval(-4, -2), interval(-0, -0)]))
+    push!(rts, Quad(interval(-1, 1), interval(-2, 2), interval(-1, 1), [interval(-∞, -1), interval(-1, -1), interval(-1, ∞)]))
+    push!(rts, Quad(interval(1, 2), interval(1, 3), interval(2, 6), [interval(-2, -1)]))
+    push!(rts, Quad(interval(1, 1), interval(2, 2), interval(1, 1), [interval(-1, -1), interval(-1, -1)]))
+
+    for i in 1:length(rts)
+            @test quadratic_roots(rts[i].a, rts[i].b, rts[i].c) == rts[i].x
+    end
+end

test/runtests.jl

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