# Bisection method in Julia

numerical-analysis root-finding juliaThe 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
```