Chapitre III : Polynômes et interpolation polynomiale
                                   Résolution des équations non linéaires

1. Opérations sur les polynômes dans MATLAB
    1.1. Multiplication des polynômes
    1.2. Division des polynômes

2. Manipulation de fonctions polynomiales dans MATLAB
    2.1. Évaluation d’un polynôme
    2.2. Interpolation au sens des moindres carrés

3. Interpolation linéaire et non linéaire

4. Interpolation de Lagrange

5. Résolution d’équations et de Systèmes d’équations non Linéaire
    5.1. Résolution d’équations non Linéaires
    5.2. Résolution de Systèmes d’équations non Linéaires

L’objectif principal de l’interpolation est d’interpoler des données connuesà partir des points discrets. Dans ce cas, la valeur de la fonction entre ces points peut être estimée. Cette méthode d’estimation peut être étendue et utilisée dans divers domaines ; à savoir la dérivation et l’intégration numérique des polynômes.
 

1. Opérations sur les polynômes dans MATLAB

Dans MATLAB, les polynômes sont représentés sous forme de vecteurs lignes dont les composantes sont données par ordre des puissances décroissantes. Un polynôme de degré n est représenté par un vecteur de taille (n+1).

Exemple :

Le polynôme : est représenté par :

>>f=[8 0 2 -3 4 -2]
f =8 0 2 -3 4 -2

D’autres fonctions dans MATLAB telles que : ‘conv’, ‘deconv’, ‘roots’, etc. peuvent être utilisées en plus des opérations propres aux vecteurs (cf. chapitre précédent).
 

1.1. Multiplication des polynômes

La fonction ‘conv’ donne le produit de convolution de deux polynômes. L’exemple suivant montre l’utilisation de cette fonction.

Soient :

Le produit de convolution :  est donné par :

>>f=[3 2 -1 4];

>>g=[2 0 -3 5 -1];

>>h=conv(f,g)
h =
6 4 -11 17 10 -19 21 -4

Ainsi, le polynôme obtenu est :

1.2. Division des polynômes

La fonction ‘deconv’ donne le rapport de convolution de deux polynômes (déconvolution des coefficients du polynôme). L’exemple suivant montre l’utilisation de cette fonction.

Soient les mêmes fonctions précédentes  et  :

La division de  par  :

est donnée par la fonction ‘deconv’ :

>>f=[3 2 -1 4];

>>g=[2 0 -3 5 -1];

>>h=deconv(g,f)
h =
0.6667 -0.4444

et le polynôme  obtenu est :

2. Manipulation de fonctions polynomiales dans MATLAB

Soit le polynôme suivant :

où n est degré du polynôme et  sont les coefficients du polynôme.

Ce polynôme peut être écrit sous la forme :

Après factorisation, on a :

où  sont les racines du polynômes .
 

Exemple :

Ce polynôme est équivalent à :

ou encore :

Un polynôme d’ordre n possède n racines qui peuvent être réelles ou complexes.

Dans MATLAB, un polynôme est représenté par un vecteur contenant les coefficient dans un ordre décroissant.

Exemple :

Le polynôme : qui est représenté dans MATLAB par :

>>P=[2 1 4 5];

a pour racines . Pour trouver ces racines, on doit exécuter la fonction ‘roots’.

D’où :

>>r=roots(P);

et le résultat donné est :

>> r
r =
0.2500 + 1.5612i
0.2500 - 1.5612i
-1.0000

Les trois racines de ce polynôme (dont 2 sont complexes) sont données sous forme d’un vecteur colonne. Quand les racines  sont connues, les coefficients peuvent être recalculés par la commande ‘poly’.

Exemple :

>>poly(r)
ans =
1.0000 0.5000 2.0000 2.5000

La fonction ‘poly’ accepte aussi une matrice comme argument dont elle retourne le polynôme caractéristique.

Exemple :

>>A=[3 1;2 4];

>>p=poly(A)
p =
1 -7 10

Ainsi, le polynôme caractéristique de la matrice A est :

