Découvrir Scilab/Calcul numérique

Un livre de Wikilivres.
Sauter à la navigation Sauter à la recherche

Table des matièresIndex



3. Calcul numérique


Sections

Vocabulaire[modifier | modifier le wikicode]

Dans le vocabulaire habituel de la programmation, une procédure est un sous-programme, c'est-à-dire un bout de programme auquel on donne un nom, et que l'on exécute à la demande en « appelant » ce nom. On peut voir cela comme une macro-définition (ou simplement « macro ») : le nom de la procédure est « remplacé » par le contenu du sous-programme. Une fonction est une procédure qui renvoit une valeur.

Dans la documentation officielle de Scilab, les termes function, procedure et macro sont équivalents ; par contre, ils sont utilisés pour les sous-programme définis par l'utilisateur. Les fonctions intégrées dans Scilab, dites built-in, sont appelées « primitives ».

Lorsqu'une primitive a pour argument une fonction, cette fonction est dite « externe » (external). C'est le cas par exemple des primitives de résolution d'équations différentielle ode() ou d'optimisation optim() : on indique à ces primitives quelle est la fonction à utiliser.

Par soucis de simplicité, nous ne feront pas la différence entre les fonctions définies par l'utilisateur et les primitives, sauf lorsque cette nuance est pertinente. Les primitives sont parfois appelées « commandes » ou « instructions » dans le présent document.

Définition de fonctions extérieures[modifier | modifier le wikicode]

Une fonction extérieure est une fonction numérique définie par l'utilisateur (notez la différence entre « externe » et « extérieure »).

Pour les cas simples, la commande deff permet de définir de nouvelles fonctions pouvant s’exprimer avec les opérateurs déjà prédéfinis (les « fonctions internes » ou « primitives »). On passe deux chaînes de caractères comme paramètre :

  • la première indique le nom de la fonction et les variables utilisées en entrée et en sortie, et
  • la deuxième indique la formule.

Par exemple, la fonction

ƒ(x) = 2·x

peut se définir par

deff("[y] = f(x)", "y = 2*x")

On peut aussi utiliser l'instruction function() sous la forme

-->function [y] = f(x)
-->--> endfunction

la formule, sous la forme y = expression, étant placée à la place des trois points. Avec l'exemple ci-dessus, cela donne :

-->function [y] = f(x)
-->y=2*x
-->endfunction

(voir aussi le chapitre Programmation).

