Chapitre 3
Interpolation et approximation
Les objectifs de ce cours sont les suivants :
1. l’acquisition d’algorithmes numériques indispensable à la résolution de problèmes usuels
et la maîtrise du traitement des données expérimentales selon la Méthode Numériques.
2. Acquérir la maîtrise du logiciel choisissons ici Pytyon pour illustrer une méthode numé-
rique à l’aide d’un programme écrit.
3. Développer des « mini-projets », qui sera évalué plus tard.
Introduction
ω Dans la pratique, une fonction f est souvent connue uniquement par ses valeurs en certains
points x0 , . . . , xn ,
ω à partir de la connaissance de f en ces points, nous souhaitons "reconstruire le mieux
possible" cette fonction i.e. trouver une "bonne" approximation de f (pour toutes les autres
valeurs de x dans un intervalle fixé).
ω On pourra imposer à l’approximation de passer exactement par les points donnés. Dans ce
cas, on parle d’interpolation.
ω Dans ce cours nous chercherons principalement à approximer une fonction par des fonctions
polynômes.
ω L’espace des polynômes possède de "bonnes" propriétés : structure d’espace vectoriel, bases
adaptées, calculs analytiques simples pour la d’erivation et l’intégration, peu coûteux à évaluer.
3.1 Position du problème
On se place ici dans un cadre réel, en dimension un. On note Pn l’espace vectoriel des
polynômes de degré inférieur ou égal à n. On sait dans le cours d’algébre que sa dimension est
n + 1. On va définir une fonction f connues sur ses n + 1 points distincts
x0 < · · · < xn , et n + 1 à valeurs associées y0 , . . . , yn ,
xi x0 x1 · · · xn→1 xn
.
yi y0 y1 · · · yy→1 yn
on considère le problème d’interpolation suivant :
Trouver Pn → Pn tel que Pn (xi ) = yi , i = 0, . . . , n. (PI)
Les abscisses xi sont appelées nœuds d’interpolation. Les données yi pourront correspondre aux
valeurs prises par une fonction f : R ↑ R. Dans ce cas, on dit que Pn interpole la fonction f
aux nœuds xi .
23
3.2 Calcul du polynôme d’interpolation par di!érentes méthodes
Nous allons, dans cette première partie, calculer l’unique polynôme d’interpolation Pn so-
lution du problème (PI), par les méthodes de Vandermonde, de Lagrange, puis de Newton.
L’objectif est de comparer ces di!érentes approches du point de vue numérique.
Chaque méthode revient à calculer Pn dans une certaine base (Qk )0,...,n de Pn . Le polynôme
Pn , décomposé dans cette base, s’écrit
n
!
Pn (x) = εk Qk (x),
k=0
où les εk sont des coe"cients réels. Le calcul du polynôme d’interpolation consiste alors, une
fois les Qk donnés, à déterminer les coe"cients εk , pour k = 0, . . . , n.
3.2.1 Base canonique : système de Vandermonde
Dans la base canonique de Pn , le polynôme Pn s’écrit
n
!
Pn (x) = εk xk .
k=0
Pour que Pn soit solution de (PI), les coe"cients εk doivent vérifier
n
!
εk xki = yi , i = 0, . . . , n.
k=0
Autrement dit, les εk doivent être calculés comme solution du système linéaire
1 x10 . . . xn0 a0 y0
1 x 1 . . . x1 a 1 y 1
1 n
.. .. . . . . = .
. . . .. .. ..
1 x1n . . . xnn an yn
La matrice de ce système est appelée matrice de Vandermonde. Si les xi sont distincts deux à
deux, son déterminant est non nul, et le système admet une solution unique.
Question 1. Ecrire une fonction V=matrice_vander(xi) qui prend en argument un vecteur xi
contenant les nœuds xi , et qui renvoie la matrice de Vandermonde V associée. Afin de gagner en
temps de calcul sous Python, on essaiera de minimiser le nombre de boucles for.
Question 2. Ecrire une fonction y=eval_poly(ak,x) qui prend en argument un vecteur ( de coef-
ficients ak et un vecteur d’abscisses x, et qui renvoie les valeurs du polynôme P (x) = nk=0 ak xk
aux points de x.
Question 3. Ecrire une fonction y=interp_vander(xi,yi,x) qui calcule, aux abscisses de x,
les valeurs du polynôme Pn interpolant les (xi , yi ).
Question 4. Tester les fonctions Python écrites sur le jeu de données
xi ↓3 ↓1 0 2 4
.
yi 2 0 2 ↓1 1
A!cher la matrice de Vandermonde obtenue, et tracer sur un même graphique le polynôme
d’interpolation Pn et les points (xi , yi ).
Le principal inconvénient de cette méthode de calcul est son coût algorithmique : la résolution
du système linéaire nécessite O(n3 ) opérations, ce qui devient vite intolérable pour les grandes
valeurs de n.
24
3.2.2 Base de Lagrange
Les n + 1 polynômes qui constituent la base de Lagrange sont définis par
)n * x↓x + (x ↓ x0 )(x ↓ x1 ) · · · (x ↓ xk→1 )(x ↓ xk+1 ) · · · (x ↓ xn )
j
Lk (x) = =
xk ↓ xj (xk ↓ x0 )(xk ↓ x1 ) · · · (xk ↓ xk→1 )(xk ↓ xk+1 ) · · · (xk ↓ xn )
j=0
j↑=k
k=0,...,n
où les xj sont les nœuds d’interpolation du problème. On sait que, dans cette nouvelle base, les
coe"cients du polynôme d’interpolation Pn sont directement donnés par les valeurs à interpoler
yk , k = 0, . . . , n. Autrement dit, le polynôme Pn s’écrit, dans la base de Lagrange,
n
!
Pn (x) = y0 L0 (x) + · · · + yn Ln (x) = yk Lk (x).
k=0
La di"culté n’est plus, comme avec la base canonique, de déterminer les coe"cients εk , mais
plutôt de construire les polynômes Lk , qui dépendent des nœuds d’interpolation considérés.
Question 5. Ecrire une fonction y=poly_lagrange(xi,k,x) qui prend en argument un vecteur
xi contenant les nœuds d’interpolation, un entier k et un vecteur d’abscisses x, et qui renvoie
les valeurs du polynôme de Lagrange Lk aux points de x.
Question 6. A!cher, sur un même graphique, les di"érents polynômes Lk construits à partir
du jeu de données de la question 4.
Question 7. Ecrire une fonction y=interp_lagrange(xi,yi,x) qui calcule les valeurs du po-
lynôme d’interpolation Pn associé aux (xi , yi ), aux points de x.
Question 8. Tracer, sur un même graphique, le polynôme Pn et les points (xi , yi ). Comparer
au résultat obtenu par la méthode de Vandermonde.
3.2.3 Base de Newton
Motivation :
En rajoutant un point xj au noeud, le calcul du polynôme de Lagrange necessite le calcul de
tous toute la base de Lagrange Lk (x). Newton a trouver un reméde à ce problème en construisant
une nouvelle base et obtenir le au final le même polynôme.
Les n + 1 polynômes de Newton ϑk associés au problème (PI) sont définis par
ϑ0 (x) = 1,
k→1
)
ϑ k (x) = (x ↓ x 0 ) · · · (x ↓ x k→1 ) = (x ↓ xj ), k = 1, . . . , n.
j=0
Supposons que les valeurs à interpoler yi soient les images d’une certaine fonction f , c’est-à-dire
yi = f (xi ), i = 0, . . . , n. Le polynôme d’interpolation de f s’écrit, dans la base de Newton,
n
!
Pn (x) = f [x0 , . . . , xk ]ϑk (x) =,
k=0
où les f [x0 , . . . , xk ] sont les di"érences divisées de Newton, calculées par la formule de récurrence
suivante :
f (xi ) si k = i,
f [xi , . . . , xk ] = f [x , . . . , x ] ↓ f [x , . . . , x ]
i+1 k i k→1
si k ↔ i.
xk ↓ xi
Remarquons que :
25
— f [x0 , . . . , xn ] est le coe"cient dominant dans la forme canonique
Pn (x) = f [x0 , . . . , xn ]xn + · · ·
— f [x0 , . . . , xn ] est invariant par des permutations sur les xi : en e!et, permuter les xi ne
change pas le polynôme d’interpolation et donc ne change pas le coe"cient du terme de
plus haut degré.
On peut représenter cette formule par un tableau, dans lequel le terme (i, j) est calculé à
partir des éléments (i ↓ 1, j ↓ 1) et (i, j ↓ 1) :
x0 f [x0 ]
x1 f [x1 ] f [x0 , x1 ]
x2 f [x2 ] f [x1 , x2 ] f [x0 , x1 , x2 ]
.. .. .. .. ..
. . . . .
xn f [xn ] f [xn→1 , xn ] f [xn→2 , xn→1 , xn ] ... f [x0 , x1 , . . . , xn→1 , xn ]
on peut créer le programme sous Python suivant :
//Fonction diff [Link]
function m = diff div(x,y)
n = length(x) ;
m = zeros(n,n) ;
m( :,1) = y;
for j = 2 :n
for i = j :n
m(i,j) = (m(i,j-1)-m(i-1,j-1))/(x(i)-x(i-j+1)) ; end
end
Question 9. Ecrire une fonction Python T=diff_div(xi,yi), qui prend en argument les points
(xi , yi ), et qui renvoie, sous la forme d’une matrice T, le tableau des di"érences divisées ci-dessus.
Question 10. Ecrire une fonction y=interp_newton(xi,yi,x) qui renvoie les valeurs, aux
abscisses de x, du polynôme Pn interpolant les (xi , yi ).
Question 11. Comme précédemment, tester les fonctions Python écrites sur le jeu de données
de la question 4.
Question 12. On souhaite ajouter le point (0.5, 2) au jeu de données précédent. Comment
obtenir le nouveau polynôme d’interpolation sans recalculer toutes les di"érences divisées ? Écrire
un programme Scilab pour le faire.
3.2.4 Quelques propriétés de l’interpolation polynomiale
Nous avons vu, dans la partie précédente, comment calculer le polynôme d’interpolation de
trois façons di!érentes. Nous allons maintenant étudier certaines propriétés de l’interpolation
polynomiale, qui sont indépendantes de la méthode utilisée pour le calcul du polynôme interpo-
lateur. On pourra donc utiliser au choix l’une des trois méthodes programmées précédemment.
26
3.2.5 Stabilité du polynôme d’interpolation
On considère la fonction f (x) = sin(2ϖx) sur l’intervalle [↓1, 1], ainsi qu’une subdivision
avec 22 nœuds équirépartis.
Question 13. A l’aide de la commande Scilab rand, générer des données perturbées yi des
valeurs yi = f (xi ), avec la condition
yi ↓ yi | ↗ 9.5 · 10→4 .
max |
i=0,...,21
Question 14. Tracer, sur un même graphique, les polynômes P21 et P21 qui interpolent respec-
tivement les points (xi , yi ) et (xi , yi ). Evaluer numériquement
sup P21 (x) ↓ P21 (x) ,
x↓[→1,1]
et comparer cette valeur à l’amplitude de la perturbation imposée.
3.2.6 Majoration de l’erreur d’interpolation.
Reste à estimer l’écart entre une fonction et son polynôme interpolateur,
f (x) = PL + erreur, PL = Pn
de plus on a le :
Théorème 6. Soit f une fonction n+1 fois dérivable sur un intervalle I = [a, b] de R, x0 , ..., xn
des réels distincts de I. Soit Pn le polynome de Lagrange donné par les xj et yj = f (xj ). Pour
tout réel x → I, il existe un réel ϱx → [a, b] (qui dépend de x) tel que :
n
f [n+1] (ϱx ) )
En (x) = f (x) ↓ Pn (x) = (x ↓ xj ) (3.1)
(n + 1)!
j=0
avec
a < min(x, x0 ) < ϱ < max(x, xn ) < b
Ainsi l’erreur commise dépend d’une majoration de la taille de la dérivée n + 1-ième sur
l’intervalle, mais aussi de la disposition des points xj par rapport à x. Par exemple
1. Si f (x) est un polyôme de p ↘ n alors f (n+1) (ϱx ) = 0 et l’erreur est nulle.
2. L’erreur diminue en prenant plus de points
3. Elle sera d’autant plus petite que x est proche d’un point d’interpolation f ne varie pas
brusquement.
Preuve du théorème : Si x est l’un des xj l’égalité est vraie. Soit
n
)
C = (f (x) ↓ P (x))/ (x ↓ xj )
j=0
on considère maintenant la fonction :
n
)
g(t) = f (t) ↓ P (t) ↓ C (t ↓ xj )
j=0
27
elle s’annule en xj pour j variant de 0 à n ainsi qu’en x suite au choix de la constante C, donc g
s’annule au moins n + 2 fois sur l’intervalle contenant les xj et x, donc g ↔ s’annule au moins n + 1
fois sur ce même intervalle, donc g ↔↔ s’annule au moins n fois, etc. et finalement g [n+1] s’annule
une fois au moins sur cet intervalle. Or
g [n+1] = f [n+1] ↓ C(n + 1)!
car P est de degré inférieur ou égal à n et nj=0 (x ↓ xj ) ↓ xn+1 est de degré inférieur ou égal à
n. Donc il existe bien un réel ϱx dans l’intervalle contenant les xj et x tel que
f [n+1] (ϱx )
C=
(n + 1)!
Du théorème on a :
(x ↓ x0 )(x ↓ x1 ) · · · (x ↓ xn )
f (x) ↓ PL (x) ↘ max f [n+1] (x)
(n + 1)! x↓[a,b]
Comme (x ↓ xi ) ↘ (b ↓ a), ≃ i = 0, · · · , n alors
n+1
(x ↓ x 0 )(x ↓ x 1 ) · · · (x ↓ x n ↘ (b ↓ a)
)
et (b ↓ a)n+1
f (x) ↓ PL (x) ↘ max f [n+1] (x)
(n + 1)! x↓[a,b]
Lemme 1. (Lien entre di"érences divisées et dérivées)
Si f est de classe C (n) [a, b], alors pour tout noeud x0 < · · · < xn dans [a, b], il existe ϱ → [a, b]
tel que
f (n) (ϱ)
f [x0 , · · · , xn ] =
n!
On a alors l’erreur par le méthode de Newton s’écrit pour f → C n+1 sur [a, b] :
n n
f (n+1) (ϱ) ) )
En (x) = (x ↓ xj ) = f [x0 , · · · , xn , x] (x ↓ xj )
(n + 1)!
j=0 j=0
n
)
f (x) ↓ Pn (x) = En (x) = f [x0 , · · · , xn , x] (x ↓ xj )
j=0
Attention, l’erreur d’interpolation peut devenir très grande aussi lorsqu’on utilise beaucoup
de points d’interpolation.
Phénomène de Runge
On cherche à interpoler la fonction
1
f (x) =
1 + x2
sur l’intervalle [↓5, 5].
Question 15. On choisit une répartition uniforme des xi . Sur un même graphique, tracer la
fonction f et les polynômes d’interpolation P5 , P10 et P15 . Que remarque-t-on ?
28