#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Script example running Bisection and Newton's Method

In Python, we're defining the functions first and putting
the code at the end  (In Matlab, we did the opposite).

"""
import math

#%% Function Definitions:

def bisection(f,a,b,tol):
 """
  Bisection Method
     Computes an approximation to f(x)=0 given that the 
     root is bracketed in [a,b] with f(a)f(b)<0.
     Will run until TOL is reached, and will output the solution xc
     
    In this version, we assume the function f is passing both the
    value of the function AND its derivative (so it lines up with 
    Newton's Method')
    
"""
 u=f(a)
 fa=u[0]
 u=f(b)
 fb=u[0]
 if fa*fb >= 0:
        print("Root may not be bracketed.")
        return None
 
 iter=0
 while (b-a)/2>tol:
        iter+=1
        c = (a + b)/2
        u = f(c)
        fc=u[0]
        if fc==0:
            print("Found exact solution.")
            return c
        elif fa*fc < 0:
            b=c
            fb=fc
        else:
            a=c
            fa=fc
        
 print('Finished after %3d iterations.\n' % iter )       
 return (a + b)/2

def NewtonMethod(F,x0,numits,tol):
    '''

    Inputs to Newton's Method:
    ----------
    f - function
        Function for which we are searching for a solution f(x)=0.
        Should return both f(x) AND f'(x) in that order
    x0 - Initial guess for a solution f(x)=0.
    numits - Maximum number of iteration.    
    tol : tolerance (like 1e-6) for stopping.

    Examples
    --------
    >>> def F(x)
            y=x**2 - x - 1
           dy=2*x-1
           return y,dy

    >>> NewtonMethod(F,1,10,1e-8)
    Found solution after 5 iterations.
    1.618033988749989
    '''
    
    for n in range(0,numits):
        g, dg = F(x0)
        if dg==0:
            print('Zero derivative. No solution found.')
            return None
        
        xnew = x0 - g/dg
        
        d=abs(xnew-x0)
        if d < tol:
            print('Found solution after',n,'iterations.')
            return x0
        x0=xnew
        
    print('Exceeded maximum iterations. No solution found.')
    return None


def testfunc01(x):
	y=math.exp(x)-3*x
	dy=math.exp(x)-3
	return y,dy

#%% Main code starts here:

x=bisection(testfunc01,0,1,1e-9)
print(x)

x=NewtonMethod(testfunc01,0,20,1e-9)
print(x)