Interpolation et Approximation par des B-splines

Pierre Pansu February 9, 2004

Interpolation

Interpoler une famille de N + 1 points Qj de Rn , cest trouver une courbe qui passe par tous ces points. En plus de prescrire des positions, on peut aussi vouloir prescrire des tangentes. Le prototype du probl`me dinterpolation consiste ` se donner a = t0 < t1 < < tN = b et e a Q0 , . . . , QN R, et chercher une fonction x : [a, b] R telle que f (ti ) = Qi pour tout i = 0, . . . , N , voir gure 1.

1.1

Polynme dinterpolation de Lagrange o

Une rponse classique ` la question consiste ` chercher f parmi les polynmes de degr au plus N . e a a o e Thor`me 1 Etant donns a = t0 < t1 < < tN = b et Q0 , . . . , QN R, il existe un unique e e e polynme x : [a, b] R de degr N tel que x(ti ) = Qi pour tout i = 0, . . . , N . On lappelle le o e polynme dinterpolation de Lagrange. o Preuve. Soit E lespace vectoriel des polynmes ` coecients rels, de degr N . Alors dimE = o a e e N + 1. On consid`re lapplication linaire e e E RN +1 , x (x(t0 ), . . . , x(tN )).

Pour montrer quelle est bijective, il sut de vrier quelle est injective. Or si x(t0 ) = = e x(tN ) = 0, le polynme x a N + 1 > deg(f ) racines distinctes, donc il est nul. o

1.2

Erreur dinterpolation

La qualit dun procd dinterpolation se mesure par une estimation derreur. Il sagit de vrier e e e e si linterpolante approche bien une fonction donne susamment rguli`re. Autrement dit, tant e e e e

Q0 Q3 Q
2

Figure 1: Interpolation en dimension 1

donne une fonction f susamment direntiable, de majorer la dierence entre f et la fonction e e x qui interpole f aux points ti , i.e. telle que x(ti ) = f (ti ) pour i = 0, . . . , N . Pour le polynme dinterpolation de Lagrange, lestimation derreur classique est la suivante. o Thor`me 2 Soit f une fonction de classe C N +1 sur [a, b] et xLag,N +1 le polynme dinterpolation e e o de Lagrange de f aux points a = t0 < < tN = b. Alors f xLag,N +1 o` qN +1 (t) = u Preuve.
N i=0 (t

1 (N + 1)!

qN +1

f (N +1)

ti ).

Posons g = f xLag,N +1 . Fixons t [a, b] distinct des ti . Posons k(u) = g(u) qN +1 (u)g(t)/qN +1 (t).

Alors k(ti ) = 0, k(t) = 0 et k (N +1) = f (N +1) (N + 1)!g(t)/qN +1 (t). Le thor`me de Rolle donne e e un zro de k dans chacun des N + 1 intervalles dlimits par les ti et t. Une nouvelle application e e e du thor`me de Rolle donne N zros pour k . On continue jusqu` trouver un zro pour k (N +1) e e e a e dans lintervalle [a, b]. Il vient |g(t)| = 1 qN +1 (t)f (N +1) () (N + 1)! 1 qN +1 f (N +1) (N + 1)!

1.3

Des polynmes aux splines o

La faiblesse de linterpolation de Lagrange, cest que lerreur dinterpolation cro avec N . Cela t se traduit exprimentalement par de grandes oscillations du polynme dinterpolation, mme si e o e f est tr`s simple. Par exemple (Runge 1901), lorsquon interpole la fonction x 1/1 + 25x2 e en des points uniformment rpartis sur lintervalle [1, 1], les polynmes de Lagrange xLag,N ne e e o convergent pas vers f . Do` lide dinterpoler par des fonctions polynmiales par morceaux, dont u e o le degr naugmente pas avec le nombre de points dinterpolation. e Par exemple, la fonction X1 qui est ane sur chaque intervalle [ti , ti+1 ] et interpole une fonction f de classe C 2 satisfait, dapr`s le thor`me 2 appliqu sur chaque intervalle, e e e e f X1

2h2 f (2)

Comme on souhaite interpoler par des fonctions plus direntiables que X1 , on cherche linterpoe lante dans lespace des fonctions de classe C k1 polynmiales de degr k sur chaque intervalle o e [ti , ti+1 ], appeles fonctions splines (la terminologie a t introduite par Schoenberg en 1946 : en e ee anglais, spline dsigne une bande de mtal souple utilise par les dessinateurs pour tracer une e e e jolie courbe entre deux points). Cest ainsi que sont nes les fonctions et les courbes B-splines. e

1.4

Interpolation par des courbes B-splines

On se limite au courbes de degr 3, pour simplier. On cherche ` faire passer une courbe B-spline e a de degr k de positions et vitesses aux extrmits prescrites par N 1 points Qi . Le probl`me se e e e e divise en deux phases. Premi`re phase : On se xe un vecteur de noeuds t et on cherche un polygone de contrle P tel e o que la courbe B-spline Xk correspondante passe par les Qi aux noeuds. Linterpolation se traduit alors par la rsolution dun syst`me linaire. e e e Deuxi`me phase : on cherche ` optimiser le choix du vecteur de noeuds. Cest typiquement e a non linaire. e

splines

Figure 2: Linterpolante est une projection orthogonale

1.5

Le probl`me linaire e e

