## Résolution équation f(x) = 0
## Armel MARTIN - Septembre 2021
## distribué sous licence CC-by-nc-sa
import numpy as np
import matplotlib.pyplot as plt
import time

def f(x): # fonction dont on cherche le zéro
    return 0.15*( 3.7/2/np.sqrt(x) -3.4 - 1.2*x**3 ) 

def fp(x): # dérivée de f(x) pour la méthode de Newton
    return 0.15*( -3.7/4*x**(-3/2) - 1.2*3*x**2 )

## Graphe de la fonction
xx = np.linspace(0.1,0.8,100)
ff = f(xx)
plt.plot(xx,ff,'b-')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
#plt.show()

## Fonctions de résolution pré-compilées SCIPY
import scipy.optimize as opt
tic = time.time()
resultat = opt.bisect(f,0.1,0.5,full_output=True)
print('\n *** Dichotomie Scipy :\n',resultat[1],time.time()-tic)
tic = time.time()
resultat = opt.newton(f,0.1,fprime=fp,full_output=True)
print('\n *** Newton Scipy :\n',resultat[1],time.time()-tic)
tic = time.time()
resultat = opt.fsolve(f,0.1)
print('\n *** fsolve Scipy : x=',resultat,' ; durée=',time.time()-tic)

## Dichotomie "MAISON"
def dichotomie(f,a,b,eps = 1e-10):
    '''Calcule la racine de f(x) = 0 sur [a,b] par dichotomie
    avec une précision de l'ordre de eps.'''
    n = 0 # indice de nombre d'itérations
    c = (a+b)/2    
    while (abs(a-b) > eps and n < 500): # critère de sortie double
        c = (a+b)/2 # milieu de l'intervalle courant
        if f(c)*f(a) < 0 : # signes différent
            b = c
        else:
            a = c
        n += 1 
    return (a+b)/2,abs(a-b),n # racine, erreur, nb itérations
               
tic = time.time()
racine,erreur,niter = dichotomie(f,0.1,0.5)
duree = time.time()-tic
print('\n *** Dichotomie MAISON: \n \t x=',racine,'\n \t erreur=',erreur, \
      '\n \t niter=',niter,'\n \t durée =',duree)

## Newton "MAISON"
def newton(f,fp,xini,eps = 1e-10):
    '''Calcule la racine de f(x) = 0 sur par Newton en partant
    de xini avec une précision de l'ordre de eps.'''
    n = 0 # indice de nombre d'itérations
    x = xini # initialisation
    erreur = 1
    while(erreur > eps and n < 500): # critère de sortie double
        xo = x # stocker la valeur précédente pour évaluer l'erreur
        x = x - f(x) / fp(x)
        erreur = abs(x-xo)
        n += 1
    return x,erreur,n # racine, erreur, nb itérations

tic = time.time()
racine,erreur,niter = newton(f,fp,0.1)
duree = time.time()-tic
print('\n *** Newton MAISON: \n \t x=',racine,'\n \t erreur=',erreur, \
      '\n \t niter=',niter,'\n \t durée =',duree)