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.0
b = 2.0
eps = 1e-10
root,success = bisect_method(a, b, eps, xsqm2)
print( root, success )