Interpolation (Cours Math 4 "Usthb")

Interpolation
Prof. Alfio Quarteroni
Semestre d’automne 2007-2008
epfl
Prof. Alfio Quarteroni Interpolation

Exemples et motivations
Example
On considère un test mécanique pour établir le lien entre
contrainte (MPa = 100N/cm2 ) et déformations relatives (cm/cm)
d’un échantillon de tissu biologique (disque intervertébral, selon P.
Komarek, Ch. 2 de Biomechanics of Clinical Aspects of
Biomedicine, 1993, J. Valenta ed., Elsevier).
disques intervertebraux
epfl

F
A
test i contrainte σ déformation ǫ
= F=A
1
2
0.00
0.06
0.00
0.08
= L=L
3 0.14 0.14
4 0.25 0.20 L L
5 0.31 0.23
6 0.47 0.25
7 0.60 0.28
8 0.70 0.29
A partir des ces données, on veut estimer la déformation

correspondante à un effort σ = 0.9 MPa.
epfl

A travers la méthode des moindres carrés on obtient que la droite
qui approche mieux ces données est p(x) = 0.3938x − 0.0629. On
peut utiliser cette droite (dite de régression) pour estimer ǫ lorsque
σ = 0.9 MPa: on trouve p(0.9) ≃ 0.4.
0.8
0.6
0.4 valeur estimée pour σ = 0.9

ǫ
0.2
−0.2
0 0.1 0.2 0.3

σ
0.4 0.5 0.6 0.7 0.8
epfl

Example
La population de la Suisse dans les années 1900 à 2000 est
recensée comme il suit (en milliers),
année 1900 1910 1920 1930 1941 1950

population 3315 3753 3880 4066 4266 4715
année 1960 1970 1980 1990 2000
population 5429 6270 6366 6874 7288
◮ Peut-on estimer le nombre d’habitants de la Suisse pendant

les années où il n’y a pas eu de recensement, par exemple en
1945 et en 1975?
◮ Peut-on envisager un modèle pour prédire la population de la epfl
Suisse en 2010?
Le polynôme de degré deux (parabole) qui approche ces données au
sens des moindres carrés est p(x) = 0.19x 2 − 710.29x + 657218.92.
Population de la Suisse au cours des années

9000
valeurs estimées par interpolation

8000
population en milliers
7000
6000
5000
4000
valeurs mesurées
3000
1900 1920 1940 1960 1980 2000
année
epfl

Example
Les points dans la figure suivante représentent les mesures du débit
du sang dans une section de l’artère carotide pendant un battement
cardiaque. La fréquence d’acquisition de données est constante et
égale à 10/T où T=1 sec. est la période du battement.
On veut déduire de ce signal discret un signal continu représenté
par une combinaison linéaire de fonctions connues (par exemple de
fonctions trigonométriques, cette représentation est adaptée a
notre problème parce que la fonction obtenue va bien approximer
un signal périodique).
40
35
30
25
20
15
10
0
epfl
−5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Soit f (t) le signal dont on connaı̂t N = 10 échantillons
[f (t0 ), . . . , f (tN−1 )], avec tj = jT /N. On cherche {ck } ∈ C,
k ∈ [0, N − 1] ⊂ N tels que:
N−1
1 X −kj
f (tj ) = c k ωN , j = 0, . . . , N − 1, (1)
N
k=0
où
ωN = e (−2πi)/N = cos(2π/N) − i sin(2π/N),
étant i l’unité imaginaire. On peut calculer le vecteur des
coefficients c = [c1 , . . . , c10 ]T par l’algorithme F.F.T. (Fast Fourier
Transform, Cooley et Tuckey, 1965), implémenté dans Octave avec
la commande fft.
epfl

A partir de l’expression (1), il y a plusieurs techniques pour
reconstruire la valeur de f dans chaque point t ∈ [0, T ]. Une fois
que l’on a reconstruit la fonction, on peut la représenter sur un
graphe comme ceci
40
35
30
25
20
15
10
−5
−10
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
epfl

Position du problème [Sec. 3.1 du livre]
Soit n ≥ 0 un nombre entier. Etant donnés n + 1 points distincts
x0 , x1 ,. . . xn et n + 1 valeurs y0 , y1 ,. . . yn , on cherche un
polynôme p de degré n, tel que
p(xj ) = yj pour 0 ≤ j ≤ n. (2)
Dans le cas affirmatif, on note p = Πn et on appelle Πn le

