Bisection method in Julia

numerical-analysis root-finding 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
        elseif fa*fc > 0
            a = c  # Root is in the right half of [a,b].
            fa = fc
            b = c  # Root is in the left half of [a,b].
    return c

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)

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()

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

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