Chapitre II

Interpolation et Approximation

Probleme de linterpolation : on recherche des fonctions simples (polynomes, polynomes par

morceaux, polynomes trigonometriques) passant par (ou proche) des points donnes

(x0 , y0 ), (x1 , y1 ), . . . , (xn , yn ), (0.1)

c.-a-d., on cherche p(x) avec p(xi ) = yi pour i = 0, 1, . . . , n. Si les valeurs de yi satisfont yi =

f (xi ) ou f (x) est une fonction donnee, il est alors interessant detudier lerreur de lapproximation

f (x) p(x) = ? (0.2)

Bibliographie de ce chapitre
J.H. Ahlberg, E.N. Nilson & J.L. Walsh (1967): The Theory of Splines and Their Applications.
Academic Press, New York. [MA 65/4]
C. de Boor (1978): A Practical Guide to Splines. Springer-Verlag. [MA 65/141]
G.D. Knott (2000): Interpolating Cubic Splines. Birkhauser. [MA 65/431]
H.J. Nussbaumer (1981): Fast Fourier Transform and Convolution Algorithms. Springer-Verlag.
H. Spath (1995): One Dimensional Spline Interpolation. AK Peters. [MA 65/362]

II.1 Differences divisees et formule de Newton

. . . tho I will not undertake to prove it to others.
(Newton, letter to Collins, Nov. 8, 1676 ; publ. Cotes 1711, p. 38)

Probleme (Newton 1676). Etant donnes les n + 1 points (0.1), chercher un polynome

p(x) = axn + bxn1 + cxn2 + . . . (1.1)

de degre n qui satisfasse

p(xi ) = yi pour i = 0, 1, . . . , n. (1.2)

Pour un exemple voir la fig. II.1.

Interpolation et Approximation 25

10
p(x)
5

0
2 4 6 8 10

F IG . II.1: Polynome dinterpolation de degre 5

Solution. En inserant les conditions (1.2) dans (1.1), le probleme se transforme en un systeme
lineaire (a matrice du type Vandermonde ; ici ecrit pour n = 2)

c + bx0 + ax20 = y0 y1 y0
b + a(x1 + x0 ) =
x1 x0
c + bx1 + ax21 = y1 soustraire et diviser (1.3)
y2 y1
c + bx2 + ax22 = y2 b + a(x2 + x1 ) =
x2 x1
et, si on soustrait et divise une deuxieme fois, on trouve
1  y2 y1 y1 y0 
a= . (1.4)
x2 x0 x2 x1 x1 x0
Le meme calcul a ete effectue pour n = 4 dans un manuscript de Newton datant de 1676 ; comme a
laccoutumee, Newton refusa de le publier (voir citation). Cotes le publia comme dernier chapitre
Methodus differentialis du livre Analysis per quantitatum series, fluxiones, ac differentias, Londini
1711 (voir fac-simile en figure II.2 1 ).

F IG . II.2: Fac-simile du calcul de Newton pour le probleme de linterpolation

Dans tous ces calculs apparaissent les differences divisees :

1
On peut observer que Newton matrise les eliminations de variables dans un systeme lineaire avec brio ; plus tard,
toute la gloire pour cette methode reviendra a Gauss.
26 Interpolation et Approximation

Definition 1.1 (differences divisees) Pour (xi , yi ) donnes (xi distincts) on definit
y[xi ] := yi
y[xj ] y[xi ]
y[xi , xj ] :=
xj xi
y[xj , xk ] y[xi , xj ]
2 y[xi , xj , xk ] :=
xk xi
y[xj , xk , xl ] 2 y[xi , xj , xk ]
2
3 y[xi , xj , xk , xl ] := etc.
xl xi
Ensuite Newton, charmant comme toujours, donne la formule suivante sans la moindre demonstra-
tion (voir citation) :

Theoreme 1.2 (formule de Newton) Le polynome dinterpolation de degre n qui passe par les
n + 1 points (x0 , y0 ), (x1 , y1), . . . , (xn , yn ), ou les xi sont distincts, est unique et donne par

p(x) = y[x0 ] + (x x0 ) y[x0, x1 ] + (x x0 )(x x1 ) 2 y[x0 , x1 , x2 ]

(1.5)
+ . . . + (x x0 )(x x1 ) . . . (x xn1 ) n y[x0 , x1 , . . . , xn ].

Demonstration. Nous utilisons deux idees :

1. On procede par recurrence. Pour n = 1, et en tenant compte des premiers deux points, nous
avons
y1 y0
p(x) = y0 + (x x0 ) . (1.6)
x1 x0
Il sagit dune formule bien connue des Geometres (voir , figure II.1.8).
Puis, pour n = 2, en rajoutant le point (x2 , y2 ), on essaie de batir la dessus un polynome de
degre 2, qui ne change plus les valeurs de y0 et de y1 . Il est donc de la forme

y1 y0
p(x) = y0 + (x x0 ) + a (x x0 )(x x1 ) (1.7)
x1 x0
0
2 4

ou le coefficient a est a determiner. Mais il sagit du coefficient de x2 de p(x) : nous savons deja
(voir (1.4)) que celui-ci est la deuxieme difference divisee 2 y[x0 , x1 , x2 ].
Pour demontrer le cas general, nous supposons que

p1 (x) = y[x0 ] + (x x0 ) y[x0, x1 ] + . . . + (x x0 ) . . . (x xn2 ) n1y[x0 , x1 , . . . , xn1 ]

soit le polynome unique de degre n 1 qui passe par (xi , yi ) pour i = 0, 1, . . . , n 1. Alors,
comme auparavant, le polynome p(x) a necessairement la forme

p(x) = p1 (x) + a (x x0 )(x x1 ) . . . (x xn1 ),

ou a est determine par p(xn ) = yn .

2. Lidee de Aitken-Neville. Pour montrer que a = n y[x0 , x1 , . . . , xn ], ce qui acheve la demonstra-
tion, nous considerons egalement le polynome de degre n 1

p2 (x) = y[x1 ] + (x x1 ) y[x1 , x2 ] + . . . + (x x1 ) . . . (x xn1 ) n1 y[x1 , x2 , . . . , xn ],

Interpolation et Approximation 27

p1 (x) p(x)
10

0
2 4 6 8 10
5

10 p2 (x)
F IG . II.3: Les polynomes p1 (t), p2 (t) et p(t) de lalgorithme dAitken-Neville

qui passe par (xi , yi) pour i = 1, . . . , n (voir figure II.3). Ensuite, on pose (Aitken - Neville, 1929,
1932 2 )
1  
p(x) = (xn x)p1 (x) + (x x0 )p2 (x) . (1.8)
xn x0
Il sagit dun polynome de degre n, qui satisfait la condition (1.2) pour le point x0 (ici, le facteur
(x x0 ) est nul), pour le point xn (ici, le facteur (x xn ) est nul), et pour les points x1 , . . . , xn1
(ici, les deux polynomes p1 et p2 sont egaux a yi ). Le polynome desire est donc trouve.
En considerant le coefficient de xn dans (1.8), nous obtenons
1  n1 
a= y[x1 , . . . , xn ] n1 y[x0 , . . . , xn1 ] = n y[x0 , . . . , xn ],
xn x0
ce qui demontre la formule (1.5).

TAB . II.1: Differences divisees pour les donnees de la fig. II.1

xi yi y 2y 3y 4y 5y
0 1
1
2 1 3/8
5/2 77/120
4 6 17/6 167/960
6 3/4 287/9600
5 0 5/3 1/8
2/3 1/4
8 2 1/6
3/2
10 5

Exemple 1.3 Pour les donnees de la fig. II.1, les differences divisees sont presentees dans le
tableau II.1. Le polynome dinterpolation est alors donne par
3 77 167
p(x) = 1 + x + x(x 2) x(x 2)(x 4) + x(x 2)(x 4)(x 5)
8 120 960
287
x(x 2)(x 4)(x 5)(x 8) .
9600
ou mieux encore pour la programmation (ou le calcul a la main)
 3  77  167 287 
p(x) = 1 + x 1 + (x 2) + (x 4) + (x 5) (x 8) .
8 120 960 9600
2
Il fallait plus de deux siecles pour avoir cette idee !...
28 Interpolation et Approximation

Remarque. Lordre des {xi } na aucune importance pour la formule de Newton (1.5). Si lon
permute les donnees (xi , yi ), on obtient evidemment le meme polynome. Pour lexemple ci-dessus
et pour les {xi } choisis dans lordre {4, 5, 2, 8, 0, 10}, on obtient ainsi
  17 3  167 287 
p(x) = 6 + (x 4) 6 + (x 5) + (x 2) + (x 8) x .
6 4 960 9600

En observant que n y[xi0 , . . . , xin ] est une fonction symetrique de ses arguments (par exemple,
2 y[x2 , x3 , x1 ] = 2 y[x1 , x2 , x3 ], voir exercices), on peut utiliser les valeurs calculees dans le
tableau II.1.
Si x est entre 4 et 5, les deux facteurs x4 et x5 dans la formule precedente sont relativement
petits, ce qui favorise la diminution des erreurs darrondi.