A nouveau, on se contente dnoncer le thor`me pour des B-splines de degr 3. Une gnralisation e e e e e e a ` tout degr impair se trouve dans [K], Theorem 8.30 page 172. e Thor`me 3 Soient Q0 , . . . , QN des points de Rn . Soient va , vb deux vecteurs de Rn . Soit t un e e vecteur de noeuds viss aux extrmits, de la forme e e e t0 = t1 = t2 = t3 = a < t4 < < tN +2 < b = tN +3 = tN +4 = tN +5 = tN +6 . Il existe un unique polygone de contrle P = (P0 , . . . , PN +2 ) tel que la courbe B-spline de degr 3 o e associe satisfasse e j = 0, . . . , N, X3 (tj+3 ) = Qj , X3 (a) = va , et X3 (b) = vb .

Preuve. Comme chaque coordonne se traite indpendamment, on peut supposer que n = 1. e e Dans ce cas, on consid`re lapplication linaire e e RN +3 RN +3 , (P0 , . . . , PN +2 (X3 (a), X3 (t3 ), . . . , X3 (tN +3 ), X3 (b)).

Pour prouver quelle est bijective, il sut de prouver quelle est injective. Cela rsulte immdiatee e ment du lemme suivant appliqu ` la fonction f = 0. En eet, si les points et vecteurs interpols ea e sont tous nuls, le lemme donne X3 0. Avec la condition initiale X3 (a) = X3 (a) = 0, cela entra ne que X3 0, et donc que les Pi sont tous nuls. Lemme 1.1 Soient f , x : [a, b] R deux fonctions de classe C 2 . On suppose que x est polynmiale de degr 3 sur chaque intervalle [ti , ti+1 ], i = 3, . . . , N + 2 ; o e f (ti ) = x(ti ) pour i = 3, . . . , N + 3 et f (a) = x (a), f (b) = x (b). Alors
a b

(f (t) x (t))2 dt =
a

f (t)2 dt

b a

x (t)2 dt.

Remarque. Ce lemme signie que la spline x qui interpole une fonction f est la projection b orthogonale de f sur le sous-espace des splines, pour le produit scalaire f g = a f (t)g (t) dt, voir gure 2. Preuve. On vrie immdiatement que e e
b a

(f (t) x (t))2 dt
a

f (t)2 dt +
a

x (t)2 dt = 2R

o` u R=
a

(f (t) x (t))x (t)dt.

On int`gre par parties sur chaque intervalle e

(f (t)
ti

x (t))x (t) dt = (f (ti+1 ) x (ti+1 ))x (ti+1 ) (f (ti ) x (ti ))x (ti )
ti+1

ti

(f (t) x (t))x (t) dt

(f (ti+1 ) x (ti+1 ))x (ti+1 ) (f (ti ) x (ti ))x (ti ) (f (ti+1 ) x(ti+1 ))x (ti+1 ) (f (ti ) x(ti ))x (ti +)
ti+1

+
ti

(f (t) x(t))x(4) (t) dt

(f (ti+1 ) x (ti+1 ))x (ti+1 ) (f (ti ) x (ti ))x (ti ),

car x(4) 0 et f (ti+1 ) x(ti+1 ) = f (ti ) x(ti ) = 0. En additionnant, il vient R = = (f (b) x (b))x (b) (f (a) x (a))x (a) 0,

car f (a) x (a) = f (b) x (b) = 0.

1.6

Estimation de lerreur dinterpolation

Thor`me 4 Soit f : [a, b] R une fonction de classe C 2 . Soit X3 la fonction B-spline de degr e e e 3 qui linterpole en N + 1 point plus les drives aux bornes, selon le thor`me 3. Alors e e e e f X3 o` h = max{|ti+1 ti |}. u Preuve. Posons g = f X3 . Le lemme 1.1 donne g
2 2 b

h3/2 2

et

f X3

h1/2 f

=
a b

g (t)2 dt f (t)2 dt
a 2 2 b

=
a

X3 (t)2 dt

Par construction, g sannule aux ti , donc, dapr`s le thor`me de Rolle, g sannule au moins une e e e fois dans chaque intervalle [ti , ti+1 ]. Tout point t [a, b] est donc ` distance au plus h dun point a t tel que g (t ) = 0. On crit e |g (t)| = = |g (t) g (t )|
t

|
t t

g (s) ds| ds)1/2 (

( h
t 1/2

g (s)2 ds)1/2

2 2,

h1/2 f 4

Q Q

Q4 Q
2

Figure 3: Paramtrisation chordale e

dapr`s lingalit de Cauchy-Schwarz. Ceci montre que e e e g

h1/2 f

Tout point t [a, b] est ` distance au plus h/2 dun point t tel que g(t ) = 0. Par consquent a e |g(t)| = |g(t) g(t )|
t

= |
t

g (s) ds|
2

|t t | g h1/2 2 f

Remarque. Dans [K], Theorem 8.34 page 178, on trouve lamlioration suivante. Sous les hye poth`ses du Thor`me 4, et si de plus f est de classe C 4 , alors e e e f X3

h4 16

f (4)

1.7

Choix du vecteur de noeuds

Etant donns les points Qi ` interpoler, quel est le meilleur choix des ti ? Pour viter que la e a e drive de la courbe interpolante soit grande, il vaut mieux que des points Qi et Qi+1 loigns e e e e soient interpols en des valeurs ti et ti+1 loignes. Autrement dit, il faut corrler les espacements e e e e ti+1 ti avec les distances Qi+1 Qi . Un choix simple consiste ` poser a ti+1 ti = Qi+1 Qi . Ce choix sappelle la paramtrisation cordale (chordal parametrization). Elle a pour eet de e produire une courbe interpolante paramtre ` vitesse a peu pr`s constante, voir gure 3. e e a ` e Lorsque les distances Qi+1 Qi varient brutalement, il y a de meilleurs choix, voir [HL], gures 4.36a et 4.36b pages 205 et 206. On gagne de la souplesse ` xer des valeurs sj o` interpoler indpendantes des noeuds et ` xer a u e a ensuite les noeuds de faon ` optimiser le conditionnement de la matrice du syst`me linaire, voir c a e e [HL].

2
2.1

Approximation
Position du probl`me e

Le probl`me dapproximation dune famille de N + 1 points Qj par une courbe consiste, au lieu e dexiger que la courbe passe exactement par les Qj , i.e. x(sj ) = Qj (ce qui ncessite autant de e 5

Q Q

Q4 Q
2

t1

t2

t3

t4

Figure 4: Approximation dun nuage de points par une courbe

degr de libert que dquations, ` savoir n(N + 1)), ` chercher ` minimiser lcart quadratique e e e a a a e x(sj ) Qj
j 2

