COURS SYS806

Fait par A. Soulaïmani, Ph.D.

INTRODUCTION

1.1 Représentation des nombres

Un nombre réel x   est représenté par un nombre en point flottant fl(x)   , c.a.d. x est
tronqué par la machine en fl(x).

Formats de Matlab :

Exemple : x= 1/7

Format long 0.14285714285714

Format short e 1.4286e-01
Format long e 1.428571428571428e-01
Format long g 0.142857142857143

x  (1) s .(0.a1a2 ...at ). e  (1) s .m. e t

s  0 ou 1
  2 la base
m un entier appelé la mantisse de longueur t
t nombre d ' entiers significatifs stockés
e un entier , l 'exp osant L  e  U

En MATLAB : (  , t , L, U )  (2,53, 1021,1024)

53 nombres significatifs en base 2 correspondent à 15 en base 10

1
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

Erreur d’arrondi
x  fl ( x) 1
 M avec  M   1t
x 2

En Matlab :  M   1t  252  2.22 1016

peut être obtenu avec la commande eps : c ' est le zero de la machine

xmin   L 1
xmax   U (1    t )

En Matlab : xmin  2.225073858507201.10308 real min

xmax  1.7976931348623158.10308 real max

Algèbre de l’ensemble  :
Commutativité pour l’addition : fl(x+y)= fl(x)+fl(y)
Commutativité pour la multiplication : fl(x y)= fl(x) fl(y)
L’associativité est violée pour l’addition dans une situation de Overflow ou Underflow :
Exemple : a= 1.0e +308; b= 1.1e+308 et c=-1.001e+308
a+(b+c) (a+b)+c

Les erreurs d’arrondi peuvent être importantes dans certaines situations :

(1  x)  1
Exemple :
x
Si x est très petit on peut avoir de sérieux problèmes numériques :
>> format short e
>> x= 1.e-15

>> ((1+x)-1)/x
ans = 1.1102e+000

2
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

(1  x)  1
=1.1102 au lieu de 1, soit 11% d’erreur!
x

Autre exemple
>> format bank
>> b= -54.32; a=1;c=0.1;x1 = (-b+sqrt(b^2-4*a*c))/(2*a)

ans x1 = 54.32 (au lieu de 54.3218)

>> b= -54.32; a=1;c=0.1;x2 = (-b-sqrt(b^2-4*a*c))/(2*a)
ans x2= 0.00
au lieu de 0.001841

L’ordre de convergence :
L’erreur numérique totale = ec = erreurs de troncature du modèle numérique + erreurs
d’arrondis

Le modèle numérique est dit convergent lorsque : ec = C hp avec h le paramètre de

discrétization, p un entier positif et C indépendant de h.

---------------------------------------------------------------------------------------------------------------------
Exemples :

1-Utiliser la fonction funtool pour représenter

f(x)= (x-1)^7
f(x)= x^7-7*x^6+21*x^5-35*x^4+35*x^3-21*x^2+7*x-1

2- La série:
z2  2