II.2 Erreur de linterpolation

Supposons que les points (xi , yi ) soient sur le graphe dune fonction f : [a, b] IR, c.-a-d.,

yi = f (xi ), i = 0, 1, . . . , n, (2.1)

etudions alors lerreur f (x) p(x) du polynome dinterpolation p(x). Deux exemples sont donnes
dans la fig. II.4. A gauche, on voit un polynome dinterpolation pour la fonction f (x) = sin x, et
a droite pour la fonction 1/(1 + x2 ). Pour mieux rendre visible lerreur, on a dessine la fonction
f (x) en une courbe pointillee.

1
1 f (x) = sin x f (x) =
1 1 + x2

0
0 2 4 6 8
0
4 2 0 2 4
1

F IG . II.4: Polynome dinterpolation pour sin x (gauche) et pour 1/(1 + x2 ) (droite)

Les resultats suivants sont dus a Cauchy (1840, Sur les fonctions interpolaires, C.R. XI, p. 775-789,
Oeuvres ser. 1, vol. V, p. 409-424). Commencons par une relation interessante entre les differences
divisees pour (2.1) et les derivees de la fonction f (x).

Lemme 2.1 Soit f (x) n-fois differentiable et yi = f (xi ) pour i = 0, 1, . . . , n (xi distincts). Alors,
il existe un (min xi , max xi ) tel que

f (n) ()
n y[x0 , x1 , . . . , xn ] = . (2.2)
n!
Demonstration. Soit p(x) le polynome dinterpolation de degre n passant par (xi , yi ) et notons
d(x) = f (x) p(x). Par definition de p(x), la difference d(x) sannule en n + 1 points distincts :

d(xi ) = 0 pour i = 0, 1, . . . , n.
Interpolation et Approximation 29

Comme d(x) est differentiable, on peut appliquer n fois le theoreme de Rolle (voir le cours
dAnalyse I) et on en deduit que
d (x) a n zeros distincts dans (min xi , max xi ).
Le meme argument applique a d (x) donne
d (x) a n 1 zeros distincts dans (min xi , max xi ),
i i

et finalement encore
d(n) (x) a 1 zero dans (min xi , max xi ).
i i
(n)
Notons ce zero de d (x) par . Alors, on a
f (n) () = p(n) () = n! n y[x0 , x1 , . . . , xn ]. (2.3)
La deuxieme identite dans (2.3) resulte du fait que n y[x0 , x1 , . . . , xn ] est le coefficient de xn dans
p(x).

Theoreme 2.2 Soit f : [a, b] IR (n + 1)-fois differentiable et soit p(x) le polynome dinterpo-
lation de degre n qui passe par (xi , f (xi )) pour i = 0, 1, . . . , n. Alors, pour x [a, b], il existe un
(min(xi , x), max(xi , x)) tel que
f (n+1) ()
f (x) p(x) = (x x0 ) . . . (x xn ) . (2.4)
(n + 1)!
Demonstration. Si x = xi pour un indice i {0, 1, . . . , n}, la formule (2.4) est verifiee car
p(xi ) = f (xi ). Fixons alors un x dans [a, b] qui soit different de xi et montrons la formule (2.4)
pour x = x.
Lidee est de considerer le polynome p(x) de degre n + 1 qui passe par (xi , f (xi )) pour i =
0, 1, . . . , n et par (x, f (x)). La formule de Newton donne
p(x) = p(x) + (x x0 ) . . . (x xn ) n+1 y[x0 , . . . , xn , x]. (2.5)
Si lon remplace la difference divisee dans (2.5) par f (n+1) ()/(n + 1)! (voir le lemme precedent)
et si lon pose x = x, on obtient le resultat (2.4) pour x = x car p(x) = f (x). Comme x est
arbitraire, la formule (2.4) est verifiee pour tout x.

Exemple 2.3 Dans la situation de la fig. II.4, on a n + 1 = 7. Comme la 7eme derivee de sin x est
bornee par 1, on a que
1
|p(x) sin x| |x(x 1.5)(x 3)(x 4.5)(x 6)(x 7.5)(x 9)| ,
7!
par exemple
|p(4) sin 4| 0.035 ou |p(1) sin 1| 0.181.
Pour le deuxieme exemple, f (x) = 1/(1 + x2 ), la 7eme derivee est donnee par
(x + 1)(x 1)x(x2 2x 1)(x2 + 2x 1)
f (7) (x) = 8! ,
(1 + x2 )8
qui est maximale pour x 0.17632698. On obtient ainsi
 1 
 4392
p(x)  |(x2 20.25)(x2 9)(x2 2.25)x| .
1 + x2 7!
Alors, lerreur peut etre 4392 fois plus grande que pour linterpolation de sin x.
30 Interpolation et Approximation

Convergence de linterpolation.
Une grande surprise en mathematiques fut la decouverte, dabord par Riemann (1854), puis
par Weierstrass (1872), de lincroyable complexite quont certaines fonctions continues,
p. ex. , de netre nulle part differentiables ;
puis la deuxieme grande surprise : toutes ces fonctions, aussi compliquees quelles puissent
etre, peuvent etre approchees, aussi pres quon le veut et uniformement, par les fonctions les
plus simples qui existent, des polynomes (Weierstrass 1885 ; voir [HW96], III.9) ;
personne ne pensait alors que les polynomes dinterpolation, si on prend seulement les points
suffisamment proches les uns des autres, ne convergeaient pas vers la fonction donnee. La
decouverte que cela nest meme pas assure pour les fonctions rationnelles, les deuxiemes
fonctions les plus simples, (voir dessin de figure II.5), a choque enormement les mathemati-
ciens vers 1900 (en particulier E. Borel).
Carl David Tolme Runge (1856-1927), premier prof de maths appliquees de lhistoire et, en tant
queleve de Weierstrass, ayant aussi une solide formation en maths pures, fut certes lhomme
ideal pour expliquer ce phenomene de maniere claire et elegante (1901, Zeitschr. Math. u. Physik
vol. 46).

1 0 1 1 0 1 1 0 1 1 0 1 1 0 1
n=2 n=4 n=6 n=8 n = 10

1 0 1 1 0 1 1 0 1 1 0 1 1 0 1
n = 12 n = 14 n = 16 n = 18 n = 20

F IG . II.5: Le phenomene de Runge pour f (x) = 1/(1 + 25x2 )

II.3 Polynomes de Chebyshev

La formule (2.4) montre que lerreur de linterpolation est un produit de la (n + 1)eme derivee de
f (x), evaluee a un point inconnu, avec lexpression (x x0 ) . . . (x xn ) qui ne depend que de
la division {x0 , . . . , xn }. Nous arrivons a la question suivante :
Chercher, pour un n donne, la division de [a, b] pour laquelle
max |(x x0 ) . . . (x xn )| est minimal. (3.1)
x[a,b]

Nous considerons lintervalle [1, 1] et avons le probleme :

Probleme. Chercher un polynome (avec coefficient principal = 1)
(x) = xn + an1 xn1 + . . . + a0 tel que L = max | (x)| est minimal. (3.2)
x[1,1]
Interpolation et Approximation 31

xi equidis. xi points Cheb.

0 0
0 1 0 1

F IG . II.6: Le produit (x x0 ) (x x1 ) . . . (x xn ) pour n = 10 et xi equidistants (gauche), xi poits

de Chebyshev (droite).

On trouve la reponse a cette question dans un travail de P.L. Chebyshev (transcription francaise :
Tchebychef, 1854, uvres I, p. 109) sur la conception optimale des tiges des pistons de loco-
motives a vapeur . . . . . . Entre-temps, les locomotives a vapeur sont au musee, et les polynomes de
Chebyshev sont encore et toujours des outils importants en mathematiques.

1 1 1 1
a = .00 L = 1.00 a = .60 L = .40 a = .75 L = .25 a = .90 L = .33

1 0 1 1 0 1 1 0 1 1 0 1

1 1 1 1

F IG . II.7: Valeurs maximales de (x) = x3 ax

Le cas n = 3. Cherchons a resoudre cette question pour n = 3, i.e., posons (symetrie oblige !...)
(x) = x3 ax ou a est un parametre a trouver. Nous voyons en figure II.7 que, pour a = 0, nous
avons L = 1 ; puis cette valeur diminue quand a crot, jusquau moment ou la courbe touche la
borne +L et L en meme temps (pour a = 3/4) ; apres, L recommence de grandir. La solution
optimale est donc
3
3 (x) = x3 x . (3.3)
4
Chebyshev a vu que pour un n quelconque, le polynone optimal possede, comme on dit au-
jourdhui, une alternante de Chebyshev :

Theoreme 3.1 Le polynome (3.2), minimisant L sur [1, 1], prend alternativement les valeurs +L
et L exactement n + 1 fois.

+L +L +L +L

1 0 1 1 0 1

L L
32 Interpolation et Approximation

Demonstration. La preuve est simple : supposons, par exemple pour n = 3, que le polynome
3 (x) = x3 + . . . prenne seulement trois fois les valeurs maximales, disons, +L, puis L, puis +L
(voir dessin precedent). Il existe alors un polynome q2 (x) de degre 2, qui est > 0, < 0, et > 0 a
ces trois points. Le polynome 3 (x) q2 (x) va donc, pour > 0, diminuer en valeur absolue a
tous les trois points. Par consequent, le polynome netait pas optimal.