polynôme d’interpolation aux points xj , j = 0, . . . , n.
(n = 4)
6 6
5 5
4 4
Πn
3 3
y0 y0
yn
2 2
yn
1 1
0 0
−1
x0 xn −1
x0 xn epfl
−1 0 1 2 3 4 5 6 −1 0 1 2 3 4 5 6

Soit f ∈ C 0 (I ) et x0 , . . . , xn ∈ I . Si comme valeurs yj on prend
yj = f (xj ), 0 ≤ j ≤ n, alors le polynôme d’interpolation Πn (x)est
noté Πn f (x) et est appelé l’interpolant de f aux points x0 ,. . . xn .
6 6
(n = 4)
5 5
4 4
3
f 3
Πn f
2 2
1 1 f
0 0
x0 xn x0 xn
−1 −1
−1 0 1 2 3 4 5 6 −1 0 1 2 3 4 5 6
epfl

Base de Lagrange [Sec. 3.1.1 du livre]
On considère les polynômes ϕk , k = 0, . . . , n de degré n tels que
ϕk (xj ) = δjk , k, j = 0, . . . , n,
où δjk = 1 si j = k et δjk = 0 si j 6= k. Explicitement, on a
n
Y (x − xj )
ϕk (x) = .
(xk − xj )
j=0,j6=k
epfl

La figure qui suit montre trois polynômes de Lagrange de degré
n = 6 relatifs aux points d’interpolations x0 = −1,
x1 = −2/3,. . . ,x5 = 2/3, et x6 = 1.
1.5
(n = 6)
ϕ0 (x) ϕ3 (x)
1
0.5
x1 x4 x6
0
x0 x2 x3 x5
epfl
−0.5
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

Example
Pour n = 2, x0 = −1, x1 = 0, x2 = 1 les polynômes de la base de
Lagrange sont
(x − x1 )(x − x2 ) 1
ϕ0 (x) = = x(x − 1),
(x0 − x1 )(x0 − x2 ) 2
(x − x0 )(x − x2 )
ϕ1 (x) = = −(x + 1)(x − 1),
(x1 − x0 )(x1 − x2 )
(x − x0 )(x − x1 ) 1
ϕ2 (x) = = x(x + 1).
(x2 − x0 )(x2 − x1 ) 2
1.2
1
ϕ0 ϕ1 ϕ2
0.8
0.6
0.4
0.2
0
epfl
−0.2
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

Le polynôme d’interpolation Πn des valeurs yj aux points xj ,
j = 0, . . . , n, s’écrit
n
X
Πn (x) = yk ϕk (x), (3)
k=0
Pn
car il vérifie Πn (xj ) = k=0 yk ϕk (xj ) = yj .
Par conséquent on aura
n
X
Πn f (x) = f (xk )ϕk (x).
k=0
epfl

On montre maintenant que le polynôme Πn en (3) est le seul
polynôme de degré n interpolant les données yi aux nœuds xi .
Soit, en effet, Qn (x) un autre polynôme d’interpolation. Alors on a
Qn (xj ) − Πn (xj ) = 0, j = 0, . . . , n.
Donc, Qn (x) − Πn (x) est un polynôme de degré n qui s’annule en

n + 1 points distincts; il en suit que Qn = Πn , d’où l’unicité du
polynôme interpolant.
epfl

Interpolation d’une fonction regulière
Théorème (Erreur d’interpolation) (voir prop. 3.2 du livre)

Soient x0 , x1 , . . ., xn , n + 1 nœuds équirépartis dans I = [a, b] et
soit f ∈ C n+1 (I ). Alors, on a
n+1
1 b−a
En (f ) = max |f (x) − Πn f (x)| ≤ max |f (n+1) (x)|.
x∈I 4(n + 1) n x∈I
(4)
On remarque que l’erreur d’interpolation dépend de la dérivée

n + 1-ième de f .
epfl

Example
Polynômes d’interpolation Πi f pour i = 1, 2, 3, 6 et
x +1
f (x) = sin(x), noeuds équirépartis sur [0, 6].
5
1
Π3 f
0.5 f Π f
6
Π2 f
0
Π1 f
−0.5
−1
−1.5
−2
−1 0 1 2 3 4 5 6 7
epfl

Interpolation numérique avec Octave
En Octave on peut calculer les polynômes d’interpolation en

