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)
        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 f1pp(x): return 2.0

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

def f3(x): return x**3.0 - 8.0*x**2.0 + 7.0*x
def f3p(x): return 3.0*x**2.0 - 16.0*x + 7.0
def f3pp(x): return 6.0*x -16.0

# f1
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'))

# f2
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'))

# f3
r1 = newton(f3, f3p, 3)
print "newton =",r1
r1_err = np.array(error)

r2 = bisect(f3, 0.5, 3)
print "bisect =",r2
r2_err = np.array(error)

r3 = regula_falsi(f3, 0.5, 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'))
plt.show()