Par contre, pour n quelconque, il est beaucoup plus difficile de trouver des formules explicites pour
ces polynomes. Lidee suivante est due a Zolotarev :
Idee. Multiplions 3 (x) de (3.3) par 4 et posons

T3 (x) = 4x3 3x ce qui ressemble a cos(3) = 4 cos3 3 cos (3.4)

(cf. lequation pour la trisection de langle ; , p. 79, ou [HW97], p. 7). On voit bien
en figure II.8, que cela represente la projection dune guirlande sur un tambour x = cos , y =
sin , z = cos(3) sur le plan (x, z). Il est maintenant facile detendre cette idee au cas general :

F IG . II.8: Les guirlandes z = cos(n) autour dun tambour x = cos , y = sin et leurs projections, les
polynomes de Chebyshev, sur le plan (x, z).

T4 (x)
T3 (x)
T1 (x)

1 0 1

T2 (x)

F IG . II.9: Les premiers 4 (a gauche) respectivement 30 (a droite) polynomes de Chebyshev

Definition 3.2 (Polynomes de Chebyshev) Pour n = 0, 1, 2, . . . et pour x [1, 1], on definit

Tn (x) = cos(n) ou x = cos . (3.5)

Interpolation et Approximation 33

Par les formules de Moivre (cf. [HW97], p. 45), on peut voir que, malgre cette definition etrange,
Tn (x) est un polynome en x. Mais il y a mieux : la formule
cos((n + 1)) + cos((n 1)) = 2 cos cos(n),
en posant cos = x, devient
T0 (x) = 1, T1 (x) = x, Tn+1 (x) = 2xTn (x) Tn1 (x). (3.6)
Par consequent, Tn (x) est un polynome de degre n dont le coefficient de xn est 2n1 , c.-a-d.,
Tn (x) = 2n1 xn + . . .. Les premiers sont T0 (x) = 1, T1 (x) = x,
T2 (x) = x2 1, T3 (x) = 4x3 3x, T4 (x) = 8x4 8x2 +1, T5 (x) = 16x5 20x3 +5x . (3.7)
Pour des dessins voir la figure II.93 . Les polynomes de Chebyshev sannulent en
  (2k + 1) 
Tn cos =0 pour k = 0, 1, . . . , n 1. (3.8)
2n

Revenons maintenant a la question de trouver une division satisfaisant (3.1).

Car le coefficient principal de Tn+1 (x) est 2n , nous voyons que
max |(x x0 ) . . . (x xn )| est minimal
x[1,1]

si et seulement si (x x0 ) . . . (x xn ) = 2n Tn+1 (x), c.-a-d., par (3.8), si

 (2k + 1) 
xk = cos , k = 0, 1, . . . , n (3.9)
2n + 2
(points de Chebyshev). Pour repondre a la question (3.1), il faut encore utiliser la translation
x 7 a+b
2
+ ba
2
x, qui envoie lintervalle [1, 1] sur [a, b]. On obtient alors

1 1

0 0
4 2 0 2 4 4 2 0 2 4

1 1
F IG . II.10: Interpolation avec des points equidistants (a gauche) et les points de Chebyshev (a droite)

Theoreme 3.3 Lexpression (3.1) est minimale parmi toutes les divisions {x0 , x1 , . . . , xn } si et
seulement si
a+b ba  (2k + 1) 
xk = + cos , k = 0, 1, . . . , n. (3.10)
2 2 2n + 2
Exemple 3.4 Comme dans la fig. II.4, nous considerons la fonction f (x) = 1/(1+x2 ) sur lintervalle
[4.5, 4.5]. Dans la fig. II.10, on compare le polynome dinterpolation base sur des points equidistants
avec celui base sur les points de Chebyshev. Tout commentaire est superflu !...
3
Pour une etude des courbes blanches dans la fig. II.9 (a droite) voir la page 209 du livre: Th.J. Rivlin, Chebyshev
Polynomials. 2nd ed., John Wiley & Sons, 1990 [MA 41/36]
34 Interpolation et Approximation

II.4 Influence des erreurs darrondi sur linterpolation

Chaque annee, les ordinateurs doublent leur performance et, aujourdhui, en un clin doeil, ils
executent des milliards doperations arithmetiques. La precision des calculs na toutefois pas
augmente ! Par consequent, en un clin doeil, un ordinateur fait aussi des milliards derreurs
darrondi !... Il est alors primordial de savoir si, apres toutes ces erreurs darrondi, les resultats
obtenus ont encore la moindre fiabilite. Nous appelons ce sujet letude de la stabilite numerique.
Encore une mauvaise surprise !... Equipe du merveilleux Theoreme de Runge, choisissons la
fonction f (x) = sin x sur lintervalle [0, 5]. Cette fonction na aucun pole fini, donc la convergence
du polynome dinterpolation est assuree pour tout x ! Et alors, que se passe-t-il en figure II.11 ?
Pour expliquer ce phenomene, interessons nous aux erreurs darrondi.

1 2 3 4 1 2 3 4 1 2 3 4

n = 30 n = 40 n = 50

F IG . II.11: Polynome dinterpolation pour f (x) = sin x sur [0, 5] a points equidistants

Representation en virgule flottante. Depuis 30 a 40 ans, la representation dun nombre reel

sur ordinateur a ete standarisee sur 32 bits binaires (64 en double precision). Prenons un lap-top
dernier cri, achete en 2005, et donnons lui les nombres 13 , 32 , 34 , 38 a digerer. Le resultat est donne
en Fig. II.12. Rappelons, quen base 2 la division 1 : 11 donne 0.0101010101 . . ., la division
10 : 11 = 0.101010101 . . ., la division 100 : 11 = 1.0101010101 . . ., etc. Nous constatons que
le premier chiffre non-nul subit une translation vers la virgule (virgule flottante); lexposant de 2
correspondant est stocke dans les bits 39. Le bit 2 indique les puissances positives, le bit 1 le
signe du nombre. Puis, dans les bits 1032 suit la mantisse, de facon a ce que le premier bit, 1, est
supprime pour gagner une place. Enfin, nous voyons que la suite reguliere des bits 1, 0, 1, 0, 1, 0...
est interrompue au bit 32, a cause dun arrondissement (correcte) des bits 33,34,35,... qui seraient
1, 0, 1, ... .
Estimation de lerreur darrondi. Pour un exposant p donne, le nombre correspondant x est
|x| 2p 1 . 0 0 . . .
9 10 11

tandis que lerreur darrondi est

|err| 2p 0 . 0 . . . 0 0 1 1 1 . . . = 2p 0 . 0 . . . 0 1
9 10 32 33 34 35 36 9 10 32 33

Estimation importante :
|err| |x| 224 = |x| eps. (4.1)
En double precision, nous disposons de 2 32 = 64 bits, dont 12 sont utilises pour lexposant. En
resume
REAL 4, eps = 224 5.96 108
REAL 8, eps = 253 1.11 1016
REAL 16, eps = 2113 9.63 1035 .
Interpolation et Approximation 35

.3333333 0 0 1 1 1 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

.6666667 0 0 1 1 1 1 1 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

1.3333334 0 0 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

2.6666667 0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

5.3333335 0 1 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

10.6666670 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

F IG . II.12: Les nombres 13 , 23 , 34 , . . . en virgule flottante en REAL*4

Theoreme 4.1 Si on denote le nombre arrondi dun x 6= 0 par arr (x), alors

|arr (x) x|
eps , (4.2)
|x|
ou
arr (x) = x(1 + ) ou || eps . (4.3)
Ces formules sont la base de toute etude derreurs darrondi.

Influence des erreurs dans yi sur le polynome dinterpolation. Supposons que les donnees yi ,
pour lesquelles on devrait calculer le polynome dinterpolation p(x), soient erronees,

ybi = yi (1 + i ) ou |i | eps, (4.4)

et quon calcule en fait avec ces valeurs un polynome p(x).

b
La difference de ces polynomes est
n
X
p(x)
b p(x) = i yii (x) ,
i=0

ou les i (x) sont les polynomes de Lagrange (voir Chap. I, formule (1.4) et Fig. II.13). On obtient
n
X
|p(x)
b p(x)| eps max |yi | |i (x)|. (4.5)
i=0,...,n
i=0
Pn
La fonction i=0 |i(x)| decrit lamplification de lerreur dans les donnees (voir figure II.14).

TAB . II.2: Constantes de Lebesgue pour noeuds equidistants

n 10 20 30 40 50 60 80 100
n 3.0 101 1.1 104 6.6 106 4.7 109 3.6 1012 3.0 1015 2.2 1021 1.8 1027

Sa valeur maximale n
X
n := max |i (x)| (4.6)
x[a,b]
i=0
36 Interpolation et Approximation

6 6 6 6 6 6

5 5 5 5 5 5

4 4 4 4 4 4

3 3 3 3 3 3

2 2 2 2 2 2

1 1 1 1 1 1