utilisant les commandes polyfit et polyval. Voyons plus en
détaille comment peut-on utiliser ces commandes.
p = polyfit(x,y,n) calcule les coefficients du polynôme de

degré n qui interpole les valeurs y aux points x.
epfl

Exemple A. On veut interpoler les valeurs
y = [3.38, 3.86, 3.85, 3.59, 3.49] aux points x = [0, 0.25, 0.5, 0.75, 1]
par un polynôme de degré 4. Il suffit d’utiliser les commandes
Octave suivantes:
>> x=[0:0.25:1]; % vecteur des points d’interpolation

>> y=[3.38 3.86 3.85 3.59 3.49]; % vecteur des valeurs
>> p1=polyfit(x,y,4)
p1 =
1.8133 -0.1600 -4.5933 3.0500 3.3800
p1 est le vecteur des coefficients du polynôme interpolant:

Π4 (x) = 1.8133x 4 − 0.16x 3 − 4.5933x 2 + 3.05x + 3.38.
Pour calculer le polynôme de degré n qui interpole une fonction f

donnée dans n + 1 points, il faut d’abord construire le vecteur y en
évaluant f dans les nœuds x. epfl

Exemple B. Considérons le cas suivant:
>> f=’cos(x)’;
>> x=[0:0.25:1];
>> y=eval(f);
>> p=polyfit(x,y,4)
p =
0.0362 0.0063 -0.5025 0.0003 1.0000
Remarque
Si la dimension m + 1 de x et y est m + 1 > n + 1 (n étant le
degré du polynôme d’interpolation), la commande
polyfit(x,y,n) retourne le polynôme interpolant de degré n au
sens des moindres carrés (voir sec. 3.4 du livre). Dans le cas où
m + 1 = n + 1, alors on trouve le polynôme d’interpolation
“standard”, puisque dans ce cas les deux coı̈ncident. epfl

y = polyval(p,x) calcule les valeurs y d’un polynôme de degré
n, dont les n+1 coefficients sont memorisés dans le vecteur p, aux
points x, c’est-à-dire: y = p(1)*x.^n + . . . + p(n)*x + p(n+1).
Exemple C. On veut évaluer le polynôme trouvé dans l’Exemple A

dans le point x = 0.4 et, ensuite, on veut tracer son graphe. On
peut utiliser les commandes suivantes:
>> x=0.4;
>> y=polyval(p1,x)
y =
3.9012
>> x=linspace(0,1,100);
>> y=polyval(p1,x);
>> plot(x,y)
epfl

Example
On veut interpoler la fonction sin(x) + x sur n = 2, 3, 4, 5, 6
nœuds. En Octave on peut utiliser la commande polyfit pour
calculer les coefficients du polynôme interpolant et polyval pour
évaluer un polynôme dont on connaı̂t les coefficients dans une
suite de points. Voilà les commandes Octave :
>> f = ’sin(x) + x’;

>> x=[0:3*pi/100:3*pi]; x_sample=x;
>> plot(x,eval(f),’b’); hold on
>> for i=2:6
x=linspace(0,3*pi,i); y=eval(f);
c=polyfit(x,y,i-1);
plot(x_sample,polyval(c,x_sample),’b--’)
end
epfl

La figure suivante montre les 5 polynômes obtenus.
12
10
8
Π4 f (x) Π1 f (x)
f (x)
4
0
Π2 f (x)
Π5 f (x)
−2
0 1 2 3 4 5 6 7 8 9 epfl

Remarque
Seul le fait que
n+1
1 b−a
lim =0
n→∞ 4(n + 1) n
n’implique pas que En (f ) tend vers zéro quand n → ∞.
Example
1
(Runge) Soit f (x) = 1+x 2 , x ∈ [−5, 5]. Si on l’interpole dans des
points équirépartis, au voisinage des extrémités de l’intervalle,
l’interpolant présente des oscillations, comme on peut le voir sur la
figure.
epfl

Fonction de Runge et oscillations des polynômes interpolants dans
des points équirépartis.
2
1.5 Π f (x)
10
1
Π f (x)
5
0.5
0
f(x)
−0.5 epfl
−5 −2.5 0 2.5 5

Remède : Interpolation par intervalles.
epfl

Interpolation par intervalles [Sec. 3.2 du livre]
Soient x0 = a < x1 < · · · < xN = b des points qui divisent

