TP2 Interpolation Polynomiale Corrige [PDF]

Zitiervorschau

Calcul Scientifique TP2 corrig?

October 6, 2019

1 1.1

TP2 : Interpolation polynomiale Objectif

L’objectif de ce TP est : • d’interpoler un nombre de points donné par un polynôme en utilisant les méthodes d’interpolation de Lagrange et Newton, • d’étudier l’évolution de l’erreur d’interpolation, d’une fonction, en fonction du nombre de points à interpoler.

1.2

Rappel sur l’interpolation

Soient n + 1 points ( x0 , y0 ), ( x1 , y1 ) . . . .. ( xn , yn ). Interpoler ces points correspond à déterminer le polynôme P ∈ Rn [ X ] passant par ces derniers : ∀i ∈ {0, · · · , n}, P( xi ) = yi . Les abscisses ( xi )0≤i≤n et les ordonnées (yi )0≤i≤n sont appelées, respectivement, les points et les valeurs d’interpolation. Pour interpoler une fonction f , on définit ses valeurs d’interpolation comme suit : yi = f ( xi ), ∀ 0 ≤ i ≤ n. Dans ce qui suit, nous présentons deux méthodes d’interpolation : par les polynômes de Lagrange et par les polynômes de Newton.

2 2.1

Pôlynomes d’interpolation de Lagrange Théorème

Soient n + 1 points de coordonnées ( xi , yi )0≤i≤n tels que xi ̸= x j , pour 0 ≤ i, j ≤ n et i ̸= j. Il existe alors un unique polynôme d’interpolation de Lagrange Pn ∈ Rn [ X ] vérifiant Pn ( xi ) = yi , ∀i ∈ {0, · · · , n}. Ce polynôme s’exprime comme : n

Pn ( x ) =

∑ y i L i ( x ),

i =0 n

où Li ( x ) =

x − xj

∏ xi − x j . j =0 j ̸ =i

1

x∈R

La famille de polynômes de Lagrange { L0 , L1 , · · · , Ln } associés aux points ( xi , yi )0≤i≤n est une base de l’ensemble Rn [ X ]. Ecrire une fonction Lagrange(t,i,x) qui évalue, au point t, le polynôme de Lagrange Li , i ∈ {0, · · · , n} associé aux points d’interpolation x = ( xi )0≤i≤n . [1]: import numpy as np import matplotlib.pyplot as plt [2]: def Lagrange(t,i,x): n=len(x) L=1 for j in np.arange(0,n): if j!=i: L*=(t-x[j])/(x[i]-x[j]) return L Tester la fonction Lagrange(t,i,x)sur les points {−1, 0, 1}. Représenter les polynômes L0 , L1 et L2 sur $[-1,1]$. [3]: x=np.arange(-1,2,1) t=np.linspace(-1,1,21) # permet dobtenir un tableau 1D allant de -1 à 1␣ ,→contenant 21 éléments. plt.figure(figsize=(20,10)) plt.plot(t,Lagrange(t,0,x),'ro--',t,Lagrange(t,1,x),'b^--',t,␣ ,→Lagrange(t,2,x),'g*--',linewidth=3,markersize=12) plt.xlabel('t',fontsize=30) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.legend(('L0','L1','L2'),fontsize=20, loc = 0) plt.grid(True) plt.text(-1,-0.1,"(-1,0)",ha="center",va="top",fontsize=30) plt.text(0,-0.1,"(0,0)",ha="center",va="top",fontsize=30) plt.text(1,-0.1,"(1,0)",ha="center",va="top",fontsize=30) plt.text(-1,1.05,"(-1,1)",ha="center",va="bottom",fontsize=30) plt.text(0,1.05,"(0,1)",ha="center",va="bottom",fontsize=30) plt.text(1,1.05,"(1,1)",ha="center",va="bottom",fontsize=30) [3]: Text(1, 1.05, '(1,1)')

2

