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(c), its value is cached in the variables
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