Understanding the number of iterations to find a solution using the Bisection method

Question

I was asked to find the root of an equation using the Bisection method and only for loops with Python 3. This thread shows how to use the method, but not with the explanation for the number in range().

As an example, I have the function

f(x) = x² - 2*x - 3

and I want to find its negative root, starting with the interval [-4, 1].

I managed to write the function with a for loop, but I don't understand which range I should use, or how to come up with it.

Here is my code, that solves the problem:

...
a = -4
b = 1
c = (a + b)/2

for i in range(1000):
    if f(c) == 0:
        break
    if f(a) * f(c) < 0:
        b = c
    elif f(a) * f(c) > 0:
        a = c
    c = (a + b) / 2

return c, f(c), i

c = -1 (found the negative root), f(c) = 0 confirming that the program works, and i = 52 which indicates that after 52 "tries" of bisection, it found the right answer.

I put a very large number in range() to make sure I found the root, but why does it need only 52 iterations?

Plus, if my interval changed to [-2, 1], I would need 53 tries. Why is that so?

Answer 1

If you print([a, b]) in the loop, you can see the range evolve:

[-4, 1]
[-1.5, 1]
[-1.5, -0.25]
[-1.5, -0.875]
[-1.1875, -0.875]
[-1.03125, -0.875]
...
...
...
[-1.0000000000000284, -0.9999999999999929]
[-1.0000000000000107, -0.9999999999999929]
[-1.0000000000000018, -0.9999999999999929]
[-1.0000000000000018, -0.9999999999999973]
[-1.0000000000000018, -0.9999999999999996]
[-1.0000000000000007, -0.9999999999999996]

The calculated mean of -1.0000000000000007 and -0.9999999999999996 is exactly -1. Why? Because you've reached the limit of the what floats can represent. Here are the exact values involved:

>>> '%.60f' % -1.0000000000000007
'-1.000000000000000666133814775093924254179000854492187500000000'

>>> '%.60f' % -0.9999999999999996
'-0.999999999999999555910790149937383830547332763671875000000000'

>>> '%.60f' % (-1.0000000000000007 + -0.9999999999999996)
'-2.000000000000000000000000000000000000000000000000000000000000'

>>> '%.60f' % ((-1.0000000000000007 + -0.9999999999999996) / 2)
'-1.000000000000000000000000000000000000000000000000000000000000'

Floats store 52 bits of fraction, meaning 52 bits after the leading 1-bit. Meaning you lose what's smaller than about 1/2⁵² of your value. After 52 steps, your initial range of size 5 has become about size 5/2⁵². That's about 1/2⁵² of your value -1. So around there, you reach a pretty good chance to stumble upon exactly -1 because of the imprecision.

It could take two or maybe three steps more, because 5/2⁵² is still a bit larger than 1/2⁵². You got lucky there. With your other initial range [-2, 1] you just got less lucky. There your range shrinks until [-1.0000000000000002, -0.9999999999999999] before you reach -1.

If you start with [-4000000, 1] instead, then you'll need 72 steps. It's 20 steps more because the initial range is a million times larger, which is about 2²⁰.

Another case: If you use function x**2 - 1000000 and initial range [999.3, 1000.3], then it takes 41 steps. Why? The final value (i.e., the root) is 1000 and the initial range has size 1. That's 1/1000, so about 1/2¹⁰. So to get to 1/2⁵² you only need about 42 bisections.