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?

Have you at least read the theory behind the bisection method, or, in other words, do you know what's happening with your code? — nbro, Mar 01 '17 at 22:51
Um, the accepted answer on the page you link to *does* use a for-loop, no? — Stefan Pochmann, Mar 01 '17 at 22:51
Who voted to close this question simply doesn't understand the topic. It's very clear what he/she is asking. — nbro, Mar 01 '17 at 23:00
By the way, in general, you should not check for this condition `f(c) == 0`, instead you should check that `f(c)` is close enough to zero, i.e. something like this `f(c) < 1e-09`. This is because when using floating-point numbers, you get "roundoff" errors, i.e. you lose precision when doing computations with floats, because of the finite nature of computers... — nbro, Mar 01 '17 at 23:09
@nbro you are totally right and normally I do that using math.isclose(), but in this case I was forced to avoid more sophisticated functions. I tried doing what you suggest (first condition f(c) < 0.000000001), but I get c = -0.25 instead, which got me even more confused! — AVLZ, Mar 01 '17 at 23:42
@AVLZ Yes, Stefan is right, I forgot to specify the absolute value, so it should be `abs(f(c)) < 1e-09` (for my explanation above). — nbro, Mar 02 '17 at 00:06

Prune · Answer 1 · 2017-03-01T23:13:29.827

With each iteration, you cut the search interval in half. At each iteration, you check to see whether f(c), the functional value at the midpoint, is 0, within the accuracy of the computer's floating-point representation.

If you were to choose the interval [-101, 99] you would get the solution in only 1 iteration, as that's when you hit c = -1. The program will stop whenever it gets close enough to the actual root that the evaluation comes up 0.000000

You started with a range of width 5. What is 5 divided by 2, 52 times? What is the accuracy of a floating-point number on your implementation? I'm going to bet that you are close to the minimum difference between two floats.

If you really want to see this in action, add one simple line inside the loop, right at the top:

print a, b, c, f(c)

This will show you the progress of finding the root.

print statements are a low-tech, effective way to trace a program.

COMMENT RESPONSE

Good point: I haven't called out the specific case hard enough.

You finished in 52 iterations because that's how many it took for the program to "stumble" across the proper value. When you changed and made it in 53 iterations with a smaller range ... the simple way is to say that you got a little lucky the first time. As I did point out, if you start with something that has -1 as its midpoint, such as [-101, 99], then you'll finish in only one iteration, despite having a much larger interval.

This doesn't really explain its special case, but it's only a general explanation of the method, which isn't what the OP is asking. — nbro, Mar 01 '17 at 23:10
Thanks, I did notice that with an interval that had zero as its midpoint I would need only 1 iteration (the first condition would be true the first time). While I understand the method on a superficial level and why it ultimately works, I still don't get why 52 is the answer in this case or 53 in the case of the other interval I mentioned. I'm guessing it has to do with the limited accuracy of floating point on Python, but how exactly? — AVLZ, Mar 01 '17 at 23:54
It's a matter of when that pure bisection gets sufficiently close to the actual root, which you see in the answer you chose. — Prune, Mar 02 '17 at 00:39

Stefan Pochmann · Accepted Answer · 2017-03-03T06:50:40.093

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.

This was exactly what I wanted to understand!! After seeing your answer I re-read the Python Tutorial (15.) and it also helps to understand this question, but I wouldn't have got there without your explanation. Many thanks! Btw your last example was very useful to understand how the denominator of powers of 2 links to the steps. I checked with [-4000, 1] and guessed the number of additional steps to be around 10 since 2^10 = 1024 ~ 4 * 1000. — AVLZ, Mar 02 '17 at 00:36
@AVLZ Yeah, I'm quite happy with that last example. Shows quite well how the different parts matter, I think. Here's a challenge: Guess how many steps it'll take for function `x**2 - 1000000` and initial range `[999.3, 1000.3]`. — Stefan Pochmann, Mar 02 '17 at 01:44
I checked that it's 41, but I don't get it. Especially because [999.6, 1000.6] would need 42 steps -> it's the same range, and the midpoint is closer to the root, yet it needs more tries. Does it have to do with the fact that the midpoint is above and not below the root? — AVLZ, Mar 02 '17 at 17:01
@AVLZ I don't think there's a meaningful reason of that 41 vs 42, I think that's just luck, like the 52 vs 53 in your original examples. But do you understand why you should expect around 42 here? — Stefan Pochmann, Mar 02 '17 at 18:06
Not really, though I've tried to apply the same rationale as in the other function's examples. Something like the range being 1/4 of length so I should expect 2^2 = 4 steps less? But it doesn't do the trick. — AVLZ, Mar 02 '17 at 22:26
@AVLZ No, the final value is 1000 and the initial range has size 1. That's 1/1000, so about 1/(2^10). So to get to 1/(2^52) you only need about 42 bisections. — Stefan Pochmann, Mar 03 '17 at 06:47

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

2 Answers2