0 0 0 0 0 0
0 0 0 0 0 0
1 1 1 1 1 1

2 2 2 2 2 2

3 3 3 3 3 3

4 4 4 4 4 4

5 5 5 5 5 5

6 0 (x) 6 1 (x) 6 2 (x) 6 3 (x) 6 4 (x) 6 5 (x)

16 16 16 16 16 16 16

12 12 12 12 12 12 12

8 8 8 8 8 8 8

4 4 4 4 4 4 4

0 0 0 0 0 0 0
0 0 0 0 0 0 0

4 4 4 4 4 4 4

8 8 8 8 8 8 8

12 12 12 12 12 12 12

16 16 16 16 16 16 16
0 (x) 1 (x) 2 (x) 3 (x) 4 (x) 5 (x) 6 (x)

F IG . II.13: Polynomes de Lagrange a points equidistants pour n = 10 et n = 12

6 6
5 5
4 4
3
7 3
2 2
1 1

0 0
0 0
1 1
2 2
3 3
7
4 4
5 5
6 6

F IG . II.14: Constante de Lebesgue pour n = 7

sappelle la constante de Lebesgue 4 associee aux points x0 , x1 , . . . , xn et a lintervalle [a, b]. Cette
4
Pour une raison obscure ; on ne trouve dans les oeuvres de Lebesgue quun article general sur lapproximation
Interpolation et Approximation 37

constante peut etre calculee numeriquement (voir Table II.2).

Calcul asymptotique pour points equidistants. Nous voyons, en figure II.13, que les i dont les
xi se trouvent au centre de lintervalle sont dangereux vers le bord de lintervalle. Essayons donc
de trouver une valeur asymptotique pour le premier maximum de m (x) ou n = 2m :

x (x1) (x2) (x(m1)) (x(m+1)) (x(n2))(x(n1))(xn)

m (x) =
m(m1) (m2) 1 1 (m2) (m1) m
(4.7)
Nous avons
1 1 1 1 1 1 
m (x) = m (x) + + + ...+ + + ...+ .
x x1 x2 xm+1 xm1 xn

Donc m (x) = 0 si

1 1 1 1 1 1
= + + ...+ + + ...+ .
x 1x 2x m1x m+1x nx
1
Car la serie harmonique diverge, on voit que x
et x 0. Ainsi, asymptotiquement,

1 1 1 1 1 1
+ + ...+ + + . . . + log n .
x 1 2 m1 m+1 n

Cela insere dans (4.7) donne

1 1 2 3 (m1)m(m+1) n n 23/2
max |m (x)| 2 (4.8)
log n m(m!)2 log n n3/2

ou la celebre formule asymptotique de Stirling (voir [HW97], II.10)

nn
n! 2n
en

a ete utilisee. Le facteur 2n dans (4.8) montre clairement que chaque augmentation de n par 1 fait
perdre 1 bit de precision !... Ainsi, avec les quelque 30 bits de mantisse en REAL*4, la limite est
atteinte vers n = 30 (comme on lobserve en figure II.11).

1 1 1 1 1 1

0 0 0 0 0 0
0 0 0 0 0 0
1 0 (x) 1 1 (x) 1 2 (x) 1 3 (x) 1 4 (x) 1 5 (x)
F IG . II.15: Polynomes de Lagrange a points de Chebyshev pour n = 10

Points de Chebyshev. Est-il encore necessaire de mentionner que, pour les points de Chebyshev
xi = a+b
2
+ ba
2
cos( (2i+1)
2n+2
), on na pas cette horrible instabilite. Voir Fig. II.15.

(Oeuvres 3, p. 256), ou une note en bas de page discute le phenomene de Runge, qui fit fureur a lepoque. Le meme
phenomene a ete independamment decouvert par E. Borel (1905, Lecons sur les fonctions de variables reelles et les
developpements en series de polynomes, Chap. IV, p. 74-82), MA 26/49.
38 Interpolation et Approximation

II.5 Transformation de Fourier discrete

Commencons par un rappel de la Transformation CONTINUE de Fourier.
Origines :
Probleme de la corde vibrante (Taylor 1713, Joh. Bernoulli 1727, Euler 1739) ;
Theorie du son (Lagrange 1759) ;
Theorie calculi sinuum dEuler (1754, Opera 14, p. 542), Euler 1777 ;
Theorie de la chaleur (Fourier 1807, 1822).
Une serie trigonometrique est une serie du type
 
a0 X
y= + ak cos kx + bk sin kx (5.1)
2 k=1

(voir Figure II.16). Elle represente une fonction f (x) periodique de periode 2

f (x + 2k) = f (x) k ZZ . (5.2)

Une telle fonction est determinee par ses valeurs dans un intervalle de reference, soit [0, 2], soit
[, ]. Les sommes partielles sappellent des polynomes trigonometriques.

0.15 + 0.4 cos x + 0.6 sin x . . . 0.5 cos 2x 0.4 sin 2x . . . + 0.25 cos 3x 0.5 sin 3x
1 1 1

0 2 0 2 0 2

1 1 1

F IG . II.16: Premiers termes dune serie trigonometrique

Un des grand problemes des mathematiques, elucide seulement apres un lutte dun bon siecle, etait
la question si une quelconque fonction puisse etre decomposee ainsi en frequences basses et
frequences hautes (harmoniques ; voir Cours dAbalyse II).
Representation complexe. Les formules les plus simples sobtiennent en passant aux complexes.
Grace a
eix + eix eix eix
eix = cos x i sin x , cos x = , sin x = (5.3)
2 2i
la serie (5.1) devient

X
ikx
ck = 12 (ak ibk ), ak = ck + ck
y= ck e ou (k 0) ou (5.4)
k= ck = 21 (ak + ibk ) bk = i(ck ck ).

Orthogonalite. La cle fondamentale permettant le calcul des series trigonometriques a ete decouverte
par Euler (1777, Opera vol. 16, Pars 1, p. 333) :
2
Z Z  0 si 6= k
2 2 1 i(k)x 
eix eikx dx = ei(k)x dx = e 

= (5.5)
0 0 i(k ) 0 2 si = k .
Interpolation et Approximation 39

Grace a cette formule, il suffit de multiplier la serie (5.4) par eix et integrer terme par terme de 0
a 2. 5 Tous les termes, sauf un, disparaissent et on obtient
Z
1 2
ck = f (x)eikx dx, (k ZZ) (5.6)
2 0

et, par (5.4),

1 Z 2 1 Z 2
ak = f (x) cos kx dx, bk = f (x) sin kx dx (k 0). (5.7)
0 0

1
Amplitude et phase. En ecrivant ck = 2
rk eik on a encore

ak cos kx + bk sin kx = ck eikx + ck eikx = rk Re eikx+ik = rk cos(kx + k ), (5.8)

i.e. les deux termes cos kx et sin kx se confondent en un terme cos kx damplitude rk = 2|ck | avec
une phase k . Ces amplitudes sont normalement representees dans un spectrogramme.

Transformation DISCRETE de Fourier.

Origines :
Interpolation de fonctions periodiques en astronomie (Clairaut 1754, Euler 1753, Opera 14,
p. 463) ;

Traitement de donnees periodiques en meteorologie (Bessel 1828).

Donnees discretes. 4
y5 y6 y7
Supposons que la fonction 2-periodique f (x) soit seule- 3 y4
ment connue pour les x de la division equidistante 2 y1 y2 y3
1y0
2j
xj = , j = 0, 1, . . . , N 1
N 0 2
et posons yj = f (xj ). Si necessaire, on peut prolonger {yj } a une suite N-periodique en posant
yj+N = yj pour tout entier j (N = 8 dans la petite figure).
Probleme. Trouver des coefficients zk (k = 0, 1, 2, . . . , N 1) tels que
N
X 1 N
X 1 N
X 1
2ikn 2i
yn = zk eikxn = zk e N = zk kn avec =eN n = 0, 1, 2, . . . , N 1 .
k=0 k=0 k=0
(5.9)

Interpolation trigonometrique. Comme

N = 1 ei(N k)xn = eikxn

les termes pour k = N 1, k = N 2, . . . correspondent aux termes k = 1, 2 . . . (voir figure

II.17). Le probleme (5.9) est alors equivalent a rechercher un polynome trigonometrique
N/2
1  X X
p(x) = zN/2 eiN x/2 + zN/2 eiN x/2 + zk eikx =:
zk eikx (5.10)
2 |k|<N/2 k=N/2

5
ce qui peut causer des problemes de legitimite ..., voir le cours dAnalyse IIA ou le cours Math. pour Info.,
Chap. VI.
40 Interpolation et Approximation

2
2 51 1 5 6 6
3 1 1 3 1 3 1 3

76543210 0 1 2 3 4 5 6 7 4 0 0 4 62 40 0 2 4 6 4 0 0 4
0 0 0 0
5 7 5 7 7 5 5 7
6 6 73 3 7 2 2

0
1
Im 1
0 Im 1
0 Im 1
0 Im
2
3 k=0 3
2
k=1 2
3 k=2 2
3 k=3
4 4 4 4
5 5 5 5
6 6 6 6
7 Re 7 Re 7 Re 7 Re

