import numpy as np
import matplotlib.pylab as plt

def newton(f, fp, x, tolerance=1.0e-12):
    global error
    error = []
    corr = 10
    while abs(corr) > tolerance:
        fx = f(x)
        fpx = fp(x)
        corr = -fx/fpx
        oldx = x
        x += corr
        error.append(abs(corr))
    return x

def bisect(f, a, b, tolerance=1.0e-12):
    global error
    error = []
    fa = f(a)
    fb = f(b)
    if fa*fb > 0:
        raise Exception("Probably no root in interval [%lf, %lf]" % (a, b))

    corr = 10
    while abs(corr) > tolerance:
        m = 0.5*(a+b)
        fm = f(m)
        corr = a-m
        if fa*fm < 0: 
            b = m
            fb = fm
        else: 
            a = m
            fa = fm
        error.append(abs(corr))

    return 0.5*(a+b)

def regula_falsi(f, a, b, tolerance=1.0e-12):
    global error
    error = []
    fa = f(a)
    fb = f(b)
    if fa*fb > 0:
        raise Exception("Probably no root in interval [%lf, %lf]" % (a, b))
    corr = 10

    while abs(corr) > tolerance:
        m = (a*fb-b*fa)/(fb-fa)
        fm = f(m)
        print "a=%lf c=%lf b=%lf" % (a, m, b)
        if fa*fm < 0: 
            corr = b-m
            b = m
            fb = fm
        else: 
            corr = a-m
            a = m
            fa = fm
        error.append(abs(corr))
    return m


def f1(x): return x**2.0 - 1.0
def f1p(x): return 2.0*x

def f2(x): return np.cos(x)
def f2p(x): return -np.sin(x)

r1 = newton(f1, f1p, 3)
print "newton =",r1
r1_err = np.array(error)

r2 = bisect(f1, 0, 3)
print "bisect =",r2
r2_err = np.array(error)

r3 = regula_falsi(f1, 0, 3)
print "regula_falsi =",r3
r3_err = np.array(error)

plt.figure()
plt.semilogy(r1_err, 'o-')
plt.semilogy(r2_err, 'o-')
plt.semilogy(r3_err, 'o-')
plt.legend(('newton', 'bisection', 'regula falsi'))

r1 = newton(f2, f2p, 2)
print "newton =",r1
r1_err = np.array(error)

r2 = bisect(f2, 0, 2)
print "bisect =",r2
r2_err = np.array(error)

r3 = regula_falsi(f2, 0, 2)
print "regula_falsi =",r3
r3_err = np.array(error)

plt.figure()
plt.semilogy(r1_err, 'o-')
plt.semilogy(r2_err, 'o-')
plt.semilogy(r3_err, 'o-')
plt.legend(('newton', 'bisection', 'regula falsi'))

# 6.2.1

def f3(x): return x**3 - 13*13
def f3p(x): return 3.0*x**2

print "13**(2/3)=", newton(f3, f3p, 13), " (via newton)"
print "13**(2/3)=", 13**(2./3), " (via python)"
    
def root(x, k):
    f = lambda y: y**k-x
    fp = lambda y: k*y**(k-1.)
    return newton(f, fp, x)

print "13**(1/3)=", root(13,3), " (via root)"
print "13**(1/3)=", 13**(1./3), " (via python)"

plt.show()