o` x varie dans un espace possdant typiquement bien moins de degrs de libert, voir gure 4. u e e e Exemple. Le probl`me dapproximation de N + 1 points Qj par une application constante, est e bien pos. e En eet, lintervalle [0, 1] est subdivis par exemple par sj = j/N , j = 0, . . . , N . Linconnue e est une application constante de [0, 1] dans Rn , i.e. un point x Rn . Il doit minimiser
N

E(P, v) =
j=0

x Qj

La solution est le centre de gravit des points Qj . Dans ce probl`me, linconnue x vit dans un e e espace de dimension n < n(N + 1).

2.2

Approximation par une courbe B-spline

Montrons que le probl`me se ram`ne ` nouveau ` la rsolution dun syst`me linaire. e e a a e e e Soit a = t0 = = tk < tk+1 < < tmk = = tm = b un vecteur de noeuds viss aux extrmits. Soit P = (P0 , . . . , Pmk2 ) Rn(mk1) un polygone e e e de contrle. La courbe B-spline de degr k correspondante est donne sur [a, b] par Xk (t) = o e e mk2 Bi,k (t)Pi . Soient a s0 < < sN b des valeurs distinctes dans [a, b] et Q0 , . . . , QN i=0 des points dans Rn . Le probl`me dapproximation consiste ` trouver P qui minimise e a
N