l’intervalle I = [a, b] dans une réunion d’intervalles Ii = [xi , xi+1 ]
de longueur H où
replacements
b−a
H= .
N
H
a xi−1 xi xi+1 b
Sur chaque sous-intervalle Ii on interpole f|Ii par un polynôme de
degré 1. Le polynôme par morceaux qu’on obtient est noté
ΠH
1 f (x), et on a:
f (xi+1 ) − f (xi )
ΠH
1 f (x) = f (xi ) + (x − xi ) pour x ∈ Ii .
xi+1 − xi
epfl

Exemple 7 (suite)
On considère les polynômes par morceaux de degré n = 1
interpolants la fonction de Runge, pour 5 et 10 sous-intervalles de
[−5, 5].
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−5 −4 −3 −2 −1 0 1 2 3 4 5
La figure montre les polynômes ΠH1

1 f et ΠH2
1 f pour H1 = 2.5 et
epfl
H2 = 1.0.

Les commandes Octave utilisées sont:
>> f = inline(’1./(1+x.^2)’,’x’);
>> H1 = 2.5; H2=1.0;
>> x1 = [-5:H1:5]; x2 = [-5:H2:5];
>> y1 = feval(f, x1); y2 = feval(f, x2);
>> x_plot = [-5:.1:5];
>> y1_plot = interp1(x1,y1,x_plot);
>> y2_plot = interp1(x2,y2,x_plot);
>> fplot(f, [-5 5]); hold on;
>> plot(x_plot, y1_plot); plot(x_plot, y2_plot);
epfl

Théorème 1 (Prop. 3.3 du livre)
Si f ∈ C 2 (I ), (I = [x0 , xN ]) alors telle que
H2
E1H (f ) = max | f (x) − ΠH
1 f (x) |≤ max |f ′′ (x)|.
x∈I 8 x∈I
Démonstration. D’après la formule (4), sur chaque intervalle Ii on a
H2
max | f (x) − ΠH
1 f (x) |≤ max | f ′′ (x) | .
x∈[xi ,xi+1 ] 4(1 + 1) x∈Ii
Remarque
On peut montrer que si l’on utilise un polynôme de degré n (≥ 1)
dans chaque sous-intervalle Ii on trouve
H n+1
EnH (f ) ≤ max |f (n+1) (x)| .
4(n + 1) x∈I
epfl

Exemple 7 (suite)
On considère la fonction de Runge f (x) sur [−5, 5] et on prend un
nombre K croissant de sous-intervalles K = 20, 40, 80, 160 et on
estime l’erreur d’interpolation commise en évaluant la différence
|f (x) − ΠH
1 f (x)| sur une grille très fine :
>> f=inline(’1./(1+x.^2)’,’x’);
>> K=[20 40 80 160]; H=10./K;
>> x_fine = [-5:0.001:5];
>> f_fine = feval(f,x_fine);
>> for i=1:4
x = [-5:H(i):5]; y = feval(f, x);
y_fine = interp1(x,y,x_fine);
err1(i) = max(abs(f_fine - y_fine));
end
>> loglog(H,err1);
epfl

La figure qui suit montre (en échelle logarithmique) l’erreur
commise en fonction de la taille H = 10/K des sous-intervalles.
On voit que l’erreur E1H f pour l’interpolation linéaire par morceaux
se comporte comme CH 2 : ce résultat est en accord avec le
théorème 1. En plus, si on calcule les rapports E1H /H 2 on peut
estimer les constantes C :
>> err1./H.^2
ans =
0.16734 0.22465 0.24330 0.24829
epfl

Erreur d’interpolation de la fonction de Runge par le polynôme
composite ΠH 1 en fonction de H.
−1
10
EH
1
−2
10
−3
10
−4
10 −2 −1 0
10 10 10
H
epfl

Approx. au sens des moindres carrés [Sec. 3.4 du livre]
Supposons que l’on dispose de n + 1 points x0 , x1 , . . . , xn et n + 1

valeurs y0 , y1 , . . . , yn . On a vu que, si le nombre de données est
grand, le polynôme interpolant peut présenter des oscillations
importantes.
Pour avoir une meilleure représentation des données, on peut
chercher un polynôme de degré m < n qui approche “au mieux”
les données.
epfl

