#!/usr/bin/env python
# coding: utf-8

# In[ ]:


#does one iteration of Newton's method
def newton(f, df, a):
    x = a - f(a)/df(a)
    return x


# In[8]:


#does Newton's method
#until the absolute error (difference between two iteration)
#is less than a given tolerance

def newton_while(f, df, a):
    not_converged = True
    tol = 1E-6
    max_iter = 100
    x = a
    iter = 0
    while not_converged:
        change = f(x)/df(x)
        not_converged = (abs(change) > tol)
        x = x - change
        iter = iter + 1
        if iter > max_iter:
            print('No convergence after maximum number of iterations')
            break
    return x


# In[9]:


f = lambda x: x**2 - 2
df = lambda x: 2*x
newton_while(f, df, 1)


# In[30]:


newton_while(f, df, -1)


# In[31]:


newton_while(f, df, 0)


# In[7]:


from numpy import sqrt
sqrt(2)


# In[16]:


def bisection(f, a, b, n):
    fa = f(a)
    fb = f(b)
    if (fa >= 0) | (fb <= 0):
        raise ValueError("f should be negative at a and positive at b")
    for i in range(n):
        c = (a + b) / 2
        fc = f(c)            
        if fc == 0:
            return c
        elif fc > 0:
            b = c
            fb = fc
        else:
            a = c
            fa = fc
    return c   


# In[17]:


bisection(f, 1, 2, 20)


# In[25]:


def secant_while(f, a, b):
    not_converged = True
    tol = 1E-6
    max_iter = 100
    fa = f(a)
    iter = 0
    while not_converged:
        fb = f(b)
        s = (f(b) - f(a)) / (b - a) 
        change = fb/s
        not_converged = (abs(change) > tol)
        a = b
        fa = fb
        b = b - change
        print([a, b])
        iter = iter + 1
        if iter > max_iter:
            print('No convergence after maximum number of iterations')
            break
    return b


# In[26]:


secant_while(f, 1, 2)


# In[28]:


#false position method
def fp(f, a, b, n):
    fa = f(a)
    fb = f(b)
    if (fa >= 0) | (fb <= 0):
        raise ValueError("f should be negative at a and positive at b")
    for i in range(n):
        s = (f(b) - f(a)) / (b - a) 
        change = fb/s
        c = b - change
        fc = f(c)            
        if fc == 0:
            return c
        elif fc > 0:
            b = c
            fb = fc
        else:
            a = c
            fa = fc
    return c   


# In[29]:


fp(f, 1, 2, 20)


# In[32]:


from scipy.optimize import root_scalar
root_scalar(f, x0=1)


# In[34]:


root_scalar(f, bracket=[1, 2])


# In[35]:


root_scalar(f, bracket=[-1, -2])


# In[37]:


#does one iteration of Newton's method for a system of equations
from numpy.linalg import solve
def newton_s(f, J, a):
    x = a - solve(J(a), f(a))
    return x


# In[38]:


from numpy import array
f = lambda x: array([x[0]**4 + x[1]**4 - 2,
               x[0] - 2*x[1]])
J = lambda x: array([[4*x[0]**3, 4*x[1]**3],
               [1, -2]])


# In[39]:


a = [1, 1]
newton_s(f, J, a)


# In[40]:


for i in range(5):
    a = newton_s(f, J, a)
    print(a)


# In[41]:


f(a)


# In[42]:


from scipy.optimize import root
root(f, [1, 1])


# In[43]:


root(f, [1, 1], jac=J)


# In[49]:


from numpy.linalg import inv
inv(J(a))


# In[ ]: