[Back to main course page]

Root Bracketing Methods

We implement the bisection and regula falsi methods following the pseudocode in the lecture notes.

Bisection

def bisection(f, a, b, ftol, xtol, maxit):
    print(f"{'n':<3} |  {'x_n':<22} |  {'f(x_n)':<10} |  bound on |x_n - r|")
    print("-" * 64)
    fa = f(a)
    x = (a + b) / 2
    fx = f(x)
    print(f"{0:<3d} | {x:< 23.16e} | {fx:< 11.4e} | {(b - a) / 2:< .4e}")
    for n in range(1, maxit + 1):
        if fa * fx <= 0:
            b = x
        else:
            a = x
            fa = fx
        x = (a + b) / 2
        fx = f(x)
        print(f"{n:<3d} | {x:< 23.16e} | {fx:< 11.4e} | {(b - a) / 2:< .4e}")
        if abs(fx) <= ftol or (b - a) / 2 <= xtol:
            return x
    print("Bisection did not converge in maximum number of iterations")
    return x

We run it on Example 2.4: $f(x) = (3x - 2)(x^2 + 1)$ on $[0, 1]$, whose only root is $r = 2/3$, and count the function evaluations.

def f(x):
    return (3 * x - 2) * (x**2 + 1)

evals = 0

def f_counted(x):
    global evals
    evals += 1
    return f(x)

evals = 0
x = bisection(f_counted, 0.0, 1.0, ftol=1e-8, xtol=1e-8, maxit=100)
print(f"\nFunction evaluations: {evals}")

n   |  x_n                    |  f(x_n)     |  bound on |x_n - r|
----------------------------------------------------------------
0   |  5.0000000000000000e-01 | -6.2500e-01 |  5.0000e-01
1   |  7.5000000000000000e-01 |  3.9062e-01 |  2.5000e-01
2   |  6.2500000000000000e-01 | -1.7383e-01 |  1.2500e-01
3   |  6.8750000000000000e-01 |  9.2041e-02 |  6.2500e-02
4   |  6.5625000000000000e-01 | -4.4708e-02 |  3.1250e-02
5   |  6.7187500000000000e-01 |  2.2678e-02 |  1.5625e-02
6   |  6.6406250000000000e-01 | -1.1258e-02 |  7.8125e-03
7   |  6.6796875000000000e-01 |  5.6491e-03 |  3.9062e-03
8   |  6.6601562500000000e-01 | -2.8195e-03 |  1.9531e-03
9   |  6.6699218750000000e-01 |  1.4110e-03 |  9.7656e-04
10  |  6.6650390625000000e-01 | -7.0519e-04 |  4.8828e-04
11  |  6.6674804687500000e-01 |  3.5267e-04 |  2.4414e-04
12  |  6.6662597656250000e-01 | -1.7632e-04 |  1.2207e-04
13  |  6.6668701171875000e-01 |  8.8164e-05 |  6.1035e-05
14  |  6.6665649414062500e-01 | -4.4081e-05 |  3.0518e-05
15  |  6.6667175292968750e-01 |  2.2041e-05 |  1.5259e-05
16  |  6.6666412353515625e-01 | -1.1020e-05 |  7.6294e-06
17  |  6.6666793823242188e-01 |  5.5101e-06 |  3.8147e-06
18  |  6.6666603088378906e-01 | -2.7551e-06 |  1.9073e-06
19  |  6.6666698455810547e-01 |  1.3775e-06 |  9.5367e-07
20  |  6.6666650772094727e-01 | -6.8876e-07 |  4.7684e-07
21  |  6.6666674613952637e-01 |  3.4438e-07 |  2.3842e-07
22  |  6.6666662693023682e-01 | -1.7219e-07 |  1.1921e-07
23  |  6.6666668653488159e-01 |  8.6096e-08 |  5.9605e-08
24  |  6.6666665673255920e-01 | -4.3048e-08 |  2.9802e-08
25  |  6.6666667163372040e-01 |  2.1524e-08 |  1.4901e-08
26  |  6.6666666418313980e-01 | -1.0762e-08 |  7.4506e-09

Function evaluations: 28

Regula falsi

def regula_falsi(f, a, b, ftol, xtol, maxit):
    print(f"{'n':<3} |  {'x_n':<22} |  {'f(x_n)':<10} |  bound on |x_n - r|")
    print("-" * 64)
    fa = f(a)
    fb = f(b)
    x = (a * fb - b * fa) / (fb - fa)
    fx = f(x)
    print(f"{0:<3d} | {x:< 23.16e} | {fx:< 11.4e} | {max(x - a, b - x):< .4e}")
    for n in range(1, maxit + 1):
        if fa * fx < 0:
            b = x
            fb = fx
        else:
            a = x
            fa = fx
        x = (a * fb - b * fa) / (fb - fa)
        fx = f(x)
        print(f"{n:<3d} | {x:< 23.16e} | {fx:< 11.4e} | {max(x - a, b - x):< .4e}")
        if abs(fx) <= ftol or max(x - a, b - x) <= xtol:
            return x
    print("Regula falsi did not converge in maximum number of iterations")
    return x