Les racines de ce polynôme sont les valeurs propres de la matrice A. ces racines peuvent être obtenues par la fonction ‘eig’ :

>>Val_prop=eig(A)
Val_prop =
2
5
 

2.1. Évaluation d’un polynôme
 

• Pour évaluer le polynôme  en un point donné, on doit utiliser la fonction ‘polyval’. On évalue ce polynôme pour x=1, par exemple :

>>polyval(P,1)
ans =
12
 

• Exemple :

>>C=[3 -7 2 1 1];

>>x=2.5;

>> y=polyval(C,x)
y =
23.8125

Si x est un vecteur contenant plusieurs valeurs, y sera aussi un vecteur qui contiendra le même nombre d’éléments que x.
 

2.2. Interpolation au sens des moindres carrés

Dans le domaine de l’analyse numérique des données, on a souvent besoin d’établir un modèle mathématique liant plusieurs séries de données expérimentales. L’interpolation polynomiale consiste à approcher la courbe liant les deux séries de mesures par un polynôme. Les coefficients optimaux de ce polynôme sont ceux qui minimisent la variance de l’erreur d’interpolation. Ce principe (voir cours) est connu sous le nom de la méthode des moindres carrés. La fonction ‘polyfit’ retourne le polynôme P de degré n permettant d’approcher la courbe  au sens des moindres carrés.

Exemple :

>> x=[1.1 2.3 3.9 5.1];

>>y=[3.887 4.276 4.651 2.117];

>>a=polyfit(x,y,length(x)-1)
a =
-0.2015 1.4385 -2.7477 5.4370

Ainsi, le polynôme d’interpolation de y (d’ordre length(x)-1=3) est :

Pour déduire l’erreur entre les valeurs expérimentales et le modèle obtenu par la fonction ‘polyfit’, on dispose de la fonction ‘polyval’ qui retourne la valeur du polynôme P pour toutes les composantes du vecteur (ou de la matrice) x.

Ainsi, cette fonction donne :

>>yi=polyval(a,x)
yi =
3.8870 4.2760 4.6510 2.1170

On remarque ici que les valeurs expérimentales de y sont bien restituées (y=yi). Dans ce cas, on a un coefficient de corrélation qui est égal à 1 (voir la définition de ce coefficient dans le cours).

Pour mieux comprendre ces fonctions prédéfinis dans MATLAB, on va simuler une courbe expérimentale par une sigmoïde à laquelle on superpose un bruit du type Gaussien. Cette courbe sera donnée par :

Le programme correspondant à l’approximation des ces données dans MATLAB (regres.m, par exemple) est le suivant :

clc; % Effacer l'écran
clear all; % Effacer des variables de l'espace de travail
x=-5:0.1:5; % Intervalle de définition et de calcul de la sigmoïde

% Fonction sigmoïde bruitée
y=1./(1+exp(-x))+0.05*randn(1,length(x));
plot (x,y); % Tracé de la sigmoïde bruitée
title('Fonction sigmoïde bruitée - Polynôme d''interpolation');
xlabel('x');ylabel('y');

% Polynôme d'interpolation d'ordre 1
P=polyfit(x,y,1);

% Valeurs du polynôme d'interpolation
Vp=polyval(P,x);

% Tracé du polynôme d'interpolation
hold on;

plot(x,Vp,'--');
% Calcul de l'erreur d'interpolation
erreur=y-Vp;
% Tracé de la courbe de l'erreur
plot(x,erreur,':')
grid

gtext('Mesures')
gtext('Erreur')
gtext('Modèle')
hold off
% Affichage du polynôme d'interpolation
disp('Polynôme d''interpolation')
P
Var_erreur=num2str(std(erreur).^2);
disp(['La variance de l''erreur d''interpolation est : ',Var_erreur])

>> regres
Polynôme d'interpolation
P =
0.1309 0.5008

La variance de l'erreur d'interpolation est : 0.011277

