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 {\em non-trivial bracket} for a root of $f$ . 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

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}

(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. The bisection method is shown in Algorithm~(1). We have put an upper limit on the number of iterations and we also specify the stopping criterion in terms of the length of the bracketing interval and/or the function value.

The way we have written the algorithm, we store the entire history of computations 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).

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 no 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)

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)
f = fun(x)
plot(x,f,'r-')
grid(True); xlabel('x'); ylabel('f');

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:
        b  = c
        fb = fc
        sb = sc
    else:
        a  = c
        fa = fc
        sa = sc
    print("%5d %16.8e %16.8e %16.8e"%(i+1,c,abs(b-a),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

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.

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)

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)
y = f2(x)
plot(x,y)
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

1Convergence rate¶

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

\alpha = \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

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

(7)