Root finding using the Bisection Method

I was motivated to write this program after joining in an Open University discussion on root finding using the bisection method (also known as the binary search method or the interval halving method). Chris has nicely described this method and the the errors in determining the root in his blog.

For completeness I will describe the method here. Suppose we wish to determine the single root of a continuous function f(x) that we know lies somewhere between x=a and x=b. That is, we want to find an approximate (or exact) value for α where f(α)=0 in the interval a<x<b.

As f(x) is continuous, the graph of y=f(x) must change sign as the graph crosses the x-axis at x=α. If we consider the value of f at the midpoint c=(a+b)/2 of the interval a<x<b, then if f(c) is the same sign as f(a), then the graph of f must cross the x-axis in the interval c<x<b whereas if f(c) is the opposite sign to f(a) then the graph of f must cross the x-axis in the interval a<x<c. Clearly if f(c)=0 then c=α and the root has been found.

Thus, even if we haven't found the root, we have at least narrowed down its location to either the interval a<x<c or c<x<b. By repeating the process iteratively we can home in on the location of the root by confining it to smaller and smaller intervals. On the first iteration our estimate of the root α is c and we know that the error in determining this root is less than (b-a)/2 since the root lies somewhere between x=a and x=b. On each iteration the maximum error will halve and so we can stop the iterations when this error is less than some tolerance δ.