Definition
On appelle polynôme aux moindres carrés de degré m f˜m (x) le
polynôme de degré m tel que
n
X n
|yi − f˜m (xi )|2 ≤
X
|yi − pm (xi )|2 ∀pm (x) ∈ Pm
i=0 i=0
Remarque
Lorsque yi = f (xi ) (f étant une fonction continue) alors f˜m est dit
l’approximation de f au sens des moindres carrés.
epfl

Autrement dit, le polynôme aux moindres carrés est le polynôme
de degré m qui, parmi tous les polynômes de degré m, minimise la
distance des données.
Si on note f˜m (x) = a0 + a1 x + a2 x 2 + . . . + am x m , et on définit la
fonction
n
yi − a0 + a1 xi + a2 xi2 + . . . + am xim 2
X
Φ(a0 , a1 , . . . , am ) =
i=0
alors les coefficients du polynôme aux moindres carrés peuvent être

déterminés par les relations
∂Φ
= 0, 0≤k ≤m (5)
∂ak
ce qui nous donne m + 1 relations linéaires entre les ak .
epfl

Utilisons cette méthode dans un cas simple. Considérons les points
x0 = 1, x1 = 3, x2 = 4 et les valeurs y0 = 0, y1 = 2, y2 = 7 et calculons
le polynôme interpolant de degré 1 au sens des moindres carrés (ligne de
régression).
Le polynôme recherché a la forme f˜1 (x) = a0 + a1 x. On définit:
P2 ∂Φ ∂Φ
Φ(a0 , a1 ) = i=0 [yi − (a0 + a1 xi )]2 et on impose ∂a 0
= 0 et ∂a 1
= 0:
2 2 2
!
∂Φ X X X
= −2 [yi − (a0 + a1 xi )] = −2 yi − 3a0 − a1 xi
∂a0
i=0 i=0 i=0
= −2(9 − 3a0 − 8a1 )
2 2 2 2
!
∂Φ X X X X
= −2 xi [yi − (a0 + a1 xi )] = −2 xi yi − a0 xi − a1 xi2
∂a1
i=0 i=0 i=0 i=0
= −2(34 − 8a0 − 26a1 )
Donc les coefficients a0 et a1 du polynôme sont les solutions du système

linéaire:
3a0 + 8a1 = 9
8a0 + 26a1 = 34 epfl

En général, on observe que si pour calculer le polynôme interpolant aux
moindres carrés f˜m (x) = a0 + a1 x + a2 x 2 + . . . + am x m on impose les
conditions d’interpolation f˜m (xi ) = yi pour i = 0, . . . , n, alors on trouve
le système linéaire Ba = ỹ, où B est la matrice de dimension
(n + 1) × (m + 1)
1 x1 . . . x1m
 
 1 x2 . . . x2m 
B= . .. 
 
 .. . 
1 xn . . . xnm
Puisque m < n le système est surdéterminé, c’est-à-dire que le nombre de

lignes est plus grand que le nombre de colonnes. Donc on ne peut pas
résoudre ce système de façon classique, mais on doit le résoudre au sens
des moindres carrés, en considérant:
B T Ba = B T ỹ.
De cette façon on trouve le système linéaire Aa = y (avec A = B T B et

y = B T ỹ), dit système d’équations normales. On peut montrer que les epfl
équations normales sont équivalentes au système (5).

Exemple 1 (suite)
On reprend l’exemple 1 proposé au début du chapitre: on va
présenter l’interpolation des données avec Octave. D’abord, il faut
définir les données expérimentales (9 couples de valeurs σ-ǫ):
>> sigma = [0.00 0.06 0.14 0.25 0.31 0.47 0.50 0.70];
>> epsilon = [0.00 0.08 0.14 0.20 0.22 0.26 0.27 0.29];
On veut extrapoler la valeur de ǫ pour σ = 0.4. On considère le

polynôme Π7 et le polynôme aux moindres carrés de degré 1 (la
ligne de régression).
On peut utiliser les commandes suivantes:
>> plot(sigma, epsilon, ’*’); hold on; % valeurs connues
>> sigma_sample = linspace(0,1.0,100);
>> p7 = polyfit(sigma, epsilon, 7);
>> pol = polyval(p7, sigma_sample); % polynome interpolant
>> plot(sigma_sample,pol, ’r’);
>> p1 = polyfit(sigma,epsilon,1);
>> pol_mc = polyval(p1, sigma_sample); % moindres carres
epfl
>> plot(sigma_sample, pol_mc, ’g’); hold off;