Ecrire une fonction Interpolation_Labgrange(t,x,y) qui évalue, au point t, le polynôme d’interpolation P de Lagrange associé aux points d’interpolation ( xi , yi )0≤i≤n , avec x = ( xi )0≤i≤n et y = (yi )0≤i≤n . Représeneter ensuite le polynôme P graphiquement qui interpole les points (−1, 8), (0, 3) et (1, 6). [4]: def Interpolation_Lagrange(t,x,y): n=len(x) P=np.zeros((len(t))) for i in np.arange(0,n): P+=y[i]*Lagrange(t,i,x) return P [5]: y=[8,3,6] P=Interpolation_Lagrange(t,x,y) plt.figure(figsize=(20,10)) plt.plot(t,P,'mo--',linewidth=3,markersize=12) plt.xlabel('t',fontsize=30) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.grid(True) plt.legend(('Polynome d\'interpolation de Lagrange',),fontsize=20, loc = 0) plt.text(-1,8.5,"(-1,8)",ha="center",va="top",fontsize=30) plt.text(0,3.5,"(0,3)",ha="center",va="top",fontsize=30) plt.text(1,6.5,"(1,6)",ha="center",va="top",fontsize=30) [5]: Text(1, 6.5, '(1,6)')

3

Exercice: 1. Interpoler, par la méthode de Lagrange, la fonction cos aux points d’interpolation : {−π, − π2 , 0, π2 , π }. 3π 2. Sur [− 3π 2 , 2 ], tracer sur un même graphe la courbe de la fonction cos et P son polynôme d’interpolation de Lagrange associé aux points {−π, − π2 , 0, π2 , π }. 3. Interpréter les résultats.

[6]: x=np.linspace(-np.pi,np.pi,5) t=np.linspace(-1.5*np.pi,1.5*np.pi,100) P=Interpolation_Lagrange(t,x,np.cos(x)) plt.figure(figsize=(20,10)) plt.plot(t,P,'r-',t,np.cos(t), 'b--',x,np.cos(x),'mo',linewidth=3,markersize=12) plt.xlabel('t',fontsize=30) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.grid(True) plt.legend(('Polynome d\'interpolation de Lagrange', 'cosinus',␣ ,→'points'),fontsize=20, loc = 0) [6]:

4

Inconvénient Un inconvénient majeur de la méthode d’interpolation par les polynômes de Lagrange réside en l’ajout d’un point ( xn+1 , yn+1 ) à l’ensemble de n points d’interpolation. Dans ce cas, il n’est numériquement pas évident de déduire Pn+1 de Pn . Tous les calculs seront refaits de zéro. D’où l’introduction de l’interpolation par les polynômes de Newton.

3

Pôlynomes d’interpolation de Newton

3.1

Théorème

Soient n + 1 points de coordonnées ( xi , yi )0≤i≤n tels que xi ̸= x j , pour 0 ≤ i, j ≤ n et i ̸= j. Il existe alors un unique polynôme d’interpolation de Newton Pn ∈ Rn [ X ] vérifiant Pn ( xi ) = yi , ∀i ∈ {0, · · · , n}. Ce polynôme s’exprime comme : n

Pn ( x ) =

∑ β i ωi ( x ),

x∈R

(1)

i =0

= β 0 .|{z} 1 + β 1 ( x − x0 ) + β 2 ( x − x0 )( x − x1 ) + .... + β n ( x − x0 )( x − x1 )...( x − xn−1 )(2) {z } | {z } | {z } | ω0

où ωi ( x ) =

ω1

ω2

ωn

i −1

∏(x − xi ), ∀1 ≤ i ≤ n et ω0 (x) = 1. j =0

La famille de polynômes de Newton {ω0 , ω1 , · · · , ωn } associés aux points ( xi , yi )0≤i≤n forme une base de l’ensemble Rn [ X ]. Les réels β i , i ∈ {0, · · · , n} correspondent aux coefficients du polynôme. Pour les déterminer, nous utilisons la méthode des différences divisées. Ecrire une fonction Newton(t,i,x) qui évalue, au point t, le polynôme de Newton ωi , i ∈ {0, · · · , n} associé aux points d’interpolation x = ( xi )0≤i≤n . 5

