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 )