evals = 0
x = regula_falsi(f_counted, 0.0, 1.0, ftol=1e-8, xtol=1e-8, maxit=100)
print(f"\nFunction evaluations: {evals}")

n   |  x_n                    |  f(x_n)     |  bound on |x_n - r|
----------------------------------------------------------------
0   |  5.0000000000000000e-01 | -6.2500e-01 |  5.0000e-01
1   |  6.1904761904761907e-01 | -1.9760e-01 |  3.8095e-01
2   |  6.5330188679245282e-01 | -5.7207e-02 |  3.4670e-01
3   |  6.6294285668753095e-01 | -1.6081e-02 |  3.3706e-01
4   |  6.6563138063691352e-01 | -4.4820e-03 |  3.3437e-01
5   |  6.6637901783862730e-01 | -1.2461e-03 |  3.3362e-01
6   |  6.6658675885355667e-01 | -3.4624e-04 |  3.3341e-01
7   |  6.6664446963809787e-01 | -9.6185e-05 |  3.3336e-01
8   |  6.6666050079346373e-01 | -2.6719e-05 |  3.3334e-01
9   |  6.6666495392164626e-01 | -7.4219e-06 |  3.3334e-01
10  |  6.6666619090397072e-01 | -2.0616e-06 |  3.3333e-01
11  |  6.6666653451034752e-01 | -5.7268e-07 |  3.3333e-01
12  |  6.6666662995657688e-01 | -1.5908e-07 |  3.3333e-01
13  |  6.6666665646941947e-01 | -4.4188e-08 |  3.3333e-01
14  |  6.6666666383409789e-01 | -1.2274e-08 |  3.3333e-01
15  |  6.6666666587984202e-01 | -3.4096e-09 |  3.3333e-01

Function evaluations: 18

Comparison: faster convergence for regula falsi

Example 2.8: $f(x) = x(10x - 1)(x - 1) + \frac{1}{2}\cos(5\pi x)$ on $[0, 1]$, whose only root in $(0, 1)$ is $r = 1/10$.

import math

def h(x):
    return x * (10 * x - 1) * (x - 1) + 0.5 * math.cos(5 * math.pi * x)

x = bisection(h, 0.0, 1.0, ftol=1e-8, xtol=1e-8, maxit=100)

n   |  x_n                    |  f(x_n)     |  bound on |x_n - r|
----------------------------------------------------------------
0   |  5.0000000000000000e-01 | -1.0000e+00 |  5.0000e-01
1   |  2.5000000000000000e-01 | -6.3480e-01 |  2.5000e-01
2   |  1.2500000000000000e-01 | -2.1869e-01 |  1.2500e-01
3   |  6.2500000000000000e-02 |  2.9976e-01 |  6.2500e-02
4   |  9.3750000000000000e-02 |  5.4319e-02 |  3.1250e-02
5   |  1.0937500000000000e-01 | -8.2498e-02 |  1.5625e-02
6   |  1.0156250000000000e-01 | -1.3696e-02 |  7.8125e-03
7   |  9.7656250000000000e-02 |  2.0469e-02 |  3.9062e-03
8   |  9.9609375000000000e-02 |  3.4183e-03 |  1.9531e-03
9   |  1.0058593750000000e-01 | -5.1320e-03 |  9.7656e-04
10  |  1.0009765625000000e-01 | -8.5496e-04 |  4.8828e-04
11  |  9.9853515625000000e-02 |  1.2821e-03 |  2.4414e-04
12  |  9.9975585937500000e-02 |  2.1372e-04 |  1.2207e-04
13  |  1.0003662109375000e-01 | -3.2059e-04 |  6.1035e-05
14  |  1.0000610351562500e-01 | -5.3430e-05 |  3.0518e-05
15  |  9.9990844726562500e-02 |  8.0144e-05 |  1.5259e-05
16  |  9.9998474121093750e-02 |  1.3357e-05 |  7.6294e-06
17  |  1.0000228881835938e-01 | -2.0036e-05 |  3.8147e-06
18  |  1.0000038146972656e-01 | -3.3394e-06 |  1.9073e-06
19  |  9.9999427795410156e-02 |  5.0091e-06 |  9.5367e-07
20  |  9.9999904632568359e-02 |  8.3484e-07 |  4.7684e-07
21  |  1.0000014305114746e-01 | -1.2523e-06 |  2.3842e-07
22  |  1.0000002384185791e-01 | -2.0871e-07 |  1.1921e-07
23  |  9.9999964237213135e-02 |  3.1307e-07 |  5.9605e-08
24  |  9.9999994039535522e-02 |  5.2178e-08 |  2.9802e-08
25  |  1.0000000894069672e-01 | -7.8267e-08 |  1.4901e-08
26  |  1.0000000149011612e-01 | -1.3044e-08 |  7.4506e-09