[7]: def Newton(t,i,x): n=len(x) if i==0: return np.ones((len(t))) else: W=1 for j in np.arange(0,i): W*=(t-x[j]) return W Tester la fonction Newton(t,i,x) sur les points {−1, 0, 1}. Représenter les polynômes ω0 , ω1 et ω2 sur $[-1,1]$. [8]: x=np.arange(-1,2,1) t=np.linspace(-1,1,21) plt.figure(figsize=(20,10)) plt.plot(t,Newton(t,0,x),'ro--',t,Newton(t,1,x),'b^--',t,␣ ,→Newton(t,2,x),'g*--',linewidth=3,markersize=12) plt.xlabel('t',fontsize=30) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.legend(('W0','W1','W2'),fontsize=20, loc = 0) plt.grid(True) plt.text(-1,-0.1,"(-1,0)",ha="center",va="top",fontsize=30) plt.text(0,-0.1,"(0,0)",ha="center",va="top",fontsize=30) plt.text(1,-0.1,"(1,0)",ha="center",va="top",fontsize=30) plt.text(-1,1.05,"(-1,1)",ha="center",va="bottom",fontsize=30) plt.text(0,1.05,"(0,1)",ha="center",va="bottom",fontsize=30) plt.text(1,1.05,"(1,1)",ha="center",va="bottom",fontsize=30) [8]: Text(1, 1.05, '(1,1)')

6

3.2

Différences divisées

3.2.1 Définition On considère (n + 1) points ( xi , yi )0≤i≤n , deux à deux distincts : 1. La différence divisée d’ordre 1 de xi−1 et xi , 0 < i ≤ n est : f [ x i −1 , x i ] =

y i − y i −1 x i − x i −1

2. La différence divisée d’ordre n des n + 1 points est définie par récurrence entre deux différences divisées d’ordre n − 1 comme suit : f [ x0 , x1 , · · · , x n ] =

f [ x 1 , · · · , x n ] − f [ x 0 , x 1 , · · · , x n −1 ] x n − x0

3.2.2 Explication de la méthode des différences divisées pour le calcul des β i , 0 ≤ i ≤ n. Le polynome d’interpolation de Newton de degré n, Pn , évalué au point x0 donne n

Pn ( x0 ) =

∑ β 0 ωi ( x0 ) = β 0 = y0 =

f [ x0 ].

i =0

On note f [ x0 ] = y0 la différence divisée d’ordre 0, correspondante à β 0 . De même, on évalue le polynome d’interpolation de Newton, Pn , au point x1 . On obtient : n

Pn ( x1 ) =

∑ β i ωi ( x1 )

(3)

i =0

= β 0 + β 1 ( x1 − x0 ) = f [ x0 ] + β 1 ( x1 − x0 ) = f [ x1 ] = y1 d’où β1 =

(4) (5) (6)

f [ x1 ] − f [ x0 ] = f [ x0 , x1 ] x1 − x0

$ f[x_0, x_1]$ correspond à la différence divisée d’ordre 1. On procède par par récurrence pour obtenir : f [ x1 , ..., xk ] − f [ x0 , ...., xk−1 ] βk = = f [ x0 , ..., xk ] x k − x0 où f [ x0 , ..., xk ] désigne la différence divisée d’ordre k. Ecrire une fonction diff_div(x,y)qui renvoit les coefficients β i , i ∈ {0, · · · , n} en utilisant la méthode des différences divisées. [9]: def diff_div(x,y): n=len(y) beta=y.copy() for i in np.arange(1,n): 7

for j in np.arange(n-1,i-1,-1): beta[j]=(beta[j]-beta[j-1])/(x[j]-x[j-i]) return beta Ecrire une fonction Interpolation_Newton(t,x,y) qui évalue, au point t, le polynôme d’interpolation de Newton associé aux points d’interpolation ( xi , yi )0≤i≤n , avec x = ( xi )0≤i≤n et y = (yi )0≤i≤n . Représeneter ensuite le polynôme P graphiquement qui interpole les points (−1, 8), (0, 3) et (1, 6). [10]: def Interpolation_Newton(t,x,y): n=len(x) P=np.zeros((len(t))) beta=diff_div(x,y) for i in np.arange(0,n): P+=beta[i]*Newton(t,i,x) return P [11]: P=Interpolation_Newton(t,x,y) plt.figure(figsize=(20,10)) plt.plot(t,P,'mo--',linewidth=3,markersize=12) plt.xlabel('t',fontsize=30) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.grid(True) plt.legend(('Polynome d\'interpolation de Newton',),fontsize=20, loc = 0) plt.text(-1,8.5,"(-1,8)",ha="center",va="top",fontsize=30) plt.text(0,3.5,"(0,3)",ha="center",va="top",fontsize=30) plt.text(1,6.5,"(1,6)",ha="center",va="top",fontsize=30) [11]: Text(1, 6.5, '(1,6)')

8

3.3

Avantages du polynôme de Newton