On remarque ici que le polynôme d’interpolation d’ordre 1 n’est pas une bonne approximation pour ces données. En effet, un polynôme d’ordre 5 donne une meilleure approximation de la sigmoïde (fig. 2). Ainsi, dans le programme précédent (sigreg.m), on change dans le ‘1’ par ‘5’ dans la fonction ‘polyfit’, et on obtient les résultats suivants :

>> regres
Polynôme d'interpolation
P =
0.0002 -0.0000 -0.0111 0.0008 0.2326 0.4844

La variance de l'erreur d'interpolation est : 0.002279

3. Interpolation linéaire et non linéaire

Une interpolation consiste à relier les points expérimentaux par une courbe sous forme de segments de droites ou de courbes polynomiales. Ceci peut être réalisé par la fonction ‘interp1’. La commande ‘interp1(x,y,x_i,’type’)’ retourne un vecteur de mêmes dimensions que x_i et dont les valeurs correspondent aux images des éléments de x_i déterminées par interpolation sur x et y. Si f est l’interpolation de y, la chaîne ‘type’ spécifie alors le type d’interpolation qui doit être parmi les suivants :

• ‘linear’ à interpolation linéaire
• ‘spline’ à interpolation par splines cubiques,
• ‘cubic’ à interpolation cubique.

Exemple 1 :

>>x = 0:10; y = sin(x); xi = 0:.25:10;
>>yi = interp1(x,y,xi,’cubic’); plot(x,y,'o',xi,yi)

Fig. 3 : Interpolation cubique

Exemple 2 :

Dans le cas suivant, on étudiera ces différents types d’interpolation sur un même exemple de valeurs discrètes de la fonction ‘cosinus’. On appellera l’algorithme : ‘interpol.m’.

clear all;
clc;
x=0:10;
y=cos(x); % Points à interpoler
z=0:0.25:10; % Le pas du vecteur z est inférieur à celui de x

% Interpolation linéaire
figure(1);
f=interp1(x,y,z);

% Tracé des valeurs réelles et de la courbe d'interpolation
plot(x,y,'*r',z,f);
grid on;
xlabel('Interpolation');

% Interpolation par splines cubiques
figure(2);
f=interp1(x,y,z,'spline');
plot(x,y,'*r',z,f);
grid on;
xlabel('Interpolation par splines cubiques');

En exécutant ce programme, on obtient les courbes suivantes :

La fonction ‘interp2’ réalise l’interpolation dans l’espace trois dimensions (3D).
 

4. Interpolation de Lagrange

Exemple :

Les masses volumiques du matériau pour différentes températures sont données par le tableau ci-dessous :
 

i	1	2	3
Température T (en °C)	94	205	371
Masse volumique R(T) : (en kg/m3)	929	902	860

• Écrire la formule d’interpolation de Lagrange qui permet d’interpoler les différents points de données précédentes.
• Trouver les masses volumiques du sodium pour T=251 °C, T=305 °C et T=800 °C en utilisant l'interpolation de Lagrange.

a) Puisque le nombre de points est égal à 3, alors le polynôme de Lagrange sera de degré 2. Ce polynôme s’écrit :

Soit :

D’où pour T=251 °C, on obtient R(251) = 890.5 kg/m³.

L'algorithme de calcul de R(251), R(305) et R(800) est : ‘polylag.m’ qui est listé ci-dessous :

La fonction ‘lag.m’ est un sous programme permettant de donner l’interpolation de Lagrange. La liste de ce sous programme est :

for i=1:n
    y=ones(size(Ti));
    for j=1:n
        if i~=j
            y=y.*(Ti-T(j))/(T(i)-T(j));
        end
        Ri=Ri+y*R(i)
    end
end
return

Il suffit maintenant d’exécuter le programme ‘polylag.m’ pour obtenir les résultats de Ri.

Ces résultats sont les suivants :
 

>>Ri
Ri =
890.5561 876.9316 742.45559
 

5. Résolution d’équations et de Systèmes d’équations non Linéaire

5.1. Résolution d’équations non Linéaires

Dans cette partie, on s’intéressera à la méthode de Newton-Raphson pour la résolution des équations non linéaires à une ou à plusieurs variables. On cherchera ensuite à trouver la valeur  qui annulera la fonction .

