Bisection method - Numerical Analysis

#%config InlineBackend.figure_format = 'svg'
from pylab import *

Let $f$ be a continuous function. Given an interval $[a,b]$ such that

\sign f(a) \ne \sign f(b), \qquad f(a) \ne 0 \ne f(b)

(1)

then $[a,b]$ is called a non-trivial bracket for a root of $f$ .

3.3.1Main idea¶

The basic idea of bisection method is: starting from a non-trivial bracket, find a new bracketing interval which is smaller in size. As the name suggests, we divide the interval into two equal and smaller intervals. First divide $[a,b]$ into $[a,c]$ and $[c,b]$ where $c = \half(a+b)$ . Then there are three possibilities.

$f(c) = 0$ , then we have found the root exactly.
$f(c) \ne 0$ and $\sign f(c) \ne \sign f(b)$ . Then $[c,b]$ is a non-trivial bracket.
$f(c) \ne 0$ and $\sign f(a) \ne \sign f(c)$ . Then $[a,c]$ is a non-trivial bracket.

We repeat the above process until the root is found or the length of the bracketing interval is reduced to a small desired size. The method may not give a unique value for the root since it can lie anywhere in the bracketing interval. We can take the mid-point of the bracketing interval as the estimate of the root, in which case the error in the root is at most equal to half the length of the bracketing interval.

The length of the bracketing interval reduces by half in each iteration, which means that the algorithm has a monotonic convergence behaviour. If

L_0 = b-a = \textrm{initial length}

(2)

then after $k$ iterations

\textrm{length of bracket interval} = L_k = \frac{1}{2^k} L_0

(3)

As a convergence criterion, we can put a tolerance on the length: stop if $L_k \le \epsilon$ for some $\epsilon > 0$ small. To achieve this tolerance, we require

k \approx \log_2 \frac{L_0}{\epsilon} \qquad \textrm{iterations}

(4)

If $L_0 = 1$ and $\epsilon = 10^{-6}$ then we need about 20 iterations. If we demand more accuracy, then the number of iterations will increase, though only logarithmically.

3.3.2Algorithm¶

The bisection method is shown in Algorithm (1). It requires some input:

initial bracketing interval $[a,b]$
tolerance on function value δ
tolerance on length of bracketing interval ε
maximum number of iterations $N$

In iterative algorithms, when we are not sure that the iterations will always converge, it is good put an upper limit on the number of iterations, since otherwise, the program may never end. But the bisection method is guranteed to satisfy the tolerance on bracketing interval, so it is safe to remove the limit on maximum iterations, but in the code below, we put an upper limit.

The way we have written the algorithm, we store the entire history of the computations in arrays a[k], b[k], c[k] which is not necessary. In actual practice, we only need to keep some of the latest results in memory and this is illustrated in Algorithm (2) which does not require any arrays.

Convergence criteria. We are checking the tolerance in terms of absolute error in the length of bracketing interval without regard to the size of the root.

If $\epsilon = 10^{-6}$ and if the root is near one, then bisection method will return roughly six accurate digits.
If $\epsilon = 10^{-6}$ and the root is $\approx 10^{-7}$ then we cannot expect any figures of accuracy.
On the other hand, if the root is $\gg 1$ and we use $\epsilon \ll 1$ , then we may never be able to achieve this tolerance or it may require far too many iterations.

It is better to use a relative error tolerance

c = \half(a+b), \qquad \textrm{If } |b - a| < \epsilon |c| \quad \textrm{return } c

(5)

which is done in the implementation below.

Example 1

First we define the function for which the root is required.

def fun(x):
    f = exp(x) - sin(x)
    return f

Let us plot and visualize the function.

x = linspace(-4,-2,100)
plot(x,fun(x),'r-')
title('f(x) = exp(x) - sin(x)')
grid(True), xlabel('x'), ylabel('f(x)');

From the figure, we see that [−4,−2] is a bracketing interval. We now implement the bisection method.

# Initial interval [a,b]
a, b = -4.0, -2.0

N     = 100    # Maximum number of iterations
eps   = 1.0e-4 # Tolerance on the interval
delta = 1.0e-4 # Tolerance on the function