2 2 73 3 7 6 6
7 5 5 7 5 7 5 7

7531 6420 0 1 2 3 4 5 6 7 4 0 0 4 62 40 0 2 4 6 4 0 0 4
0 0 0 0
1 3 1 3 3 1 1 3
6 6 51 1 5 2 2

1
0 Im 1
0 Im 1
0 Im 1
0 Im
3
2
k=4 2
3 k=5 2
3 k=6 3
2
k=7
4 4 4 4
5 5 5 5
6 6 6 6
7 Re 7 Re 7 Re 7 Re

F IG . II.17: Fonctions de base pour la transformation discrete de Fourier

satisfaisant p(xn ) = yn pour n = 0, 1, . . . , N 1. De plus, si les donnees yn sont reelles, nous

aurons zN k = zk = zk . Ainsi, notre polynome devient avec (5.4)
N/21  ak = 2 Rezk , aN/2 = RezN/2 ,
a0 X
p(x) = + ak cos kx + bk sin kx +aN/2 cos(N/2)x
2 k=1 bk = 2 Imzk
(5.11)
et nous voyons clairement la similitude avec la transformee continue de Fourier.
Orthogonalite. Pour la solution du probleme, i.e., le calcul des coefficients zk , nous attendons un
miracle semblable a celui de la condition dorthogonalite ci-dessus (5.5). En fait, nous avons

1( k )N 1( N )k 11k
N
X 1
ixn ikxn
N
X 1
(k)n

1 k
= 1 k
= 1 k
= 0 si 6= k
e e = = (5.12)


n=0 n=0 1 + 1 +...+ 1 = N si = k .

En parfaite analogie avec la preuve ci-dessus pour (5.6), on multiplie la somme (5.9) par eixn et
on additionne termes a termes de 0 a N 1.6 Tous les termes, sauf un, disparaissent et on obtient :

Theoreme 5.1
1 NX
1
1 NX
1 N
X 1
zk = yn eikxn = yn kn yn = zk nk
N n=0 N n=0 k=0 (5.13)
z = FN (y) y = FN1 (z) = N F N (z) .

Relation avec lalgebre lineaire. Les deux formules (5.13) representent un systeme lineaire et
son inverse. La matrice est orthogonale (ou hermitienne comme on dit dans le cas complexe) et
jouit de toutes les belles proprietes quon connat en algebre et en geometrie. En particulier, ces
systemes sont bien conditionnes (details au chapitre IV) ; ici, il ny a pas de mechants phenomenes,
tels que celui de Runge ou lexplosion des erreurs darrondi (de la section precedente).
6
sans probleme de legitimite cette fois-ci.
Interpolation et Approximation 41

200 data eeee

100

0
0 100 200 300 400 500 600 700 800 900 1000
101
spectre periode 1024

100

101
0 50 100 150
200 data periode 930

100

0
0 100 200 300 400 500 600 700 800 900 1000
101
spectre periode 930

100

101
0 50 100 150
101
compression

100

101
0 50 100 150
200 decodage

100

0
0 100 200 300 400 500 600 700 800 900 1000
101
lissage

100

101
0 50 100 150
200 data lissees

100

0
0 100 200 300 400 500 600 700 800 900 1000
F IG . II.18: Le spectrogramme pour un son eee a 118 Hz

Exemple. Dans la premiere image de la fig. II.18, on voit la digitalisation dun son eeee prononce
par un celebre mathematicien 7 . Nous partons dune fenetre de N = 1024 donnees pour un signal
enregistre a 22000 impulsions par seconde ; nous montrons en deuxieme image le spectre |zk |.
Celui-ci est chaotique. Point de miracle : on applique une theorie pour donnees periodiques a
des donnees qui ne le sont apparemment pas, cela etant du a la coupure arbitraire apres 1024
7
les historiens de la science ne sont pas unanimes quant a son identite ; il sagit soit de Martin Hairer, soit de
G. Wanner. Une chose est certaine : tous deux etaient presents lors de lenregistrement...
42 Interpolation et Approximation

