scientific-programming-exercises/ex_46.py

import numpy as np

def bisec(f, a, b, eps, nmax):
	"""
	Takes f:[a,b] -> |R and tries to compute the null point 
	using bisection.
	"""

	if(f(a)*f(b) >= 0):
		raise ValueError("f(a)*f(b) >= 0")
	
	x_minus = a
	x_plus = b

	if(a > b):
		x_minus = b
		x_plus = a

	for i in range(nmax):
		x_minus, x_plus, err = bisec_one(x_minus, x_plus, f, eps)
		if(err < eps):
			break
	if(err < eps):
		return x_minus, err
	raise ValueError("bisection hit nmax")

	
def bisec_one(x_minus, x_plus, f, eps):
	a = x_minus
	b = x_plus

	x_m = (b + a) / 2
	y_m = f(x_m)

	if(y_m < 0):
		if(y_m > -1*eps):
			return x_m, x_m, -y_m

		if(f(x_m)*f(b) >= 0):
			return a, x_m, -y_m
		return x_m, b, -y_m
	
	if(y_m < eps):
		return x_m, x_m, y_m
	
	if(f(a)*f(x_m) >= 0):
		return x_m, b, y_m
	return a, x_m, y_m


if( __name__ == "__main__"):
	f1 = lambda x: x
	f2 = lambda x: x**3
	f3 = lambda x: -x + 1
	f4 = lambda x: -x**3 + 1
	f4 = lambda x: (x - 2)*np.exp(-x**2)

	fs = [f1, f2, f3, f4]
	for f in fs:
		print(bisec(f, -12, 10, 0.0000001, 100))

print(bisec(f4, 1.2, 2.4, 0.001, 100))