fa, fb = fun(a), fun(b)
sa, sb = sign(fa), sign(fb)

for i in range(N):
    e = b-a
    c = a + 0.5*e
    if abs(e) < eps*abs(c):
        print("Interval size is below tolerance\n")
        break
    fc = fun(c)
    if abs(fc) < delta:
        print("Function value is below tolerance\n")
        break
    sc = sign(fc)
    if sa != sc: # new interval is [a,c]
        b  = c
        fb = fc
        sb = sc
    else:        # new interval is [c,b]
        a  = c
        fa = fc
        sa = sc
    print("%5d %16.8e %16.8e %16.8e"%(i+1,c,abs(b-a),abs(fc)))

print("c = %16.8e, |f(c)| = %16.8e" % (c,abs(fc)))

    1  -3.00000000e+00   1.00000000e+00   1.90907076e-01
    2  -3.50000000e+00   5.00000000e-01   3.20585844e-01
    3  -3.25000000e+00   2.50000000e-01   6.94209267e-02
    4  -3.12500000e+00   1.25000000e-01   6.05288259e-02
    5  -3.18750000e+00   6.25000000e-02   4.61629389e-03
    6  -3.15625000e+00   3.12500000e-02   2.79283147e-02
    7  -3.17187500e+00   1.56250000e-02   1.16471966e-02
    8  -3.17968750e+00   7.81250000e-03   3.51301957e-03
    9  -3.18359375e+00   3.90625000e-03   5.52273640e-04
   10  -3.18164062e+00   1.95312500e-03   1.48021741e-03
   11  -3.18261719e+00   9.76562500e-04   4.63932552e-04
Function value is below tolerance

c =  -3.18310547e+00, |f(c)| =   4.41804335e-05

Example 2

Let us define three different functions for which we will find the root.

def f1(x):
    f = exp(x) - sin(x)
    return f

def f2(x):
    f = x**2 - 4.0*x*sin(x) + (2.0*sin(x))**2
    return f

def f3(x):
    f = x**2 - 4.0*x*sin(x) + (2.0*sin(x))**2 - 0.5
    return f

The following function implements the bisection method. Note that it takes some default arguments.

# Returns (root,status)
def bisect(fun,a,b,N=100,eps=1.0e-4,delta=1.0e-4,debug=False):
    fa, fb = fun(a), fun(b)
    sa, sb = sign(fa), sign(fb)

    if abs(fa) < delta:
        return (a,0)

    if abs(fb) < delta:
        return (b,0)

    # check if interval is admissible
    if fa*fb > 0.0:
        if debug:
            print("Interval is not admissible\n")
        return (0,1)

    for i in range(N):
        e = b-a
        c = a + 0.5*e
        if abs(e) < eps*abs(c):
            if debug:
                print("Interval size is below tolerance\n")
            return (c,0)
        fc = fun(c)
        if abs(fc) < delta:
            if debug:
                print("Function value is below tolerance\n")
            return (c,0)
        sc = sign(fc)
        if sa != sc:
            b = c
            fb= fc
            sb= sc
        else:
            a = c
            fa= fc
            sa= sc
        if debug:
            print("%5d %16.8e %16.8e %16.8e"%(i+1,c,abs(b-a),abs(fc)))

    # If we reached here, then there is no convergence
    print("No convergence in %d iterations !!!" % N)
    return (0,2)

You call this as

root, status = bisect(...)

and check the value of status == 0 to make sure the computation was successfull.

First function

M     = 100        # Maximum number of iterations
eps   = 1.0e-4     # Tolerance on the interval
delta = 1.0e-4     # Tolerance on the function
a, b  = -4.0, -2.0 # Initial interval
r,status = bisect(f1,a,b,M,eps,delta,True)

    1  -3.00000000e+00   1.00000000e+00   1.90907076e-01
    2  -3.50000000e+00   5.00000000e-01   3.20585844e-01
    3  -3.25000000e+00   2.50000000e-01   6.94209267e-02
    4  -3.12500000e+00   1.25000000e-01   6.05288259e-02
    5  -3.18750000e+00   6.25000000e-02   4.61629389e-03
    6  -3.15625000e+00   3.12500000e-02   2.79283147e-02
    7  -3.17187500e+00   1.56250000e-02   1.16471966e-02
    8  -3.17968750e+00   7.81250000e-03   3.51301957e-03
    9  -3.18359375e+00   3.90625000e-03   5.52273640e-04
   10  -3.18164062e+00   1.95312500e-03   1.48021741e-03
   11  -3.18261719e+00   9.76562500e-04   4.63932552e-04
