The bisection method

Neither Newton's method nor the secant method can guarantee that an existing solution will be found (see Exercise 72: Understand why Newton's method can fail and Exercise 73: 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.m):

function bisection_method()
    f = @(x) x^2 - 9;
    eps = 1e-6;
    a = 0;   b = 1000;
    [solution, no_iterations] = bisection(f, a, b, eps);
    if solution <= b   % Solution found
        fprintf('Number of function calls: %d\n', 1+2*no_iterations);
        fprintf('A solution is: %f\n', solution);
    else
        fprintf('Abort execution.\n');
    end
end

function [result1, result2] = bisection(f, x_L, x_R, eps)
    if f(x_L)*f(x_R) > 0
        fprintf('Error! Function does not have opposite\n');
        fprintf('signs at interval endpoints!')
        exit(1)
    end
    x_M = (x_L + x_R)/2.0;
    f_M = f(x_M);
    iteration_counter = 1;
    while abs(f_M) > eps
        left_f = f(x_L);
        right_f = f(x_R);
        if left_f*f_M > 0   % i.e., same sign
            x_L = x_M;
        else
            x_R = x_M;
        end
        x_M = (x_L + x_R)/2;
        f_M = f(x_M);
        iteration_counter = iteration_counter + 2;
    end
    result1 = x_M;
    result2 = iteration_counter;
end

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.m 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 stored as a separate file bisection.m for easy use by other programs.