======Bisection root-finding in R======

=====Iterative=====
The following code was adapted from an original version by Ravi Varadhan posted on the R-help mailing list 20/9/2010.

<code>
bisectRavi <- function(fn, lo, hi, tol = 1e-7, ...) {

    flo <- fn(lo, ...) 
    fhi <- fn(hi, ...) 

    if(flo * fhi > 0) 
        stop("root is not bracketed by lo and hi")

    chg <- hi - lo
    while (abs(chg) > tol) {
        mid <- (lo + hi) / 2
        fmid <- fn(mid, ...)
        if (abs(fmid) <= tol) break
        if (flo * fmid < 0) hi <- mid 
        if (fhi * fmid < 0) lo <- mid
        chg <- hi - lo
    }

    return(mid)
}
</code>

=====Recursive=====
The following code was posted to R-help by myself (BioStatMatt) as a recursive alternative to the iterative method above.

<code>
bisectMatt <- function(fn, lo, hi, tol = 1e-7, ...) {

    flo <- fn(lo, ...)
    fhi <- fn(hi, ...)

    if(flo * fhi > 0)
        stop("root is not bracketed by lo and hi")

    mid <- (lo + hi) / 2
    fmid <- fn(mid, ...)
    if(abs(fmid) <= tol || abs(hi-lo) <= tol)
        return(mid)

    if(fmid * fhi > 0)
        return(bisectMatt(fn, lo, mid, tol, ...))

    return(bisectMatt(fn, mid, hi, tol, ...))
} 
</code>

=====Illustration=====

These functions find the root of a smooth function, provided two 'bracketing' starting values:

<code>
> testFn <- function(x, a = 1) exp(-x) - a*x
> curve(testFn, from = -1, to = 3, xlab = "x", ylab = "testFn(x)")
> # Plot the 'bracketing points'
> points(c(0, 2),testFn(c(0,2)), pch = "x", col= "red" )
> # Compute and plot the root
> abline(h = 0, lty = 2, col = "red")
> abline(v = bisectMatt(testFn, 0, 2), lty = 2, col = "red")
</code>

This code produces the following graphic:

{{:bisect.png|Bisection}}

=====Timing=====
The iterative function is roughly twice faster than the recursive function:
<code>
> testFn <- function(x, a) exp(-x) - a*x 
> system.time(replicate(1000, bisectMatt(testFn, 0, 2, a=1)))
   user  system elapsed 
  0.656   0.000   0.658 
> system.time(replicate(1000, bisectRavi(testFn, 0, 2, a=1)))
   user  system elapsed 
  0.384   0.000   0.393 
</code>