Bisection method in Julia

The bisection method is a simple root-finding method. Methods for finding roots are iterative and try to find an approximate root \(x\) that fulfills \(|f(x)| \leq \epsilon\), where \(\epsilon\) is a small number referred later as tolerance.

Being \(f\) a nonlinear equation and \([a,b]\) the interval where the root has to be found, in each iteration, the bisection method algorithm evaluates \(f\) at the midpoint \(c\) of the interval \([a,b]\) and compares it against the evaluation at \(a\) and \(b\). If \(f(c)=0\) the problem is solved and the root is \(c\), otherwise the next iteration will continue with the interval \([a,c]\) if \(f(x)\) changes the sign in the left half interval or with \([c,b]\) if it changes in the right half interval. This procedure will continue till find the root in the iteration or, most commonly, by reaching a minimum interval, defined by the tolerance, where the solution is found, get an approximation of the root as the last midpoint \(c\).

The Julia implementation may be written as follows:

function bisection(f::Function, a::Number, b::Number;
                   tol::AbstractFloat=1e-5, maxiter::Integer=100)
    fa = f(a)
    fa*f(b) <= 0 || error("No real root in [a,b]")
    i = 0
    local c
    while b-a > tol
        i += 1
        i != maxiter || error("Max iteration exceeded")
        c = (a+b)/2
        fc = f(c)
        if fc == 0
            break
        elseif fa*fc > 0
            a = c  # Root is in the right half of [a,b].
            fa = fc
        else
            b = c  # Root is in the left half of [a,b].
        end
    end
    return c
end

For the equation \(f(x) = x^2-2,\ (x \in [0,2])\), the solution is \(x = +\sqrt 2 = 1.414213 \dots\)

julia> bisection(x -> x^2-2,0,2)
1.4142074584960938

In order to avoid multiple evaluations of f(a) and f(c), its value is cached in the variables fa and fc and used both in the conditionals and the new designation fa = fc. Also, the comparison of the signs of the two halves of the interval is done by fa*fc > 0 instead of sign(fa) == sign(fc), since, in terms of performance, the former uses simliar, but less memory allocation for float values:

julia> x, y = rand(), -rand()
(0.48656317905335866,-0.0722375704738607)

julia> @time sign(x) == sign(y)
  0.000002 seconds (6 allocations: 192 bytes)
false

julia> @time x*y > 0
  0.000002 seconds (5 allocations: 176 bytes)
false