La méthode de Newton-Raphson permet d’approcher par itérations la valeur  au moyen de la relation suivante :

Dans toutes les méthodes itératives, il est nécessaire pour éviter une divergence de la solution, de bien choisir la valeur initiale . Celle-ci peut être obtenue graphiquement.
 

On se propose d’appliquer cette méthode pour la recherche des racines de la fonction non linéaire suivante :

Dans un premier temps, on se propose de tracer la courbe représentative de cette fonction en utilisant le programme ci-dessous ‘pnum1.m’:

x=-1:0.1:1;
f=exp(x)-2*cos(x);
plot(x,f); grid on;
title('Fonction : f(x)=exp(x)-2.cos(x)');

La fonction dérivée  a pour expression : .

Pour chercher la solution de , on peut rajouter au programme précédent ‘pnum1.m’ quelques lignes :

clear all;
clc;
x(1)=input('Donner la valeur initiale x(1): \n');
e=1e-10;
n=5000;
for i=2:n
    f=exp(x(i-1))-2*cos(x(i-1));
    diff=exp(x(i-1))+2*sin(x(i-1));
    x(i)=x(i-1)-f/diff;
    if abs(x(i)-x(i-1))<=e
        xp=x(i);
        fprintf('xp=%f\n',x(i));
        break;
    end
end

j=1:i;
figure(2);
plot(j,x(j),'*r',j,x(j));
xlabel('Nombre d''itérations');
title('Convergence de la solution : Méth. de Newt.-Raph.');
disp('Les valeurs successives de x(i) sont :');
x'

>>pnum1
Donner la valeur initiale x(1):
0.5
xp=0.539785
Les valeurs successives de x(i) sont :
ans =
0.5000
0.5408
0.5398
0.5398
0.5398

Fig. 7 : Évolution de la solution x(i) en fonction du nombre d’itérations

Dans MATLAB, on peut obtenir la même solution avec la fonction ‘fzero’. Pour cela, il faut créer le fichier.m MATLAB (‘f.m’, par exemple) dans lequel sera programmé .

Pour obtenir la solution de  au voisinage de , on exécute la commande ‘fzéro(‘f’,0.5)’. Cette fonction affiche les valeurs obtenues de  à chaque itération pour donner à la fin la solution recherchée .

>>d=fzero('f',0.5)
f =
-0.1064
-0.1430
-0.0692
-0.1579
-0.0536
-0.1788
-0.0313
-0.2080
5.8950e-004
-3.0774e-005
-4.1340e-009
2.2204e-016
-8.8818e-016

d =
0.5398

5.2. Résolution de systèmes d’équations non Linéaires

Nous allons appliquer cette méthode à un cas très simple : résolution d’un système de deux équations non linéaires à 2 inconnues, par exemple :

l’algorithme de Newton-Raphson devient :

Exemple :

Les fonctions  et  admettent comme dérivées partielles :

Afin de déterminer si le système admet des solutions sur , on tracera les surfaces représentatives des fonctions  et  dans le but de montrer l’existence éventuelle d’un couple  pour lequel elles s’annulent. L’algorithme ‘pnum2.m’ ci-dessous permet de tracer en 3 dimensions ces deux fonctions sur le même graphe.
 

x=-1:0.1:1;
y=x;
[x,y]=meshgrid(x,y);
f=x.^2+y.^2-x;
g=x.^2-y.^2-y;
figure(1);
mesh(f);
grid on;
hold on;
mesh(g);
title('Courbes f(x,y) et g(x,y)');
xlabel('x'); ylabel('y');zlabel('f(x,y) et g(x,y)');
hold off;

On peut remarquer sur le graphe obtenu (en 3D) que les 2 fonctions  et  s’annulent simultanément en différents points. Grâce aux fonctionnalités de MATLAB, on peut rechercher les solutions éventuelles ou les conditions initiales pour l’algorithme de Newton-Raphson en faisant varier simultanément les variables  et  à l’aide de la fonction ‘meshgrid’. On obtient alors les courbes de niveaux des fonctions  et .
 
 

