import numpy as np
import matplotlib.pyplot as plt


print('------ 3.16 ---------------')
print("Détermination d'un polynôme d'interpolation.")
ax = np.array([-2, -1, 0, 1, 2], np.float64) # Donnée
ay = np.sin(ax) # Donnée

# Création de la matrice
matA = np.zeros((ax.shape[0], ax.shape[0]), np.float64)
matA[:, 0] = np.ones(ax.shape[0]) # Première colonne de la matrice

# Remplissage des autres colonnes de la matrice
for jj in range(1, ax.shape[0]):
    matA[:, jj] = ax ** jj

print("A =\n", matA)

# Détermination des coefficients du polynôme
coefs = np.linalg.lstsq(matA, ay, rcond=None)[0] # Coefficients du polynôme
print("Solution = ", coefs)

# Graphique, les points, le sinus et le polynôme.
xx = np.linspace(-2.5, 2.5, 21,endpoint=True) # Abscisses des points pour tracer le polynôme
yy = coefs[0] + coefs[1]*xx + coefs[2]*xx**2 + coefs[3]*xx**3 + coefs[4]*xx**4
plt.plot(ax, ay, 'o') # Affiche les points
plt.plot(xx, yy, '-') # Trace le polynôme
plt.plot(xx, np.sin(xx), '-') # Trace le sinus
plt.title('Polynôme passant par les points.')
plt.xlabel('x')
plt.ylabel('y = sin(x)')
plt.legend(['points', 'y=poly(x)', 'y=sin(x)'], loc='upper left')
plt.axis([-2.5, 2.5, -1.5, 1.5]) # Limites des axes
plt.show()

print('------ 3.16 bis ---------------')
ax = np.array([-2, -1, 0, 1, 2], np.float64) # Donnée
ay = np.sin(ax) # Donnée

coefs = np.polyfit(ax, ay, 4) # Coefficients du polynôme de degré 4
print("Solution = ", coefs)
coefs = np.flip(coefs) # Change l'ordre des nombres dans le vecteur
print("Solution = ", coefs)

# Graphique, les points, le polynôme et le sinus.
xx = np.linspace(-2.5, 2.5, 81, endpoint=True) # Abscisses des points pour tracer le polynôme
yy = coefs[0] + coefs[1]*xx + coefs[2]*xx**2 + coefs[3]*xx**3 + coefs[4]*xx**4
plt.plot(ax, ay, 'o') # Affiche les points
plt.plot(xx, yy, '-') # Trace le polynôme
plt.plot(xx, np.sin(xx), '-') # Trace le sinus
plt.title('Polynôme passant par les points.')
plt.xlabel('x')
plt.ylabel('y = sin(x)')
plt.legend(['points', 'y=poly(x)', 'y=sin(x)'], loc='upper left')
plt.axis([-2.5, 2.5, -1.5, 1.5]) # Limites des axes
plt.show()