DZone Snippets is a public source code repository. Easily build up your personal collection of code snippets, categorize them with tags / keywords, and share them with the world

Snippets has posted 5883 posts at DZone. View Full User Profile

Bisection, Rope And Newton Methods

06.15.2008
| 8388 views |
  • submit to reddit
        
from math import log

def bisect(f, a, b, e):
	""" Determines zero between a and b using Bisection. """
	n = 0
	fa = f(a)
	if fa == 0.0: return (a, n)
	fb = f(b)
	if fb == 0.0: return (b, n)
		
	while (abs(a-b) > e):
		c = 0.5*(a+b)
		fc = f(c)
		
		if fc == 0.0: return (c, n)
		n = n + 1
		if fb*fc < 0.0:
			a = c
			fa = fc
		
		else:
			b = c
			fb = fc
	
		
	if fa < fb:
		return (a, n)
	else:
		return (b, n)

def rope(f, a, b, e):
	""" Determines zero between a and b using the Rope methode. """
	n = 0
	fa = f(a)
	if fa == 0.0: return (a, n)
	fb = f(b)
	if fb == 0.0: return (b, n)
	
	while (abs(a-b) > e):
		c = (a*fb - b*fa) / (fb - fa)
		fc = f(c)
		if fc == 0.0: return (c, n)
		n = n + 1
		if fb*fc < 0:
			a = c
			fa = fc
		
		else:
			b = c
			fb = fc
		
	
	if fa < fb:
		return (a, n)
	else:
		return (b, n)

# Note: must verify that for the function f and guess c
#		the method will _converge_.
def newton(f, df, c, t):
	""" Determines zero between a and b using Newton """
	n = 0
	fc = f(c)
	if fc == 0.0: return (c, n)
	
	
	while (True):
		fc = f(c)
		dfc = df(c)
		if dfc == 0:
			print "dfc is 0"
			return (0, -1)
		
		dc = -fc/dfc
		
		c = c + dc
		n = n + 1
		print c
		if abs(dc) < t: return (c, n)
	
		

##Tests
#def f(x): return -log(x)+4.0
#def df(x): return -1.0/x
#x= bisect(f, 1, 70, 0.00000001)
#print x
#x = rope(f, 1, 70, 0.00000001)
#print x
#x = newton(f, df, 0.1, 0.0001)
#print x

def f(x): return (x+2)*log(2*x**3)+3.0
def df(x): return log(2*x**3)+(3*(x+2))/x

x,y = newton(f, df, 1.0, 0.0001)
print x,y
    

Comments

Snippets Manager replied on Wed, 2011/01/19 - 11:20am

Good bisect() works with python3. Thank you.