import numpy as np np.set_printoptions(suppress=True) # (pour mieux arrondir à 0 lors des print) import matplotlib.pyplot as plt # On active le "mode interactif" de pyplot. Ainsi, les plt.show() ne sont plus nécessaires. plt.ion() # -------------------------------------------------------------------------------------- # Fonction de "fin temporaire" pendant l'avancée du TP. Attend une frappe Entrée, puis quitte le script. def stop_ici(): plt.pause(0.1) # nécessaire pour que matplotlib ait le temps d'afficher les figures input() exit() # -------------------------------------------------------------------------------------- # Tous les passages indiqués "TODO()" sont à remplacer par vos soins def TODO(): print("à vous!") stop_ici() ################################################################################################# ### ### EXERCICE 2 : Interpolation polynômiale en 1D ### ################################################################################################# ### Mise en jambes : afficher le graphe d'un polynôme ##################################################### # Aux TPs précédents, vous avez vu comment Numpy permet facilement d'appliquer une fonction directement # à l'ensemble d'un tableau de valeurs print('''On considère le polynôme P(X)=3-X+X^2. Utilisez Numpy pour calculer P(t) aux points t=0, 0.5, 1, 1.5, ..., 4''') t = np.linspace(0,4,9) print("t=\n",t) Pt = 3-t+t**2 print("P(t)=\n",Pt) # stop_ici() # -------------- Supprimez cette ligne pour passer à la suite ------------------ # -------------------------------------------------------------------------------------- print("-"*80) # Toutefois, dans ce TP, on va calculer les valeurs P(t) avec une autre méthode, encore plus puissante, # basée directement sur les coefficients du polynôme. Le principe est le suivant: # * Un polynôme P(X) de degré N-1 peut être représenté par une liste "p" de taille N, contenant ses N coefficients. # Les coefficients doivent être stockés en commençant par le plus petit ordre (= le terme constant). # Par exemple, le polynôme 1+3X-7X^2 est représenté par la liste p=[1,3,-7] # * Le module numpy.polynomial.polynomial possède ensuite une fonction intitulée "polyval" # qui permet d'évaluer directement un polynôme de coefficients "p" aux valeurs de la variable données par "t" from numpy.polynomial.polynomial import polyval print('''Définissez une liste "p" représentant le polynôme P(X) = 3 -X +X^2 ''') p = [3, -1, 1] Pt = polyval(t,p) # La ligne importante ! print("p (coefficients)=\n",p) print("t=\n",t) print("P(t)=\n",Pt) # Vérification : dans P(t) vous devez retrouver exactement les mêmes valeurs qu'à la question précédente ! # Si c'est bien le cas, vous pouvez passer à la suite. # stop_ici() # -------------- Supprimez cette ligne pour passer à la suite ------------------ # -------------------------------------------------------------------------------------- print("-"*80) print('''En utilisant la fonction polyval, représentez le graphe du polynôme P(X) = 2 -3X -4X^2 +X^3 sur l'intervalle [-5,5]''') # Coefficients du polynôme P(X) p = [2, -3, -4, 1] # Discrétisation de l'intervalle [-5,5] en 200 points t = np.linspace(-5,5,200) # Calcul de P(t) pour chacune des valeurs dans le tableau t Pt = polyval(t,p) # Création d'une figure, et représentation de la courbe (t,P(t)) sur l'intervalle [-5,5] plt.figure() plt.plot(t,Pt) plt.grid() # stop_ici() # --------------- Supprimez cette ligne pour passer à la suite ------------- # -------------------------------------------------------------------------------------- print("-"*80) ### Interpolation polynômiale quadratique (= de degré 2) ######################################################## # # On arrive maintenant dans le vif du sujet : l'INTERPOLATION polynômiale # (c'est-à-dire le calcul d'un polynôme qui passe par un ensemble de points donnés). # ### On définit 3 points "cible" # Abscisses (t) des 3 cibles ts = 1,2,3 # Ordonnées (y) des 3 cibles ys = 4,-1,0 # Représentation graphique des points cible plt.figure() plt.plot(ts,ys,linestyle='none',marker='o') plt.grid() plt.title("Problème d'interpolation") # stop_ici() # --------------- Supprimez cette ligne pour passer à la suite ------------- # -------------------------------------------------------------------------------------- print("-"*80) print('''(Question papier) On cherche un polynôme P(X) de degré 2 P(X) = u + v.X + w.X^2 qui INTERPOLE ces 3 points, c'est-à-dire qui vérifie P(t)=y pour chacun des 3 points cible. Soit, ici : P(1) = 4 P(2) = -1 P(3) = 0 Montrez que ces 3 équations constituent un *système linéaire* sur les 3 coefficients inconnus (u,v,w). Précisez la matrice de ce système (appelons-la A), et le vecteur du second membre (appelons-le b). ''') A = np.array([[1,1,1],[1,2,4],[1,3,9]]) b = np.array([4,-1,0]) # stop_ici() # -------------- Supprimez cette ligne pour passer à la suite ------------------ # -------------------------------------------------------------------------------------- print("-"*80) print("Résolvez ce système linéaire, c'est-à-dire trouvez la valeur des coefficients (u,v,w) recherchés.") p = np.linalg.solve(A,b) # stop_ici() # -------------- Supprimez cette ligne pour passer à la suite ------------------ # -------------------------------------------------------------------------------------- print("-"*80) print("Dans la même figure que précédemment, représentez le graphe du polynôme P. Vérifiez qu'il INTERPOLE bien les 3 points imposés.") plt.figure() plt.plot(ts,ys,linestyle='none',marker='o') t = np.linspace(0,4,200) Pt = polyval(t,p) plt.plot(t,Pt) plt.grid() plt.title("Problème d'interpolation") # stop_ici() # --------------- Supprimez cette ligne pour passer à la suite ------------- # -------------------------------------------------------------------------------------- print("-"*80) ### Interpolation polynômiale quelconque (= de degré N) ######################################################## # On considère maintenant les abscisses # t=1,2,3,...,N # et on se donne N valeurs 'cible' en ordonnées : # y1,y2,...,yN ys = np.array([ 3,1,2,-2,-4 ]) N = len(ys) # (Ici, N=5, mais votre code doit être robuste si cette valeur change) plt.figure() plt.plot(range(1,N+1),ys,linestyle='none',marker='o') plt.grid() plt.title("Problème d'interpolation") # stop_ici() # --------------- Supprimez cette ligne pour passer à la suite ------------- # -------------------------------------------------------------------------------------- print("-"*80) # On peut montrer qu'il existe toujours un UNIQUE polynôme INTERPOLATEUR de degré N-1 qui vérifie #P(1) = y1 #P(2) = y2 #... #P(N) = yN # Nota : en réalité l'interpolation marcherait pour n'importe quel choix de N nombres "ti" en abscisses (tant qu'ils sont deux à deux différents). Dans ce TP, on choisit que les abscisses "ti" soient 1,2,...,N pour simplifier. print('''(Question papier) On cherche un polynôme P(X) de degré 4 P(X) = u + v.X + w.X^2 + z.X^3 + r.X^4 qui INTERPOLE les 5 points yi donnés en entrée, c'est-à-dire qui vérifie P(i)=yi pour i allant de 1 à 5. (Soit dans notre exemple, concrètement: P(1)=3 P(2)=1 P(3)=2 P(4)=-2 P(5)=-4) Montrez que ces 5 équations constituent un *système linéaire* sur les 5 coefficients inconnus (u,v,w,z,r). Précisez la matrice de ce système (appelons-la A), et le vecteur du second membre (appelons-le b). ''') # stop_ici() # -------------- Supprimez cette ligne pour passer à la suite ------------------ # -------------------------------------------------------------------------------------- print("-"*80) print('''Complétez le code suivant, définissant la matrice A. (Remarque: l'intérêt de cette écriture est qu'elle reste valide lorsqu'on change le nombre de points N à interpoler.''') A = np.array( [ [ i**j for j in range(N) ] for i in range(1,N+1) ] ) print("A=\n",A) # stop_ici() # -------------- Supprimez cette ligne pour passer à la suite ------------------ # -------------------------------------------------------------------------------------- print("-"*80) print('''Trouvez les coefficients du polynôme interpolateur P recherché. Puis tracez son graphe sur l'intervalle [1,N]. Vérifiez qu'il interpole bien les N ordonnées imposées.''') p = np.linalg.solve(A,ys) plt.figure() plt.plot(range(1,N+1),ys,linestyle='none',marker='o') t = np.linspace(1,N,200) Pt = polyval(t,p) plt.plot(t,Pt) plt.grid() plt.title("Problème d'interpolation") # stop_ici() # --------------- Supprimez cette ligne pour passer à la suite ------------- # -------------------------------------------------------------------------------------- print("-"*80) ### Une fonction d'interpolation polynômiale ############################################ print('''Encapsulez votre code des questions précédentes dans une fonction, qui: - prend en entrée une liste d'ordonnées ys=[y1,y2,...,yN] (le nombre pouvant varier) - renvoie l'unique polynôme interpolateur de degré N-1 tel que P(1) = y1 ... P(N) = yN ''') def interpol(ys): N = len(ys) A = np.array( [ [ i**j for j in range(N) ] for i in range(1,N+1) ] ) p = np.linalg.solve(A,ys) return p # Test : ys = np.array( [ 6,-1,4,3,0,2,1,1,1] ) p = interpol(ys) plt.figure() plt.plot(range(1,len(ys)+1),ys,linestyle='none',marker='o') t = np.linspace(1,len(ys),200) plt.plot(t,polyval(t,p)) plt.grid() plt.title("Problème d'interpolation") # Remarque : copiez votre fonction "interpol", elle ressert dans l'exo 3 stop_ici()