The bisection method

Neither Newton's method nor the secant method can guarantee that an existing solution will be found (see Exercise 6.1: Understand why Newton's method can fail and Exercise 6.2: See if the secant method fails). The bisection method, however, does that. However, if there are several solutions present, it finds only one of them, just as Newton's method and the secant method. The bisection method is slower than the other two methods, so reliability comes with a cost of speed.

To solve $ x^2 - 9 = 0 $, $ x \in \left[0, 1000\right] $, with the bisection method, we reason as follows. The first key idea is that if $ f(x) = x^2 - 9 $ is continuous on the interval and the function values for the interval endpoints ($ x_L = 0 $, $ x_R =1000 $) have opposite signs, $ f(x) $ must cross the $ x $ axis at least once on the interval. That is, we know there is at least one solution.

The second key idea comes from dividing the interval in two equal parts, one to the left and one to the right of the midpoint $ x_M = 500 $. By evaluating the sign of $ f(x_M) $, we will immediately know whether a solution must exist to the left or right of $ x_M $. This is so, since if $ f(x_M) \ge 0 $, we know that $ f(x) $ has to cross the $ x $ axis between $ x_L $ and $ x_M $ at least once (using the same argument as for the original interval). Likewise, if instead $ f(x_M) \le 0 $, we know that $ f(x) $ has to cross the $ x $ axis between $ x_M $ and $ x_R $ at least once.

In any case, we may proceed with half the interval only. The exception is if $ f(x_M) \approx 0 $, in which case a solution is found. Such interval halving can be continued until a solution is found. A "solution" in this case, is when $ |f(x_M)| $ is sufficiently close to zero, more precisely (as before): $ |f(x_M)| < \epsilon $, where $ \epsilon $ is a small number specified by the user.

The sketched strategy seems reasonable, so let us write a reusable function that can solve a general algebraic equation $ f(x)=0 $ (bisection_method.py):

def bisection(f, x_L, x_R, eps, return_x_list=False):
    f_L = f(x_L)
    if f_L*f(x_R) > 0:
        print "Error! Function does not have opposite \
                 signs at interval endpoints!"
        sys.exit(1)
    x_M = float(x_L + x_R)/2.0
    f_M = f(x_M)
    iteration_counter = 1
    if return_x_list:
        x_list = []

    while abs(f_M) > eps:
        if f_L*f_M > 0:   # i.e. same sign
            x_L = x_M
            f_L = f_M
        else:
            x_R = x_M
        x_M = float(x_L + x_R)/2
        f_M = f(x_M)
        iteration_counter += 1
        if return_x_list:
            x_list.append(x_M)
    if return_x_list:
        return x_list, iteration_counter
    else:
        return x_M, iteration_counter

def f(x):
    return x**2 - 9

a = 0;   b = 1000

solution, no_iterations = bisection(f, a, b, eps=1.0e-6)

print "Number of function calls: %d" % (1 + 2*no_iterations)
print "A solution is: %f" % (solution)

Note that we first check if $ f $ changes sign in $ [a,b] $, because that is a requirement for the algorithm to work. The algorithm also relies on a continuous $ f(x) $ function, but this is very challenging for a computer code to check.

We get the following printout to the screen when bisection_method.py is run:

Number of function calls: 61
A solution is: 3.000000

We notice that the number of function calls is much higher than with the previous methods.

Required work in the bisection method. If the starting interval of the bisection method is bounded by $ a $ and $ b $, and the solution at step $ n $ is taken to be the middle value, the error is bounded as $$ \begin{equation} \frac{|b-a|}{2^n}, \tag{6.4} \end{equation} $$ because the initial interval has been halved $ n $ times. Therefore, to meet a tolerance $ \epsilon $, we need $ n $ iterations such that the length of the current interval equals $ \epsilon $: $$ \frac{|b-a|}{2^n}=\epsilon\quad\Rightarrow\quad n = \frac{\ln ((b-a)/\epsilon)}{\ln 2}\thinspace . $$ This is a great advantage of the bisection method: we know beforehand how many iterations $ n $ it takes to meet a certain accuracy $ \epsilon $ in the solution.

As with the two previous methods, the function bisection is placed in the file nonlinear_solvers.py for easy import and use.