x = regula_falsi(h, 0.0, 1.0, ftol=1e-8, xtol=1e-8, maxit=100)

n   |  x_n                    |  f(x_n)     |  bound on |x_n - r|
----------------------------------------------------------------
0   |  5.0000000000000000e-01 | -1.0000e+00 |  5.0000e-01
1   |  1.6666666666666666e-01 | -5.2561e-01 |  3.3333e-01
2   |  8.1252830676146068e-02 |  1.5912e-01 |  8.5414e-02
3   |  1.0110135377158658e-01 | -9.6505e-03 |  6.5565e-02
4   |  9.9966365457377207e-02 |  2.9443e-04 |  1.8714e-02
5   |  9.9999967681573834e-02 |  2.8291e-07 |  1.1014e-03
6   |  9.9999999968992559e-02 |  2.7144e-10 |  1.1014e-03

Comparison: slower convergence for regula falsi

Example 2.9: $f(x) = (7x - 1)^2 (3x - 2)$ on $[0, 1]$, whose roots are $r = 1/7$ and $r = 2/3$.

def g(x):
    return (7 * x - 1)**2 * (3 * x - 2)

x = bisection(g, 0.0, 1.0, ftol=1e-8, xtol=1e-8, maxit=100)

n   |  x_n                    |  f(x_n)     |  bound on |x_n - r|
----------------------------------------------------------------
0   |  5.0000000000000000e-01 | -3.1250e+00 |  5.0000e-01
1   |  7.5000000000000000e-01 |  4.5156e+00 |  2.5000e-01
2   |  6.2500000000000000e-01 | -1.4238e+00 |  1.2500e-01
3   |  6.8750000000000000e-01 |  9.0845e-01 |  6.2500e-02
4   |  6.5625000000000000e-01 | -4.0359e-01 |  3.1250e-02
5   |  6.7187500000000000e-01 |  2.1427e-01 |  1.5625e-02
6   |  6.6406250000000000e-01 | -1.0399e-01 |  7.8125e-03
7   |  6.6796875000000000e-01 |  5.2779e-02 |  3.9062e-03
8   |  6.6601562500000000e-01 | -2.6193e-02 |  1.9531e-03
9   |  6.6699218750000000e-01 |  1.3146e-02 |  9.7656e-04
10  |  6.6650390625000000e-01 | -6.5606e-03 |  4.8828e-04
11  |  6.6674804687500000e-01 |  3.2834e-03 |  2.4414e-04
12  |  6.6662597656250000e-01 | -1.6409e-03 |  1.2207e-04
13  |  6.6668701171875000e-01 |  8.2065e-04 |  6.1035e-05
14  |  6.6665649414062500e-01 | -4.1028e-04 |  3.0518e-05
15  |  6.6667175292968750e-01 |  2.0515e-04 |  1.5259e-05
16  |  6.6666412353515625e-01 | -1.0257e-04 |  7.6294e-06
17  |  6.6666793823242188e-01 |  5.1287e-05 |  3.8147e-06
18  |  6.6666603088378906e-01 | -2.5643e-05 |  1.9073e-06
19  |  6.6666698455810547e-01 |  1.2822e-05 |  9.5367e-07
20  |  6.6666650772094727e-01 | -6.4108e-06 |  4.7684e-07
21  |  6.6666674613952637e-01 |  3.2054e-06 |  2.3842e-07
22  |  6.6666662693023682e-01 | -1.6027e-06 |  1.1921e-07
23  |  6.6666668653488159e-01 |  8.0135e-07 |  5.9605e-08
24  |  6.6666665673255920e-01 | -4.0068e-07 |  2.9802e-08
25  |  6.6666667163372040e-01 |  2.0034e-07 |  1.4901e-08
26  |  6.6666666418313980e-01 | -1.0017e-07 |  7.4506e-09

x = regula_falsi(g, 0.0, 1.0, ftol=1e-8, xtol=1e-8, maxit=100)

