bisect_method.py

import numpy as np
def bisect_method (a, b, eps, f):
#  Accompanying program for the text
#
#     Classical and Modern Numerical Analysis:
#     Theory, Methods and Practice
#  by Azmy S. Ackleh, Edward J. Allen,
#     R. Baker Kearfott, and Padmanabhan Seshaiyer
#
#  [root, success]  =  bisect_method (a, b, eps, func) returns the result
#  of the method of bisection, with starting interval [a, b],
#  tolerance eps, and with function defined by y  =  f(x).  For example,
#  suppose an m-file xsqm2.m is available in Matlab's working
#  directory, with the following contents:
#      function [y]  =  xsqm2(x)
#      y  =  x^2-2;
#      return
#  Then, issuing
#      [root,success]  =  bisect (1, 2, 1e-10, 'xsqm2')
#  from Matlab's command window will cause an approximation to
#  the square root of 2 that, in the absence of excessive roundoff
#  error, has absolute error of at most 10^{-16}
#  to be returned in the variable root, and success to be set to
#  'true'.
#
#  success is set to 'False' if f(a) and f(b) do not have the same
#  sign.  success is also set to 'False' if the tolerance cannot be met.
#  In either case, a message is printed, and the midpoint of the present
#  interval is returned in the variable root.

#  Algorithm 2.1 is basically used.  However,  Theorem 2.2 is used to
#  compute the number of iterations N required for bisection to achieve
#  that tolerance.  This is to avoid an infinite loop that might be
#  caused by an excessively small tolerance (in which case the midpoint of
#  the interval is repeatedly rounded to one of the end points and the
#  tolerance is never met).

        error = b-a
        fa = f(a)
        fb = f(b)

        #  First, handle incorrect arguments --
        if fa*fb > 0:
                print('Error: f(a)*f(b)>0')
                return a + (b-a)/2, False

        if eps <= 0:
                print('Error: eps is less than or equal to 0')
                return a + (b-a)/2, False

        if b < a:
                print('Error: b < a')
                return (a+b)/2, False


        #  Set N to be the smallest integer such that N iterations of bisection
        #  suffices to meet the tolerance --
        N  =  int(np.ceil( np.log2((b-a)/eps) - 1 ))

        #  This is where we actually do Algorithm 2.1 --
        print(' -----------------------------')
        print('   Error           Estimate   ')
        print(' -----------------------------')
        for i in range(N):
                x =  a + (b-a)/2
                fx = f(x)
                if fx*fa > 0:   # JR This is bad. Will underflow prematurely. Should use np.sign(fx) == np.sign(fa) instead.
                    a = x
                else:
                    b = x
                error = b-a
                print(' {:12.4e} {:12.4e}'.format(error, x))

        #  Finally, check to see if the tolerance was actually met.  (With
        #  additional analysis of the minimum possible relative error as in Theorem
        #  1.6 and Definition 1.5, unreasonable values of epsilon for the particular
        #  floating point system can be determined before the loop on i, avoiding
        #  unnecessary work.)
        error  =  (b-a)/2
        root  =  a + (b-a)/2
        if error > eps:
                print('Error: epsilon is too small for tolerance to be met')
                return root, False

        return root, True

if __name__ == '__main__':

        #  run_bisect_method.m
        #  Accompanying program for the text
        #
        #     Classical and Modern Numerical Analysis:
        #     Theory, Methods and Practice
        #  by Azmy S. Ackleh, Edward J. Allen,
        #     R. Baker Kearfott, and Padmanabhan Seshaiyer
        #
        #  This Matlab script is an example of how to run bisect.m

        from xsqm2 import xsqm2
        a = 1
        b = 2
        eps = 1e-10
        root,success = bisect_method(a, b, eps, xsqm2)
        print( root, success )