points. On constate que yN +n yn est plutot minimale pres de N = 930. On calcule donc 1024
nouvelles valeurs du signal sur une periode de 930 (par interpolation lineaire) ; leur transformee de
Fourier discrete dans la troiseme image est nettement plus claire et montre une frequence principale
(5 22000/930 118Hz), ainsi que des harmoniques (multiples de la frequence principale).
Compression des donnees. Lidee est de supprimer tous les termes de la serie de Fourier dont les
coefficients sont en dessous dun certain seuil (p.ex. 3% du coefficient maximal). Ainsi, la vraie
information contenue dans le signal ne contient que 14 nombres (au lieu de 1024 ; 4eme image de
la figure II.18) et peut toujours etre decompressee par F 1 (5eme image).
Autres applications. La transformation de Fourier permet de nombreuses autres operations sur les
signaux (filtering, denoising, repair, crossing ...) et un excellent code Amadeus II est vivement
conseille (http://www.hairersoft.com/).

II.6 Transformation de Fourier rapide (FFT)

10% dinspiration et 90% de transpiration.
(Un grand scientifique (Einstein?) sur le secret de ses inspirations)
En 1965, John Tukey, celebre statisticien a Princeton et aux Bell Labs, demande au jeune program-
meur J.W. Cooley dinterrompre son travail pendant deux semaines pour laider a programmer une
petite idee et a tester son utilite pratique. Les deux semaines sont devenues 20 annees et la pe-
tite idee un des algorithmes les plus importants du 20eme siecle. Comme souvent dans pareil cas,
on a trouve plus tard des references anterieures dans la litterature, en particulier un expose clair de
la methode dans les cours de Runge (voir Runge-Konig, 1924, 66).
Idee. Cherchons a calculer FN1 (z) par (5.13). Commencons par la transpiration, i.e., ecrivons
cette formule dans tout ses details pour N = 4 :

0 0 0 0 z0 y0

0 1 2 3 z1 y1

=

0 2 0 2
z2


y2

0 3 2 1 z3 y3

Un calcul brut de ce produit necessiterait N 2 multiplications et additions. Essayons detre plus

elegants : lidee est de permuter les donnees z (et les colonnes de la matrice), en prenant dabord
les indices pairs, puis les indices impairs :

0 0 0 0 z0 y0


0 2 1 3
z2


y1

=

0 0 2 2
z1


y2

0 2 3 1 z3 y3
Nous voyons que la matrice se decompose en quatre blocs, dont les deux blocs de gauche sont la
matrice de la transformee de Fourier pour 7 2 , i.e., celle avec N 7 N/2. Comme 2 = 1,
nous avons encore (en ecrivant u pour les donnees paires et v pour les donnees impaires)

0 0 0 0 u0 y0


0 2 1 3
u1


y1

=

0 0 0 0

v0


y2

0 2 1 3 v1 y3
Interpolation et Approximation 43

Maintenant, la structure est parfaitement claire et notre inspiration montre que

1 1
yk = (FN1z)k = (FN/2 u)k + k (FN/2 v)k
(6.1)
yk+N/2 = (FN1 z)k+N/2 = (FN/2
1
u)k k (FN/2
1
v)k

(k = 0, 1, . . . N/21). Pour la transformee de Fourier directe, nous devons tenir compte du facteur
1/N et remplacer par 1 . Ainsi, nous arrivons au resultat :
1 
(FN y)k = (FN/2 u)k + k (FN/2 v)k
2
(6.2)
1 
(FN y)k+N/2 = (FN/2 u)k k (FN/2 v)k .
2
Cest formidable : le meme procede peut etre applique recursivement aux suites u et v, aussi
longtemps que leurs longueurs restent paires. Si lon suppose que N = 2m , lalgorithme developpe
toute sa puissance : on peut pousser les simplifications jusquau bout. Par exemple, pour N = 8 =
23 lalgorithme est presente dans le schema suivant
FN/8 y0 = y0
/
\ FN/8 y4 = y4
 
y0
y0 . (N/4) FN/4
y0
y4

y1 , (N/2) FN/2
y2
/
FN/8 y2 = y2
  /
y2 y 4 y2
(N/4) FN/4 \ FN/8 y6 = y6
y y6 y6
3 -
N FN
y4   (6.3)
y1 . y1 / FN/8 y1 = y1
y
5 y (N/4) FN/4
y5

y6 (N/2) FN/2 3 / \
y5   FN/8 y5 = y5
y7 y y3
7 (N/4) FN/4
y7 /
FN/8 y3 = y3
\
FN/8 y7 = y7
La programmation de cet algorithme se fait en deux etapes. Dabord, on met les {yi} dans lordre
exige par lalgorithme (6.3), c.-a-d. quil faut inverser les bits dans la representation binaire des
indices: . .
0 = (0, 0, 0) 0 = (0, 0, 0)
. .
1 = (0, 0, 1) 4 = (1, 0, 0)
. .
2 = (0, 1, 0) 2 = (0, 1, 0)
. .
3 = (0, 1, 1) 6 = (1, 1, 0)
. .
4 = (1, 0, 0) 1 = (0, 0, 1)
. .
5 = (1, 0, 1) 5 = (1, 0, 1)
. .
6 = (1, 1, 0) 3 = (0, 1, 1)
. .
7 = (1, 1, 1) 7 = (1, 1, 1)
Apres, on effectue les operations de (6.2) comme indiquees dans le schema (6.3). Une explication
detaillee du code en FORTRAN ancien est donnee dans le livre Numerical Recipies.8 Les
programmeurs modernes profiteront de la recursivite (voir TP).
8
W.H. Press, B.R. Flannery, S.A. Teukolsky & W.T. Vetterling (1989): Numerical Recipies. The Art of Scientific
Computing (FORTRAN Version). Cambridge University Press.
44 Interpolation et Approximation

TAB . II.3: Comparaison de nombres doperations

N N2 N log2 N quotient
25 = 32 103 160 6.4
210 103 106 104 100
220 106 1012 2 107 5 104

Pour passer dune colonne a une autre (dans le schema (6.3)), on a besoin de N/2 multi-
plications complexes et de N additions (ou soustractions). Comme m = log2 N passages sont
necessaires, on a
N
Theoreme 6.1 Pour N = 2m , le calcul de FN y peut etre effectue en log2 N multiplica-
2
tions complexes et N log2 N additions complexes.

Pour mieux illustrer limportance de cet algorithme, comparons dans le tableau II.3 le nombre
doperations necessaires pour le calcul de FN y avec ou sans FFT.

II.7 Transformee cosinus discrete (DCT) et JPEG

Lalgorithme du paragraphe precedent marche tres bien si les donnees ont une certaine periodicite.
Si elles ne sont pas periodiques (par exemple pour la compression dimages), on utilise souvent
une variante de la transformee de Fourier discrete. Cette variante a en plus lavantage deviter le
calcul avec des nombres complexes.

Transformee de Fourier en cosinus. Soit f (x) une fonction continue, definie sur lintervalle
[0, ]. On la prolonge en une fonction paire par f (x) = f (x) et en une fonction 2-periodique
par f (x + 2) = f (x). La serie de Fourier (5.1) dune telle fonction contient seulement les termes
en cosinus et secrit
a0 X
f (x) = + ak cos kx. (7.1)
2 k=1

Sur lintervalle [0, ], les fonctions cos kx (pour k 0) verifient la relation dorthogonalite

Z 1
Z   0 si k 6=
cos x cos kx dx = cos( + k)x + cos( k)x dx = /2 si k = 6= 0 (7.2)
0 2 0
si k = = 0.
La demarche habituelle (multiplier (7.1) par cos x et integrer terme par terme de 0 a ) nous
donne
2Z
ak = f (x) cos kx dx. (7.3)
0

Transformee cosinus discrete (DCT). Comme f (0) 6= y1 y2 y3

3
f () en general, nous considerons les N points au milieu
2 y0
des sous-intervalles, c.-a-d. les points
1
(2j + 1) 0
xj = , j = 0, 1, . . . , N 1
2N 0 2
Interpolation et Approximation 45

et nous posons yj = f (xj ) (N = 4 dans la petite figure). Par analogie avec (7.1), nous exprimons
cette suite par (voir la fig. II.19 pour les fonctions de base cos kxj )

z0 NX1
yj = + zk cos kxj . (7.4)
2 k=1

Avec la relation dorthogonalite discrete (pour 0 k, N 1)

N
X 1
1 NX
1  0 si k 6=
cos xj cos kxj = cos( + k)xj + cos( k)xj = N/2 si k = 6= 0 (7.5)
2 j=0
j=0 N si k = = 0.
nous trouvons (multiplier lequation (7.4) par cos xj et additionner de j = 0 a j = N 1) la
transformee cosinus discrete (DCT)

2 NX
1
zk = yj cos kxj . (7.6)
N j=0

La valeur zk de (7.6) est le resultat de la regle du point milieu appliquee a lintegrale dans (7.3).
Transformee cosinus discrete en dimension 2. Une image digitale est donnee par un tableau
{Yi,j }, ou i parcourt les pixels verticalement et j horizontalement. La valeur de Yi,j represente le
niveau de gris (ou la couleur) du pixel (i, j). Motive par lalgorithme precedent, nous essayons
dexprimer ces donnees par
N
X 1 N
X 1
Yi,j = Zek, cos kxi cos xj (7.7)
k=0 =0

ou xj = (2j +1)/2N comme pour la transformee cosinus discrete. Pour compenser le facteur 1/2
dans (7.4), nous utilisons la notation Ze0,0 = Z0,0 /4, Zek,0 = Zk,0/2, Ze0, = Z0, /2 et Zek, = Zk, .
La relation dorthogonalite (7.5) nous permet de calculer les Zk, par la formule

4 NX 1 N
X 1
Zk, = Yi,j cos kxi cos xj . (7.8)
N 2 i=0 j=0

Les fonctions de base cos kxi cos xj (k, = 0, . . . , N 1) sont illustrees dans la fig. II.20.

JPEG (Joint Photographic Experts Group). Lutilisation de cette base dans la compression
dimages est due a Ahmed, Natarajan et Rao (IEEE 1974) 9 . On decompose limage entiere en
blocs de 8 8 pixels et on calcule pour chaque bloc la transformee cosinus discrete (7.8). Les
coefficients Zk, sont alors quantifies et ceux qui sont en-dessous dun seuil sont remplaces par
zero. Pour des images contenant des parties uniformes (par exemple : ciel bleu), seulement deux
ou trois coefficients vont rester ce qui donne un important facteur de compression.
9
Pour des explications plus detaillees et references voir JPEG still image data compression standard par W. B. Pen-
nebaker et J. L. Mitchell, New York, 1992. (Merci a Kyle Granger pour ces renseignements).

F IG . II.19: Fonctions de base pour la transformee cosinus discrete (N = 8)

46 Interpolation et Approximation

F IG . II.20: Fonctions de base pour la transformee cosinus discrete en 2 dimensions (8 8 pixels)

La reconstruction de limage a droite (un conglomerat entre Asterix, Mickey

Mouse et Donald Duck, dessine par un grand artiste) est illustree dans la fig. II.21.
Le premier dessin (en haut a gauche) correspond au premier terme de la somme
(7.7). En rajoutant un terme apres lautre, et en suivant un chemin astucieux en
zig-zag, la reconstruction est demontree.

II.8 Interpolation par fonctions spline

Le mot spline (anglais) signifie languette elastique. On sinteresse a la courbe decrite par
une languette forcee de passer par un nombre fini de points donnes (disons par (xi , yi ) pour i =
0, 1, . . . , n). La fig. II.22 montre le spline passant par les memes donnees que pour la fig. II.1 (pour
pouvoir le comparer avec le polynome dinterpolation).
Interpolation et Approximation 47

F IG . II.21: Une image toujours mieux reconstituee en zig-zag par JPEG

La theorie des splines a ete developpee dans les annees 1950 et 60 par I.J. Schoenberg pour servir
au calcul scientifique (approximations, integrales, equations differentielles) ; de nos jours elle est
costamment appliquee pour la representation de courbes (voir figure II.23) et surfaces en Computer
Graphics.
Mathematiquement, ce probleme peut etre formule comme suit: on cherche une fonction s :
[a, b] IR (a = x0 , b = xn ) satisfaisant

(S1) s(xi ) = yi pour i = 0, 1, . . . , n;

(S2) s(x) est 2 fois continument differentiable;
Z b
(S3) (s (x))2 dx min.
a

Lintegrale dans (S3) represente lenergie de la languette deformee qui, par le principe de Mauper-
tius, est supposee minimale.
48 Interpolation et Approximation

10

s(x)
5

0
2 4 6 8 10

F IG . II.22: Spline cubique (a comparer avec fig. II.1)

F IG . II.23: Un dessin en zig-zag (a gauche) et en splines (a droite)

Theoreme 8.1 Soit a = x0 < x1 < . . . < xn = b une division donnee, s : [a, b] IR et
f : [a, b] IR deux fonctions qui verifient (S1) et (S2). Supposons que
   
s (b) f (b) s (b) = s (a) f (a) s (a) (8.1)

et que s(x) soit un polynome de degre 3 sur chaque sous-intervalle [xi1 , xi ]. Alors, on a
Z b Z b
(s (x))2 dx (f (x))2 dx. (8.2)
a a

Demonstration. Cherchons des conditions suffisantes pour s(x) afin de satisfaire (S3). Pour ceci
nous considerons des fonctions voisines f (x) = s(x) + h(x) ou IR et h(x) est de classe C 2 et
verifie
h(xi ) = 0 pour i = 0, 1, . . . , n. (8.3)
Chaque fonction f (x) satisfaisant (S1) et (S2) peut etre obtenue de cette maniere. La condition
(S3) secrit alors
Z b Z b
2

(s (x)) dx (s (x) + h (x))2 dx
a a
Z b Z b Z b
2
=
(s (x)) dx + 2
s (x)h (x) dx + 2
(h (x))2 dx
a a a

(pour tout IR), ce qui est equivalent a

Z b
s (x)h (x) dx = 0 (8.4)
a
Interpolation et Approximation 49

pout tout h C 2 verifiant (8.3). En supposant s(x) 3 fois differentiable sur chaque sous-intervalle
de la division, on obtient par integration par parties que
b Z b

s (x)h (x)  s (x)h (x) dx = 0. (8.5)
a a

Lhypothese (8.1) implique que la premiere expression de (8.5) est nulle. Comme s (x) est con-
stant sur (xi1 , xi ), disons egal a i , la deuxieme expression de (8.5) devient
Z b n
X Z xi n
X
s (x)h (x) dx = i h (x) dx = i (h(xi ) h(xi1 )) = 0
a i=1 xi1 i=0

par (8.3). Ainsi, (8.4) et par consequent (8.2) aussi sont verifies.

Definition 8.2 Soit a = x0 < x1 < . . . < xn = b une division de [a, b]. Une fonction s C 2 [a, b]
sappelle spline (cubique) si, sur chaque intervalle [xi1 , xi ], elle est un polynome de degre 3.

Pour satisfaire la condition (8.1), on les possibilites suivantes sont les plus utiles:
spline naturel: on suppose que

s (a) = 0 et s (b) = 0. (8.6)

spline periodique: on suppose que

s (a) = s (b) et s (a) = s (b). (8.7)

Calcul du spline interpolant.

Nous supposons les points xi equidistants de distance 1. Pour le cas plus general, les idees sont les
memes, mais les calculs sont plus compliques. Nous utilisons ici une idee de Schoenberg (1946) et
L.L. Schumaker (1969). Une autre approche (utilisant linterpolation dHermite) sera lobjet dun
exercice.
1ere idee: Si nous developpons un spline s(x) dans le voisinage de xi en puissances de (x xi ):

s(x) = a + b(x xi ) + c(x xi )2 + d(x xi )3 ,

les coefficients a, b, c representent les derivees dordre 0, 1, et 2; elles sont donc les memes a droite
et a gauche de xi . La difference dun spline a droite et a gauche de xi est donc un multiple de
(x xi )3 . Si nous introduisons la fonction
(
(x xi )3 si x xi ,
(x xi )3+ = (8.8)
0 si x xi ;

alors un spline aura comme base les fonctions (voir premiere colonne de Fig. II.24)

(x x0 )3+ , (x x1 )3+ , (x x2 )3+ , . . .

Inutile de dire que cette base serait tres mal conditionnee.

Idee des B-splines. Il vaut donc mieux de prendre les differences de ces fonctions. Et puisque les
4emes differences dun polynome de degre 3 sont nulles, nous avons la bonne surprise que les 4
50 Interpolation et Approximation
0

(x 2)3+
1

2 1 0 1 2 2
(x 1)3+
2 1 0 1 2
1 3
6

2 1 0 1 2
1
2 1 0
1
1 2 4
6
5 B(x)
(x)3+ 7
1 12 1
2 1 0 1 2 2 1 0 1 2 4
6 c+ (x)
8
1 19 1 1
1
2 1 0 1 2 2 1 0 1 2 2 1 0 1 2
18 6 6
+ (x) 5
(x + 1)3+ 7
1 12 1
27 2 1 0 1 2 2 1 0 1 2
6
8
1 19 1
2 1 0 1 2 2 1 0 1 2

(x + 2)3+ 7
1
27 2 1 0 1 2

8
1
2 1 0 1 2

F IG . II.24: Differences successives des fonctions (x xi )3+ .

de (xxi )3+ sont a support compact (voir Fig. II.24). Notre base est donc compose des translations
de

(x + 2)3 si 2 x 1,



(x + 2)3 4(x + 1)3 si 1 x 0,


B(x) = 4 (x 2)3+ = (x + 2)3 4(x + 1)3 + 6x3 si 0 x 1, (8.9)





(x + 2)3 4(x + 1)3 + 6x3 4(x 1)3 si 1 x 2,


0 sinon

Petite difficulte. Les B-splines de (8.9) ne reproduisent pas, pour les splines naturels, les polynomes
de degre 1 en dehors des points 0, 1, . . . , n. Nous allons donc rajouter a notre base quelques
deuxiemes et troisiemes differences (en gris dans la Fig. II.24; veut dire lineaire, c veut
Interpolation et Approximation 51

dire constant en dehors des xi ; c (x) = c+ (x) et (x) = + (x) sobtiennent par symetrie)
s(x) = 0 (x) + 1 c (x1) + 2 B(x2) + . . .
+n2 B(x(n2)) + n1 c+ (x(n1)) + n + (xn)
+
c c+ (8.10)
B B B
6 6
5 5
4 4 4

1 1 1 1 1
0 1 2 3 4 5 6
Calcul du spline interpolant. Les conditions s(i) = yi, en introduisant les valeurs des B, c, et ,
deviennent des equations lineaires pour les coefficients i , que voici (pour n = 6)

6 6 0 y0
1 5 1 y
1 1



1 4 1 2

y2

1 4 1 3 = y3 (8.11)


1 4 1 4 y4

1 5 1 5 y5
6 6 6 y6

Remarque. Pour le spline periodique, il ny a que des B(x) tout au tour et la matrice dans
(8.11) na que des 1, 4, 1 avec deux 1 supplementaires dans les coins n, 0 et 0, n.

0
5 4 3 2 1 0 1 2 3 4 5

F IG . II.25: Spline de Lagrange compare au polynome de Lagrange.

Phenomene de Runge, erreurs darrondi, splines de Lagrange. Nous allons proceder a

lelimination des elements sous la diagonale du systeme (8.11):
           
6 6 6 6 4 1 4 1 3.75 1 3.75 1
7 6 7 1 7 1
1 5 0 5 6
1 4 0 4 4
1 4 0 4 3.75
Nous voyons que les elements de la diagonale convergent vite vers une valeur qui satisfait
1
=4 donc =2+ 3 = 3.73205...

Si maintenant les premiers yi sont toutes nulles, nous obtenons pour les i approximativement
1
i + i+1 = 0 ou i = i+1 = 0.268 i+1

 52 Interpolation et Approximation

i.e., en seloignant de la position ou yj = 1, les i tendent vers zero comme une suite geometrique
a rapport 0.268 (une decouverte de M. Powell 1966 ; voir illustration en Fig. II.25). Pour les
splines, il ny a donc ni Phenomene de Runge, ni influence catastrophique des erreurs darrondi.

II.9 Ondelettes (One Lecture on Wavelets)

Origines.
Base de Haar (Math. Annalen vol. 69, 1910) ;

Morlet 1982/83 (Analyse dondes sismiques) ;

MorletGrossmann 1984, Y. Meyer 1985 (Wavelets de Meyer ; 1984 : lannee de nais-

sance des ondelettes . . .) ;

Mallat 1987 (Analyse multiresolution ; 1987 : lannus mirabilis) ;

Daubechies 1987 (Wavelets a support compact) ;

LemarieBattle 1987 (Wavelets spline) ;

Holschneider 1988, Daubechies 1992 (Fast Wavelet Transform).

Litterature.
1. I. Daubechies, Ten Lectures on Wavelets 1992 ;

2. Y. Meyer, Ondelettes 1990 ;

3. P.G. Lemarie, Les Ondelettes en 1989, (Lecture Notes 1438, de nombreuses applications) ;

4. LouisMaassRieder, Wavelets (en allemand) 1994.

Motivation. Lanalyse de Fourier est basee sur les fonctions

eikx =

Elle determine les frequences, mais ne localise pas bien les sons. Prenons pour exemple trois
notes de Beethoven (sol - sol - mibemol) en figure II.26. Lanalyse de Fourier montre bien les deux
frequences dominantes (en rapport 5 sur 4), mais la localisation de ces notes est brouillee dans le
chaos des hautes frequences.

101
200 data spectre
100
100

101
0
0 100 200 300 400 500 600 700 800 900 1000 0 50 100 150
F IG . II.26: Trois notes de la 5eme de Beethoven, la transformee de Fourier et le resultat dune trop forte
compression.

Exemple : la base de Haar. Apres le desastre des theoremes de convergence pour series de
Fourier et fonctions continues (fausse preuve de Cauchy, correction de Dirichlet, contre-exemples
Interpolation et Approximation 53

de Dubois-Reymond, H.A. Schwarz, Lebesgue et Fejer, voir cours dAnalyse II), Hilbert pose a
son etudiant A. Haar le probleme suivant : trouver (enfin) une base de fonctions orthogonales,
ou la convergence est assuree pour toute fonction continue. Le resultat de ces recherches est la
base de Haar (1910, voir figure II.27) ; la convergence (meme uniforme) nest pas trop difficile
a demontrer (voir exercices).
Mais : la base de Haar est composee de fonctions discontinues et la convergence est lente. Donc,
nous ne sommes pas au bout de nos peines. Avant de progresser, nous constatons :
2 2 2 2 2 2 2 2
1 2 3 4 5 6 7 8

1 1 1 1 1 1 1 1

0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1

1 1 1 1 1 1 1 1

2 2 2 2 2 2 2 2

F IG . II.27: La base de Haar

Lemme 9.1 Soit (x) = 1 pour 0 < x < 12 et = 1 pour 21 < x < 1 la fonction 2 de la base de
Haar, alors toutes les autres fonctions de cette base (a lexception de 1 ) sont de la forme

m,n (x) = 2m/2 (2m x n) (m = 1, 2, 3, . . . , n = 0, 1, 2, . . .). (9.1)

Une telle fonction sappelle ondelette mere. Les valeurs de n representent des translations en
x et servent a localiser un evenement, les valeurs de m correspondent aux frequences comme
precedemment.

GRAND PROBLEME. Trouver une ondelette mere (x) (continue, differentiable, a support
compact . . .), pour laquelle toutes les fonctions dans (9.1) (pour tout n et pour tout m) forment
une base orthogonale.

Les diverses reponses a cette question (ondelettes de Meyer, ondelettes de Daubechies . . .) neces-
sitent des moyens lourds danalyse fonctionnelle (transformee continue de Fourier), seules les
ondelettes spline de LemarieBattle sont, pour nous, relativement faciles a comprendre.
Comme cest souvent le cas, il est conseille de sattaquer dabord a un probleme plus facile :
Probleme plus facile. Trouver une fonction (x) (quon appellera le pere de londelette) orthog-
onale a tous ses fonctions translatees
Z
(x n) (x m) dx = 0 pour tout n, m ZZ, n 6= m. (9.2)

54 Interpolation et Approximation

0 = 0.328294
4 1 = 0.112072
1.0 2 = 0.044784
3 = 0.019753
3
4 = 0.009199
5 = 0.004420
.5
(x) 6 = 0.002166
2
7 = 0.001076

1
.0
4 2 0 2 4

2 1 0 1 2

F IG . II.28: B-spline (a gauche), pere orthogonal de londelette de BattleLemarie (a droite).

Calcul du pere spline de londelette. Nous partons dun spline a support compact, le B-spline.
Au lieu de la condition (9.2), cette fonction satisfait

604


d0 = 35
si n = m,



1191



d1 = 140
si |n m| = 1,
Z
6
B(x n) B(x m) dx = d2 = 7
si |n m| = 2, (9.3)



1

d3 = 140
si |n m| = 3,



0 sinon
Nous posons alors
X
(x) = 0 B(x)+1 (B(x1)+B(x+1))+2 (B(x2)+B(x+2))+. . . = k B(xk) (9.4)
k

ou nous supposons pour linstant 0 = 1. Les conditions dorthogonalite (9.2) deviennent alors,
grace a (9.3),
d0 S + d1 (S1 + S+1 ) + d2 (S2 + S+2 ) + d3 (S3 + S+3 ) = 0 (9.5)
ou
1X
S = k k = + S
2 k
ou S est la meme somme sans les termes k = 0 et k = . Nous resolvons alors les conditions (9.5)
par calcul numerique brutal en posant iterativement
d1 d2 d3
:= S (S1 + S+1 ) (S2 + S+2 ) (S3 + S+3 ) = 1, 2, 3, . . . .
d0 d0 d0
ApresRla convergence de cet algorithme en une trentaine diterations, nous normalisons les pour

avoir ((x))2 dx = 1 et nous arrivons aux valeurs presentees en figure II.28 a droite et a la
fonction recherchee.
Analyse multi-echelles.
Cette analyse (Mallat 1987) est la clef pour trouver finalement londelette mere. Lidee est de
raffiner pas a pas la resolution par un facteur 2 en remplacant x 7 2x. Commencons par le lemme
suivant :
 Interpolation et Approximation 55

Lemme 9.2 Lespace V 0 engendre par

. . . B(x + 3), B(x + 2), B(x + 1), B(x), B(x 1), B(x 2), B(x 3), . . .

est un sous-espace de V1 engendre par

. . . B(2x + 3), B(2x + 2), B(2x + 1), B(2x), B(2x 1), B(2x 2), B(2x 3), . . .

et nous avons
X2
1 1 3 1 1
B(x) = B(2x + 2) + B(2x + 1) + B(2x) + B(2x 1) + B(2x 2) = ck B(2x k).
8 2 4 2 8 k=2
(9.6)

4 4 4

3 3 3

2 2 2

1 1 1

2 1 0 1 2 2 1 0 1 2 2 1 0 1 2
B(x) 81 (B(2x + 2) + B(2x 2)) 12 (B(2x + 1) + B(2x 1)) . . . = 43 B(2x)
F IG . II.29: Decomposition de B(x) dans la base des B(2x n).

La preuve est indiquee en figure II.29. Vers lextremite gauche B(x) = (x + 2)3 , par contre
B(2x + 2) = (2x + 2 + 2)3 = 8(x + 2)3 . Donc B(x) 18 (B(2x + 2) + B(2x 2)) a le support
reduit des deux cotes. Puis on reduit une deuxieme fois en soustrayant 12 (B(2x + 1) + B(2x 1)).
Il reste un spline en 2x dont le support est [1, +1]. Celui doit donc etre un multiple de B(2x).
Vers une base fan orthogonale de V 1 . Les fonctions (xn) forment une base orthogonale de
V 0 . Les fonctions (2x n) forment une base orthogonale de V 1 , mais elle nest pas orthogonale
a V 0 . Essayons de trouver une formule du type

X
(x) = 2 hn (2x n) , (9.7)
n=

formule qui correspond a (9.6) si on remplace B par . Nous developpons les deux cotes de cette
formule dans la base des B(2x m) :
X X XX 
(x) = B(x ) = ck B(2x 2 k) = cm2 B(2x m)
,k m

pour lun, et
X X X X 
2 hn (2x n) = 2 hn B(2x n) = 2 hm B(2x m)
n n, m
56 Interpolation et Approximation

pour lautre. Nous egalisons les deux termes et inserons les ck de (9.6). Cela nous donne un grand
systeme lineaire pour les hm que nous ecrivons sous la forme

1 3 1
X 1 8 m2 1 + 4 m2 + 8 m2 +1 si m est pair
0 hm = hm +
6=0 2 12 m1 + 21 m+1 si m est impair.
2 2

Vu que 0 est le terme dominant, nous resolvons ce systeme par iterations et obtenons, en une
dizaine de cycles, les valeurs
h0 = 0.7661300 h1 = 0.4339226 h2 = 0.0502017 h3 = 0.1100370
h4 = 0.0320809 h5 = 0.0420684 h6 = 0.0171763 h7 = 0.0179823
h8 = 0.0086853 h9 = 0.0082015 h10 = 0.0043538 h11 = 0.0038824 .
h12 = 0.0021867 h13 = 0.0018821 h14 = 0.0011037 h15 = 0.0009272
h16 = 0.0005599 h17 = 0.0004621 h18 = 0.0002853 h19 = 0.0002324
Theoreme 9.3 Si lon pose, avec les valeurs hn de (9.7),
X
(x) = 2 gn (2x n) avec gn = (1)n h1n , (9.8)
n=

alors les fonctions

. . . (x + 3), (x + 2), (x + 1), (x), (x 1), (x 2), (x 3), . . .
sont orthogonales entre elles et orthogonales aux fonctions (x n). En dautres mots, lespace
W 0 engendre pas ces s forme avec V 0 une decomposition orthogonale de V 1 :
V 1 = V 0 W0 .
Une illustration est donnee en figure II.30 a gauche : Supposons que seulement h0 et h1 soient 6= 0.
Alors g0 = h1 et g1 = h0 et le vecteur (x) est visiblement orthogonal a (x).
Le miracle est maintenant quon peut continuer en posant a nouveau x 7 2x et on aura succes-
sivement des espaces V m , engendres par B(2m x n), avec
V 2 = V 1 W1 = V 0 W0 W1 , V 3 = V 0 W0 W1 W2
etcetera. Ainsi toutes les (2m x n) sont orthogonales. Notre ondelette mere est enfin trouvee et
est fierement presentee dans toute sa beaute en figure II.30 a droite.
Il reste a faire la preuve des orthogonalites : on calcule
Z X Z
(x j)(x k) dx = 2 hn hm (2x 2j n)(2x 2k m) dx
m,n
Z X Z
(x j)(x k) dx = 2 hn gm (2x 2j n)(2x 2k m) dx
m,n
Z X Z
(x j)(x k) dx = 2 gn gm (2x 2j n)(2x 2k m) dx
m,n

A cause de lorthogonalite des (2x n), les trois sommes doubles se reduisent a trois sommes
simples X X X
hn hn2 hn gn2 gn gn2 .
n n n
Celle du milieu est zero par definition des gn . Car on sait que les (x n) sont orthogonaux, la
premiere somme doit etre zero. Donc la troisieme somme, qui est la meme par definition des gn ,
est zero aussi.
Interpolation et Approximation 57

(2x 1) g-1 = 0.050202 g0 = 0.433923

1
g-2 = 0.110037 g1 = 0.766130
V0
g-3 = 0.032081 g2 = 0.433923
h1 (x) g-4 = 0.042068 g3 = 0.050202
V1 g4 = 0.110037
h1 (x) g5 = 0.032081
h0 (2x)
2 1 0 1 2 3

h0 (x)

W0

F IG . II.30: Analyse multi-echelles (a gauche) ; Wavelet spline de BattleLemarie (a droite).