E(P) =
j=0

Xk (sj ) Qj

On remarque que E est un polynme du second degr sur Rn(mk1) . Si on reprsente les coo e e ordonnes de tous les Pi simultanment par une colonne V ` n(m k 1) coecients, alors E e e a scrit e E(P) = T V AV + 2 T BV + c o` A est une matrice symtrique positive, B est la colonne des coordonnes des vecteurs u e e Bi,k (sj )Qj ,
j

et c=
j

Qj

La fonction E atteint donc son minimum aux vecteurs V solutions de lquation vectorielle e AV + B = 0. Il sagit dun syst`me linaire ` n(m k 1) quations et n(m k 1) inconnues. Noter que la e e a e taille ne dpend pas du nombre de points ` approximer. e a Typiquement, lapproximation est utilise pour allger des donnes : remplacer une courbe e e e composite ou comportant un grand nombre de points de contrle par une seule B-spline ayant peu o de points de contrle. o

Rsolution numrique dun probl`me dinterpolation e e e

On se donne le vecteur de noeuds

dans lintervalle [0, N ]. Il sagit de trouver le polygone de contrle (` N +3 sommets) de la B-spline o a qui passe par le point Qi en ti+3 et a pour drives v0 (resp. vN ) aux extrmits. On peut supposer e e e e que n = 1. On utilise les valeurs des fonctions B-splines relatives au vecteur de noeuds ci-dessus (en exercice du cours de gomtrie direntielle). e e e X3 (0) = P0 , X3 (0) = 3(P1 P0 ), 7 1 1 P1 + P2 + P3 , X3 (t4 ) = 4 12 6 1 2 1 X3 (ti+3 ) = Pi + Pi+1 + Pi+2 6 3 6 (1) (2) (3) (4)

pour i 2. Il sagit de rsoudre le syst`me AP = Q pour Q = (Q0 , v0 , Q1 , . . . , QN 1 , vN , QN ) et P RN +3 e e et 1 0 0 0 0 3 3 0 0 0 1 7 1 A = 0 4 12 6 0 . 0 0 1 2 1 6 3 6 . . . . . .. . . . . . . . . . . . La matrice A nest pas tout ` fait symtrique. Comme on prf`re les matrices symtriques, on a e ee e rsoud le syst`me quivalent e e e T AAP = T AQ. La matrice T AA est coercive, 1 P 2. 3 Cela indique que la rsolution devrait tre ecace. On applique la mthode du gradient conjugu. e e e e Ap Ap

3.1

Rappels sur la mthode du gradient conjugu e e

Soit E un espace de Hilbert et q une forme quadratique dnie positive sur E, q(x) = Ax x. e Alors x est solution de lquation Ax = y si et seulement si x est le minimum de la fonction e 1 e a z f (z) = 2 Az z y z. La mthode consiste ` chercher les minima xr de la fonction f sur les sous-espaces vectoriels croissants Er = vect{y, Ay, . . . , Ar y}. Alors xr xr1 est donn par une formule simple au moyen de produits scalaires. La suite xr e converge vers la solution vritable x . En particulier, si E est de dimension N , alors xN = x . e Lexprience prouve que xN/5 est souvent dj` tr`s proche de x . e ea e Proposition 3.1 On dnit trois suites xr , vr et er par x0 = 0, v0 = y, les relations de rcurrence e e xr+1 = xr vr+1 = er+1 et er = Axr y.
1 Alors la fonction f : x 2 Ax x y x atteint son minimum sur Er en xr . La suite xr converge vers lunique solution x de lquation Ax = y. De plus, on a lestimation a posteriori e

er vr vr , Avr vr Avr er+1 vr Avr vr

A(x xr ) (x xr ) = er (xr+1 xr ). Si A est coercive, i.e. x

C Ax x, on a lestimation a posteriori x xr C Axr y .

Preuve. Dnissons xr comme le point de Er o` la restriction de f ` Er atteint son minimum. e u a Il existe et est unique par stricte convexit. Pour la mme raison, si E est ladhrence de la e e e runion des Er , la fonction uniformment convexe f atteint son minimum sur E en un unique e e point x . On utilise le produit scalaire associ ` la forme quadratique q, ea a A b = Aa b. Relativement ` cette structure euclidienne, les ensembles de niveau de f sont des sph`res centres a e e en x . Les ensembles de niveau de la restriction de f ` Er sont des sph`res centres en xr . Soit a e e s r 1. Au point xs o` f atteint son minimum sur Es cette sph`re est tangente ` Es , donc u e a xr xs A Es . En particulier x xs A Es . Posons er = Axr y. Alors er Er+1 . Or er = A(xr x ) Er , donc en particulier er er1 = 0. Considrons la suite vr dnie par v0 = y et par la relation de rcurrence e e e vr+1 = er+1 vr A er+1 vr . vr A vr

Montrons par rcurrence sur r que vr Er+1 , vr A Er et e xr+1 = xr er vr vr . Avr vr

Lorsque r = 0, v0 = y E1 par dnition et E0 = {0}. e

Supposons connu que vr1 Er , vr1 A Er1 et er1 vr1 xr = xr1 vr1 . Avr1 vr1 On remarque que er AEr Er+1 donc vr Er+1 . Dautre part, x xr A Er , donc er = A(x xr ) A A1 Er . Or AEr1 Er donc er A Er1 . Comme vr est un vecteur du plan engendr par er et vr1 , il est orthogonal ` Er1 . e a Par construction, vr A vr1 . Si vr1 = 0, alors dapr`s lhypoth`se de rcurrence, xr = xr1 , donc er = er1 . Comme e e e er er1 = O, ncessairement er1 = 0, donc x = xr1 . On conclut que xr+1 = xr , vr = 0 et la e relation de rcurrence pour xr+1 est satisfaite. e On peut donc supposer que vr1 = 0. Alors Er est engendr par Er1 et vr1 , donc vr A Er . e Si vr = 0, alors er est colinaire ` vr1 . Or er Er1 . Ncessairement er = 0, donc xr+1 = xr e a e et la relation de rcurrence pour xr+1 est ` nouveau satisfaite. e a On peut donc supposer que vr = 0. Alors vr est une base de lorthogonal de Er1 dans Er , donc il existe un rel tel que xr+1 xr = vr . Comme cette valeur de minimise f (xr + vr ), e ncessairement e xr A vr y vr = , vr A vr comme annonc. e Par construction, la distance (mesure au moyen de x AX x) de x au sous-espace Er est e atteinte en xr . Comme la runion des Er est dense dans E , cette distance, gale ` A(x xr ) e e a (x xr ), tend vers 0. Remarquer que les vecteurs xr+1 xr sont deux ` deux orthogonaux (ils a sont colinaires aux vecteurs obtenus en orthonormalisant la base y, Ay,...,Ar y...). On a donc e

A(xr+1 xr ) (xr+1 xr ) = Ax x .
r=0

Enn, xr x
2

C A(xr x ) (xr x ) = (Axr y) (xr x ) Axr y xr x

do` u x xr C Axr y .

3.2

Mthode de Jacobi e

On peut aussi utiliser la mthode de rsolution de Jacobi (voir [K], Theorem 4.2 page 55). Soit D e e est la diagonale de A. Alors Ax = y x = D1 (A D)x + D1 y et on it`re la transformation ane x D1 (A D)x + D1 y. Celle-ci est contractante lorsque e Rn+1 est muni de la norme 1 . En eet, q = D1 (A D)
1 1

= max{|Ajj |1
i=j

|Aji | ; j = 0, . . . , m}

qui dans le cas prsent vaut q = 1/2. Posant x0 = O et e x+1 = D1 (A D)x + D1 y, on peut armer que la solution x satisfait x x et dans le cas prsent e q = 21 . 1q 9
1

q 1q

D1 y

References
[HL] J. HOSCHEK, D. LASSER, Grundlagen der geometrischen Datenverarbeitung. Teubner, Stuttgart (1989). English translation : Fundamentals of computer aided geometric design. A. K. Peters, Wellesley, Mass. (1993). [K] R. KRESS, Numerical Analysis. Graduate Texts in Math. 181. Springer, Heidelberg (1998). [R] J.-J. RISLER, Mthodes mathmatiques pour la CAO. Masson, Paris (1991). e e

10