n   |  x_n                    |  f(x_n)     |  bound on |x_n - r|
----------------------------------------------------------------
0   |  5.2631578947368418e-02 | -7.3480e-01 |  9.4737e-01
1   |  7.1581654522074586e-02 | -4.4440e-01 |  9.2842e-01
2   |  8.2902780573149426e-02 | -3.0846e-01 |  9.1710e-01
3   |  9.0693968963268076e-02 | -2.3038e-01 |  9.0931e-01
4   |  9.6476053216862132e-02 | -1.8031e-01 |  9.0352e-01
5   |  1.0097889491904012e-01 | -1.4584e-01 |  8.9902e-01
6   |  1.0460618861254856e-01 | -1.2089e-01 |  8.9539e-01
7   |  1.0760287017041165e-01 | -1.0214e-01 |  8.9240e-01
8   |  1.1012767274103073e-01 | -8.7638e-02 |  8.8987e-01
9   |  1.1228869942499750e-01 | -7.6150e-02 |  8.8771e-01
10  |  1.1416249198502471e-01 | -6.6874e-02 |  8.8584e-01
11  |  1.1580497299078098e-01 | -5.9260e-02 |  8.8420e-01
12  |  1.1725807255759692e-01 | -5.2925e-02 |  8.8274e-01
13  |  1.1855392161890105e-01 | -4.7590e-02 |  8.8145e-01
14  |  1.1971760392893346e-01 | -4.3050e-02 |  8.8028e-01
15  |  1.2076901876626434e-01 | -3.9151e-02 |  8.7923e-01
16  |  1.2172417552370361e-01 | -3.5776e-02 |  8.7828e-01
17  |  1.2259611436380773e-01 | -3.2832e-02 |  8.7740e-01
18  |  1.2339557426948575e-01 | -3.0248e-02 |  8.7660e-01
19  |  1.2413148651094148e-01 | -2.7965e-02 |  8.7587e-01
20  |  1.2481134498630427e-01 | -2.5939e-02 |  8.7519e-01
21  |  1.2544148814438080e-01 | -2.4131e-02 |  8.7456e-01
22  |  1.2602731637544504e-01 | -2.2510e-02 |  8.7397e-01
23  |  1.2657346160769459e-01 | -2.1052e-02 |  8.7343e-01
24  |  1.2708392103124283e-01 | -1.9734e-02 |  8.7292e-01
25  |  1.2756216356852018e-01 | -1.8539e-02 |  8.7244e-01
26  |  1.2801121540726085e-01 | -1.7452e-02 |  8.7199e-01
27  |  1.2843372928279009e-01 | -1.6460e-02 |  8.7157e-01
28  |  1.2883204102736071e-01 | -1.5552e-02 |  8.7117e-01
29  |  1.2920821605502114e-01 | -1.4718e-02 |  8.7079e-01
30  |  1.2956408782625223e-01 | -1.3952e-02 |  8.7044e-01
31  |  1.2990128987276187e-01 | -1.3244e-02 |  8.7010e-01
32  |  1.3022128261467741e-01 | -1.2591e-02 |  8.6978e-01
33  |  1.3052537593859251e-01 | -1.1985e-02 |  8.6947e-01
34  |  1.3081474830330445e-01 | -1.1423e-02 |  8.6919e-01
35  |  1.3109046298469240e-01 | -1.0901e-02 |  8.6891e-01
36  |  1.3135348195050575e-01 | -1.0414e-02 |  8.6865e-01
37  |  1.3160467776141868e-01 | -9.9590e-03 |  8.6840e-01
38  |  1.3184484382033548e-01 | -9.5342e-03 |  8.6816e-01
39  |  1.3207470323296627e-01 | -9.1364e-03 |  8.6793e-01
40  |  1.3229491649565403e-01 | -8.7634e-03 |  8.6771e-01
41  |  1.3250608818869283e-01 | -8.4131e-03 |  8.6749e-01
42  |  1.3270877282292609e-01 | -8.0838e-03 |  8.6729e-01
43  |  1.3290347996271620e-01 | -7.7738e-03 |  8.6710e-01
44  |  1.3309067872824806e-01 | -7.4815e-03 |  8.6691e-01
45  |  1.3327080176364081e-01 | -7.2056e-03 |  8.6673e-01
46  |  1.3344424874377972e-01 | -6.9450e-03 |  8.6656e-01
47  |  1.3361138948157314e-01 | -6.6984e-03 |  8.6639e-01
48  |  1.3377256668804094e-01 | -6.4650e-03 |  8.6623e-01
49  |  1.3392809842989581e-01 | -6.2437e-03 |  8.6607e-01
50  |  1.3407828032280270e-01 | -6.0337e-03 |  8.6592e-01
51  |  1.3422338749306703e-01 | -5.8343e-03 |  8.6578e-01
52  |  1.3436367633592594e-01 | -5.6448e-03 |  8.6564e-01
53  |  1.3449938609474985e-01 | -5.4644e-03 |  8.6550e-01
54  |  1.3463074028218380e-01 | -5.2927e-03 |  8.6537e-01
55  |  1.3475794796147123e-01 | -5.1291e-03 |  8.6524e-01
56  |  1.3488120490382552e-01 | -4.9730e-03 |  8.6512e-01
57  |  1.3500069463568290e-01 | -4.8240e-03 |  8.6500e-01
58  |  1.3511658938792484e-01 | -4.6817e-03 |  8.6488e-01
59  |  1.3522905095766102e-01 | -4.5457e-03 |  8.6477e-01
60  |  1.3533823149187033e-01 | -4.4156e-03 |  8.6466e-01
61  |  1.3544427420107921e-01 | -4.2911e-03 |  8.6456e-01
62  |  1.3554731401029113e-01 | -4.1718e-03 |  8.6445e-01
63  |  1.3564747815353770e-01 | -4.0575e-03 |  8.6435e-01
64  |  1.3574488671769269e-01 | -3.9479e-03 |  8.6426e-01
65  |  1.3583965314055038e-01 | -3.8427e-03 |  8.6416e-01
66  |  1.3593188466761288e-01 | -3.7417e-03 |  8.6407e-01
67  |  1.3602168277154242e-01 | -3.6447e-03 |  8.6398e-01
68  |  1.3610914353780484e-01 | -3.5514e-03 |  8.6389e-01
69  |  1.3619435801965532e-01 | -3.4617e-03 |  8.6381e-01
70  |  1.3627741256528320e-01 | -3.3754e-03 |  8.6372e-01
71  |  1.3635838911964090e-01 | -3.2923e-03 |  8.6364e-01
72  |  1.3643736550322280e-01 | -3.2123e-03 |  8.6356e-01
73  |  1.3651441566982930e-01 | -3.1352e-03 |  8.6349e-01
74  |  1.3658960994514885e-01 | -3.0609e-03 |  8.6341e-01
75  |  1.3666301524780927e-01 | -2.9892e-03 |  8.6334e-01
76  |  1.3673469529438836e-01 | -2.9200e-03 |  8.6327e-01
77  |  1.3680471078973139e-01 | -2.8533e-03 |  8.6320e-01
78  |  1.3687311960379306e-01 | -2.7888e-03 |  8.6313e-01
79  |  1.3693997693610921e-01 | -2.7264e-03 |  8.6306e-01
80  |  1.3700533546889923e-01 | -2.6662e-03 |  8.6299e-01
81  |  1.3706924550970900e-01 | -2.6080e-03 |  8.6293e-01
82  |  1.3713175512442208e-01 | -2.5517e-03 |  8.6287e-01
83  |  1.3719291026139141e-01 | -2.4971e-03 |  8.6281e-01
84  |  1.3725275486737840e-01 | -2.4444e-03 |  8.6275e-01
85  |  1.3731133099592552e-01 | -2.3933e-03 |  8.6269e-01
86  |  1.3736867890873425e-01 | -2.3438e-03 |  8.6263e-01
87  |  1.3742483717057102e-01 | -2.2958e-03 |  8.6258e-01
88  |  1.3747984273818004e-01 | -2.2493e-03 |  8.6252e-01
89  |  1.3753373104364150e-01 | -2.2043e-03 |  8.6247e-01
90  |  1.3758653607257676e-01 | -2.1605e-03 |  8.6241e-01
91  |  1.3763829043756992e-01 | -2.1181e-03 |  8.6236e-01
92  |  1.3768902544714479e-01 | -2.0769e-03 |  8.6231e-01
93  |  1.3773877117060879e-01 | -2.0369e-03 |  8.6226e-01
94  |  1.3778755649905117e-01 | -1.9981e-03 |  8.6221e-01
95  |  1.3783540920275950e-01 | -1.9604e-03 |  8.6216e-01
96  |  1.3788235598529841e-01 | -1.9237e-03 |  8.6212e-01
97  |  1.3792842253447543e-01 | -1.8881e-03 |  8.6207e-01
98  |  1.3797363357040171e-01 | -1.8535e-03 |  8.6203e-01
99  |  1.3801801289083931e-01 | -1.8198e-03 |  8.6198e-01
100 |  1.3806158341401323e-01 | -1.7870e-03 |  8.6194e-01
Regula falsi did not converge in maximum number of iterations