Un des avantages de la méthode de Newton pour l’interpolation des points ( xi , yi )0≤i≤n , deux à deux distincts, est le suivant : si on note par Pk le polynôme d’interplation tronqué : le polynôme de degré inférieur ou égal à k, 0 ≤ k ≤ n, qui interpole que les points ( xi , yi )0≤i≤k , exprimé dans la base de polynômes de Newton {ω1 , · · · , ωk }, comme suit : Pk ( x ) = β 0 .|{z} 1 + β 1 ( x − x0 ) + β 2 ( x − x0 )( x − x1 ) + .... + β k ( x − x0 )( x − x1 )...( x − xk−1 ), | {z } | | {z } {z } ω0

ω1

ω2

x ∈ R,

ωk

alors Pk+1 , 0 ≤ k < n, le polynôme tronqué de degré inférieur ou égal à k + 1 interpolant les points ( xi , yi )0≤i≤k+1 , sera exprimé en fonction de Pk comme suit : Pk+1 ( x ) = Pk ( x ) + β k+1 ( x − x0 )( x − x1 )..( x − xk ). | {z } ω k +1

Par conséquent, si l’on connait le polynôme Pn , interpolant les points ( xi , yi )0≤i≤n , et que l’on rajoute un point d’interpolation ( xn+1 , yn+1 ), alors le polynôme Pn+1 interpolant les n + 2 points sera déduit de celui interpolant les anciens n + 1 points comme suit : Pn+1 ( x ) = Pn ( x ) + β n+1 ( x − x0 )( x − x1 )..( x − xn ). | {z } ω n +1

Les coefficients β 0 , ..., β n resteront les mêmes, il suffit de calculer le coefficient β n+1 et de déduire le polynôme de Newton associé ωn+1 = ωn ( x − xn ). Exercice 1. Ecrire une fonction Interpolation_Newton_opt(t,x,y) qui optimise le calcul de la détermination du polynôme d’interpolation par la méthode de Newton. 2. Représeneter ensuite le polynôme P graphiquement qui interpole les points (−1, 8), (0, 3) et (1, 6). 3. Comparer les deux fonctions Interpolation_Newton(t,x,y) et Interpolation_Newton_opt(t,x,y) en terme de temps d’exécution. [12]: def Newton_opt(t,i,x): n=len(x) if i==0: return np.ones((len(t))) elif i==1: return t-x[0] else: return (t-x[i-1])*Newton_opt(t,i-1,x) [13]: def Interpolation_Newton_opt(t,x,y,beta=diff_div(x,y)): n=len(y) if len(x)==2: return beta[0]+beta[1]*Newton_opt(t,1,x[0:1]) else: return Interpolation_Newton_opt(t,x[0:n-1],y[0: ,→n-1])+beta[len(x)-1]*Newton_opt(t,len(x)-1,x)

9

[14]: P=Interpolation_Newton_opt(t,x,y) plt.figure(figsize=(20,10)) plt.plot(t,P,'mo--',linewidth=3,markersize=12) plt.xlabel('t',fontsize=30) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.grid(True) plt.legend(('Polynome d\'interpolation de Newton',),fontsize=20, loc = 0) plt.text(-1,8.5,"(-1,8)",ha="center",va="top",fontsize=30) plt.text(0,3.5,"(0,3)",ha="center",va="top",fontsize=30) plt.text(1,6.5,"(1,6)",ha="center",va="top",fontsize=30) [14]: Text(1, 6.5, '(1,6)')

[15]: %timeit Interpolation_Newton(t,x,y) %timeit Interpolation_Newton_opt(t,x,y) %timeit Interpolation_Lagrange(t,x,y) 55.7 ţs ś 8.02 ţs per loop (mean ś std. dev. of 7 runs, 10000 loops each) 15.7 ţs ś 1.09 ţs per loop (mean ś std. dev. of 7 runs, 100000 loops each) 68.4 ţs ś 7.54 ţs per loop (mean ś std. dev. of 7 runs, 10000 loops each)

10

4 4.1

Erreur d’interpolation Théorème

Soient f une fonction de classe C n+1 ([ x0 , xn ]) et Pn son polynôme d’interpolation aux points x0 < x1 < · · · < xn . L’erreur d’approximation de f par Pn est donnée par : ∀ x ∈ [ x0 , xn ], max | f (n+1) (t) |