La figure montre le polynôme interpolant Π8 (en rouge) et le
polynôme d aux moindres carrés de degré 1 (ligne de régression, en
bleu).
8
4
ε
−1
0 0.2 0.4 0.6 0.8 1
σ
On remarque que pour σ > 0.7 MPa les comportements des

fonctions qui approchent les données sont tout à fait différents. epfl

En particulier, pour σ = 0.9 les valeurs de ǫ(σ) extrapolées avec le
polynôme interpolant Π7 et le polynôme d aux moindres carrés de
degré 1, sont obtenus par:
>> polyval(p7, 0.9)
ans =
1.7221
>> polyval(p1, 0.9)

ans =
0.4173
◮ La déformation obtenue par Π7 (172.21%) est sûrement trop

grande pour être considérée physiologique.
◮ Par contre, la valeur obtenue grâce à la ligne de régression
peut être utilisée pour estimer la déformation pour σ = 0.9
Mpa.
epfl

On revient à l’exemple sur la population suisse.
Exemple 2 (suite) On veut estimer la population de la Suisse en

2010 en utilisant les valeurs connues en 1900-2000. A ce propos,
on utilise un polynôme au moindres carrés p2 (x) de degré 2.
>> annee = [1900, 1910, 1920, 1930, 1941, 1950,...
1960, 1970, 1980, 1990, 2000];
>> population = [3315, 3753, 3880, 4066, 4266, 4715,...
5429, 6270, 6366, 6874, 7288];
>> p2 = polyfit(annee, population, 2);
>> annee_sample = [1900:2:2010];
>> vp2 = polyval(p2, annee_sample);
>> plot(annee, population, ’*r’, ...
annee_sample, vp2, ’b’);
epfl

9000
8000
7000
population
6000
5000
4000
3000
1900 1920 1940 1960 1980 2000
année
La valeur éstimées par la méthode pour l’année 2010 est
>> polyval(p2, 2010)

ans = 8084 epfl

Finalement, on montre comment il est possible d’interpoler des
données périodiques à travers la fonction interpft.
Exemple 3 (suite) On définit les 10 couples de données (temps,

débit) par les vecteurs t et deb; puis, à l’aide de la commande
interpft, on interpole les données dans 1000 nœuds équirépartis
dans l’intervalle de périodicité [0,1]:
>> t = [0 .1 .2 .3 .4 .5 .6 .7 .8 .9];
>> deb = [0 35 .15 5 0 5 .6 .3 .15 0];
>> f = interpft(deb, 1000);
>> plot(linspace(0, 1, 1000), f); hold on;
>> plot(t, deb, ’r*’)
Le résultat est montré par la figure qui suit. epfl

40
35
30
25
20
15
10
−5
−10
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
epfl

Octave: traitement des polynômes
Dans Octave il y a des commandes spécifiques pour faire des
calcules avec polynômes. Soit x un vecteur de abscisses, y un
vecteur d’ordonnées et p (respectivement pi ) le vecteur des
coefficients d’un polynôme P(x) (respectivement Pi ); alors on a les
commandes suivantes:
commande action
y=polyval(p,x) y = valeurs de P(x)
p=polyfit(x,y,n) p = coefficients du polynôme d’interpolation Πn
z=roots(p) z = zéros de P tels que P(z) = 0
p=conv(p1 ,p2 ) p = coefficients du polynôme P1 P2
[q,r]=deconv(p1 ,p2 ) q = coefficients de Q, r = coefficients de R
tels que P1 = QP2 + R
y=polyderiv(p) y = coefficients de RP ′ (x)
y=polyinteg(p) y = coefficients de P(x) dx
epfl

Interpolation par fonctions splines [Sec. 3.3 du livre]
Soient a = x0 < x1 < · · · < xn = b des points qui divisent

l’intervalle I = [a, b] dans une réunion d’intervalles Ii = [xi , xi+1 ].
Definition
On appelle spline cubique interpolant f une fonction s3 qui
satisfait
1. s3 |I ∈ P3 pour tout i = 0, . . . n − 1, P3 étant l’ensemble de
i
polynômes de degré 3,
2. s3 (xi ) = f (xi ) pour tout i = 0, . . . n,
3. s3 ∈ C 2 ([a, b]).
epfl

Cela revient à vérifier les conditions suivantes (on indique par
s3 (xi− ) la limite à gauche de s3 au point xi et par s3 (xi+ ) la limite à
droite) :
s3 (xi− ) = f (xi ) pour tout 1 ≤ i ≤ n − 1,

