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

# In[1]:


#function to do n iterations of Newton's method
def newton(f, df, x0, n):
    x = x0
    for i in range(n):
        x = x - f(x)/df(x)
        print(x)
    return x


# In[3]:


f = lambda x: x**3 - 3*x + 1
x0 = 2
df = lambda x: 3 * x**2 - 3
n = 10
newton(f, df, x0, n)


# In[4]:


#plotting a function is a good starting point for solving equations that involve it
from matplotlib.pyplot import plot, figure, grid, title, xlabel, ylabel, legend, show
from numpy import linspace
x = linspace(-10, 10, 200)
y = f(x)
plot(x, y)


# In[6]:


x = linspace(-4, 4, 200)
y = f(x)
plot(x, y)
grid()


# In[8]:


#function to do n iterations of bisection
#assumes that f(a) < 0 and f(b) > 0
def bisection(f, a, b, n):
    for i in range(n):
        c = (a + b) / 2 #midpoint
        fc = f(c)
        if fc == 0:
            return c
        elif fc < 0:
            a = c
        elif fc > 0:
            b = c
        print([a, b])
    return c


# In[9]:


bisection(f, 1, 2, 10)


# In[10]:


#function to do n iterations of secant method
def secant(f, a, b, n):
    for i in range(n):
        fa = f(a)
        fb = f(b)
        s = (fb - fa) / (b - a)
        a = b
        b = b - fb/s
        print([a, b])
    return b


# In[11]:


secant(f, 1, 2, 10)


# In[13]:


#function to do n iterations of false position method
#assumes that f(a) < 0 and f(b) > 0
def fp(f, a, b, n):
    for i in range(n):
        fa = f(a)
        fb = f(b)
        s = (fb - fa) / (b - a)
        c = b - fb/s
        fc = f(c)
        if fc == 0:
            return c
        elif fc < 0:
            a = c
        elif fc > 0:
            b = c        
        print([a, b])
    return c


# In[14]:


fp(f, 1, 2, 10)


# In[16]:


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


# In[20]:


sol


# In[21]:


root_scalar(f, x0=2, fprime=df, method='newton')


# In[22]:


root_scalar(f, x0=0, fprime=df, method='newton') #a different starting point can lead to another root


# In[23]:


#find root of system of equations
f = lambda x: [x[0] + x[1] + x[0]*x[1], 
              x[0]**2 - x[1] + 3]


# In[25]:


f([1, 1])


# In[26]:


J = lambda x: [[1 + x[1], 1 + x[0]], [2*x[0], -1]] #Jacobian of f


# In[27]:


J([0, 0])


# In[31]:


#function to do n iterations of Newton's method (for a system of equations)
def newton_md(f, J, x0, n):
    x = x0
    from numpy.linalg import solve
    for i in range(n):
        x = x - solve(J(x), f(x))
        print(x)
    return x


# In[34]:


x = newton_md(f, J, [0, 0], 10)


# In[36]:


f(x)


# In[37]:


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


# In[38]:


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


# In[41]:


root(f, [1, 0], jac=J) #other starting points seem to still give the same root as the [0, 0] starting point


# In[ ]: