import numpy as np
import matplotlib.pyplot as plt


# Soit on définit les fonctions dans des fichiers à part, que l'on appelle grâce à la commande import (plus propre) soit on les définit directement dans le fichier principal (plus rapide), c'est ce que j'ai choisi de faire ici.

##1.1.1.a
def lagrange(f, a, b, n):
    x = np.linspace(a, b, n+1);
    X = np.poly1d([1, 0]);
    P = 0;
    for i in range(n+1):
        L = 1
        for j in range(0, i):
            L = L * (X - x[j]) / (x[i] - x[j])
        for j in range(i+1, n+1):
            L = L * (X - x[j]) / (x[i] - x[j])

        P = P + L * f(x[i]);
    
    return P
    
##1.1.1.b
def f(x):
    return np.sin(x)


##1.1.1.c
N = 100
a = 0
b = 3 * np.pi
nmax = 5

xp = np.linspace(a, b, N)
##1.1.1.d
Psin = [lagrange(f, a, b, n) for n in range(1, nmax+1)]


# bien sûr on aurait aussi pu faire :
"""
Psin=[]
for n in range (nmax):
    Psin[n]=lagrange(f, a, b, n)
"""
## 1.1.1.e
plt.figure(1) #inutile sur spyder. 
plt.plot(xp, f(xp))

for n in range(nmax):
    
    plt.plot(xp, np.polyval(Psin[n], xp))
    
plt.title("Interpolation de Lagrange de la fonction sinus")
leg = ["f"]
for n in range(1, nmax+1):
    #leg.append("n = " + str(n))
    leg=leg+["n = " + str(n)]
plt.legend(leg)

plt.show(block=False) #selon le logiciel utilisé, cette commande est nécessaire ou non. Avec Spyder, la commande
#plt.plot trace le graph mais par exemple avec Pyzo, plt.plot ne fait que définir ce qui doit apparaitre sur le
#graph mais celui ci n'est affiché que lorsque l'on utilise plt.show()


## 1.1.2.a et 1.1.2.b
plt.figure(2) #inutile sur spyder. Sur Pyzo, cette commande per met d'ouvrir une nouvelle fenêtre graphique
Err = []
Errmax = []
for n in range(nmax):
    Err=Err+[(np.amax(np.abs(f(xp) - np.polyval(Psin[n], xp))))]
    Errmax=Errmax+[(1./(4 * (n + 2)) * ((b - a) / (n + 1))**(n+2))]
    
plt.plot(range(1, nmax + 1), Err)
plt.plot(range(1, nmax + 1), Errmax)
leg2=["Erreur", "Maj. theorique"]
plt.legend(leg2)
plt.title("Comparaison erreur et majoration theorique")

plt.show()


#Bonus :  Pour être plus clair, on peut afficher les graphs précédents dans des sous-figures :
plt.subplot(1,2,1)
plt.plot(xp, f(xp))

for n in range(nmax):    
    plt.plot(xp, np.polyval(Psin[n], xp))
plt.title("Interpolation de Lagrange de la fonction sinus")
leg = ["f"]
for n in range(1, nmax+1):
    #leg.append("n = " + str(n))
    leg=leg+["n = " + str(n)]
plt.legend(leg)
    

plt.subplot(1,2,2)
plt.plot(range(1, nmax + 1), Err)
plt.plot(range(1, nmax + 1), Errmax)
leg2=["Erreur", "Maj. theorique"]
plt.legend(leg2)
plt.title("Comparaison erreur et majoration theorique")

plt.show()


##1.2.1.a
def f(x):
    return 1.0 / (1 + x**2)
    

N = 100
a = -5.
b = 5.

xp = np.linspace(a, b, N)

##1.2.1.b

def lagrange(f, a, b, n):
    x = np.linspace(a, b, n+1)
    X = np.poly1d([1, 0])
    P = 0;
    for i in range(n+1):
        L = 1
        for j in range(0, i):
            L = L * (X - x[j]) / (x[i] - x[j])
        for j in range(i+1, n+1):
            L = L * (X - x[j]) / (x[i] - x[j])

        P = P + L * f(x[i])
    
    return P
    
nvals = [2, 4, 8, 12]
P = [lagrange(f, a, b, n) for n in nvals]

plt.plot(xp, f(xp))
for n in range(len(nvals)):
    plt.plot(xp, np.polyval(P[n], xp))
    
plt.title("Interpolation de Lagrange de la fonction f (points equidistants)")
leg = ["f"]
for n in nvals:
    leg.append("n = " + str(n))
plt.legend(leg)

plt.show()


##1.2.1.c

Err = []
for n in range(4):
    Err=Err+[np.amax(np.abs(f(xp) - np.polyval(P[n], xp)))]
    
plt.plot(nvals, Err)
plt.title("Erreur d'interpolation En'")

plt.show()
# On constate que l'erreur augmente avec le nombre de points d'interpolation, ce qui peut 
#paraître paradoxal, mais on voit bien sur la figure précédente que les polynômes de haut degrés oscillent et s'éloignent de la courbe d'origine
# C'est le phénomène de Runge, qui sera corrigé en modifiant la répartition des points xi dans la question suivante

## 1.2.2.a

def lagrangeTcheb(f, a, b, n):
    x = (a +b) / 2 + (b - a) / 2 * np.cos((2 * np.arange(n+1) + 1) * np.pi / (2 * (n + 1)))
    X = np.poly1d([1, 0])
    P = 0;
    for i in range(n+1):
        L = 1
        for j in range(0, i):
            L = L * (X - x[j]) / (x[i] - x[j])
        for j in range(i+1, n+1):
            L = L * (X - x[j]) / (x[i] - x[j])

        P = P + L * f(x[i])
    
    return P


## 1.2.2.b



P = [lagrangeTcheb(f, a, b, n) for n in nvals]

plt.plot(xp, f(xp))
for n in range(len(nvals)):
    plt.plot(xp, np.polyval(P[n], xp))
    
plt.title("Interpolation de Lagrange de la fonction f (point de Tchebytchev)")
leg = ["f"]
for n in nvals:
    leg.append("n = " + str(n))
plt.legend(leg)

plt.show()


ErrT = []
for n in range(4):
    ErrT=ErrT+[(np.amax(np.abs(f(xp) - np.polyval(P[n], xp))))]
    
plt.plot(nvals, ErrT)
plt.title("Erreur d'interpolation En'")

plt.show()