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

As with the two previous methods, the function `bisection`

is
placed in the file `nonlinear_solvers.py`

for easy import and use.