s3 (xi+ ) = f (xi ) pour tout 1 ≤ i ≤ n − 1,
s3 (x0 ) = f (x0 ),
s3 (xn ) = f (xn ),
s3′ (xi− ) = s3′ (xi+ ) pour tout 1 ≤ i ≤ n − 1,
s3′′ (xi− ) = s3′′ (xi+ ) pour tout 1 ≤ i ≤ n − 1,
c’est-à-dire 2(n − 1) + 2 + 2(n − 1) = 4n − 2 conditions.
epfl

Il faut trouver 4n inconnues (qui sont les 4 coefficients de chacune
des n restrictions s3 |I , i = 0, . . . n − 1) et on dispose de 4n − 2
i
relations.
On rajoute alors 2 conditions supplémentaires à vérifier.
◮ Si l’on impose
s3′′ (x0+ ) = 0 et s3′′ (xn− ) = 0, (6)
alors la spline s3 est complètement déterminée et s’appelle

spline naturelle.
◮ Une autre possibilité est d’imposer la continuité des dérivées
troisièmes dans le nœuds x2 et xn−1 , c’est à dire:
s3′′′ (x2− ) = s3′′′ (x2+ ) et −

s3′′′ (xn−1 +
) = s3′′′ (xn−1 ). (7)
Les conditions (7) sont dites not-a-knot.

epfl

Spline cubique naturelle interpolant la fonction
f (x) = (x−0.3)12 +0.01 + (x−0.9)12 +0.04 − 6, (nœuds équirépartis
xj = −1 + j/4, j = 0, . . . , 16)
100
90
f
80
70
60
50
40
30
s3
20
10
−10
−1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5
epfl

Interpolation (Cours Math 4 "Usthb")

Transféré par

Droits d'auteur :

Formats disponibles

Interpolation (Cours Math 4 "Usthb")

Transféré par

Informations du document

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Interpolation (Cours Math 4 "Usthb")

Transféré par

Droits d'auteur :

Formats disponibles

Interpolation

Prof. Alfio Quarteroni

Semestre d’automne 2007-2008

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

A partir des ces données, on veut estimer la déformation

Prof. Alfio Quarteroni Interpolation

0.4 valeur estimée pour σ = 0.9

0 0.1 0.2 0.3

Prof. Alfio Quarteroni Interpolation

année 1900 1910 1920 1930 1941 1950

◮ Peut-on estimer le nombre d’habitants de la Suisse pendant

Population de la Suisse au cours des années

valeurs estimées par interpolation

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

p(xj ) = yj pour 0 ≤ j ≤ n. (2)

Dans le cas affirmatif, on note p = Πn et on appelle Πn le

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

On considère les polynômes ϕk , k = 0, . . . , n de degré n tels que

où δjk = 1 si j = k et δjk = 0 si j 6= k. Explicitement, on a

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

Par conséquent on aura

Prof. Alfio Quarteroni Interpolation

Soit, en effet, Qn (x) un autre polynôme d’interpolation. Alors on a

Donc, Qn (x) − Πn (x) est un polynôme de degré n qui s’annule en

Prof. Alfio Quarteroni Interpolation

Théorème (Erreur d’interpolation) (voir prop. 3.2 du livre)

On remarque que l’erreur d’interpolation dépend de la dérivée

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

En Octave on peut calculer les polynômes d’interpolation en

p = polyfit(x,y,n) calcule les coefficients du polynôme de

Prof. Alfio Quarteroni Interpolation

>> x=[0:0.25:1]; % vecteur des points d’interpolation

p1 est le vecteur des coefficients du polynôme interpolant:

Pour calculer le polynôme de degré n qui interpole une fonction f

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

Exemple C. On veut évaluer le polynôme trouvé dans l’Exemple A

Prof. Alfio Quarteroni Interpolation

>> f = ’sin(x) + x’;

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

n’implique pas que En (f ) tend vers zéro quand n → ∞.

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

Soient x0 = a < x1 < · · · < xN = b des points qui divisent

Prof. Alfio Quarteroni Interpolation

La figure montre les polynômes ΠH1

Prof. Alfio Quarteroni Interpolation

Prof. Alfio Quarteroni Interpolation

Démonstration. D’après la formule (4), sur chaque intervalle Ii on a