« Découvrir Scilab/Interpolation, extrapolation et lissage » : différence entre les versions
Contenu supprimé Contenu ajouté
ébauche |
import de Calcul numérique |
||
Ligne 5 :
<span style="font-size:25px;">5. Interpolation, extrapolation et lissage</span>
----
== Interpolation dans le plan ==
=== Présentation du problème ===
On dispose de points expérimentaux (''x<sub>i</sub>'', ''y<sub>i</sub>''), qui sont censés suivre une loi ''y'' = ƒ(''x'') inconnue, et l'on désire avoir une valeur ''y'' correspondant à une valeur ''x'' arbitraire, mais contenue dans l'intervalle [min{''x<sub>i</sub>''}, max{''x<sub>i</sub>''}]. Une manière de faire consiste à utiliser une interpolation ; Scilab propose deux démarches :
* interpolation linéaire : si ''x''<sub>''i''</sub> ≤ ''x'' ≤ ''x''<sub>''i'' + 1</sub>, alors on prend ''y'' sur le segment de droite [(''x<sub>i</sub>'', ''y<sub>i</sub>'') ; ((''x''<sub>''i'' + 1</sub>, ''y''<sub>''i'' + 1</sub>)] :<br /><math>y = y_i + \frac{y_{i + 1} - y_i}{x_{i + 1} - x_i}(x - x_i)</math> ;
* interpolation par un polynôme de degré trois, en utilisant des ''splines'' cubiques d'Hermite :<br /> ''y'' = ''p''(''x'')<br />la fonction ''p'' étant un polynôme de degré trois par parties, les portions de polynôme étant raccordées aux points expérimentaux (même valeur, même pente).
{{loupe|wikipedia:fr:Interpolation linéaire|wikipedia:fr:Spline cubique d'Hermite}}
=== Interpolation linéaire ===
Le nuage de ''n'' points (''x<sub>i</sub>'', ''y<sub>i</sub>'')<sub>1 ≤ ''i'' ≤ ''n''</sub> est décrit par deux vecteurs : <code>x</code>, qui contient les ''x<sub>i</sub>'' classés par ordre croissant, et <code>y</code>, qui contient les ''y<sub>i</sub>'' correspondant.
Pour un interpolation linéaire, on considère une valeur <code>xp</code> arbitraire, ou un vecteur de valeurs, les valeurs de <code>yp</code> correspondantes sont obtenues par
<source lang="scilab">
yp = interp1(x, y, xp)
</source>
Plutôt que d'avoir les valeurs interpolées, on peut préférer les valeurs les plus proches ; pour cela, on utilise l'option <code>"nearest"</code> :
<source lang="scilab">
yp = interp1(x, y, xp, "nearest")
</source>
On peut aussi utiliser la commande d'interpolation à ''n'' dimensions <code>linear_interpn()</code> :
<source lang="scilab">
yp = linear_interpn(xp, x, y)
</source>
Si l'on veut faire une ''extrapolation'', c'est-à-dire avoir une valeur hors de l'intervalle [''x''<sub>1</sub>, ''x<sub>n</sub>''], on peut ajouter une option :
* <code>"by_zero"</code> : ''y'' = 0 hors intervalle ;
* <code>"C0"</code> : ''y'' = ''y''<sub>1</sub> si ''x'' < ''x''<sub>1</sub>, ''y'' = ''y<sub>n</sub>'' si ''x'' > ''x<sub>n</sub>'' ;
* <code>"natural"</code> : si la fonction est supposée monotone hors de l'intervalle de définition, on utilise le segment de droite le plus proche ;
* <code>"periodic"</code> : si la fonction est supposée périodique, on répète le motif obtenu.
Par exemple :
<source lang="scilab">
yp = linear_interpn(xp, x, y, "natural")
</source>
=== Interpolation par une ''spline'' cubique ===
Pour une interpolation par ''spline'' cubique d'Hermite, la fonction d'interpolation est décrite par les coordonnées des points et un vecteur de coefficients <code>d</code> (les dérivées), soit trois vecteurs ( <code>x</code>, <code>y</code>, <code>d</code>). Le vecteur <code>d</code> est obtenu par la commande <code>splin()</code> :
<source lang="scilab">
d = splin(x, y)
</source>
si l'on a maintenant une valeur <code>xp</code> arbitraire, ou un vecteur de valeurs, les valeurs de <code>yp</code> correspondantes sont obtenues par
<source lang="scilab">
yp = interp(xp, x, y, d)
</source>
=== Lissage ===
L'interpolation considère que les données sont peu bruitées ; les points mesurés sont « fiables », les valeurs des <code>x</code> et <code>y</code> sont « exactes ». Si le signal est bruité, il faut alors effectuer un lissage.
[[Fichier:Lissage sinus bruite spline cubique.svg|vignette|Lissage par une ''spline'' cubique.]]
Scilab propose la commande <code>lsq_splin()</code> qui utilise des ''splines'' cubique comme modèle, déterminées par la méthode des moindres carrés ''(least square)''. Si l'on a deux vecteurs ''x'' et ''y'' de ''m'' valeurs chacun, on choisit ''n'' abscisses ''xs'' pour déterminer la ''spline'', avec dans l'idéal ''n'' << ''m''. Typiquement :
<source lang="scilab">
// ***** Génération des données *****
m = 200; // nombre de points mesurés
bruit = 0.2; // amplitude du bruit
x = linspace(0, 2*%pi, m);
y = sin(x) + 0.2*rand(x, "normal"); // valeurs bruitées
// ***** Détermination de la spline *****
n = 6; // nombre de points de contrôle de la spline
xs = linspace(x(1), x($), n); // abscisse des points de contrôle
// répartition uniforme
[ys, d] = lsq_splin(x, y, xs); // détermination de la spline
// ***** Calcul des valeurs lissées de y *****
yliss = interp(x, xs, ys, d); // calcul des données lissées
// ***** Tracé *****
plot(x, y) // données originales
plot(x, yliss, "r") // données lissées
</source>
Ainsi, Scilab a construit une ''spline'' cubique passant par les abscisses ''xs'' (en l'occurrence 0, 2π/5, 4π/5, 6π/5, 8π/5 et 2π), et a déterminé les dérivées ''d'' en ces ''n'' points.
Si l'on veut éviter des effets de bord, on peut choisir une répartition différente des points de contrôle. On utilise fréquemment une répartition sinusoïdale pour avoir plus de points aux bords qu'au centre (comme pour les [[w:fr:Polynôme de Tchebychev|polynômes de Tchebychev]]).
<source lang="scilab">
theta = linspace(-%pi, 0, n); // angles uniformément répartis
xs = mean([x(1), x($)]) + 0.5*(x($) - x(1))*cos(theta);
// abscisse des points de contrôle
</source>
Ceci est en particulier important si l'on veut extrapoler, c'est-à-dire calculer des valeurs de ''y'' pour des ''x'' hors des points de mesure.
Cela ne présente d'intérêt que pour ''n'' ≥ 4.
[[Fichier:Savitzky-golay pic gaussien bruite.svg|vignette|upright=2|Lissage et détermination des dérivée et dérivée seconde par l'algorithme de Savitzky-Golay.]]
L'utilisateur peut évidemment coder un autre algorithme. L'image ci-contre montre le calcul de la dérivée et de la dérivée seconde par l'algorithme de Savitzky-Golay ; cliquer sur l'image pour avoir accès à la page de description indiquant le code utilisé.
On peut aussi effectuer un simple lissage par la méthode des moyennes glissantes en utilisant la fonction <code>convol()</code> (voir [[#Produit de convolution|ci-après]]). Si l'on appelle ''p'' le nombre de points de l'intervalle glissant :
<source lang="scilab">
p = 10; // largeur de l'intervalle
h = 1/p*ones(1, p); // fonction porte pour le lissage
yliss = conv(h, y);
plot(x, y) // données originales
plot(x, yliss(1+p/2:$-p/2+1), "r") // données lissées
</source>
== Interpolation d'une surface ==
=== Présentation du problème ===
On dispose de valeurs ''z<sub>i</sub>'' à des positions déterminées (''x<sub>i</sub>'', ''y<sub>i</sub>''), qui sont censés suivre une loi ''z'' = ƒ(''x'', ''y'') inconnue, ce qui définit donc une surface dans l'espace, et l'on désire avoir une valeur ''z'' correspondant à un point (''x'', ''y'') arbitraire.
=== Interpolation linéaire sur un quadrillage ===
Les données sont décrites par :
* un vecteur <code>x</code> de dimension ''m'' et un vecteur <code>y</code> de dimension ''n'', strictement croissants, qui définissent un quadrillage du plan, donc un ensemble de ''m''×''n'' points {(''x''<sub>1</sub>, ''y''<sub>1</sub>) ; (''x''<sub>1</sub>, ''y''<sub>2</sub>) ; … ; (''x''<sub>1</sub>, ''y<sub>n</sub>'') ; (''x''<sub>2</sub>, ''y''<sub>1</sub>) ; … ; (''x<sub>m</sub>'', ''y<sub>n</sub>'') ;
* une matrice <code>z</code> de dimension ''m''×''n'', contenant les valeurs correspondante : ''z''<sub>''i'', ''j''</sub> = ƒ(''x<sub>i</sub>'', ''y<sub>j</sub>'').
L'interpolation linéaire se fait de la même manière que précédemment : on définit un nouveau couple de valeurs <code>xp</code> et <code>yp</code> (ou un couple de vecteurs défnissant un nouveau quadrillage), et l'on a :
<source lang="scilab">
zp = linear_interp(xp, yp, x, y, z)
</source>
c'est-à-dire qu'à partir des points (''x'', ''y'', ''z''), on définit les points interpolés (''xp'', ''yp'', ''zp'').
=== Interpolation cubique sur un quadrillage ===
Comme précédemment, les données sont décrites par :
* un vecteur <code>x</code> de dimension ''m'' et un vecteur <code>y</code> de dimension ''n'', strictement croissants, qui définissent un quadrillage du plan, donc un ensemble de ''m''×''n'' points {(''x''<sub>1</sub>, ''y''<sub>1</sub>) ; (''x''<sub>1</sub>, ''y''<sub>2</sub>) ; … ; (''x''<sub>1</sub>, ''y<sub>n</sub>'') ; (''x''<sub>2</sub>, ''y''<sub>1</sub>) ; … ; (''x<sub>m</sub>'', ''y<sub>n</sub>'') ;
* une matrice <code>z</code> de dimension ''m''×''n'', contenant les valeurs correspondante : ''z''<sub>''i'', ''j''</sub> = ƒ(''x<sub>i</sub>'', ''y<sub>j</sub>'').
Et comme dans le plan, la fonction d'interpolation est décrite par (<code>x</code>, <code>y</code>, <code>C</code>), les valeurs coordonnées des points du quadrillage et une matrice de coefficients <code>C</code> obtenue par
<source lang="scilab">
C = splin2d(x, y, z)
</source>
Si l'on considère un point du plan (<code>xp</code>, <code>yp</code>) (ou un couple de vecteurs défnissant un nouveau quadrillage), les valeurs de ''z'' sur ces nouveaux points sont :
<source lang="scilab">
zp = interp2d(xp, yp, x, y, C)
</source>
=== Interpolation cubique sur un nuage de points dispersés ===
On ne dispose pas toujours de relevés de valeurs sur des coordonnées régulières formant un quadrillage. Dans ce cas-là, on utilise la méthode de Shepard (pondération des plus proches voisins). Les ''n'' points expérimentaux (''x<sub>i</sub>'', ''y<sub>i</sub>'', ''z<sub>i</sub>'') sont décrits par une matrice M de trois colonnes et de ''n'' lignes
: <math>\mathrm{M} = \begin{pmatrix}
x_1 & y_1 & z_1 \\
\vdots & \vdots & \vdots \\
x_n & y_n & z_n
\end{pmatrix}</math>
La fonction d'interpolation est définie par une [[../Structures de données avancées#Structures|liste typée]] T obtenue par
<source lang="scilab">
T = cshep2d(M)
</source>
Si l'on considère un point du plan (<code>xp</code>, <code>yp</code>) (ou un couple de vecteurs défnissant un nouveau quadrillage), les valeurs de ''z'' sur ces nouveaux points sont :
<source lang="scilab">
zp = eval_cshep2d(xp, yp, T)
</source>
== Notes ==
| |||