In [1]:
from numpy import *
In [13]:
def mybisect(g,a,b,tol):
    ga = g(a)
    gb = g(b)
    while abs(a-b)>tol:
        #print(a,b)
        c = (a+b)/2
        gc = g(c)
        if gc==0: return c
        if sign(gc)==sign(ga):
            a  = c
            ga = gc
        else:
            b = c
            gb = gc
    return (a+b)/2
In [15]:
def g1(x): 
    global ncalls
    ncalls += 1
    return cos(x)-x
ncalls = 0
r = mybisect(g1,0,pi/2,1.e-12)
print(r,g1(r),ncalls)
0.7390851332149163 4.08895139969e-13 44
In [2]:
sign(4.3)
Out[2]:
1.0
In [3]:
sign(-18)
Out[3]:
-1
In [4]:
sign(0)
Out[4]:
0

Newton's method

In [16]:
def g1(x): return cos(x)-x
def g1prime(x): return -sin(x)-1
In [17]:
def newton(g,gprime,x,tol):
    while True:
        s = - g(x)/gprime(x)
        x += s
        if abs(s)<tol: return x
In [18]:
r = newton(g1,g1prime,1,1.e-12)
print(r)
0.739085133215
In [23]:
def newton1(g,gprime,x,tol):
    global ncalls
    while True:
        print(x)
        ncalls += 2
        s = - g(x)/gprime(x)
        x += s
        if abs(s)<tol: return x
        
ncalls = 0
r = newton1(g1,g1prime,1,1.e-12)
print(r,g1(r),ncalls)        
1
0.75036386784
0.739112890911
0.739085133385
0.739085133215
0.739085133215 0.0 10