zn 1  2n 1/ 2 (1  1  41 n zn2 n  2,3,...

converge vers  . Représentez avec Matlab l’erreur relative obtenue numériquement.

Représentez la convergence du processus numérique.

3
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

1.2 Quelques commandes de MATLAB

fplot : fplot(fun,lims,tol,n,’LineSpec’)

» fun='1/(1+x^2)'; lims=[-5,5];fplot(fun,lims,9,'o')
>> grid
» x=5;y=eval(fun)

fzero : fzero(fun,x0)

>> syms x
>> fun='x^2-1+exp(x)';
fzero(fun,1)

>> fun='x^7-7*x^6+21*x^5-35*x^4+35*x^3-21*x^2+7*x-1';
>> fzero(fun,1)

ans =

>> fzero(fun,0)

ans =

1.0092

>> fun= 'x^7+21*x^5+35*x^3+7*x-1-7*x^6-35*x^4-21*x^2';

4
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

>> fzero(fun,0)

ans =

1.0087

polyval: y= polyval(p,x)

>> p= [1 0 0 0 0 3 0 -1 ]; x= [-1:0.25:1];
>> y=polyval(p,x)

y =

Columns 1 through 8

1.0000 0.5540 -0.2578 -0.8126 -1.0000 -0.8124 -

0.2422 0.8210

Column 9

3.0000

roots : donne les zéros du polynôme dont les coefficients sont

dans p

>> roots(p)

ans =

5
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

-0.3739 + 1.2305i
-0.3739 - 1.2305i
0.9698 + 0.7716i
0.9698 - 0.7716i
-1.1793
-0.5840
0.5716

À cause des erreurs d’arrondis on peut obtenir de faux résultats!

Exemple : fun= (x-1)^7 = x^7-7*x^6+21*x^5-35*x^4+35*x^3-
21*x^2+7*x-1

>> p =[1 -7 21 -35 35 -21 7 -1];

>> roots(p)

ans =

1.0101
1.0063 + 0.0079i
1.0063 - 0.0079i
0.9977 + 0.0099i
0.9977 - 0.0099i
0.9909 + 0.0044i
0.9909 - 0.0044i

conv : p= conv(p1,p2) polynôme formé par le produit p1 x p2

deconv : [q,r]= deconv(p1,p2) p1= q p2+ r

>> p1= [1 0 0 0 -1];

>> p2= [1 0 0 -1];

6
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

>> p= conv(p1,p2)

p =

1 0 0 -1 -1 0 0 1

>> [q,r]= deconv(p1,p2)

q =

1 0

r =

0 0 0 1 -1

diff : donne f’(x)

» syms x
»» diff(f)

x
int : donne l’intégrale F(x)= 
a
f (t )dt

>> syms x;
>> f= 'x';
>> int(f,1,x)

7
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

ans =

1/2*x^2-1/2

taylor : taylor(f,x,n+1) donne le développement en série de Taylor jusqu'à

l’ordre n autour de x0= 0 :
( x  x0 ) n ( n )
Tn ( x)  f ( x0 )  ( x  x0 ) f ( x0 )  .... 
'
f ( x0 )
n!

>> syms x;
>> f= 'exp(x)';
>> taylor(f,x,5)

ans =

1+x+1/2*x^2+1/6*x^3+1/24*x^4

taylortool

8
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

INTERPOLATION

On a n+1 couples de données ( xi , yi  f ( xi )) i = 0 ,…,n ; on cherche à construire l’interpolation


f ( x) telle que : f ( x )  y . L’expression de la fonction f(x) est en général inconnue, on cherche
i i

plutôt une approximation f ( x) qu’on peut expliciter.

Plusieurs choix sont possibles, par exemple :


- Interpolation polynomiale : f ( xi )   n f ( x)  a0  a1 x  ...  an x n


- Interpolation trigonométrique : f ( xi )  a M e  iMx  ...  a0  ...  aM e  iMx

avec M= n/2 si n est pair ou M= (n-1)/2 s’il est impair; i2=-1

 a0  a1 x  ...  ak x k
- Interpolation en fraction rationnelle : f ( xi ) 
ak 1  ak x  ...  an x n
- Autres types…

2.1 Les polynômes de Lagrange en une dimension:

On montre qu’il existe un et un seul polynôme d’interpolation noté  n f ( x) de degré n qui

passe par les (n+1) points distincts xi.

n (x  x j )
Construction : posons k ( x)  
j  0, j  k ( xk  x j )
, k  0,..., n

alors  n f ( x )   yk  k ( x )
k 0

9
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

polyfit : c= polyfit(x,y,n) donne les coefficients du polynôme de

Lagrange, comme suit : c(1) = an, …, c(n+1)= a0.

Estimation de l’erreur d’interpolation :

Il existe une valeur ξ dans l’intervalle  x0 , xn  tel que l’écart entre la fonction exacte

(inconnue) et son interpolation est donné par :

f ( n 1) ( ) n
En f ( x)  f ( x)   n f ( x) 
(n  1)! i 0
 ( x  xi )

avec f ( n 1) ( ) la n+1 dérivée de f(x) au point ξ.

En corollaire : si les points des données sont équidistants : xi+1=xi+ h, avec h le pas
d’espacement constant, alors

max( f ( n 1) ( x)) n 1

max( En f ( x))  h
4(n  1)

Cette formule montre que l’augmentation de l’ordre n ne garantit pas nécessairement la

convergence. Une divergence se produit lorsque la (n+1) dérivée n’est pas bornée ou lorsque
l’augmentation de n produit un polynôme très oscillant (phénomène de Gibbs/Runge).

On peut aussi montrer que l’erreur maximale sur la première dérivée est bornée selon la
relation :

max( f ( x)(1)  ( n f n )(1) ( x) )  C h n max( f ( n 1) ( x) )

10
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

Exemple 1 : déterminer empiriquement l’erreur maximale sur l’interpolation

clc;
clear all;
%interpolation de la fonction f(x) = 1/(x^2+1)
n= 10;
a= -5;
b= 5;
h= (b-a)/n;

%liste des abscisses dans l’intervalle (a,b)

for k=1:n+1
x(k)= a+ (k-1)*h;
end

% les ordonnées des points de données

for i=1:n+1
y(i) = 1/(1+x(i)^2);
end

%construction du polynome de Lagrange

c=polyfit(x,y,n);

%vecteur auxiliaire servant à la visualisation et le calcul de

%l erreur sur 1000 points
xv= linspace(a,b,1000);

%evaluer la fonction d’interpolation sur les 1000 points

p= polyval(c,xv);

for i=1:1000
fx(i) = 1/(1+xv(i)^2);
end
p;
fx;

% erreur maximale empirique

err= max(abs(p-fx))
plot(x,fx,'-',x,p,'o')
legend('fx= exacte','p=interpol. n=10')

Le graphique montre des oscillations proches des bornes a et b.

Le calcul d’erreur donne:
n= 5 » err = 0.4327

11
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

n= 10 » err = 1.9156
n= 15 » err = 2.1069
n= 20 » err = 59.7684
Ce qui montre une divergence, c’est le phénomène de Gibbs/Runge; il est causé par le fait
que les dérivées d’ordre élevés ont un maximum divergent. Ce phénomène se produit aussi
en techniques d’apprentissage-automatique (Machine Learning) et est appelé ‘overfitting’.

2
fx= exacte
p=interpol. n=10

1.5

0.5

-0.5
-5 -4 -3 -2 -1 0 1 2 3 4 5

12
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

10

-10

fx= exacte
-20 p=interpol. n=20

-30

-40

-50

-60
-5 -4 -3 -2 -1 0 1 2 3 4 5

2.2 Interpolation de Chebyshev:

Le phénomène d’oscillations peut être évité si on a la possibilité de choisir la distribution des

nœuds. On considère l’intervalle  a  x0 , b  xn  pour lequel on choisit les nœuds :

ab ba ^ ^ i
xi   xi , xi   cos( ), i  0,..., n
2 2 n

Si f est continuellement différentiable, on peut montrer qu’avec cette distribution  n f

converge vers f lorsque n   . A noter que les noeuds de Chebyshev ne sont pas
équidistants et sont resserrés proche des extrémités de l’intervalle.

Exemple 2:
clc;
clear all;
n= 10;
a= -5;

13
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

b= 5;
xc= -cos(pi*[0:n]/n);
x= (a+b)/2+ (b-a)*xc/2;
for i=1:n+1
y(i) = 1/(1+x(i)^2);
end
c=polyfit(x,y,n);
xv= linspace(-5,5,1000);
p= polyval(c,xv);
for i=1:1000
fx(i) = 1/(1+xv(i)^2);
end
p;
fx;
err= max(abs(p-fx))
plot(x,fx,x,p,'o')
legend('fx= exacte','p=interpol. n=10')

1
fx= exacte
0.9 p=interpol. n=10

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

14
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

1
fx= exacte
0.9 p=interpol. n=20

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

Le calcul d’erreur donne:

n= 5 » err = 0.6386
n= 10 » err = 0.1322
n= 15 » err = 0.0993
n= 20 » err = 0.0177
Ce qui montre une convergence.

2.3 Interpolation trigonométrique et FFT (à titre de complément)

On désire approximer une fonction périodique f sur l’intervalle :  0, 2  , i.e. f (0)  f (2 ) , par

une fonction f ( x) telle que :


f ( xi )  f ( xi ) pour les nœuds : x j  2 j /(n  1) j  0,..., n

15
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

L’interpolation trigonométrique est obtenue par une combinaison linéaire de cosinus et sinus. En
particulier, pour n pair :
^ a0 M
f ( x)     ak cos(kx)  bk sin(kx) , M  n / 2
2 k 1
Et si n est impair :
^ a0 M
f ( x)     ak cos(kx)  bk sin(kx)   aM 1 cos(( M  1) x), M  (n  1) / 2
2 k 1

Nous pouvons réécrire ces deux relations sous la forme :

^ M 
f ( x)   ce
k  ( M  )
k
ikx

  0 si n est pair et M  n / 2,
  1 si n est impair et M  (n  1) / 2

A cause de l’analogie avec la série de Fourier, l’interpolation est nommée ‘la transformée
discrète réelle de Fourier’.

Si on impose les conditions d’interpolation aux nœuds:

x j  jh avec h  2 /(n  1) j  0,..., n ; on obtient les coefficients:

1 n
cm  
n  1 j 0
f ( x j )e imjh

m  ( M   ),...., ( M   )

Le calcul de ces coefficients peut se faire de manière optimale (avec le nombre total
d’opérations de l’ordre n log 2 n ) en utilisant la transformée rapide de Fourier (commande fft).

On peut aussi obtenir les valeurs  f ( xi )  à partir des coefficients  cm  en utilisant la

transformée inverse (commande ifft).

16
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

Exemple :
On a la fonction : f ( x)  x( x  2 )e  x , pour x   0, 2  . On choisit n= 9

clc;
clear all;
n= 9;
x = 2*pi/(n+1)*[0:n];
y= x.*(x-2*pi).*exp(-x);
Y = fft(y); ‘donne les coefficients  cm 
%calcul de la tronsformée inverse
Z= ifft(Y);
plot(x,y,'*', x,Z,'o')
err= max(abs(Z-y))
legend('y= données','ifft. n=9')

-0.2

-0.4

-0.6
y= données
-0.8 ifft. n=9

-1

-1.2

-1.4

-1.6

-1.8

-2
0 1 2 3 4 5 6

interpft : donne l’interpolation trigonométrique pour un ensemble de données.

Elle a besoin comme entrées de : un entier N et un vecteur de valeurs prises par une fonction
périodique de période p aux noeuds x j  jp / M j  1,..., M  1

17
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

La commande retourne les N valeurs de l’interpolant trigonométrique, obtenues par la
transformée de Fourier , aux noeuds t j  jp / N j  0,..., N  1 .

Exemple :
On considère f ( x)  x( x  2 )e  x , pour x   0, 2  .

On prend ses valeurs aux 10 points équidistants : x j  j 2 /10 j  0,...,9 (M= 10).

On obtient son interpolée trigonométrique par la FFT, et on l’évalue aux nœuds

t j  j 2 /100 j  0,...,99 (N=100)

clc;
clear all;
M= 10;
x1 = 2*pi/M*[0:M-1];
y1= x1.*(x1-2*pi).*exp(-x1);
z = interpft(y1,100);

N= 100;
x = 2*pi/N*[0:N-1];
y= x.*(x-2*pi).*exp(-x);

plot(x,z,'*',x,y,'o')

legend('z= = interpol trigo. n=9','y= donnees')

18
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

0.5
z= = interpol trigo. n=9
y= donnees
0

-0.5

-1

-1.5

-2

-2.5
0 1 2 3 4 5 6 7

Dans certains cas, l’interpolée trigonométrique peut être très érronnée.

Exemple :

on considère f ( x)  sin( x)  sin(5 x), pour x   0, 2  et M=9 N=100

>> x=2*pi/100*[0:99];y=sin(x)+sin(5.*x); plot(x,y)

19
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

1.5

0.5

-0.5

-1

-1.5

-2
0 1 2 3 4 5 6 7

Son interpolée trigonométrique est donnée par :

>> x=2*pi/9*[0:8];y=sin(x)+sin(5.*x);z=interpft(y,100);zz= 2*pi/99*[0:99];
plot(zz,z)

1.5

0.5

-0.5

-1

-1.5

-2
0 1 2 3 4 5 6 7

En superposant ces deux figures :

20
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

2
z= = interpol trigo. n=9
1.5 y= donnees

0.5

-0.5

-1

-1.5

-2
0 1 2 3 4 5 6 7

Si le nombre de points M n’est pas suffisant pour représenter les hautes fréquences on peut
avoir d’importantes erreurs (phénomène d’aliaising).
clc;
clear all;
M= 10;
x1 = 2*pi/(M+1)*[0:M];
y1=sin(x1)+sin(5.*x1);
z = interpft(y1,100);

N= 99;

x = 2*pi/(N+1)*[0:N];
y=sin(x)+sin(5.*x);
plot(x,z,'*',x,y,'o')
legend('z= = interpol trigo. n=10','y= donnees')

21
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

2
z= = interpol trigo. n=10
1.5 y= donnees

0.5

-0.5

-1

-1.5

-2

-2.5
0 1 2 3 4 5 6 7

2.4 Interpolation polynomiale par morceaux

L’interpolation de Chebychev permet d’obtenir une bonne approximation de fonctions

‘régulières’ dont l’expression est connue. Il est plus simple d’utiliser une autre interpolation
qui est polynomiale (exemple linéaire) par morceaux et qui converge.

Considérons une liste de points : x0  x1......  xn . Sur l’intervalle Ii=  xi , xi 1  on construit

une fonction linéaire (c.-à-d. le polynôme de Lagrange d’ordre 1) qui passe par les valeurs
yi  f ( xi ) et yi 1  f ( xi 1 ) :

yi 1  yi
1h f ( x)  yi  ( x  xi )  yi 0 ( x)  yi 11 ( x)
xi 1  x i
xi 1  x x  xi
avec : 0 ( x)  ; 1 ( x) 
xi 1  x i xi 1  x i

avec i=0,…,n-1 et l’indice h= max(xi+1 - xi)

22
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

L’étude théorique de la convergence donne:
h2
max ( E1 f ( x)) xI  max( f (2) ( x))
8
avec I=  x0 , xn  .

Par conséquent, pour tout x  I , 1h f converge vers f (quadratiquement ou à l’ordre 2) à

condition que la dérivée seconde soit bornée dans l’intervalle I.

La commande interp1 de Matlab (et son équivalent en Python) permet de construire cette
interpolation.

De même, on a le résultat concernant l’erreur sur l’approximation de la première dérivée :

max ( E1 f (1) ( x)) xI  C h max( f (2) ( x))
avec C une constante.

Exemple 3:
Revenons à l’interpolation de la fonction f(x) = 1/(x2+1). Avec l’interpolation
linéaire par morceaux, on trouve une approximation qui converge vers la solution exacte.

clc;
clear all;
%interpolation de la fonction f(x) = 1/(x^2+1)
n= 20;
a= -5;
b= 5;
h= (b-a)/n;
for k=1:n+1
x1(k)= a+ (k-1)*h;
end
for i=1:n+1
y(i) = 1/(1+x1(i)^2);

23
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

end
s1= interp1(x1,y,x1) % polynome d interpolation

x= linspace(-5,5,1000); % points servant aux predictions et au traçage des courbes

for i=1:1000
fx(i) = 1/(1+x(i)^2);
end
s= interp1(x1,y,x) % construire le polynôme avec les points x1 et faire les prédictions sur x
err= max(abs(s-fx))
plot(x,fx,'-',x1,s1,'o')
legend('fx= exacte','s1=interpol. n=20')

1
fx= exacte
0.9 s1=interpol. n=20

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

2.5 Interpolation par des fonctions splines

L’interpolation linéaire par morceaux permet de construire une fonction continue qui passe
par les valeurs nodales (on dit qu’on a une fonction interpolante de classe C0), mais il n’y a
aucune raison pour qu’elle soit aussi de dérivée première continue (si c’est le cas on dit qu’on
a une fonction de classe C1). Dans certaines applications, comme en DAO, on désire former
des interpolations plus régulières que les fonctions C0.

Les fonctions cubiques spline s3 ( x) :

24
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

On peut alors construire sur chaque sous intervalle Ii=  xi , xi 1  un polynôme de degré 3

(cubique) :
s3 ( x)  c0  c1 ( x  xi )  c2 ( x  xi ) 2  c3 ( x  xi )3

sous les 4 conditions :

s3 ( x)  f ( xi ),
(s3 )' ( xi )  (s3 )' ( xi ), (c.a.d continuité de la première derivée en chaque point x i )
et
(s3 )'' ( xi )  (s3 )'' ( xi ) (c.a.d continuité de la deuxième derivée en chaque point x i )
ce qui donne 4 n  2 conditions pour 4 n coefficients à trouver .

On impose en plus deux autres conditions comme :

s3'' ( x0 )  0 et s3'' ( xn )  0

Exercice :
Construisez un programme en Matlab pour implémenter cette interpolation. Appliquez le
pour la fonction f(x)= 1/(1+x2) dans l’intervalle  5,5 et analyser la convergence de

l’erreur.

Exemple 4 :
On a des données discrètes et on veut construire une interpolation de classe C2
clc;
clear all;
x= [-55:10:65];
y = [-3.25 -3.37 -3.35 -3.2 -3.12 -3.02 -3.03 ...
-3.07 -3.17 -3.32 -3.3 -3.22 -3.1];
plot(x,y,'o')

z= [-55:1:65];

s= spline(x,y,z);

hold on
plot(z,s)

25
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

3- Régression par combinaison linéaire de fonctions de base

La régression (approximation ou lissage) par la méthode des moindres carrés est très utilisée
pour construire des fonctions d’approximation lorsqu’on a une grande quantité de données.
Considérons la fonction de régression f ( x)

f ( x)  c  ( x)  c  ( x)    c  ( x), (E)
1 1 2 2 m m

où les fonctions de base  j ( x), j 1,  , m sont indépendantes des coefficients c j .


Ces fonctions ne sont pas nécessairement des polynômes. La fonction f ( x) est une combinaison
linéaire des  j ( x) et c j sont des paramètres inconnus. L’objectif est de déterminer les valeurs de

c1 ,  , cm de façon à minimiser l’écart entre les données et l’approximation.

Pour être plus précis, supposons d’abord qu’on cherche les m paramètres de telle façon que

f ( x) passe exactement par les n données, cela entraîne les n équations :

26
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

c1 1 ( x1 )  c2 2 ( x1 )  ... cm m ( x1 )  y1 ,
c1 1 ( x2 )  c2 2 ( x2 )  ... cm m ( x2 )  y2 ,

c1 1 ( xn )  c2 2 ( xn )  ... cm m ( xn )  yn .

qui se mettent sous forme d’un système à m inconnues,

Bc  Y
où:

1 ( x1 ) 2 ( x1 )... m ( x1 )  c1   y1 

 ( x )  ( x )...  ( x )  c  y 
B  1 2 2 2 m 2 
, C  2 
et Y   2  .
      
     
1 ( xn ) 2 ( xn )... m ( xn )  cm   yn 

Dim(B)= (n,m) et dim (Y)=n alors que dim(C)= m

Dans le cas où (n  m) le système est surdéterminé et il n’admet pas de solution. La méthode des
moindres carrés permet de relaxer le problème en cherchant un compromis. On cherche les m
paramètres optimaux qui minimisent l’erreur:
n
 (C )   ( yi  f ( xi )) 2  Y  BC
2

i 0

Cela entraîne la résolution du système:

( B T B ) C  BT Y

La matrice BT B est symétrique et positive de dimension (m, m). Dans le cas où m=n et si les
lignes de la matrice B sont indépendantes (cela dépend du choix des fonctions de base et des
points de données), alors ce système est équivalent à :
B C Y

Fonctions de base

27
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

L’équation (E) exprime la forme générale de la fonction de lissage. Les fonctions de base sont
par exemple :
1, x , x 2 , x 2 / 3 , sin x , e x , xe 4 x , cos (ln 25 x )

Par contre, les fonctions : sin ( c1 x ) , e c3 x , x c2 ne sont pas valides (elles sont non-linéaires)
puisqu’elles sont dépendantes des paramètres c j .

Dans le cas de fonctions de base polynomiales, les commandes polyfit et polyval de MATLAB
permettent d’effectuer le lissage.

Algorithme de la régression par moindres carrés

1. Choisir les fonctions de base:  j ( x), j 1,  , m.

2. Lire les données, ( xi , yi ), i 1,  , n.

3. Calculer la matrice B.
4. Résoudre le système dont le vecteur inconnu est C.

Exemple : Considérons le lissage de paires données ( x , y ) sous la forme :

 c1
f  c2 x.
x

avec
x 0.955 1.380 1.854 2.093 2.674 3.006 3.255 3.940 4.060
y 5.722 4.812 4.727 4.850 5.011 5.253 5.617 6.282 6.255

La fonction de lissage est obtenue par les instructions suivantes:

>> B= [1./x x]; % matrice rectangulaire

>> C = (B’*B)\(B’*y) % résoudre le système
C=

28
COURS SYS806

Fait par A. Soulaïmani, Ph.D.

4.2596
1.3008

Par la suite, on évalue la fonction d’interpolation en d’autres points dans l’intervalle des x.

>> xf = linspace(min(x),max(x)); % 100 points dans l’intervalle des x

>> yf = c(1)./xf + c(2)*xf; % Evaluer l’interpolation en xf


Comme x f est un vecteur colonne, on peut écrite y f  f ( x f ) comme

 y f ,1  1 / x f , 1   x f ,1 
     
 y f ,2   c 1 / x f , 2   c  x f ,2  .
   1    2   
     
 y f , p  1 / x f , p   x f , p 

y f  Bf c ,

avec

1/ x f ,1 x f ,1 
 
 1/ x f ,2 x f , 2 
Bf   .
  
 
1/ x f , p x f , p 

>> xf = linspace(min(x),max(s));
>> Bf = [1./xf xf];
>> yf = Bf*c:
>> plot(x,y,’o’,xf,yf,’-‘);

29