La fonction s'utilise par la suite de la même manière qu'une fonction interne (sur l'exemple ci-dessus) :

-->f(2)
 ans  =

    4.

-->feval(1,f) // calcule f(1)
 ans  =

    2.

Graphe d'une fonction[modifier | modifier le wikicode]

Il est souvent utile de pouvoir afficher des fonctions, ainsi voici un exemple de la syntaxe Scilab pour afficher la représentation graphique y = ƒ(x) d'une fonction ƒ.

t=[1:0.1:30]; // t est un vecteur des valeurs de 1 à 30 par pas de 0,1

deff("[y] = f(x)", "y = cos(x)./x")  // définition de la fonction f

plot(t, f(t)) // tracé de la fonction


Pour plus de détails voir : Découvrir Scilab/Graphiques et sons#Tracé de fonction.

Statistiques[modifier | modifier le wikicode]

Pour l'analyse statistique de données, on peut se reporter au document

(en) Anna Bassi et Giovanni Borzi, « Data Analysis and Statistics: a Scilab data mining tutorial » sur Openeering.

Tendance centrale[modifier | modifier le wikicode]

Considérons une matrice X. On peut calculer les moyennes suivantes :

  • mean(X) : calcule la moyenne de tous les éléments de X ;
  • mean(X, "r") : donne un vecteur ligne (row) contenant les moyennes colonne par colonne ;
  • mean(X, "c") : donne un vecteur colonne contenant les moyennes ligne par ligne.

On peut calculer la moyenne géométrique avec geomean(), et la moyenne harmonique avec harmean(), avec la même syntaxe.

De la même manière, on peut calculer les médianes avec median(X), median(X, "r") et median(X, "c").

La commande trimmean(X, p) permet de calculer la moyenne en éliminant les p pourcents extrêmes (les p/2 % valeurs les plus élevées, et les p/2 % valeurs les plus faibles). On peut également faire le calcul par ligne ou par colonne — trimmean(X, p, "r") ou trimmean(X, p, "c"). notons que l'on a :

  • mean(X) = trimmean(X, 0) ;
  • median(X) = trimmean(X, 100).

Dispersion[modifier | modifier le wikicode]

La commande

strange(X)

calcule l'écart (range) entre la valeur maximale et la valeur minimale de la matrice. On peut calculer les valeurs pour chaque colonne ou chaque ligne avec strange(X, "r") et strange(X, "c"). La commande

mad(X)

calcule la déviation absolue moyenne (mean absolute deviation)

on peut de même la calculer pour chaque colonne ou ligne.

La commande

iqr(X)

calcule l'écart interquartile (interquartile range).

La commande stdev calcule l'écart type (standard deviation) et la commande variance la variance des éléments de la matrice (entière, par ligne ou par colonne). Par défaut, la variance est sans biais : normalisée par m×n - 1 pour une matrice de m×n éléments, par m - 1 par colonne et par n - 1 par ligne. Si l'on calcule la variance par ligne ou par colonne, on peut ajouter l'option 1 pour calculer le moment d'ordre 2 (normalisation par m ou n), par exemple

mom2 = variance(X, "r", 1) // pour une matrice non-vecteur
mom2 = variance(X, "*", 1) // pour un vecteur

L'écart type est toujours sans biais. Si l'on veut l'estimateur normalisé par m ou n, il faut utiliser

estim_ecart_type = sqrt(variance(X, "*", 1)) // pour une matrice non-vecteur

Normalisation[modifier | modifier le wikicode]

La commande center(X) permet de centrer la matrice (équivalent de X - mean(X)). La commande wcenter(X) donne la matrice centrée réduite (équivalent de 1/stdev(X)*(X - mean(X))).

Fréquences et quartiles[modifier | modifier le wikicode]

La commande tabul(X) calcule les fréquences d'apparition des valeurs dans la matrice X. Le résultat est classé par défaut par ordre décroissant des valeurs ; pour avoir un classement croissant, on ajoute le paramètre "i" (increasing order) :

tabul(X, "i")

La commande quart(X) calcule les quartiles et les renvoie dans un vecteur de dimension 3. On peut calculer les quartiles pour chaque colonne ou chaque ligne, avec quart(X, "r") et quart(X, "c").

On peut aussi travailler avec des valeurs et leur fréquence. Si la matrice X contient les valeurs et la matrice F leur fréquence — F(i, j) est la fréquence de X(i, j) —, alors

  • meanf(X, F) donne la moyenne ;
  • stdevf(X, F) donne l'écart type ;
  • variancef(X, F) donne la variance.

Notons que F n'est pas nécessairement normalisé à 1 ou à 100 ; il peut par exemple s'agir d'un nombre d'occurrences.

Exemple

X = int(10*rand(100,1, "normal")); // génère 100 valeurs entières avec une loi normale centrée d'écart type 10
F = tabul(X); // calcule les fréquences d'apparition
x = F(:, 1); // extraction de la liste des valeurs
f = F(:, 2); // extraction de la liste des fréquences

quartiles = quart(X); // calcule les quartiles de la population

moyenne = meanf(x, f); // calcul de la moyenne pondérée par les fréquences

Considérons maintenant deux vecteurs X et Y de même dimension, et la matrice F des fréquences de X et Y : F(i, j) est la fréquence de X(i) & Y(j). Alors :

  • covar(X, Y, F) donne la covariance ;
  • correl(X, Y, F) donne le coefficient de corrélation.

Exemple

n = 10; // nombre de valeurs
X = 1:n; // abscisses
Y = 2*X + 5 + 3*rand(X); // ordonnées : droite "bruitée"
plot(X, Y); // tracé des données
F = eye(n, n); // couplage (x, y) : chaque couple est présent une fois, (X(i) & Y(i)) = 1, (X(i) & Y(j)) = 0

covariance1 = covar(X, Y, F);
covariance2 = 1/n*sum(X.*Y) - 1/n^2*sum(X)*sum(Y); // calcul "à la main" pour comparer

correlation1 = correl(X, Y, F);
correlation2 = covariance2/sqrt(variance(X, "*", 1))/sqrt(variance(Y, "*", 1)); // calcul "à la main"

Fonctions de distribution et de répartition de lois statistiques[modifier | modifier le wikicode]

Par défaut, Scilab propose une loi de distribution, pour la loi binomiale, et onze lois de répartition (probabilités cumulatives).

Pour la loi binomiale, la fonction de répartition pour n épreuves (nombre entier) avec une probabilité p (nombre réel entre 0 et 1) est donnée par

f = binomial(p, n)

La variable f est un vecteur de n + 1 composantes, avec

La fonction de répartition s'obtient avec cdfbin() :

P =  = cdfbin("PQ", k, n, p, 1-p)

donne la probabilité P d'avoir entre 0 et k succès pour n tirages avec une probabilité de réussite p. On peut d'ailleurs le comparer avec

f = binomial(p, n)
c = sum(f(1:k+1))

De manière générale, les fonctions de répartition commencent par cdf… (cumulative distribution function), par exemple :

  • P = cdfchi("PQ", x, k) : loi du χ2 à k degrés de liberté ;
  • P = cdfnor("PQ", x, mu, sigma) : loi normale d'espérance mu et d'écart type sigma ;
  • P = cdfpoi("PQ", k, lambda) : loi de Poisson de paramètre lambda ;
  • P = cdft("PQ", t, k) : loi t de Student à k degrés de liberté ;
  • … voir Scilab help >> Statistics > Cumulated Distribution Functions

Notons que pour chaque cas, on peut aussi avoir la valeur complémentaire Q = 1 - P avec la syntaxe

[P, Q] = cdf…("PQ", …)

Exemple

Ci-dessous, nous traçons la fonction de répartition de la loi normale centrée réduite (μ = 0, σ = 1). Notons que pour des calculs vectoriels, l'espérance et l'écart type sont des vecteurs de même dimension que X.

X = -3:0.1:3;
mu = zeros(X);
sigma = ones(X);
[P, Q] = cdfnor("PQ", X, mu, sigma);
plot(X', [P', Q'])

Les mêmes fonctions permettent de rechercher les quantiles ou fractiles :

  • k = cdfbin("S", n, P, 1-P, p, 1-p) : nombre de succès k sur n épreuves ayant pour paramètre p, donnant une probabilité cumulée P ;
  • x = cdfchi("X", k, P, 1-P) : valeur de x pour laquelle on a une probabilité cumulée P avec une loi du χ2 à k degrés de liberté ;
  • x = cdfnor("X", mu, sigma, P, 1-P) : valeur de x pour laquelle on a une probabilité cumulée P avec une loi normale d'espérance mu et d'écart type sigma ;
  • k = cdfpoi("S", lambda, P, 1-P) : valeur de k pour laquelle on a une probabilité cumulée P avec une loi de Poisson de paramètre lambda ;
  • t = cdft("T", k, P, 1-P) : valeur de t pour laquelle on a une probabilité cumulée P avec une loi t de Student à k degrés de liberté ;

Exemple

Ci dessous, nous déterminons le quantile pour un risque α bilatéral de 5 %, c'est-à-dire que l'on rejette 0,25 % des cas les plus hauts et les plus bas, pour une loi normale centrée réduite puis pour une loi t de Student à 30 degrés de liberté

-->x = cdfnor("X", 0, 1, 0.975, 0.025)
 x  =
 
    1.959964
 
-->t = cdft("T", 30, 0.975, 0.025)
 t  =
 
    2.0422725

Enfin, ces fonctions permettent de déterminer les paramètres des loi pour atteindre un quantile visé :

  • loi binomiale :
    • n = cdfbin("Xn", p, 1-p, P, 1-P, k) : nombre d'épreuves ayant pour paramètre p afin d'avoir un quantile P valant k,
    • [p, q] = cdfbin("PrOmpr", P, 1-P, k, n) : paramètre p de chaque épreuve afin d'avoir un quantile P pour k succès parmi n ;
  • loi du χ2 :
    k = cdfchi("Df", P, 1-P, x) : valeur du nombre de degrés de liberté (degrees of freedom) k pour avoir un quantile P valant x ;
  • loi normale :
    • mu = cdfnor("Mean", sigma, P, 1-P, x) : valeur de l'espérance (moyenne) pour avoir un quantile P valant x si l'on a un écart type sigma,
    • sigma = cdfnor("Std", P, 1-P, X, mu) : valeur de l'écart type (standard deviation) pour avoir un quantile P valant x si l'on a une espérance de mu ;
  • loi de Poisson :
    lambda = cdfpoi("Xlam", P, 1-P, k) : valeur du paramètre λ pour avoir un quantile P valant k ;
  • loi t de Student :
    k = cdft("Df", P, 1-P, t) : nombre de degrés de liberté k pour avoir un quantile P valant t.

Générateurs pseudo-aléatoires[modifier | modifier le wikicode]

Les générateurs de nombres aléatoires permettent de tester des programmes ou bien de faire des simulations par la méthode de Monte Carlo.

Scilab possède une fonction rand() qui par défaut génère un nombre entre 0 et 1 avec une loi de distribution uniforme. Par exemple, pour simuler un jet de dé à six faces, on peut utiliser :

de = 1 + int(6*rand())

La fonction peut générer une matrice m×n avec les syntaxes

rand(m, n)
rand(X)

où X est une matrice m×n. Notons que si l'on veut générer n valeurs, il faut bien taper rand(n, 1) ou rand(1, n).

On peut lui faire générer une valeur selon une loi normale centrée réduite avec la clef "normal" :

rand(m, n, "normal")

Exemple

Supposons par exemple que l'on a une poutre faite dans un acier dont la résistance intrinsèque, la limite d'élasticité Re, est garantie au minimum à 235 MPa. En réalité, Re suit une loi normale d'espérance μRe = 269 MPa et d'écart type σRe = 17 MPa. Dans les conditions d'utilisation, l'acier de cette poutre doit résister à une contrainte S (stress) de 214 MPa ; en fait, les conditions de chargement sont variables et suivent une loi normale de paramètres μS = 214 MPa et σS = 15 MPa.

On veut déterminer statistiquement dans combien de cas la contrainte S est supérieure à la résistance Re.

n = 1e6; // nombre d'événements
R = 269 + 17*rand(n, 1, "normal"); // résistance
S = 214 + 14*rand(n, 1, "normal"); // contrainte
Mbool = S>R; // comparaison des valeurs (matrice de booléens)
moyenne = sum(Mbool)/n; // les "vrai" comptent pour 1 et les "faux" pour 0
disp(string(100*moyenne)+" %") // affichage du résultat

La commande grand() (get random number) permet d'avoir des nombres aléatoires selon de nombreuses lois, par exemple :

  • Y = grand(m, n, "bin", N, p) : matrice m×n de nombres aléatoires tirés selon une loi binomiale de N essais avec une probabilité p chacun ;
  • Y = grand(m, n, "nor", mu, sigma) : pour une loi normale d'espérance mu et d'écart type sigma ;
  • Y = grand(m, n, "unf", a, b) : pour une loi uniforme sur [a' ; b[ (b exclus) ;
  • Y = grand(m, n, "uin", a, b) : pour une loi uniforme sur [a' ; b] (b inclus) ;

Dans l'exemple précédent, on aurait donc :

R = grand(n, 1, "nor", 369, 17); // résistance
S = grand(n, 1, "nor", 314, 14); // contrainte

Autres lois statistiques[modifier | modifier le wikicode]

Des modules complémentaires Atoms permettent d'avoir accès à d'autres loi.

Par exemple le module CASCI donne accès à de nombreuses fonctions de distribution et de répartition, notamment des lois log-normale (pdflognormal(), cdflognormal()), exponentielle (pdfexponential(), cdfexponential()) ou de Weibull (pdfweibull(), cdfweibull()). Ainsi :

p = pdfnormal(x, mu, sigma); // densité d'une loi normale en x
y = cdflognormal(x, lambda); // répartition d'une loi log-normale d'espérance mu et d'écart type sigma
y = cdfexponential(x, mu, sigma); // répartition d'une loi exponentielle de paramètre lambda
y = cdfweibull(x, k, lambda, theta); // répartition d'une loi de Weibull
   // de paramètre de forme k, de paramètre d'échelle lambda et de paramètre de position thêta

Le module fournit également le calcul des quantiles avec les fonctions idf… (inverse cumulative distribution functions), et des générateurs de nombres aléatoires avec les fonctions rnd…. Par exemple :

x = idfweibull(y, k, lambda, theta); // quantile y d'une loi de Weibull
   // de paramètre de forme k, de paramètre d'échelle lambda et de paramètre de position thêta
X = rndweibull(n, k, lambda, theta); // génère n nombres aléatoires selon la loi de Weibull

Le module Stixbox définit également un certan nombre de fonctions statistiques utiles.

Interpolation, extrapolation et lissage[modifier | modifier le wikicode]

Interpolation dans le plan[modifier | modifier le wikicode]

Présentation du problème[modifier | modifier le wikicode]

On dispose de points expérimentaux (xi, yi), 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{xi}, max{xi}]. Une manière de faire consiste à utiliser une interpolation ; Scilab propose deux démarches :

  • interpolation linéaire : si xixxi + 1, alors on prend y sur le segment de droite [(xi, yi) ; ((xi + 1, yi + 1)] :
     ;
  • interpolation par un polynôme de degré trois, en utilisant des splines cubiques d'Hermite :
    y = p(x)
    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).


Pour plus de détails voir : wikipedia:fr:Interpolation linéaire et wikipedia:fr:Spline cubique d'Hermite.

Interpolation linéaire[modifier | modifier le wikicode]

Le nuage de n points (xi, yi)1 ≤ in est décrit par deux vecteurs : x, qui contient les xi classés par ordre croissant, et y, qui contient les yi correspondant.

Pour un interpolation linéaire, on considère une valeur xp arbitraire, ou un vecteur de valeurs, les valeurs de yp correspondantes sont obtenues par

yp = interp1(x, y, xp)

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 "nearest" :

yp = interp1(x, y, xp, "nearest")

On peut aussi utiliser la commande d'interpolation à n dimensions linear_interpn() :

yp = linear_interpn(xp, x, y)

Si l'on veut faire une extrapolation, c'est-à-dire avoir une valeur hors de l'intervalle [x1, xn], on peut ajouter une option :

  • "by_zero" : y = 0 hors intervalle ;
  • "C0" : y = y1 si x < x1, y = yn si x > xn ;
  • "natural" : si la fonction est supposée monotone hors de l'intervalle de définition, on utilise le segment de droite le plus proche ;
  • "periodic" : si la fonction est supposée périodique, on répète le motif obtenu.

Par exemple :

yp = linear_interpn(xp, x, y, "natural")

Interpolation par une spline cubique[modifier | modifier le wikicode]

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 d (les dérivées), soit trois vecteurs ( x, y, d). Le vecteur d est obtenu par la commande splin() :

d = splin(x, y)

si l'on a maintenant une valeur xp arbitraire, ou un vecteur de valeurs, les valeurs de yp correspondantes sont obtenues par

yp = interp(xp, x, y, d)

Lissage[modifier | modifier le wikicode]

L'interpolation considère que les données sont peu bruitées ; les points mesurés sont « fiables », les valeurs des x et y sont « exactes ». Si le signal est bruité, il faut alors effectuer un lissage.

Lissage par une spline cubique.

Scilab propose la commande lsq_splin() 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 :

// ***** 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

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 polynômes de Tchebychev).

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

Ceci est en particulier important si l'on veut extrapoler, c'est-à-dire calculer des valeurs de y pour des x hors des points d emesure.

Cela ne présente d'intérêt que pour n ≥ 4.

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 convol() (voir ci-après). Si l'on appelle p le nombre de points de l'intervalle glissant :

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

Interpolation d'une surface[modifier | modifier le wikicode]

Présentation du problème[modifier | modifier le wikicode]

On dispose de valeurs zi à des positions déterminées (xi, yi), 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[modifier | modifier le wikicode]

Les données sont décrites par :

  • un vecteur x de dimension m et un vecteur y de dimension n, strictement croissants, qui définissent un quadrillage du plan, donc un ensemble de m×n points {(x1, y1) ; (x1, y2) ; … ; (x1, yn) ; (x2, y1) ; … ; (xm, yn) ;
  • une matrice z de dimension m×n, contenant les valeurs correspondante : zi, j = ƒ(xi, yj).

L'interpolation linéaire se fait de la même manière que précédemment : on définit un nouveau couple de valeurs xp et yp (ou un couple de vecteurs défnissant un nouveau quadrillage), et l'on a :

zp = linear_interp(xp, yp, x, y, z)

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[modifier | modifier le wikicode]

Comme précédemment, les données sont décrites par :

  • un vecteur x de dimension m et un vecteur y de dimension n, strictement croissants, qui définissent un quadrillage du plan, donc un ensemble de m×n points {(x1, y1) ; (x1, y2) ; … ; (x1, yn) ; (x2, y1) ; … ; (xm, yn) ;
  • une matrice z de dimension m×n, contenant les valeurs correspondante : zi, j = ƒ(xi, yj).

Et comme dans le plan, la fonction d'interpolation est décrite par (x, y, C), les valeurs coordonnées des points du quadrillage et une matrice de coefficients C obtenue par

C = splin2d(x, y, z)

Si l'on considère un point du plan (xp, yp) (ou un couple de vecteurs défnissant un nouveau quadrillage), les valeurs de z sur ces nouveaux points sont :

zp = interp2d(xp, yp, x, y, C)

Interpolation cubique sur un nuage de points dispersés[modifier | modifier le wikicode]

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 (xi, yi, zi) sont décrits par une matrice M de trois colonnes et de n lignes

La fonction d'interpolation est définie par une liste typée T obtenue par

T = cshep2d(M)

Si l'on considère un point du plan (xp, yp) (ou un couple de vecteurs défnissant un nouveau quadrillage), les valeurs de z sur ces nouveaux points sont :

zp = eval_cshep2d(xp, yp, T)

Régression[modifier | modifier le wikicode]

Présentation du problème

Le problème a des similarités avec celui de l'interpolation : on a des points expérimentaux (xi, yi), et l'on détermine une fonction décrivant ces points. Les différences majeures sont :

  • la fonction est déterminée sur l'ensemble de définition complet, à partir de tous les points, et non pas par morceaux entre deux points ;
  • la fonctions ne passe en général pas par les points expérimentaux mais près de ces points.

Par ailleurs, la forme de la fonction est souvent déterminée par le type de problème auquel on s'attaque, elle a un « sens physique », là où l'interpolation se contente de prendre des fonctions simples.

Considérons un nuage de n points « expérimentaux » (X, Y) : X et Y sont des vecteurs

  • X = [X1 ; X2 ; … ; Xn] ;
  • Y = [Y1 ; Y2 ; … ; Yn].

Considérons une fonction « modèle » y = ƒ(A, x) où A est un vecteur ; les valeurs de A sont des paramètres de la fonction ƒ, par exemple :

  • ƒ(A, x) = A(1)·exp((x - A(2))2/A(3)) : une fonction gaussienne ;
  • ƒ(A, x) = A(1)·x3 + A(2)·x2 + A(3)·x + A(4) : un polynôme de degré 3.

La régression consiste à touver le vecteur A tel que la fonction ƒ « colle le mieux » au nuage de points. On utilise pour cela la méthode des moindres carrés (least square) :

  • on définit la fonction « résidus » r(A, x, y) entre le modèle et les points expérimentaux :
    r(A, x, y) = y - ƒ(A, x) ;
  • on définit la somme S des carrés des résidus comme étant :
    S = ∑i r2(A, Xi, Yi) ;

la régression par la méthode des moindres carrés consiste donc à trouver le vecteur A donnant la valeur minimale de S.

Régression linéaire[modifier | modifier le wikicode]

Régression linéaire simple[modifier | modifier le wikicode]

Dans le cas de la régression linéaire, la fonction modèle ƒ est une fonction affine :

ƒ(A, x) = A(1)·x + A(2).

Comme il n'y a que deux paramètres dans le modèle, on note ceux-ci a et b :

ƒ(a, b, x) = a·x + b.

D'une certaine manière, la régression consiste à résoudre un système d'équations surdéterminé (chaque point expérimental étant une équation, il y a plus d'équations que d'inconnues) :

ce qui peut s'écrire sous forme matricielle

A×X = Y

avec

A = (a b)
Y = (y1 y2yn )

Si l'on demande à Scilab de résoudre le système avec la division matricielle alors que la matrice X n'est pas carrée, il effectue automatiquement une régression linéaire des moindres carrés :

\\ génération des données de test
x = 1:10
y = x + 0.3*rand(x);

\\ régression
X = [x ; ones(x)];
A = y/X

Si l'on veut forcer l'ordonnée à l'origine à 0, il suffit de ne pas mettre la ligne de 1 :

A = y/x

On peut aussi travailler avec des vecteurs colonne, et utiliser l'autre diviseur matriciel A = x\y.

Mais la méthode de la division matricielle ne donne pas d'autre information que les coefficients, et en particulier ne donne pas le résidus.

Pour cela, on peut utiliser la fonction reglin(), avec la syntaxe :

// définition du nuage de points
X = [];
Y = [];

// régression
[a, b, sigma] = reglin(X, Y);
 les matrices X et Y doivent être des vecteurs ligne.

On peut aussi utiliser la commande regress(). Elle ne permet pas d'obtenir l'écart type, et ne fonctionne qu'avec deux variables, mais par contre, les vecteurs peuvent être indifféremment ligne ou colonne ; ils peuvent même être de types différents (un vecteur ligne et l'autre colonne).

coefs = regress(X, Y)

Exemple

Dans l'exemple suivant, on génère des point dispersés aléatoirement autour d'une droite y = pente*x + ordonnee, puis on fait la régression linéaire de ce nuage de points.

Pour clarifier le code source, le programme est scindé en trois fichiers situés dans le même répertoire (dossier) appelé ici monrepertoire/ :

  • le fichier fonction.sce qui contient la fonction affine ;
  • le fichier generation.sce qui crée un nuage de points et le met dans le fichier donnees.txt ;
  • le fichier regression.sce qui lit le fichier donnees.txt et effectue la régression.
Fichier fonction.sce
chdir("monrepertoire/");

// **********
// Constantes
// **********

// **********
// Fonctions
// **********

deff("[y] = affine(a, b, x)", "y = a*x + b")
Fichier generation.sce
chdir("monrepertoire/");

// **********
// Constantes
// **********

// paramètres de la loi affine
pente = 1;
ordonnee = 10;

// paramètres de la loi normale
var = 1; // variance

// nombre de points
nbpts = 20;

// **********
// Initialisation
// **********

X = zeros(nbpts)';
Y = zeros(nbpts)';

// **********
// Fonctions
// **********

exec("fonction.sce", -1) // charge le fichier (-1 : sans commentaire)

// **********
// Programme principal
// **********

// nuage de points
for i = 1:nbpts
    X(1,i) = i + var*rand(1, "normal");
    Y(1,i) = affine(pente, ordonnee, i) + var*rand(1, "normal");
end

// enregistrement

write("donnees.txt", [X,Y]);
Fichier regression.sce
clf;
chdir("monrepertoire/");

// **********
// Programme principal
// **********

// lecture des données
donnees = read("donnees.txt", -1, 2)
Xexp = donnees(:,1);
Yexp = donnees(:,2);


// regression
[aopt, bopt, sigma] = reglin(Xexp, Yexp)

// loi modèle
Yopt = affine(aopt, bopt, X);

// affichage
plot(Xexp, Yexp, "+b")
plot(Xexp, Yopt, "-r")
xstring(0, ordonnee, "a = "+string(aopt)+" pour "+string(pente)...
+" ; b = "+string(bopt)+" pour "+string(ordonnee)...
+" ; sigma = "+string(sigma))

Comme la forme de la fonction ƒ est connue, il est inutile de la définir. La fonction reglin() retourne un vecteur de trois composantes : les paramètres a et b du modèle, et l'écart type σ du résidu (qui est la racine carrée de la somme des carrés des résidus, σ = √S).

Régression multilinéaire[modifier | modifier le wikicode]

Notons que a, b, x et y peuvent être des vecteurs. Si par exemple on a un modèle linéaire à plusieurs variables :

ƒ(a1, a2, …, am, b, x1, x2, …xn) = a1·x1 + a2·x2 + … + an·xn + b

alors il suffit de définir X comme étant une matrice m×n

X = [X11, X12, … X1n ;
     X21, X22, … X2n ;
     …
     Xm1, Xm2, … Xmn]

où Xij est la valeur de xj pour le point de mesure i ; et de définir Y comme étant un vecteur ligne de n valeurs

Y = [Y1, Y2, … , Yn]

et le paramètre a retourné par la division matricielle ou par la fonction reglin() est le vecteur [a1, a2, … , an].

Régression polynomiale[modifier | modifier le wikicode]

La régression polynomiale est un cas particulier de régression multilinéaire.


Pour plus de détails voir : wikipedia:fr:Régression polynomiale.

Voici par exemple un programme générant n points pour x allant de 0 à 5 (valeurs exactes), et y suivant le polynôme x3 - 2·x2 + 3·x - 4, bruité :

// n points aléatoires approchant un  polynôme de degré 3

// chdir("mon_chemin_d_acces"); // répertoire où enregistrer le fichier

n = 10;
x = linspace(0, 5, n);
P = %s^3 - 2*%s^2 - 3*%s + 4;

y = horner(P, x); // valeurs exactes du polynôme

clf;
plot(x, y); // tracé du polynôme

y = y + 3*rand(y, "normal"); // valeurs bruitées

plot(x, y, "+"); // tracé des points bruités

write("polynome_bruite.txt", [x' , y']);

Et voici un programme permettant d'effectuer la régression sur ces données :

// Régression polynomiale de degré 3

// chdir("mon_chemin_d_acces"); // répertoire où chercher le fichier


M = read("polynome_bruite.txt", -1, 2); // lecture des données

x = M(:, 1);
y = M(:, 2);

X = [x.^3, x.^2, x, ones(x)];

A = X\y; // régression linéaire par les moindres carrés

P = [%s^3, %s^2, %s, 1]*A; // définition du polynôme

yth = horner(P, x); // valeurs correspondant au polynôme

// tracé du résultat
clf;
plot(x, yth)
plot(x, y, "+")

Régression non-linéaire[modifier | modifier le wikicode]

Dans le cas de la régression non-linéaire, ƒ est une fonction quelconque. Pour effectuer la régression, on utilise la fonction leastsq() avec la syntaxe suivante :

// fonction modèle
function y = f(A, x)
   ...
endfunction 

// fonction résidus
function e = r(A, x, y)
   e = f(A, x) - y
endfunction

// initialisation des paramètres du modèle
Ainit = [];

// définition du nuage de points
X = [];
Y = [];

[S, Aopt] = leastsq(list(r, X, Y), Ainit)

Avec cette syntaxe, la fonction leastsq() retourne donc un vecteur :

  • le premier élément, S, est la valeur de la somme des carrés des écarts ;
  • le second élément, Aopt, est le vecteur contenant la valeur « optimisée » des paramètres du modèle.
 les matrices X et Y doivent être ici des vecteurs colonne. Par ailleurs, la fonction de calcul du résidus doit s'appliquer aux vecteurs X et Y, pour renvoyer un vecteur e.

Exemple 1

Nous reprenons l'exemple de la régression linéaire, mais en utilisant cette fois-ci la fonction leastsq(). Lorsque l'on enregistre les données générées, on utilise :

write("donnees.txt", [X', Y']);

afin d'avoir des vecteurs colonne.

clf;
chdir("monrepertoire/");

// **********
// Initialisation
// **********

Ainit=[1;1];

// **********
// Fonctions
// **********

deff("[y] = f(a, b, x)", "y = a*x + b")

deff("[e] = res(A, x, y)", "e = f(A(1),A(2),x) - y")

// **********
// Programme principal
// **********
 
// lecture des données
donnees = read("donnees.txt", -1, 2)
Xexp = donnees(:, 1);
Yexp = donnees(:, 2);
 
// regression
[S, Aopt] = leastsq(list(res, Xexp, Yexp),Ainit)
 
// loi modèle
Ymod = f(Aopt(1), Aopt(2), Xexp);
 
// affichage
plot(Xexp, Ymod, "-r")
plot(Xexp, Yexp, "+b")
xstring(0, ordonnee, "a = "+string(Aopt(1))+" pour "+string(pente)...
+" ; b = "+string(Aopt(2))+" pour "+string(ordonnee)...
+" ; variance = "+string(S))
Résultat de la régression.

Exemple 2

Nous travaillons ici avec un pic généré par la somme d'une fonction gaussienne et d'une fonction lorentzienne. Ce genre de profil est classiquement utilisé pour décrire la forme des pics obtenus par diffractométrie X. On fait ce que l'on appelle une désommation de pics (on parle souvent de « déconvolution », mais ce terme est abusif puisqu'il ne s'agit pas de produit de convolution).

clf

// **********
// Constantes
// **********

paramlorentz(1) = 5; // hauteur de la courbe lorentzienne
paramlorentz(2) = 2; // largeur de la courbe lorentzienne
paramgauss(1) = 10; // hauteur de la courbe gaussienne
paramgauss(2) = 3; // largeur de la courbe gaussienne
var=0.02; // variance de la loi normale
nbpts = 100 // nombre de points
demielargeur = max(3*paramgauss(2), 3*paramlorentz(2)) // pour intervalle x
pas = 2*demielargeur/nbpts;
Ainit = [1;1;1;1]; // paramètres initiaux du modèle

// **********
// Initialisation
// **********

X = zeros(nbpts,1);
Y = zeros(nbpts,1);
Yopt = zeros(nbpts,1);

// **********
// Fonctions
// **********

// Fonction gaussienne centrée sur 0

function [y] = gauss(A, x)
    // A(1) : hauteur de pic
    // A(2) : "largeur" du pic
    y = A(1)*exp(-x^2/A(2));
endfunction

// Fonction lorentzienne centrée sur 0

function [y] = lorentz(A, x)
    // A(1) : hauteur de pic
    // A(2) : "largeur" du pic
    y = A(1)./(1 + (2.*x/A(2))^2);
endfunction

function [y] = profil(A, x)
    // A(1) : hauteur de pic lorentzien
    // A(2) : "largeur" du pic lorentzien
    // A(3) : hauteur de pic gaussien
    // A(4) : "largeur" du pic gaussien
    L = A(1:2);
    G = A(3:4);
    y = lorentz(L, x) + gauss(G, x);
endfunction

// Résidus

function [e] = residus(A, x, y)
    e = profil(A, x) - y;

endfunction

// **********
// Nuage de points
// **********

for i=1:nbpts
    x = pas*i - demielargeur;
    X(i) = x;
    Y(i) = profil([paramlorentz;paramgauss], x) + rand(var, "normal");
end

// **********
// Programme principal
// **********

[S, Aopt] = leastsq(list(residus, X, Y), Ainit)
Yopt = profil(Aopt, X);
YLopt = lorentz(Aopt(1:2),X);
YGopt = gauss(Aopt(3:4),X);

// affichage

plot(X, Yopt, "-r")
plot(X, YLopt, "-c")
plot(X, YGopt, "-g")
plot(X, Y, "+b")

hauteur = paramlorentz(1)+paramgauss(1);

xstring(-demielargeur, hauteur,...
 "lorentzienne : Al = "+string(Aopt(1))+" pour "+string(paramlorentz(1))...
+" ; Hl = "+string(Aopt(2))+" pour "+string(paramlorentz(2)))
xstring(-demielargeur, 3*hauteur/4,...
 "gaussienne : Ag = "+string(Aopt(3))+" pour "+string(paramgauss(1))...
+" ; Hg = "+string(Aopt(4))+" pour "+string(paramgauss(2)))
xstring(-demielargeur, hauteur/2,...
 "variance : S = "+string(S))

L'algorithme par défaut de leastsq() est l'algorithme de pseudo-Newton. On peut aussi utiliser l'algorithme du gradient conjugué en ajoutant "gc" aux paramètres, ou bien algorithme foncitonnant avec les fonctions non-différentiables avec "nd". On peut enfin utiliser la commande lsqrsolve(), qui utilise l'algorithme de Levenberg-Marquardt.

Dérivation[modifier | modifier le wikicode]

La dérivation peut s'aborder de différentes manières selon les données disponibles.

Pour une fonction définie[modifier | modifier le wikicode]

S'il s'agit d'un polynôme ou d'une fraction rationnelle, la commande derivat() effectue la dérivée formelle, par exemple

-->q = (1 + %s^2)/(1+%s)
 q  =
 
         2  
    1 + s   
    -----   
    1 + s   
 
-->derivat(q)
 ans  =
 
              2  
  - 1 + 2s + s   
    ----------   
              2  
    1 + 2s + s

Si l'on dispose d'une fonction extérieure (définie par function() … endfunction ou bien par deff()), alors on peut utiliser la fonction numderivative()[1] pour obtenir l'approximation de la matrice jacobienne et hessienne :

-->deff("y = f(x)", "y = x.^2")

-->numderivative(f, 2) // f'(x) = 2*x ; J = f'(2) = 4
 ans  =
 
    4.

-->[J, H] = numderivative(f, 2) // f''(x) = 2 ; H = f''(2) = 2
 H  =
 
    2.  
 J  =
 
    4.  

-->x=(1:3)';

-->numderivative(f, x)
 ans  =
 
    2.    0.    0.   
    0.    4.    0.  
    0.    0.    6.   

-->[J, H] = numderivative(f, x)
 H  =
 
    2.    0.    0.    0.    0.    0.    0.    0.    0.  
    0.    0.    0.    0.    2.    0.    0.    0.    0.  
    0.    0.    0.    0.    0.    0.    0.    0.    2.  
 J  =
 
    2.    0.    0.  
    0.    4.    0.  
    0.    0.    6.

Pour une liste de valeurs[modifier | modifier le wikicode]

Si l'on dispose d'une liste de valeurs, typiquement un vecteur x et un vecteur y de même dimension, il faut alors distinguer plusieurs cas.

Valeurs peu bruitées[modifier | modifier le wikicode]

Si l'on suppose que le bruit est faible, donc que les valeurs sont « exactes », et que de plus la dérivée seconde est faible, alors on peut se contenter de déterminer les pentes entre les points, c'est-à-dire de calculer le taux de variation, faire une interpolation linéaire. Par exemple :

x = [0 2 5 20 40 60 80];   
y = [0.0956820 0.0480457 0.0277857 0.0036214 0.0002543 0.0002543 0.0001265];

yprime = diff(y)./diff(x);

subplot(2, 1, 1)
plot(x, y)

subplot(2, 1, 2)
plot(x(1:$-1), yprime)

Si la dérivée seconde n'est pas négligeable, on peut alors utiliser une interpolation par un polynôme de degré trois (spline cubique) :

yprime1 = splin(x, y);

subplot(2, 1, 2)
plot(x, yprime1, "r")

Valeurs bruitées : régression et lissage[modifier | modifier le wikicode]

Les méthodes ci-dessus sont très sensibles au bruit. Si les valeurs sont bruitées, alors il faut lisser la courbe d'une manière ou d'une autre.

Si l'on dispose d'un modèle analytique pour les données, on peut déterminer les paramètre de ce modèle par régression, voir la section ci-dessus. Sinon, l'utilisation d'une spline par la commande lsq_splin() présentée ci-dessus permet également d'obtenir la dérivée : il suffit d'évaluer la spline avec la syntaxe suivante :

[yi, yprime2] = interp(x, xs, ys, d); // valeurs de la spline aux abscisses x
// et valeurs de la dérivée
Détermination des dérivée et dérivée seconde par l'algorithme de Savitzky-Golay.

L'algorithme de Savitzky-Golay présenté ci-dessus permet également d'obtenir la dérivée : comme on détermine localement un polynôme, on a les dérivées successives de ce polynôme.

Intégration[modifier | modifier le wikicode]

Intégrale numérique[modifier | modifier le wikicode]

Scilab peut calculer l'intégrale numérique d'une fonction entre deux bornes :

integrate("expression", "x", x0, x1)

calcule l’intégrale de la fonction décrite par expression ; expression est une chaîne de caractères interprétable par Scilab, comme par exemple "sin(x)". Le paramètre x est la variable d’intégration, et les bornes sont x0 et x1.

Par exemple :

// fonction carré
function [y] = f(x)
    y = x.^2
endfunction

// intégrale entre 0 et 4
Y = integrate("f(x)", "x", 0, 4)

// vérification
Y1 = 4^3/3

On peut aussi utiliser la fonction intg(a, b, f)a et b sont les bornes d'intégration. Avec l'exemple précédent :

-->Y = intg(0, 4, f)
 Y  =
 
    21.333333

Il peut également effectuer une intégration numérique à partir d'une suite de points (x(i), y(i)), par la méthode des trapèzes, avec la commande inttrap() :

x = linspace(0, 4, 50);
y = f(x);
inttrap(x, y)

Produit de convolution[modifier | modifier le wikicode]

Considérons deux suite de valeurs (vecteurs ligne) (hi)1 ≤ im et (xi)1 ≤ in, avec m < n. Le produit de convolution discret est défini par :

Il est obtenu par la commande conv(h, x) ou bien la commande convol(h, x). Par exemple, pour convoluer une fonction sinus avec une fonction porte de largeur 5 points :

t = linspace(0, %pi, 50);
x = sin(t);
h = 1/5*ones(1, 5);
y = conv(h, x);
plot(x)
plot(y, "r")

Les commandes conv et convol n'utilisent pas le même algorithme. La fonction convol utilise la transformée de Fourier rapide, et est mieux adaptée lorsqu'il y a beaucoup de valeurs.

Résolution d'équations[modifier | modifier le wikicode]

Scilab dispose de plusieurs primitives permettant la résolution d'équations. Nous présentons ci-après quatre fonctions internes, mais il en existe d'autres pour des systèmes spécifiques.

Équation polynomiale[modifier | modifier le wikicode]

Rappelons que la commande roots(p) détermine les racines du polynôme p, c'est-à-dire que cette fonction résout l'équation

p(x) = 0.

L'équation

p(x) = y

se résout donc avec

solution = roots(p - y)

Équation linéaire matricielle[modifier | modifier le wikicode]

La commande kernel() permet de déterminer le noyau de la matrice, c'est-à-dire, pour une matrice M, de résoudre l'équation matricielle

MX = 0

la syntaxe est simplement :

X = kernel(M)

La division de matrice consiste à résoudre une équation linéaire :

  • à droite : X = A/B est la solution de XB = A,
  • à gauche : X = A\B est la solution de AX = B,
  • on a B/A = (A' \ B')' ;

Notons que l'on peut représenter un système d'équations linéaires par une matrice ; la division matricielle est donc un moyen de résoudre un tel système.

Si la matrice A est creuse, il est recommandé d'utiliser umfpack() ; pour résoudre, X = A\B :

X = umfpack(A, "\" , B)

Par ailleurs, les résolutions numériques d'équations matricielles font souvent intervenir des décompositions d'une matrice :

  • lu(M)  : factorisation LU d'une matrice carrée, décomposition en un produit d'une matrice triangulaire inférieure L et d'une matrice triangulaire supérieure U,
  • qr(M)  : factorisation QR, décomposition en un produit d'une matrice orthogonale Q et une matrice triangulaire supérieure R,
  • svd(M) : décomposition en valeurs singulières (singular value decomposition).

Système d'équations linéaires[modifier | modifier le wikicode]

Scilab permet la résolution numérique d'un système d'équations linéaires, avec la fonction linsolve

[x, kerA] = linsolve(A, b)

où A est la matrice réelle des coefficients du système d'équations, et b un vecteur de constantes. Les solutions de

x + b = 0

sont de la forme

x0 + w×kerA,

w étant un réel arbitraire.

Exemple

Le système d'équation
est décrit par
donc
-->A = [3 1 ; 4 -1];

-->b = [-5 ; -9];

-->linsolve(A, b)
 ans  =

    2.
  - 1.
le vecteur
représente la solution x = 2 et y = -1.
Si maintenant on tape
-->[x0, kerA] = linsolve(A,b)
 kerA  = 

     []
 x0  =

    2.
  - 1.
indiquant que la solution est unique.

Exemple

S'il y a plus d'inconnues que d'équations :
-->A = [3 1] ; b = -5;

-->linsolve(A, b)
 ans  =

    1.5
    0.5

-->[x0, kerA] = linsolve(A,b)
 kerA  =

  - 0.3162278
    0.9486833
 x0  =

    1.5
    0.5
la première commande renvoie une seule solution ; la deuxième commande indique que les solutions sont de la forme
notons qu'il s'agit d'une solution approchée, puisque la vraie solution est du type :

Si le système n'admet pas de solution, Scilab affiche le message d'erreur : WARNING:Conflicting linear constraints!.

Scilab permet également la résolution symbolique d'un système linéaire avec la fonction solve

w = solve(A, c)

où A est une matrice triangulaire supérieure de réels (les coefficients du système d'équation), c est un vecteur de chaînes de caractères (la partie à droite des signes égal), et w est la matrice résultat de

w = c.

Exemple

Le système d'équation
est décrit par
-->A = ["1", "1" ; "0", "1"] ; c = ["a1" ; "a2"];

-->solve(A,c)
 ans  =

!a1-a2  !
!       !
!a2     !
indiquant que la solution est
et .

Scilab propose aussi la méthode du gradient conjugué avec la fonction conjgrad(). La syntaxe élémentaire est :

X = conjgrad(A, b)

Système d'équations quelconque[modifier | modifier le wikicode]

Scilab permet la résolution numérique d'un système d'équations quelconques, avec la fonction fsolve.

Considérons un système de n équations à n inconnues (x1, …, xn), de la forme

On définit, en tant que fonction externe (voir Programmation > Fonctions) la fonction ƒ qui associe à un vecteur x(x1, …, xn) un vecteur v(vx1, …, vn) :

Le système d'équations consiste donc à résoudre : ƒ(x) = 0. La syntaxe utilisée est :

[x] = fsolve(x0, f)

f est la fonction externe telle que définie ci-dessus et x0 une estimation de la solution. Dans le cas d'un système d'équations, l'entrée x est donc un vecteur de dimension n et la fonction f renvoie un vecteur de dimension n.

Scilab utilise un algorithme du type méthode hybride de Powell. La fonction peut donner le vecteur v correspondant à la solution x estimée (v devant être proche de 0), avec la syntaxe :

[x, v] = fsolve(x0, f)

Il est recommandé de fournir la jacobienne de ƒ, sous la forme d'une fonction externe ƒJ :

où J est la matrice jacobienne :

[x] = fsolve(x0, f, fJ)

Exemple

En reprenant le système d'équations ci-dessus
on écrit donc
function [y] = f(x)
    y = [3, 1 ; 4, -1]*x - [5 ; 9];
endfunction

x0 = [0 ; 0];
solution = fsolve(x0, f)
ce qui donne bien (2 ; -1).

Si la fonction est paramétrée, alors

  • dans la définition de la fonction, la variable doit venir avant les paramètres ;
  • lors de l'appel de fsolve(), la fonction doit être encapsulée avec ses paramètres dans une liste.

Exemple

Considérons le système
que l'on veut résoudre pour a = 1, b = 2, c = 3, d = 4. On écrit donc
function [Y] = f(X, A, b)
    Y = A*X + b;
endfunction

param1 = [1, 0 ; 0, 3];
param2 = [2 ; 4];

X0 = [0 ; 0]

solution = fsolve(X0, list(f, A = param1, b = param2));
// alternative : solution = fsolve(X0, list(f, A = param1, b = param2));
disp(solution);
ce qui donne comme résultat (-2 ; -1,333 333 3).

Équation différentielle ordinaire[modifier | modifier le wikicode]

Une d'équation différentielle ordinaire (ordinary differential equation), en abrégé EDO (ODE), peut être résolue de manière numérique avec la fonction ode. Soit une fonction y(t) ; si l'équation différentielle est

dy/dt = ƒ(t,y),

alors ƒ ayant été définie (fonction externe), la syntaxe pour déterminer y(t) est

y = ode(y0,t0,t,f)

  • y0 et t0 sont les valeurs initiales du système, y(t0) = y0,
  • t est un vecteur de valeurs pour lesquelles on calcule les solutions, et
  • y est le vecteur de solutions (plot(t,y) permet de tracer le graphique des solutions).

La fonction interne ode admet des arguments permettant de résoudre des situations spécifiques :

  • "roots" pour rajouter une équation g(t,y) = 0,
  • "discrete" pour calculer de manière récursive y(k+1) = ƒ(k,y(k)) à partir d'un état initial y(k0).
Exemple
On se place sur un segment [0,5;1] ; on veut déterminer les valeurs numériques de la fonction y vérifiant
et
on a donc en fait ƒ(t,y) = t , t0 = y0 = 0
-->t=[0.5:0.1:1]; // points de [0,5;1] espacés de 0,1

-->deff("[u] = f(t,y)", "u = t") // définition de la fonction f

-->ode(0,0,t,f)
 ans  =
    0.125    0.18    0.245    0.32    0.405    0.5
qui sont bien les valeurs de 1/2×t²
Note
La fonction ƒ doit être continue et différentiable en tout point. Si la fonction présente des singularité, alors il faut la résoudre sur la première partie D1 du domaine, puis sur la seconde partie en prenant comme conditions limites les résultats du premier calcul.
Pour aller plus loin
Voir G. Sallet, « Ordinary Differential Equations with Scilab – WATS Lectures, Provisional notes (Université de Saint-Louis) » sur Université de Metz, 2004. Consulté le 2018-02-23

Exemple complet[modifier | modifier le wikicode]

Levage d'un objet avec un pivot fixe en A et un vérin poussant en B.

Considérons un système mécanique simple : une poutre pivote autour d'un point fixe A(xA, yA) ; un vérin fixe pousse sur cette poutre selon un contact ponctuel au point B(xB, yB). La hauteur de la poutre vaut e = 200 mm, le pivot est implanté à mi-hauteur. Le vérin est implanté à une distance horizontale L = 1 000 mm. On a donc xBxA = L, et l'altitude yB varie. Le sujet présente une symétrie plane, on se place donc dans le plan de symétrie (A, x, y). Dans ce repère, on a donc :

On désire connaître l'inclinaison de la poutre en fonction du déplacement du vérin. Pour cela, on détermine l'équation de la droite génératrice de la poutre, située sur sa face inférieure. L'équation générale de cette droite (D) est :

(D) : y + ax + b = 0 soit (D) : y = –axb

et l'on veut donc déterminer les paramètres a et b. Les conditions géométriques sont[2] :

la droite passe par le point B : yB + axB + b = 0
la distance de A à (D) vaut r = e/2 :

On doit donc résoudre un système de deux équations, qui n'est pas linéaire.

On désire développer du code que l'on pourra copier-coller. Donc, au niveau des notations :

  • les inconnues à résoudre sont mises dans un vecteur X :  ;
  • les paramètres du problème sont mis dans un vecteur A : .

Le système d'équations devient donc :

Or ce problème présente deux solutions, puisqu'il existe deux droites passant par B et tangentes au cercle de centre A et de rayon r. Il faut donc :

  • choisir un point de départ X0 proche de la solution afin que l'algorithme converge vers le minimum local adapté ; on choisit la droite (A'B), avec
    A'(xA, yAr) ;
    (A'B) : y = yA + (yByA')/(xBxA)⋅(xxA), donc
    X0(1) = –(yByA')/(xBxA) = –(A(4) – A(2) + A(5))/(A(3) – A(1)) ;
    X0(2) = (yByA')/(xBxA)⋅xAyA = (A(4) – A(2) + A(5))/(A(3) – A(1)) × A(1) – A(2) ;
  • vérifier que la solution trouvée est bonne : on remarque que les deux droites possibles sont symétriques par rapport à la droite (AB), il suffit donc de vérifier que la pente de la droite trouvée est supérieure à celle de la droite (AB) :
    pente(AB) = (yByA)/(xBxA) = (A(4) – A(2))/(A(3) – A(1)).

Pour accélérer les calculs, lorsqu'une expression revient plusieurs fois, on la met dans une variable ; on rappelle que seul yB = A(4) est variable. On identifie donc :

facteur1 = 1/(A(3) – A(1)) ;
facteur2 = (–A(2) + A(5)) × facteur1 ;
facteur3 = A(4) × facteur1 + facteur2 ;

soit

X0(1) = –facteur3 ;
X0(2) = facteur3 × A(1) – A(2) ;
pente(AB) = (A(4) – A(2)) × facteur1.

Donc,

  • on connaît A(xA, yA, xB, yB, r)', la valeur yB prenant plusieurs valeurs successives ; on cherche X(a, b)' ;
  • on calcule les valeurs dérivées facteur1 = ƒ(A(1), A(3)) et facteur2 = ƒ(A(2), A(5)) qui sont fixes, et la valeur facteur3 = ƒ(facteur1, facteur2, A(4)) qui dépend de yB ;
  • pour chaque valeur de yB (et donc de facteur3) :
    • on détermine la valeur de départ X0,
    • on résout le système d'équations, donc on détermine a et b,
    • on vérifie que a est compatible avec la réalité.

On a donc deux vecteurs de valeurs a et b, les valeurs a(i) et b(i) correspondant au déplacement du piston yB(i), ce qui permet de déterminer l'inclinaison de la poutre pour chaque valeur yB(i). La valeur b(i) ne nous sert à rien (b est juste un intermédiaire de calcul permettant de résoudre le système), mais nous la conservons tout de même car elle permet de tracer rapidement la droite solution si cela nous intéresse.

Isolement de la poutre, et représentation des caractéristiques connues des forces : on ne connaît que le point d'application de FA (flèche brisée), et l'on connaît la direction de FB (double flèche avec trait d'axe).

Une fois ceci fait, on désire déterminer les effort aux liaisons, c'est-à-dire la force FA que l'axe exerce sur la poutre en A, et la force FB que le vérin exerce sur la poutre en B. La force FA est totalement inconnue :

et on connaît la direction de la force FB : c'est le vecteur normal (–a, 1). Le vecteur peut donc s'écrire sous la forme :

.

Le système a une m, le centre des masses étant en G(xG, yG) ; prenons arbitrairement m = 100 kg, xG = 700 mm et yG = 200 mm. L poids P s'écrit p(0, –mg) avec g = 9,81 m/s2. Les équations de la statique s'écrivent

équilibre des forces :
équilibre des moments par rapport à A :

soit

équilibre des forces :
équilibre des moments par rapport à A :.

Si l'on définit le vecteur des inconnues

et les vecteurs de paramètre

les équations de la statique s'écrivent sous la forme d'une équation matricielle :

AFXF + bF = 0.

En résolvant ce système, on peut en déduire les intensités des forces :

 ;
.

Concrètement, le programme développé est :

// paramètres du problème géométrique
xA = 0;
yA = 0;
e = 50; // en mm
L = 1000; // en mm
xB = L;
course = 200; // course du vérin en mm
yBmin = -e/2;
yBmax = yBmin + course;

// paramètres du problème de statique
xG = 700; // en mm
yG = 200; // en mm
m = 100; // en kg
g = 9.81; // en m/s^2

n = 11; // nombre de points calculés

// paramètres dérivés
r = e/2;
yB = linspace(yBmin, yBmax, n); // vecteur de valeurs de y pour B

// fonction
function [Y] = contraintes_geo(X, A)
    // droite d'équation y + az + b = 0
    // a et b à déterminer
    // X(1) = a ; X(2) = b
    // A(1) = xA ; A(2) = yA ; A(3) = xB ; A(4) = yB ; A(5) = r ;
    // 1. B est sur la droite
    //    yB + axB + b = 0
    Y(1) = A(4) + A(3)*X(1) + X(2);
    // 2. La distance de A à la droite vaut r
    //    (yA + axA + b)^2/(1 + a^2) - r^2 = 0
    Y(2) = (A(2) + A(1)*X(1) + X(2))^2/(1 + X(1)^2) - A(5)^2;
endfunction

// ***********************
// * Programme principal *
// ***********************

// ********** Problème géométrique **********

// initialisation
param = [xA ; yA ; xB ; yBmin ; r]; // correspond à A dans la fonction
a = zeros(n);
b = a;

// simplification du calcul
facteur1 = 1/(param(3) - param(1));
facteur2 = (param(5) - param(2))*facteur1;

// recherche des valeurs
for i = 1:n
    param(4) = yB(i);
    facteur3 = param(4)*facteur1 + facteur2;
    X0 = facteur3*[-1 ; param(1)] + [0 ; param(2)]; // point de départ
    solution = fsolve(X0, list(contraintes_geo, param));
    pente = -solution(1);
    a(i) = -pente;
    b(i) = solution(2);
    // vérification que l'on a pris la bonne solution
    penteAB = (param(4) - param(2))*facteur1;
    if pente < penteAB then // si c'est la mauvaise solution
        a(i) = %nan;
        b(i) = %nan;
    end
end

// affichage de l'angle
angle = atan(-a); // en radians

disp("position (mm) ; angle (°)");
disp([yB', 0.1*round((1800/%pi)*angle)]); // angle en °, arrondi

// ********** Problème de statique **********

// initialisation
FA = zeros(n, 2); 
C = zeros(n, 1);
iFA = C; // intensités des forces
iFB = C;

AF = [1, 0, -a(1)
0 1 1
0 0 L];
bF = -m*g*[0 ; 1 ; xG];

// Résolution des équations de la statique

for i=1:n
    AF(1, 3) = a(i);
    XF = linsolve(AF, bF);
    FA(i, :) = XF(1:2);
    C(i) = XF(3);
    iFA(i) = norm(XF(1:2));
    iFB(i) = C(i)*norm([a(i) ; 1]);
end

// affichage des intensités
disp("position (mm) ; FA (N) ; FB(N)");
disp([yB', round(iFA), round(iFB)]); // angle en °, arrondi

On obtient donc :

position (mm), angle (°)

 -25.    0.       
  -5.    1.1
  15.    2.3
  35.    3.4 
  55.    4.6
  75.    5.7
  95.    6.9
 115.    8.
 135.    9.1
 155.   10.2
 175.   11.3

 position (mm) ; FA (N) ; FB(N)

  -25.    294.   687.
  -5.     295.   687.
   15.    296.   687.
   35.    297.   688.
   55.    299.   689.
   75.    302.   690.
   95.    306.   692.
   115.   310.   693.
   135.   314.   695.
   155.   319.   698.
   175.   325.   700.

Optimisation d'une fonction[modifier | modifier le wikicode]

Soit une fonction ƒ de . L'optimisation consiste à trouver le vecteur x (vecteur à n composantes) donnant la valeur minimale de ƒ.

Optimisation linéaire[modifier | modifier le wikicode]

La fonction ƒ est une application linéaire ; elle peut donc s'écrire :

c est un vecteur de et ct est sa transposée. On peut par ailleurs restreindre le domaine de recherche à un polyèdre convexe décrit par m inéquations :

  • A est une matrice m×n ;
  • b est un vecteur de m dimensions.

Scilab dispose de la commande karmarkar() qui permet de résoudre le problème avec l'algorithme de Karmarkar.

La syntaxe la plus simple permet de résoudre le problème sur la frontière du polyèdre, donc avec des égalités partout :

xopt = karmarkar(Ae, be, c)

résout

Exemple

Si l'on veut minimiser la fonction ƒ(x1, x2, x3) = –x1x2 avec pour contraintes x1x2 = 0 et x1 + x2 + x2 = 2, on a :

La solution s'obtient donc par

Aeq = [1, -1, 0 ; 1, 1,1];
beq = [0 ; 2];
c = [-1 ; -1 ; 0];
xopt = karmarkar(Aeq, beq, c)

et le résultat est :

xopt  = 

   0.9999949
   0.9999949
   0.0000102

On peut y ajouter un système d'inéquations[3] :

xopt = karmarkar(Ae, be, c, [], [], [], [], [], Ai, bi)

résout

et si l'on n'a que des inéquations :

xopt = karmarkar([], [], c, [], [], [], [], [], Ai, bi)

résout

On peut indiquer un vecteur de départ x0 :

xopt = karmarkar(A, b, c, x0)
xopt = karmarkar(A, b, c, x0, [], [], [], [], Ai, bi)
xopt = karmarkar([], [], c, x0, [], [], [], [], Ai, bi)

et l'on peut demander de calculer la valeur de ƒ(xopt) :

[xopt,fopt] = karmarkar()
Voir aussi

Optimisation non linéaire[modifier | modifier le wikicode]

Ici, la fonction ƒ de est non linéaire. On utilise pour cela la fonction optim :

[fopt, xopt] = optim(costf, x0)

  • costf est la « fonction de coût » de ƒ ; c'est une fonction qui renvoit la valeur de la fonction ƒ en x et le gradient de ƒ en x, défini sous la forme
    function [f, g, ind] = costf(x, ind)
    f désigne ƒ(x), g est le gradient et ind est un index, un entier permettant de modifier le comportement de costf ;
  • x0 est une estimation de la solution ;
  • xopt est le vecteur trouvé ;
  • fopt = ƒ(xopt), valeur estimée du minimum.

On peut indiquer des bornes inférieures et supérieure de x, sous al forme de vecteurs xinf et xsup :

[fopt, xopt] = optim(f, x0, "b", xinf, xsup)

La fonction optim() est une méthode de quasi-Newton utilisant les critères de Wolfe.

Optimisation quadratique[modifier | modifier le wikicode]

Soit ƒ une fonction quadratique de n variables (xi)1 ≤ in, c'est-à-dire une combinaison linéaire de xixj. Cette fonction peut se décrire par deux matrices, Q, définie positive de dimension n×n, et p, matrice colonne de dimension n :

si l'on appelle x la matrice colonne [x1 ; x2 ; … ; xn]
ƒ(x1, …, xn) = ½txQx + tpx
tM est la transposée de la matrice M.

L'optimisation quadratique consiste à trouver le vecteur x donnant la valeur minimum de ƒ, en imposant de plus que la solution se trouve dans un espace convexe, ce qui se traduit par m contraintes linéaires : me conditions d'égalité

1 ≤ jnCijxj = bi, 1 ≤ ime

et m - me conditions d'inégalité

1 ≤ jnDijxjdi, 1 ≤ im - me

soit sous forme matricielle

C1x = b1
C2xb2

  • C1 est une matrice me×n ;
  • b1 est une matrice colonne de dimension me ;
  • C2 est une matrice (m - men ;
  • b2 est une matrice colonne de dimension (m - me).

On rassemble les matrices :

  • C = [C1 ; C2] ;
  • b = [b1 ; b2].

Pour résoudre un tel système, Scilab propose deux commandes.

La commande qpsolve, sous la forme :

[x] = qpsolve(Q, p, C, b, ci, cs, me)

utilise la fonction Fortran qpgen1.f (également appelée QP.solve.f), developpée par Berwin A. Turlach selon l'algorithme de Goldfarb/Idnani.

Les paramètres ci et cs sont des vecteurs colonne de dimension n contenant les limites inférieures et supérieures aux valeurs des xi

cixcs.

Si l'on n'impose pas de limite, on utilise des matrices vides [].

La commande qld, sous la forme :

[x, lagr] = qld(Q, p, C, b, ci, cs, me)

utilise la méthode de Powell modifiée par Schittkowski.

La variable lagr est un vecteur de multiplicateurs lagrangiens.

Le problème de la précision et de l'exactitude[modifier | modifier le wikicode]

Scilab est destiné à faire du calcul. La précision, l'exactitude des résultats est donc une préoccupation primordiale.

La première chose est bien évidemment d'utiliser des méthodes, des algorithmes corrects, validés, documentés. Mais la mise en code de ces algorithmes doit également être correcte. Bien évidemment, la méthode choisie doit être retranscrite fidèlement et sans erreur, mais il faut aussi prendre en compte certains « effets de bord » dus à la représentation des nombres.

En effet, Scilab calcule par défaut avec des nombre décimaux codés en virgule flottante à double précision selon la norme IEEE754. Cela permet de représenter des nombres dont la valeur absolue est comprise entre 10-307 et 10308 ; si un calcul donne une valeur inférieure, cela retourne la valeur 0 (erreur de soupassement), et s'il donne une valeur supérieure, cela retourne la valeur Inf ou -Inf (erreur de dépassement). Le nombre de chiffres significatifs est de l'ordre de 16 ; ainsi, si deux nombres diffèrent au 17e chiffre significatif, ils seront considérés comme égaux (erreur d'arrondi), leur différence sera nulle, leur rapport sera de 1.

Il faut donc s'assurer d'une part que le valeur finale est représentable, mais aussi que les résultats des calculs intermédiaires sont eux aussi représentables. Par exemple le calcul de abs(-10^-200) donne bien le résultat 10-200, alors que sqrt((-10^-200)^2) donne 0, puisque le résultat intermédiaire (-10200)2 n'est pas représentable (soupassement). Dans certains cas, un algorithme itératif, cherchant une valeur approchée, peut donner un meilleur résultat qu'un calcul direct utilisant une formule exacte… Par exemple, pour calculer l’hypoténuse d'un triangle rectangle (norme euclidienne), on peut utiliser la méthode de Moler et Morrison (1983)[4], ce que fait d'ailleurs la commande Scilab norm().

Si l'on doit faire des calculs avec des cas « extrêmes » — nombres très proches, ou bien ayant une valeur absolue élevée ou faible —, on s'intéressera à des algorithmes « robustes ». Ceux-ci font fréquemment intervenir une normalisation des données ou une factorisation pertinente. Dans la mesure du possible, si l'on a le choix entre une fonction « fait maison  » et une fonction Scilab faisant le même travail, on favorisera la fonction Scilab, qui est censée être optimisée.

Si des calculs intermédiaires font intervenir des nombres positifs très grands ou très petits, on peut utiliser la fonction logarithme, qui a initialement été inventée pour cela. Par exemple, pour calculer

avec

on peut utiliser la fonction logarithme de la fonction gamma, gammaln()[5] :

1-exp(gammaln(365 + 1)-(gammaln(365 - 23 + 1) + 23*log(365)))

On peut également retravailler la formule en

et utiliser

1 - prod( (343/365):(1/365):(364/365) )

Quoi qu'il en soit, la mise en œuvre d'un algorithme doit toujours s'accompagner de tests avec des valeurs pour lesquelles le résultat est connu. Les résultats de référence peuvent être calculés à la main, obtenus avec un logiciel de référence, ou bien les données peuvent être générées en fonction du résultat attendu.

Par ailleurs, un résultat de qualité fait figurer la précision du calcul.

On pourra se reporter à Baudin 2010a.

Voir aussi[modifier | modifier le wikicode]

Dans Wikipédia
Liens externes

Notes[modifier | modifier le wikicode]

  1. les fonctions derivative() et numdiff() sont déclarées obsolètes et ont été supprimées à partir de la version 6
  2. voir Propriétés métriques des droites et des plans sur Wikipédia.
  3. les cinq paramètres vides [] sont des paramètres que nous ne présentons pas ici à l'exception du premier, x0 ci-dessous. Ils sont décrits dans la page d'aide : (en) karmarkar: Solves a linear optimization problem. sur help.scilab.org. Consulté le 2017-02-13
  4. Baudin 2010b, p. 44–45
  5. voir (en) Computing with very large numbers in Scilab sur liste de discussion Scilab-users et gammaln sur Aide de Scilab

Calculs élémentaires < > Graphiques et sons