Function value is below tolerance

Second function

M     = 100        # Maximum number of iterations
eps   = 1.0e-4     # Tolerance on the interval
delta = 1.0e-4     # Tolerance on the function
a, b  = -4.0, -2.0 # Initial interval
r,status = bisect(f2,a,b,M,eps,delta,True)

Interval is not admissible

Lets visualize this function.

x = linspace(-3,3,500)
plot(x,f2(x))
grid(True)

It looks like there are double roots; bisection method cannot compute such roots since the function value does not change around the root !!!

Third function

M     = 100        # Maximum number of iterations
eps   = 1.0e-4     # Tolerance on the interval
delta = 1.0e-4     # Tolerance on the function
a, b  = -3.0, +2.0 # Initial interval
r,status = bisect(f3,a,b,M,eps,delta,True)

    1  -5.00000000e-01   2.50000000e+00   2.89455689e-01
    2  -1.75000000e+00   1.25000000e+00   4.52488254e-01
    3  -2.37500000e+00   6.25000000e-01   4.75412891e-01
    4  -2.06250000e+00   3.12500000e-01   4.10335430e-01
    5  -2.21875000e+00   1.56250000e-01   1.10488056e-01
    6  -2.29687500e+00   7.81250000e-02   1.42093933e-01
    7  -2.25781250e+00   3.90625000e-02   6.27310346e-03
    8  -2.23828125e+00   1.95312500e-02   5.44180845e-02
    9  -2.24804688e+00   9.76562500e-03   2.46591833e-02
   10  -2.25292969e+00   4.88281250e-03   9.34083833e-03
   11  -2.25537109e+00   2.44140625e-03   1.57095736e-03
   12  -2.25659180e+00   1.22070312e-03   2.34178305e-03
   13  -2.25598145e+00   6.10351562e-04   3.83092531e-04
   14  -2.25567627e+00   3.05175781e-04   5.94512219e-04
   15  -2.25582886e+00   1.52587891e-04   1.05854829e-04
Interval size is below tolerance

We dont need to specify all arguments to the function as some of them have default values.

r,status = bisect(f3,a,b,debug=True)

    1  -5.00000000e-01   2.50000000e+00   2.89455689e-01
    2  -1.75000000e+00   1.25000000e+00   4.52488254e-01
    3  -2.37500000e+00   6.25000000e-01   4.75412891e-01
    4  -2.06250000e+00   3.12500000e-01   4.10335430e-01
    5  -2.21875000e+00   1.56250000e-01   1.10488056e-01
    6  -2.29687500e+00   7.81250000e-02   1.42093933e-01
    7  -2.25781250e+00   3.90625000e-02   6.27310346e-03
    8  -2.23828125e+00   1.95312500e-02   5.44180845e-02
    9  -2.24804688e+00   9.76562500e-03   2.46591833e-02
   10  -2.25292969e+00   4.88281250e-03   9.34083833e-03
   11  -2.25537109e+00   2.44140625e-03   1.57095736e-03
   12  -2.25659180e+00   1.22070312e-03   2.34178305e-03
   13  -2.25598145e+00   6.10351562e-04   3.83092531e-04
   14  -2.25567627e+00   3.05175781e-04   5.94512219e-04
   15  -2.25582886e+00   1.52587891e-04   1.05854829e-04
Interval size is below tolerance

3.3.3Convergence rate¶

The bisection method always converges and the root $r$ is reached in the limit

r = \lim_{k \to \infty} c_k

(6)

After $k$ iterations, if we take the mid-point of the bracketing interval as the estimate of the root, the error is bounded by

|r - c_k| \le \frac{L_k}{2} = \frac{b-a}{2^{k+1}}

(7)

Numerical Analysis

Root finding