#!/usr/bin/env python2
# -*- coding: utf-8 -*-



import numpy as np
import matplotlib.pyplot as plt


#Partie 1: ordre des methodes de quadratures

#fonction donnant les 4 methodes de quadrature
def integration_approchee(f, a, b, N):
    X,h=np.linspace(a,b,N+1,retstep=True)
    I=np.empty(4)
    I[0]=h/6*np.sum(f(X[0:N])+4*f((X[0:N]+X[1:N+1])/2)+f(X[1:N+1])) #Simpson
    I[1]=h/2*np.sum(f(X[0:N])+f(X[1:N+1])) #trapezes
    I[2]=h*np.sum(f(X[0:N])) #rectangles a gauche
    I[3]=h*np.sum(f(X[1:N+1])) #rectangles a droite
    return I

#fonction f pour tester
def f(x):
    return np.sin(x)

    
a=0
b=np.pi/2

#fonction pour determiner l'erreur de chaque methode en fonction de n
def erreur_integration(f,a,b,valTh,nmax):
    err=np.empty((nmax,4)) #nmax: n va prendre les valeurs 1,...,nmax
    for i in range(nmax):
        I=integration_approchee(f,a,b,i+1)
        err[i,:]=np.abs(valTh-I)
    
    n=np.arange(1,nmax+1)
    plt.plot(np.log(n),np.log(err[:,0]),"r")
    plt.plot(np.log(n),np.log(err[:,1]),"b")
    plt.plot(np.log(n),np.log(err[:,2]),"m")
    plt.plot(np.log(n),np.log(err[:,3]),"k")
    plt.xlabel("log(n)")
    plt.ylabel("log(err)")
    plt.title("Ordre de convergence des methodes de quadrature")
    plt.legend(["Simpson","Trapezes","Rectangles a gauche","Rectangles a droite"])
    plt.show()
    return err


#application
nmax=50
err=erreur_integration(f,a,b,1,nmax) #la valeur theorique est 1


#on peut trouver l'ordre de convergence de chaque methode
n=np.arange(1,nmax+1)
print("Exemple f(x)=sin(x)")
print("Ordre de convergence de la methode de Simpson: ",round(-np.polyfit(np.log(n),np.log(err[:,0]),1)[0]))
print("Ordre de convergence de la methode des trapezes: ",round(-np.polyfit(np.log(n),np.log(err[:,1]),1)[0]))
print("Ordre de convergence de la methode des rectangles a gauche: ",round(-np.polyfit(np.log(n),np.log(err[:,2]),1)[0]))
print("Ordre de convergence de la methode des rectangles a droite: ",round(-np.polyfit(np.log(n),np.log(err[:,3]),1)[0]))




#Bonus: un deuxieme exemple
def g(x):
    return np.sqrt(x)

    
a=0
b=1

#application
nmax=50
err=erreur_integration(g,a,b,2./3.,nmax) #la valeur theorique est 2/3


#on peut trouver l'ordre de convergence de chaque methode dans cet exemple
n=np.arange(1,nmax+1)
print("Exemple bonus f(x)=sqrt(x)")
print("Ordre de convergence de la methode de Simpson: ",round(-np.polyfit(np.log(n),np.log(err[:,0]),1)[0]))
print("Ordre de convergence de la methode des trapezes: ",round(-np.polyfit(np.log(n),np.log(err[:,1]),1)[0]))
print("Ordre de convergence de la methode des rectangles a gauche: ",round(-np.polyfit(np.log(n),np.log(err[:,2]),1)[0]))
print("Ordre de convergence de la methode des rectangles a droite: ",round(-np.polyfit(np.log(n),np.log(err[:,3]),1)[0]))
print("Ordre 1 car la fonction g n'est pas assez reguliere")





#Partie 2: methode de Newton

#fonction f 
def f(x):
    return x+np.exp(x)+10/(1+x**2)-5

tol=1e-8 #on se donne un seuil de tolerance, on arretera si l'ecart entre 2 iteres successifs est inferieur au seuil

y=np.linspace(-2,2,100)
plt.plot(y,f(y))
plt.legend("fonction f")
plt.show()

#fonction derivee
def df(x):
    return 1+np.exp(x)-20*x/((1+x**2)**2)

#calcul approche du zero de f par une methode deja implementee
#import scipy.optimize as scopt
#zerof=scopt.fsolve(f,-0.1)
#print("approximation du zero de f=",zerof)


#fonction pour appliquer la methode de Newton
def newton(f,df,x0,tol,nmax):
    valeurs=[x0]
    x=x0
    err=tol+1
    
    niter=0
    while err>tol and niter<nmax: #arret si l'un des criteres est verifie
        x=x-f(x)/df(x)
        valeurs.append(x)
        err=np.abs(valeurs[niter+1]-valeurs[niter])
        niter=niter+1
        
    if niter==nmax and err>tol:
        print("Arret de la methode avant convergence")
        
    return valeurs, niter
#valeurs contient la suite des iteres, niter le nombre d'iterations
#au plus nmax iterations, mais l'algorithme peut terminer avant

#application a l'exemple
suite,niter= newton(f,df,-0.1,tol,100)

#suite des iteres en fonction de n
plt.plot(np.arange(niter+1),suite)
plt.legend(["suite des iteres de la methode de Newton"])
plt.show()
zerof=-0.9
#erreur et determination de l'ordre
I0=np.abs(suite[0:niter-1]-zerof)
I1=np.abs(suite[1:niter]-zerof) 
plt.plot(np.log(I0),np.log(I1))
plt.legend(['ordre de convergence de la methode de Newton'])
plt.show()


print("Ordre :", np.polyfit(np.log(I0), np.log(I1), 1)[0])



#essai avec une deuxieme condition initiale
x0=1.
zerof=scopt.fsolve(f,x0)
print(zerof)

zero,niter= newton(f,df,x0,tol,200)

plt.plot(np.arange(niter+1),zero)
plt.legend("suite des iteres avec une autre condition initiale")
plt.show()
zero
niter


#plt.plot(np.log(np.abs(zerof-zero[0:niter-1])),np.log(np.abs(zerof-zero[1:niter])))
#np.polyfit(np.log(np.abs(zerof-zero[0:niter-1])),np.log(np.abs(zerof-zero[1:niter])),1)[0]
#plt.show()
#zero
#