En ( x ) =| f ( x ) − Pn ( x ) |≤

5

t∈[ x0 ,xn ]

( n + 1) !

n

∏ | x − xi | . i =0

Phénomène de Runge

Le phénomène de Runge désigne le problème qui peut survenir lorsqu’on tente d’approcher une fonction à l’aide de polynômes d’interpolation. Contrairement à ce que l’on aurait pu penser, lorsque le degré du polynôme interpolateur augmente, les oscillations de ce dernier au bord du domaine s’intensifient. Il y a malgré tout une convergence du polynôme vers la fonction sur un intervalle réduit par rapport au domaine de définition initial.

6

Explication du phénomène de Runge sur un exemple

1 définie sur [−1, 1] et la suite de points de la subdivision 1 + 8x2 uniforme de l’intervalle [−1, 1] suivante : Considérons la fonction f ( x ) =

(n)

xk

2 = −1 + k , n

k ∈ {0, 1, ..., n}. (n)

(n)

Nous notons par Pn le polynôme d’interpolation de degré n satisfaisant Pn ( xk ) = f ( xk ) pour k ∈ {0, 1, ..., n}. Le but est d’étudier la convergence de Pn vers f lorsque l’on augmente n. On peut observer une divergence de l’interpolation pour cette fonction, aux bords de l’intervalle [−1, 1]. Ce résultat pourrait être expliqué par deux faits : 1. La dérivée d’ordre n + 1 de la fonction f croit plus rapidement que (n + 1)! lorsque n augmente. 2. Le choix des points d’interpolation équidistants qui n’est pas adéquat pour ce type de fonctions. Exercice Soit f ( x ) = 1+18x2 définie sur [−1, 1]: Rédiger un script produisant l’affichage, dans une même fenêtre graphique, de l’erreur En = | Pn ( x ) − f ( x )| de la fonction f et de ses polynômes d’interpolation Pn pour n ∈ {5, 10, 15, 20, 25, 30}. [16]: f=lambda x: 1/(1+8*x**2) N=np.arange(5,31,5) t= np.linspace(-1,1,1000) Erreur=np.zeros((len(N),len(t))) for i in np.arange(0,len(N)): x=np.linspace(-1,1,N[i]) Erreur[i,:]=np.abs(f(t)-Interpolation_Lagrange(t,x,f(x))) [17]: plt.figure(figsize=(20,10)) plt.subplot(3,2,1) plt.plot(t,Erreur[0,:],linewidth=3) 11

plt.title('Erreur pour n=5',fontsize=20) plt.grid(True) plt.subplot(3,2,2) plt.plot(t,Erreur[1,:],linewidth=3) plt.title('Erreur pour n=10',fontsize=20) plt.grid(True) plt.subplot(3,2,3) plt.plot(t,Erreur[2,:],linewidth=3) plt.title('Erreur pour n=15',fontsize=20) plt.grid(True) plt.subplot(3,2,4) plt.plot(t,Erreur[3,:],linewidth=3) plt.title('Erreur pour n=20',fontsize=20) plt.grid(True) plt.subplot(3,2,5) plt.plot(t,Erreur[4,:],linewidth=3) plt.title('Erreur pour n=25',fontsize=20) plt.grid(True) plt.xlabel('t',fontsize=30) plt.subplot(3,2,6) plt.plot(t,Erreur[5,:],linewidth=3) plt.title('Erreur pour n=30',fontsize=20) plt.grid(True) plt.xlabel('t',fontsize=30) plt.xticks(fontsize=20) plt.yticks(fontsize=20) [17]: (array([-0.5, 0. , 0.5, 1. , 1.5, 2. , )

12

2.5,

3. ]),

7

Application

En relevant toutes les 10 secondes, la vitesse d’écoulement de l’eau dans une conduite cylindrique, on a obtenu |temps (seconde)|0 | 10 |20|30| |:|:-|:—|:-|:-| |vitesse |2.00 | 1.89 |1.72 | 1.44| (a) Trouver une approximation de la vitesse pour t = 15 secondes via un polynôme interpolant de degré 2. (b) Répéter l’opération avec un polynôme de degré 3.

8

Références

[1] Kiusalaas, J. (2013). Numerical methods in engineering with Python 3. Cambridge university press. [2] Numpy Package [3] Mathplotlib Package [4] Jupyter markdowns

13