def bisect_method(f, a, b, eps):
try:
f = f._fast_float_(f.variables()[0])
except AttributeError:
pass
intervals = [(a,b)]
two = float(2); eps = float(eps)
while True:
c = (a+b)/two
fa = f(a); fb = f(b); fc = f(c)
if abs(fc) < eps: return c, intervals
if fa*fc < 0:
a, b = a, c
elif fc*fb < 0:
a, b = c, b
else:
raise ValueError, "f must have a sign change in the interval (%s,%s)"%(a,b)
intervals.append((a,b))
html("<h1>Double Precision Root Finding Using Bisection</h1>")
@interact
def _(f = cos(x) - x, a = float(0), b = float(1), eps=(-3,(-16..-1))):
eps = 10^eps
print "eps = %s"%float(eps)
try:
time c, intervals = bisect_method(f, a, b, eps)
except ValueError:
print "f must have opposite sign at the endpoints of the interval"
show(plot(f, a, b, color='red'), xmin=a, xmax=b)
else:
print "root =", c
print "f(c) = %r"%f(c)
print "iterations =", len(intervals)
P = plot(f, a, b, color='red')
h = (P.ymax() - P.ymin())/ (1.5*len(intervals))
L = sum(line([(c,h*i), (d,h*i)]) for i, (c,d) in enumerate(intervals) )
L += sum(line([(c,h*i-h/4), (c,h*i+h/4)]) for i, (c,d) in enumerate(intervals) )
L += sum(line([(d,h*i-h/4), (d,h*i+h/4)]) for i, (c,d) in enumerate(intervals) )
show(P + L, xmin=a, xmax=b)
|
Click to the left again to hide and once more to show the dynamic interactive window
|