On peut développer davantage l’algorithme ‘pnum2.m’ pour chercher, par exemple, les conditions initiales et les solutions évidentes.

Les solutions évidentes s’obtiennent par :

>> [i,j]=find(f==0 & g==0)
i =
11
j =
11

>>[x(11),y(11)]
ans =
-1 0

De cette façon, on peut retrouver plusieurs solutions évidentes si le pas d’échantillonnage est assez faible. Avec le pas de 0,1 choisi, on ne retrouve pas la solution évidente (0,0).

Graphiquement, en faisant un zoom autour de l’abscisse ‘x=11’ avec la fonction ‘axis’, on peut cerner cette solution évidente.

>>axis([10.5 11.5 -0.02 0.02]);

Fig. 10 : zoom autour de l’abscisse 

Ainsi, l’algorithme de Newton-Raphson appliqué au système précédent, est représenté par la suite du programme MATLAB (pnum2.m) :

clf;
clear all;
clc;
x=-1:0.1:1;
y=x;
[x,y]=meshgrid(x,y);
f=x.^2+y.^2-x;
g=x.^2-y.^2-y;
figure(1);
mesh(f);
grid on;
hold on;
mesh(g);
title('Courbes f(x,y) et g(x,y)');
xlabel('x'); ylabel('y');zlabel('f(x,y) et g(x,y)');
hold off;
figure(2);
plot(f);
hold on;plot(g);grid on;
title('Intersection de f et g');
xlabel('x');ylabel('y');
axis([0 20 -0.5 0.5]);
gtext('f(x,y)');
gtext('g(x,y)');
hold off;
clear x y;
% Conditions initiales
x(1)=0.85; y(1)=0.35;

% Algorithme de Newton-Raphson
for i=2:10
    diff=inv([2*x(i-1)-1 2*y(i-1);2*x(i-1) -2*y(i-1)-1]);
    f=x(i-1)^2-x(i-1)+y(i-1)^2;
    g=x(i-1)^2-y(i-1)^2-y(i-1);
    xy=[x(i-1) y(i-1)]'-diff*[f,g]';
    x(i)=xy(1);
    y(i)=xy(2);
end

Le choix des conditions initiales influe sur la convergence de l’algorithme. Dans ce qui suit (suite de ‘pnum2.m’), nous étudierons l’évolution des solutions par 2 jeux de conditions initiales.

% Tracé des courbes d'évolution des solutions
figure(3);
plot(1:10,x,'*r',1:10,x);
hold on;
plot(1:10,y,'og',1:10,y);
title('Evolution des solutions : ***-->x et ooo-->y');
grid on;
hold off;

L’évolution des solutions x et y (pour les conditions initiales : et ) est donnée par la figure 11 ci-après. Ces solutions sont les suivantes :

>>x
x =
0.8500 0.7800 0.7719 0.7718 0.7718 0.7718 0.7718 0.7718 0.7718 0.7718

>>y
y =
0.3500 0.4271 0.4197 0.4196 0.4196 0.4196 0.4196 0.4196 0.4196 0.4196

Les 2 solutions peuvent être obtenues au bout de 2 itérations grâce au bon choix des conditions initiales. Avec les conditions initiales x(1)=0,2 et y(1)=0.1, il y a convergence vers la solution (0,0) au bout de 3 itérations.

Pour les conditions initiales x(1)=0,2 et y(1)=0,1, l’évolution des solutions x et y est :

>>x
x =
0.2000 -0.1031 -0.0112 -0.0002 -0.0000 -0.0000 -0.0000 0 0 0

>>y
y =
0.1000 -0.0594 -0.0054 -0.0001 -0.0000 -0.0000 -0.0000 0 0 0

La boite à outil ‘Optimisation Toolbox’ propose les fonctions ‘fsolve’ et ‘fsolve2’ qui permettent respectivement, la recherche des racines d’une fonction dans un intervalle donné et les solutions d’équations non linéaires et de systèmes d’équations non linéaires.