Découvrir Scilab/Version imprimable

Un livre de Wikilivres.
Sauter à la navigation Sauter à la recherche
Nuvola-inspired File Icons for MediaWiki-fileicon-ps.png

Ceci est la version imprimable de Découvrir Scilab.

  • Si vous imprimez cette page, choisissez « Aperçu avant impression » dans votre navigateur, ou cliquez sur le lien Version imprimable dans la boîte à outils, vous verrez cette page sans ce message, ni éléments de navigation sur la gauche ou en haut.
  • Cliquez sur Rafraîchir cette page pour obtenir la dernière version du wikilivre.
  • Pour plus d'informations sur les version imprimables, y compris la manière d'obtenir une version PDF, vous pouvez lire l'article Versions imprimables.


Découvrir Scilab

Une version à jour et éditable de ce livre est disponible sur Wikilivres,
une bibliothèque de livres pédagogiques, à l'URL :
https://fr.wikibooks.org/wiki/D%C3%A9couvrir_Scilab

Vous avez la permission de copier, distribuer et/ou modifier ce document selon les termes de la Licence de documentation libre GNU, version 1.2 ou plus récente publiée par la Free Software Foundation ; sans sections inaltérables, sans texte de première page de couverture et sans Texte de dernière page de couverture. Une copie de cette licence est incluse dans l'annexe nommée « Licence de documentation libre GNU ».

Interface

Table des matièresIndex



1. Interface utilisateur


Saisie des commandes et affichage numérique[modifier | modifier le wikicode]

Scilab est un logiciel qui s'utilise en ligne de commande. L'écran se compose d'une zone de saisie dans laquelle on va taper des commandes.

L’invite de la ligne de commande (prompt) est constituée d’une « flèche » : deux tirets et un signe supérieur -->. L’instruction est tapée puis validée avec la touche de retour chariot ([↵], [Enter] ou [Return]). Le résultat est affiché à la suite, sauf si la ligne se termine par un point-virgule auquel cas le résultat est caché.

Par exemple 
  • sans point-virgule
-->a=1[↵]
 a  =

    1.

-->A=2[↵]
 A  =

    2.

-->a+A[↵]
ans  =

   3.
  • avec points-virgules
-->a=1;[↵]

-->A=2;[↵]

-->a+A[↵]
ans  =

   3.

Lorsque l'on écrit une ligne, on peut utiliser les touches classiques d'édition (variable selon le système et le type du clavier) :

  • [→] : avancer d'une lettre ;
  • [←] : reculer d'une lettre ;
  • touches de suppression ([Del], [Suppr], []) ;
  • touche de retour chariot [↵], [Enter] ou [Return] : validation et exécution de la commande ;
  • [↑] et [↓] : permettent de faire défiler les commandes tapées précédemment (historique des commandes).

L'endroit où se place le caractère (lettre, chiffre, signe…) tapé est signalé par le curseur, un trait de soulignement clignotant (Windows XP) ou un rectangle (MacOS X). L'édition se fait en insertion, c'est-à-dire que si le curseur est au milieu d'une ligne, les caractères suivants sont décalés d'une position à droite lorsque l'on entre un caractère.

On peut, si on le désire, fractionner l'édition sur plusieurs ligne. Pour cela, il faut terminer la ligne par trois points « ... ».

Exemple
-->a=...[↵]
-->5[↵]
 a  =

    5.

Ceci est bien sûr plutôt intéressant lorsque la ligne à taper est longue.

À l'inverse, on peut mettre plusieurs instructions sur une même ligne en les séparant d'une virgule, ou bien d'un point-virgule si l'on ne veut pas que le résultat soit affiché.


Note
Par la suite, le retour chariot [↵] sera omis, un retour à la ligne après une ligne commençant par l'invite --> l'indiquera implicitement.

Chaînes de caractère[modifier | modifier le wikicode]

Par défaut, une suite de caractères est interprété comme une commande ou comme une variable. Pour indiquer que c'est une chaîne de caractères, il faut la mettre entre guillemets simples (apostrophes, single quote) « ' » ou entre guillemets doubles (double quote) « " ».

Lorsque la chaîne doit contenir une apostrophe, on redouble l'apostrophe. De même, lorsque la chaîne doit contenir des guillemets, on redouble ces guillemets.

Exemple
-->l'amertume de la "bière"
  !--error 4
undefined variable : l


-->'l''amertume de la ""bière""'
 ans  =

 l'amertume de la "bière"

-->"l''amertume de la ""bière"""
 ans  =

 l'amertume de la "bière"

Les chaînes de caractères sont utilisées pour définir certains paramètres de fonctions, ou bien pour afficher des étiquettes et légendes sur les graphiques.

Éditeur SciNotes[modifier | modifier le wikicode]

Scilab dispose d'un éditeur de programme, SciNotes, qui peut être lancé en tapant edit dans la ligne de commande. On peut aussi taper scinotes, mais cette solution est sensible à un changement du nom de l'éditeur ; par exemple, avant la version 5.3.0 (2010), l'éditeur était SciPad[1].

Il s'agit d'un éditeur multiligne : lorsque l'on appuie sur le retour chariot, cela ne valide pas la commande mais passe simplement à la ligne suivante. Les touches [↑] et [↓] permettent de changer de ligne, d'aller éditer une autre ligne.

Ce qui est tapé dans cet éditeur peut ensuite être sauvegardé dans un fichier .SCI ou bien directement exécuté dans Scilab (l'ancienne extension était .SCE).

Affichage graphique[modifier | modifier le wikicode]

L'affichage graphique (tracé d'une courbe, dessin…) se fait dans une fenêtre dédiée que nous appelerons « fenêtre graphique ». Cette fenêtre s'affiche automatiquement dès la première instruction graphique. Par exemple, taper clf puis valider : cela affiche une fenêtre graphique vide.

Aide en ligne[modifier | modifier le wikicode]

Scilab dispose d'une aide en ligne. Pour cela, taper help.

Cela affiche une fenêtre Browse Help (« naviguer dans l'aide ») avec une arborescence à gauche. Lorsque l'on clique sur un dossier de l'arborescence, une branche se développe, faisant apparaître les articles. Chaque article correspond à une commande ou à une fonction prédéfinie ; il suffit de cliquer dessus pour que le contenu de l'aide — disponible en français, mais quelques articles restent en anglais — s'affiche dans la partie droite de la fenêtre.

Lorsque l'on clique sur le bouton Rechercher, un champ de recherche apparaît dans lequel on peut entrer un mot-clef ; lorsque l'on appuie sur le retour chariot [↵], le programme effectue une recherche et affiche les article correspondants dans la zone de texte située dans la partie gauche de la fenêtre. il suffit de cliquer dessus pour que le contenu de l'article s'affiche dans la partie droite de la fenêtre. On revient à une navigation par arborescence en cliquant sur le bouton Arborescence.

Fermer le programme[modifier | modifier le wikicode]

Scilab peut être fermé de la manière habituelle selon les systèmes d'exploitation : en général bouton en forme de croix ou menu Fichier | Quitter. On peut aussi

  • taper exit en ligne de commande ;
  • lorsque l'on est dans l'environnement initial, on peut également taper quit (voir aussi Environnement).

Spécificités selon les systèmes d'exploitation[modifier | modifier le wikicode]

Sous Microsoft Windows[modifier | modifier le wikicode]

Éditeur de ligne de commande[modifier | modifier le wikicode]

Voici les effets des touches suivantes dans l'éditeur de ligne de commande :

  • [↖] ou [Home] : aller en début de ligne ;
  • [Fin] ou [End] : aller en fin de ligne ;
  • [Échap] ou [Esc] efface la ligne courante.

On peut utiliser les actions classiques de la souris, par exemple sélectionner une ligne de texte et faire du copier-coller :

  • en utilisant le menu contextuel accessible en cliquant avec le bouton droit [2],
  • avec les boutons d'édition situés en haut de la fenêtre, ou
  • avec les raccourcis clavier, par exemple [Ctrl]+[C] et [Ctrl]+[V].

Attention
l'utilisation de [Ctrl]+[C] lorsqu'aucun texte n'a été sélectionné avec la souris a un effet totalement différent, voir le chapitre Environnement. Il est donc préférable d'éviter le raccourci clavier pour la fonction copier.

Le texte ne peut être sélectionné qu'avec la souris ; la combinaison de touches [⇑]+[→] ou [⇑]+[←], utilisable dans de nombreuses autres applications, ne fonctionne pas[3]. Il est possible de copier du texte dans une autre application (par exemple un éditeur de texte ou bien du texte affiché par un navigateur Internet) et de le coller dans l'éditeur de ligne de commande de Scilab. On peut ainsi utiliser des instructions écrites au préalable, éventuellement par d'autres personnes (voir aussi Programmation > Chargement d'une fonction).

Notons que le raccourci Scilab pointe vers un programme nommé WScilex.EXE situé dans ~\bin\, le tilde ~ indiquant ici le répertoire (dossier) d'installation (par défaut C:\Program Files\Scilab-4.0\ pour la version 4.0). Il existe un autre programme dans le même répertoire, Scilex.EXE, qui est un éditeur en ligne de commande Scilab mais sans « l'environnement habituel » Microsoft Windows : la fenêtre n'a pas de menu ni de bouton, seule la ligne de commande est active, et on ne peut pas copier le texte de la fenêtre (on en peut pas le sélectionner avec la souris), ni coller du texte copié d'une autre application Microsoft Windows. On peut tout de même lancer l'aide en tapant help.

Éditeur de programme SciNotes[modifier | modifier le wikicode]

Concernant l'éditeur de programme SciNotes, il est possible de faire du copier-coller depuis ou vers une autre application Microsoft Windows, et les combinaisons de touches [⇑]+[→] et [⇑]+[←] permettent de sélectionner du texte.

Aide en ligne[modifier | modifier le wikicode]

On peut lire directement le fichier d'aide au format CHM (compiled HTML, un format d'aide propre à Microsoft Windows). Ce fichier contient les mêmes informations que l'aide en ligne accessible par help, seule change l'interface ; le fichier s'appelle scilab_fr_FR_help.chm (la nomenclature était différente dans les versions précédentes, le fichier s'est appelé sciman-fr-3.1.1.CHM, man-fr-scilab-4.0.CHM)[4]. Il se trouve dans le dossier C:\Program Files\scilab-5.4.0\modules\helptools\chm.

On lance cette aide avec un double-clic sur l'icône dans l'Explorateur Windows. La fenêtre comporte deux parties :

  • à gauche, un panneau de navigation avec trois onglets :
    • sommaire : c'est une arborescence, similaire à l'aide help ;
    • index : permet une recherche par mots-clefs (des mots-clefs sont associés aux articles) ;
    • rechercher : permet une recherche par chaîne de caractère (lorsque l'on cherche un mot qui n'a pas été entré comme mot-clef) ;
  • à droite, un panneau où s'affiche l'article sélectionné dans le panneau de navigation.

Sous MacOS X[modifier | modifier le wikicode]

Scilab fonctionne dans l'environnement X11 ; on ne peut pas copier un texte dans une autre application puis le coller dans l'éditeur en ligne de Scilab. Par contre, c'est possible dans SciPad. SciPad peut être lancé en cliquant sur le bouton Editor en haut de la fenêtre.

Il est possible de copier du texte d'une application extérieure (y compris sous Aqua[5]) de manière classique (par exemple avec touche pomme+[C]), puis de le coller dans SciPad avec le menu Edit | Paste ou bien le raccourci clavier [Ctrl]+[V].


Note
Les raccourcis pour le copier-coller sont différents entre les applications sous Aqua (touche touche pomme) et SciPad (touche [Ctrl]).

Les touches de fonction classiques fonctionnent normalement dans l'éditeur de ligne de commande :

  • [↖] : aller en début de ligne ;
  • [↘] : aller en fin de ligne.

Voir aussi[modifier | modifier le wikicode]

Dans Wikipédia

Notes[modifier | modifier le wikicode]

  1. SciPad est toujours disponible en tant que produit à part, sur http://sourceforge.net/p/scipad/
  2. pour les droitiers, ou le bouton gauche si la souris est configurée pour les gauchers
  3. [⇑] désigne la touche permettant d'accéder aux capitales, souvent notée [Maj] pour « majuscules » bien que le terme soit impropre (voir capitale et majuscule sur Wikipédia Article sur Wikipédia) ; cette touche est appelée « shift » en anglais
  4. le terme « man » fait référence à la fonction Unix permettant d'afficher les pages de manuel.
  5. Aqua est l'environnement graphique « naturel » de MacOS X

Introduction < > Calculs élémentaires


Calculs élémentaires

Table des matièresIndex



2. Calculs élémentaires


Scilab peut s'utiliser comme une calculatrice scientifique, de type calculatrice graphique. Il suffit de taper une formule dans une syntaxe relativement standard, et d'appuyer sur le retour chariot pour avoir le résultat.

Opérateurs et fonctions classiques[modifier | modifier le wikicode]

Scilab utilise les fonctions et opérateurs classiques :

  • + (addition), - (soustraction), * (produit), / (division), ^ ou ** (élévation à une puissance),
  • modulo(a, b) : reste entier de a/b, pmodulo(a, b) donne un résultat toujours positif,
  • sqrt() (racine carrée), nthroot(x, n) (racine n-ième de x)
  • cos(), sin(), tan() pour les fonctions trigonométriques de base, acos(), asin(), atan() pour leurs réciproques, cotg() pour la cotangente ;
  • les fonctions trigonométriques hyperboliques cosh(), sinh(), tanh(), acosh(), asinh(), atanh()
  • log() pour le logarithme népérien (ln), log10() pour le logarithme en base 10 (log10), log2() pour le logarithme en base 2 (log2),
  • exp() pour la fonction exponentielle,
  • floor() pour la partie entière, int() pour la troncature (arrondit à l'entier vers 0, différent de floor pour les nombres négatifs), round() pour l’arrondi au plus proche (0,5 → 1), abs() pour la valeur absolue,
  • erf() pour la fonction d'erreur de Gauss,
  • rand pour avoir un nombre aléatoire entre 0 et 1, suivant une loi uniforme ; rand(1, "normal") pour avoir un nombre aléatoire suivant une loi normale centrée réduite (moyenne 0, variance 1) ;
  • la factorielle de n peut s'obtenir par prod(1:n) ou bien par gamma(n+1) ;

Le séparateur décimal est le point.

Les puissances de 10 se notent avec la lettre « e » en minuscule ou en capitale : 5,2·103 s'écrit 5.2e3 ou 5.2E3, 5,2·10-3 s'écrit 5.2e-3 ou 5.2E-3. On peut aussi utiliser le « d » minuscule ou capital.

Par défaut, Scilab affiche les résultats avec dix caractères (signe et point décimal inclus). On peut modifier le nombre de caractères avec la fonction format(n)n est le nombre de caractères.

Mentionnons également la fonction clear qui permet d'effacer toutes les variables (dont les fonctions définies par l'utilisateur). Cela peut être utile lorsque l'on a procédé par essai-erreur et que l'on veut remettre les « compteurs à zéro ».


Vocabulaire
Dans Scilab, une « fonction interne[1] » est une fonction prédéfinie. Par opposition, les fonctions définies par l'utilisateur sont appelées « fonctions extérieures ».

Constantes particulières[modifier | modifier le wikicode]

Pour entrer l’imaginaire i, il faut utiliser %i ; il figure simplement sous la forme « i » dans les résultats.

Pour entrer l’infini ∞, il faut utiliser %inf ; il figure simplement sous la forme « inf » dans les résultats.

La valeur de π s’obtient par %pi, et la constante de Neper e par %e.

Lorsqu'une opération invalide se produit, certains programmes retournent un nombre particulier appelé « NaN » pour not a number (littéralement « ceci n'est pas un nombre »)[2]. Dans Scilab, cette valeur est obtenue par %nan et s'affiche sous la forme Nan.

Variables[modifier | modifier le wikicode]

On peut définir des variables simplement sous la forme NomVariable=valeur. NomVariable est une chaîne de caractères commençant par une lettre ; Scilab est sensible à la casse, c'est-à-dire que les variables « a » et « A » sont deux variables différentes. Le type de la variable (entier, réel, imaginaire) est géré automatiquement.

La variable ans contient le dernier résultat (pour answer, « réponse » en anglais).

La fonction who affiche les variables déclarées (who? signifie « qui ? » en anglais). On remarque que dès l'ouverture, Scilab a défini de nombreuses variables, notamment des variables dites « d'environnement », c'est-à-dire des paramétrages de fonctionnement de Scilab.

Exemples[modifier | modifier le wikicode]

Exemple 1[modifier | modifier le wikicode]

Vérification de l'équation (= sin(45°) ≈ sin(0,78) ≈ 0,71) :

-->asin(sqrt(2)/2)
 ans  =

    0.7853982

-->ans/%pi*180
 ans  =

    45

Exemple 2[modifier | modifier le wikicode]

-->%e^(%i*%pi/2)
 ans  =

    6.1273D-17 + i

Notez que dans l'absolu, on devrait obtenir i.

Exemple 3[modifier | modifier le wikicode]

Calcul de la pression, en Pa, de 2,5 moles de gaz parfait dans un volume de un litre (10-3 m3) à une température de 300 K (26,85 °C)

-->n = 2.5; R = 8.314472; T = 300; V = 1e-3;

-->n*R*T/V
 ans  =

    6235854

soit environ 6,2·106 Pa.

Complexes[modifier | modifier le wikicode]

Les nombres complexes sont entrés sous la forme a + b*%i. On peut utiliser la commande complex(a, b) pour créer ce nombre. Les lignes suivantes sont équivalentes :

z = 1 + 2*%i
z = complex(1,2)

Les commandes real(z) et imag(z) donnent respectivement la partie réelle et la partie imaginaire du nombre complexe. La conjugaison complexe z s'obtient par conj(z).

Le module d'un nombre complexe s'obtient avec la commande abs(z). L'argument complexe de z = a + bi s'obtient par

phi = atan(b, a)
phi = atan(imag(z), real(z))

La commande isreal(z) teste si un nombre est réel, ou bien si une matrice ne contient que des réels.

Pour multiplier par l'imaginaire pur i, il vaut mieux utiliser la commande imult(a) que %i*a.

Booléens[modifier | modifier le wikicode]

Lorsque l’on assigne une valeur booléenne, on utilise %t pour « vrai » (true) et %f pour « faux » (false). Le résultat de l’opération affiché est respectivement T ou F.

L’opérateur ou (OR) est noté « | » (tube), et (AND) est noté « & » (esperluette), non (NOT) est noté « ~ » (tilde). Le ou exclusif (XOR) s'obtient avec ~=.

Exemple

-->(%t & %f) | %f
ans =

F

On peut également utiliser les fonctions and() et or() qui appliquent les opérateurs respectifs et et ou à tous les éléments d'une matrice. L'exemple ci-dessus devient :

-->or([and([%t, %f]), %f])
ans =

F

Un booléen s’obtient en général en comparant deux valeurs, avec les relations d’égalité == et de différence <> ou ~=, et les relations d'ordre <, <=, > et >=.

Exemple

-->a = 3
 a  =

    3.

-->a == 3
 ans  =

  T

-->a < 2
 ans  =

  F

Sur des matrices de même dimension, ou bien entre un scalaire et une matrice, l'égalité s'évalue élément par élément. Si l'on savoir si deux matrices sont égales, on utilise l'opérateur and(). Si l'on veut savoir si un scalaire est dans un ensemble, on peut écrire l'ensemble sous la forme d'un vecteur ou d'une matrice et utiliser l'opérateur or().

Exemple

[1 2 ; 3 4] == [1 3 ; 2 4]
 ans  =

  T F
  F T

--> and([1 2 ; 3 4] == [1 3 ; 2 4])
 ans  =

  F

--> 2 == [1 2 3 4]
 ans  =

  F T F F

--> or(2 == [1 2 3 4]) // 2 ∈ {1 ; 2 ; 3 ; 4} ?
 ans  =

  T

Il existe également un certain nombre de commandes de test renvoyant un booléen, par exemple :

  • isempty(a) : renvoie « vrai » si c'est la variable vide, « faux » sinon[3] ;
  • isnan(a) : vérifie si la variable a la valeur NaN (not a number, résultat d'une opération invalide) et renvoie « vrai » dans ce cas ;
  • isinf(a) : vérifie si la variable a une valeur infinie et renvoie « vrai » dans ce cas ;
  • isequal(a, b, …) : vérifie si les variables a, b, … sont toutes égales ;
  • isreal(a) : vérifie si a est une matrice de réels.

Si la variable x est un booléen, la commande bool2s(x) (boolean to standard) transforme la valeur « faux » en « 0 » et la valeur « vrai » en « 1 ». Si c'est un nombre, la commande transforme toute valeur réelle non nulle en « 1 », et une valeur nulle en « 0 ». Les valeurs NaN (%nan</code) et ∞ (%inf) ne sont pas nulles et donc sont transformées en « 1 ».

À l'inverse, les opérateur booléens — &, |, and() et or() — s'appliquent aux réels et aux complexes : toute valeur non nulle est considérée comme vraie, la valeur nulle (0 ou 0 + 0⋅i) est considérée comme fausse.

Polynômes et fractions rationnelles[modifier | modifier le wikicode]

Scilab manipule les polynômes et fractions rationnelles formels, c'est-à-dire non pas en tant que fonction, mais en tant qu'objets mathématiques, éléments du corps des polynômes ou des fractions rationnelles.

La première opération pour définir un polynôme est de définir le monôme x. Ceci se fait par la commande x = poly(0, "x"). Cette commande indique que la variable x est un polynôme (commande poly()) dont l'indéterminée est notée « x » (le "x" situé à l'intérieur de la parenthèse) et la racine est 0.

Puis, le polynôme ou la fraction rationnelle est définie par une formule contenant ce x.

Exemple de définition de polynôme

-->x = poly(0, "x"); p = x ^ 2 + 2 * x + 1
 p  =

              2
    1 + 2x + x

Exemple de définition de fraction rationnelle

-->x = poly(0, "x"); q = (1 + x) / (1 - x)
 q  = 

    1 + x
    -----
    1 - x

On peut aussi utiliser directement les constantes spéciales %s et %z, qui sont l'équivalent de %s = poly(0, "s") et %z = poly(0, "z") :

-->p = %s^ 2 + 2*%s + 1
 p  =

              2
    1 + 2s + s

Les polynôme et fractions rationnelles se manipulent ensuite de manière formelle comme les autres variables. Voici des exemples à partir des objets définis ci-dessus.

Exemple

-->1/q
 ans  =

    1 - x
    -----
    1 + x

--> p*q
 ans  =

               2    3
    1 + 3x + 3x  + x
    ----------------
         1 - x

La fonction derivat() fait la dérivée formelle du polynôme ou de la fonction rationnelle.

De manière interne, Scilab représente les polynômes et les fractions rationnelles comme des listes typées. Ainsi, avec les exemples précédents,

  • p(1) va nous renvoyer 1 + 2x + x2 ;
  • q(1) renvoit la liste des noms des champs de la liste ;
  • q(2) (ou q("num"), q.num) donne 1 + x ;
  • q(3) (ou q("den"), q.den) donne 1 - x ;

et p(a) donne une erreur quelle que soit la valeur de a (sauf 1), q(a) donne une erreur quelle que soit la valeur de a (sauf 1, 2, 3 et 4).

On ne peut donc pas utiliser un polynôme ou une fraction rationnelle comme une fonction classique. L'évaluation, c'est-à-dire l'utilisation du polynôme ou de la fraction rationnelle comme fonction numérique, se fait avec la fonction horner :

horner(f,val) calcule ƒ(val).

Exemple

-->horner(p, 2)

 ans  =

    9.

Fonctions spécifiques aux polynômes[modifier | modifier le wikicode]

La fonction roots(p) donne les racines du polynôme p. La fonction factors(p) donne la décomposition en facteurs irréductibles.

La fonction coeff(p) donne une matrice dont les coefficients sont les coefficients du polynôme p. ; on peut extraire le coefficient de degré n avec coeff(p, n).

La fonction varn(p) renvoie le nom de l’indéterminée de p (ici, x).

En fait, la commande x = poly(0, "x") définit que x est le polynôme le plus simple dont l’indéterminée est le caractère « x » et dont la racine est 0, c’est-à-dire le monôme x. La commande A = poly(2, "x") définit de même le polynôme le plus simple ayant pour racine 2, c’est-à-dire x - 2.


Note
Si l'on écrit x = poly(2, "x"), alors la variable x sera le polynôme « x - 2 ». Il faut bien distinguer la notion d'indéterminée (chaîne de caractère attachée au polynôme) et la variable x. Si l'on écrit A = poly(2, "x") puis x = 1, A sera toujours le polynôme « x - 2 », l'indéterminée ne sera pas remplacée par la valeur de la variable x.

Cette fonction permet de définir un polynôme ayant plusieurs racines, en utilisant une matrice ligne au lieu d’un simple nombre. Par exemple, A = poly([1 2], "x") définit le polynôme le plus simple ayant pour racines 1 et 2, soit le polynôme x2 - 3x + 2.

On peut aussi utiliser cette fonction pour générer un polynôme dont les coefficients sont contenus dans une matrice ligne, en ajoutant le caractère "c" à la fin des arguments. Par exemple, A = poly([a0, a1, …, an], "x", "c") définit le polynôme a0 + a1·x + … + an·x n.

Fonctions spécifiques aux fractions rationnelles[modifier | modifier le wikicode]

La fraction rationnelle q est le rapport de deux polynômes, notés q.num (numérateur) et q.den (dénominateur), qui peuvent se manipuler comme des polynômes normaux.

Exemple

-->q.num
 ans  =

    1 + x

-->q.num^2
 ans  =

              2
    1 - 2x + x

La fonction simp() fait la simplification de la fraction rationnelle.

Chaînes de caractères[modifier | modifier le wikicode]

Créer une chaîne[modifier | modifier le wikicode]

Comme indiqué précédemment, une chaîne de caractères est simplement délimitée par des guillemets simples ou doubles :

'a', "a"
Nous recommandons l'utilisation des guillemets doubles : en effet, le guillemet simple sert aussi à la transposition des matrices, ce qui peut porter à confusion.

Rappelons que pour faire figurer une apostrophe ou un guillemet dans la chaîne de caractères, il faut le doubler, '' ou "" (voir Interface utilisateur > Chaînes de caractère).

Mais une chaîne de caractères peut aussi être obtenues à partir des codes ASCII[4]. La fonction ascii() renvoit :

  • si on lui donne un caractère ou une chaîne : le code ASCII du caractère ou le vecteur contenant les codes ASCII de chaque caractère de la chaîne ;
  • si on lui donne un entier ou un vecteur d'entiers : le caractère dont le code ASCII est l'entier, ou la chaîne de caractères.
Exemple
--> ascii("a")
 ans  =
 
    97.
 
--> ascii(65)
 ans  =
 
 A

La fonction char() crée un vecteur colonne de caractères à partir d'une matrice de nombres, par exemple

-->char([64, 65])
 ans  =
 
 @A

Notons que @A est une chaîne de caractères.

-->char([64, 65 ; 66, 67])
 ans  =
 
!@A  !
!    !
!BC  !

Par contre, appliquée à une matrice de chaînes de caractères, cette fonction concatène les lignes.

-->char(["a", "b" ; "c", "d"])
 ans  =
 
!ab  !
!    !
!cd  !

Rappelons quelques caractères ASCII spéciaux utiles :

  • 9 : tabulation (horizontal tabulation, HT) ;
  • 10 : saut de ligne (line feed, LF) ;
  • 13 : retour de chariot (carriage return, CR) ;
  • les caractères affichables/imprimables vont du code 32 au code 126 :
    • les chiffres vont de 48 (pour 0) à 57 (pour 9),
    • les lettres capitales vont de 65 (A) à 90 (Z),
    • les lettres minuscules vont de 97 (a) à 122 (z).


Pour plus de détails voir : w:fr:American Standard Code for Information Interchange#Table des 128 caractères ASCII .

Rechercher, remplacer, découper une chaîne[modifier | modifier le wikicode]

La commande strsubst() fait un « rechercher/remplacer » dans une chaîne. Par exemple, pour remplacer une virgule par un point dans la chaîne de caractères X :

strsubst(X, ",", ".")

La commande length(X) indique la longueur de la chaîne X.

La commande part permet d'extraire des caractères d'une chaîne. Les positions des caractères sont décrites dans un vecteur. Par exemple, pour extraire les cinq premiers caractères d'une chaîne X :

part(X, 1:5)

D'autres commandes sont présentées dans le chapitre Programmation > Analyse et construction de chaînes de caractères.

Utiliser une chaîne[modifier | modifier le wikicode]

L'opérateur « + » permet la concaténation de deux chaînes de caractère. La fonction strcat() concatène une matrice de chaînes.

-->a = strcat(["a", "b" ; "c", "d"])
 a  =
 
 acbd

La fonction string permet de transformer un nombre, un booléen, un vecteur, une matrice, … en une chaîne de caractères, par exemple pour permettre son utilisation par une fonction d'affichage (avec la fonction xstring, input, …) ou bien pour fabriquer un nom de fichier.

La fonction evstr() permet de transformer une chaîne en une expression Scilab (on parle « d'évaluation de la chaîne »). La fonction execstr() permet d'exécuter le contenu de la chaîne.

Exemple
-->a = "1";
 
-->evstr(a)*2
 ans  =
 
    2.
 
-->b = "x"; c = "=5"; execstr(b+c)
 
-->x^2
 ans  =
 
    25.
la chaîne b+c est "x=5", cette chaîne est exécutée, ce qui signifie que la valeur 5 est attribuée à la variable x, ce que confirme le calcul fait par la suite.

La commande evstr() permet donc de transformer une chaîne de caractères — ou une matrice de chaînes — ne contenant que des chiffres en un nombre — ou une matrice de nombres. Si la chaîne de caractères mélange des chiffres et d'autres caractères, on peut alors utiliser la commande strtod() (string to double).

Variable chaîne comme matrice[modifier | modifier le wikicode]

Une variable contenant une chaîne de caractères peut être utilisée pour générer une matrice. Pour cela, on fait comme si la variable était elle-même une matrice acceptant une matrice rectangulaire ne contenant que des uns (ones()) en indice. Par exemple :

-->chaine = "a"
 chaine  =
 
 a   
 
-->chaine([1, 1 ; 1, 1])
 ans  =
 
!a  a  !
!      !
!a  a  !
 
-->chaine(ones(1, 5))
 ans  =
 
!a  a  a  a  a  !

La concaténation d'un vecteur ligne avec une chaîne permet alors d'intercaler des caractères.

-->strcat(chaine(ones(1, 5)), "b")
 ans  =
 
 ababababa

Cela permet de fabriquer des chaînes particulières de manière rapide.

Matrices[modifier | modifier le wikicode]

Scilab a été conçu pour le calcul matriciel. Les éléments des matrices peuvent être de tout type (nombre réel, nombre complexe, booléen, polynôme, fraction rationnelle, chaîne de caractères…).

Définir une matrice[modifier | modifier le wikicode]

Pour définir une matrice à partir de ses coefficients, on les place entre deux crochets […]. Les éléments d’une ligne sont séparés d’un espace ou d’une virgule, les lignes sont séparées retour à la ligne ou d’un point-virgule. À l’affichage, une matrice de nombres est représentée comme un tableau sans encadrement ; une matrice de chaînes de caractères est encadré par des points d’exclamation.

Exemple

-->[1,0,0;0,1,0;0,0,1]
ans =
 
1.    0.    0.
0.    1.    0.
0.    0.    1.

La matrice vide est notée par [].

L’expression M(i,j) désigne l’élément (i, j) de la matrice M (ligne i, colonne j).

Le caractère : (deux-points) signifie « tous les indices », par exemple M(1, :) est la première ligne de la matrice (c’est un vecteur ligne). Le caractère $ (dollar) désigne le dernier indice (ligne ou colonne) d’une matrice. L'expression i1:i2 désigne tous les indices entre i1 et i2. Par exemple, si l'on veut effectuer l'opération M(i + 1, j) - M(i, j) pour tous les indices, on peut écrire simplement

difference = M(2:$, :) - M(1:$-1, :)

(note que la fonction diff(M) donne le même résultat).

Par ailleurs, l’expression N1:N2 permet aussi de générer une matrice-ligne dont le premier coefficient est N1, le dernier est inférieur ou égal à N2, et le pas entre les coefficients est 1. L’expression N1:pas:N2 permet de générer une matrice-ligne en choisissant le pas. Par exemple

-->1.1:5.2
ans = 

1.1    2.1    3.1    4.1    5.1

--> 1:2:5
ans =

1    3    5

-->"a":"d"
ans =

abcd

On notera que dans le dernier exemple, "a":"d" ne crée pas une matrice mais une chaîne de caractères.

On peut aussi utiliser la fonction

linspace(x1, x2, n)

pour générer une matrice ligne de n valeurs espacées régulièrement entre x1 et x2.

La fonction

  • zeros(m,n) crée une matrice m×n remplie de 0 ;
  • emptystr(m, n)</codes> crée une matrice m×n de chaînes de caractère vides ;
  • ones(m,n) crée une matrice m×n remplie de 1 ;
  • eye(n,n) crée une matrice unité n×n ;
  • diag(V), ou V est une matrice ligne ou une matrice colonne, crée la matrice diagonale : matrice carrée dont les termes de la diagonale sont les éléments de V et les autres termes sont nuls.

On peut aussi passer une matrice M en paramètre de ces fonctions ; elles créent alors une matrice de même dimension que la matrice M. Par exemple, M = zeros(M) met tous les coefficients de M à zéro.

Par défaut, une matrice numérique est de type « double précision ». Si l'on n'a pas besoin d'une telle précision, c'est un gaspillage de mémoire. Pour créer une matrice m×n d'entiers non signés codés sur 8 bits (soit des valeurs allant de 0 à 255), on peut utiliser la fonction uint8() (unsigned integer) sur une matrice existante :

uint8(zeros(m, n))

Mais cela revient à créer une énorme matrice (m × n × 8 octets) pour ensuite la réduire (à m × n × 1 octets). Or, nous supposons ici que la taille de la matrice m × n est important.

L'idéal est donc de :

  • partir du nombre en double précision 0 (zéro) et le transformer en entier non signé, avec la fonction uint8() ;
  • redimensionner cette matrice 1×1 avec la fonction resize_matrix()

ce qui initialise la matrice avec des zéros :

resize_matrix(uint8(0), m, n)

Si ƒ est une fonction extérieure (c’est-à-dire par définie par deff ou par function, voir Calcul numérique > Définition de fonctions extérieures), et que x et y sont des vecteurs, alors la fonction feval permet de bâtir une matrice z = ƒ(x, y)

z = feval(x, y, f) : on a z(i, j ) = ƒ(x(i ), y(j ))

Ajouter ou supprimer des lignes, des colonnes[modifier | modifier le wikicode]

Considérons une matrice M. Comme indiqué précédemment, M(i, :) crée le vecteur ligne composé des coefficients de la ligne i, et M(:, j) crée le vecteur colonne composé des coefficients de la colonne j.

Si la matrice M est une matrice par blocs composée de quatre blocs A, B, C et D (eux-même des matrices), on peut créer M simplement par

M = [A, B ; C, D]

Si l'on veut insérer un vecteur ligne L après la ligne i de la matrice, on peut écrire

M = [M(1:i, :) ; L ; M(i+1:$, :)]

et si l'on veut supprimer la ligne i, en faisant remonter les lignes suivantes, il suffit d'y mettre la matrice vide :

M(i, :) = []

On peut procéder de même avec les colonnes.

Fonctions et opérateurs sur les matrices[modifier | modifier le wikicode]

Toutes les fonctions numériques s’appliquent à des matrices ; par exemple si M est une matrice, alors cos(M) sera la matrice dont les coefficients sont les cosinii des coefficients de M.

La fonction exp(M) calcule l'exponentielle de chaque terme. L'exponentielle de la matrice (∑(1/k!)Mk) s'obtient avec la fonction expm(M).

La fonction size(M) renvoie la taille de la matrice sous la forme d’une matrice 2×1 contenant le nombre de lignes puis le nombre de colonnes. La fonction size(M, "*") renvoie le produit des deux valeurs, soit le nombre total d'éléments dans la matrice. On peut avoir le nombre de lignes avec size(M, "r") (row), et le nombre de colonnes avec size(M, "c").

Les opérations spécifiques aux matrices classiques sont :

  • la transposition : il suffit de mettre un point suivi d'une apostrophe .’ après la matrice ; l'utilisation de l'apostrophe seule, ', fait la transposition de la matrice conjuguée (opérateur adjoint) ;
  • somme, différence et produit de matrices : +, - et * ;
  • produit, rapport et élévation à la puissance élément par élément : .*, ./ et .^ ;
  • le produit tensoriel .*. ;
  • le déterminant d’une matrice carrée M : det(M) (determ() pour une matrice de polynômes et detr() pour une matrice de fractions rationnelles) ;
  • la trace d’une matrice carrée M : trace(M) ;
  • l’inverse d’une matrice inversible M : inv(M) ;

Vecteur[modifier | modifier le wikicode]

Un vecteur est une matrice ligne ou colonne (matrice 1 × n ou n × 1) ; l’expression V(i) désigne la composante i du vecteur V.

Si l'on veut avoir la taille d'un vecteur v sous la forme d'un entier, on peut utiliser

vecteur ligne : size(v, "r")
vecteur colonne : size(v, "c")
dans les deux cas : max(size(v)) ou bien size(v, "*")

La norme d'un vecteur v s'obtient simplement par

nv = norm(v);

Si V1 et V2 sont des vecteurs colonnes, alors le produit scalaire est

V1' * V2;

si ce sont des vecteurs lignes, le produit scalaire est

V1 * V2'

Il n'existe pas de fonction « toute faite » pour calculer le produit vectoriel (cross product).

Vecteurs ou matrices utilisés comme suite de nombres[modifier | modifier le wikicode]

Un vecteur peut être utilisé pour stocker une suite de nombre. cela peut être les valeurs successives d'une suite, au sens mathématique du terme, ou bien des données recueillies expérimentalement. Quelques fonctions sont relatives à ces applications.

On peut faire :

  • la somme de tous les éléments de la matrice : sum(M) ;
  • le produit de tous les éléments de la matrice : prod(M).

La fonction max(M) renvoie la plus grand élément de la matrice, et la fonction min(M) renvoie le plus petit.

Si l'on tape [m,k]=max(A), alors Scilab met la valeur du maximum dans la variable m et sa position dans la matrice (sous la forme d'un vecteur ligne [i,j]) dans la variable k — on peut utiliser d'autres noms de variable que m et k. De même pour [m,k]=min(A).

On peut ajouter un paramètre après la matrice, sous la forme d'une caractère :

  • max(A, "c") et min(A, "c") renvoient les extrema de chaque colonne de la matrice ;
  • max(A, "r") et min(A, "r") renvoient les extrema de chaque ligne de la matrice (« r » pour row).

On peut mettre une suite de matrices de même taille dans min et max. Alors, Scilab renvoit une matrice m de taille identique dont chaque élément m(i, j ) est l'extremum des éléments (i, j ) des matrice, et éventuellement une matrice k dont chaque élément k(i, j ) contient le numéro d'ordre de la matrice (dans la liste) qui contient ce maximum.

Exemple

--> A = [1,2;3,4], B = [4,3;2,1]
  A  =

     1.    2.
     3.    4.
  B  =

     4.    3.
     2.    1. 

-->[m,k] = max(A,B)
 k  =

    2.    2.
    1.    1.
 m  =

    4.    3.
    3.    4.

La matrice m(1,1) contient « 4 » qui est le maximum entre A(1,1) et B(1,1), et k(1,1) contient « 2 » ce qui indique que c'est la seconde matrice (donc B) qui contient ce maximum.

La fonction cumsum(A) fait la somme cumulée des valeurs de la matrice A. Si A est un vecteur de dimension n, cela correspond à la somme de la série, la fonction renvoie un vecteur de dimension n contenant

.

Si A est une matrice rectangulaire m×n, la somme se fait par colonne, de haut en bas, le résultat est une matrice de même dimension contenant

La fonction cumprod() fait le produit cumulé, selon la même démarche.

La fonction gsort() effectue un tri :

gsort(A) # tri par ordre décroissant
gsort(A, 'i') # tri par ordre croissant
Pour plus de détails voir : Découvrir Scilab/Programmation#Tri, recherche et sélection.

Si l'on effectue des comparaisons sur des matrices de même dimensions, cela donne une matrice de booléens de même dimensions, correspondant à la comparaison terme à terme. Par exemple

-->A = ones(2,2)
 A  =
 
    1.    1.  
    1.    1.  
 
-->B = [1, 2 ; 3, 4]
 B  =
 
    1.    2.  
    3.    4.  
 
-->A == B
 ans  =
 
  T F  
  F F  
 
-->A < B
 ans  =
 
  F T  
  T T

On peut ainsi utiliser les opérateurs ~=, <>, <, <=, > et >=.

Les fonctions and() et or() appliquent les opérateurs respectivement ET et OU à la totalité des éléments de la matrice. Par exemple, si l'on veut déterminer si deux matrices sont différentes, on peut utiliser and(A ~= B) ou and(A <> B).

Importance des matrices dans Scilab[modifier | modifier le wikicode]

Les matrices sont très utilisées dans Scilab pour définir les paramètres. Ainsi, pour tracer une courbe, il utilise un nuage de point sous forme de matrice, ou encore, pour effectuer une boucle itérative, il prend les indices dans une matrice ligne ou colonne.

La syntaxe propre aux matrices se retrouve dans de nombreux endroits.

Date et heure[modifier | modifier le wikicode]

Pour certaines applications, il est nécessaire de manipuler les dates et heures. Scila fournit les commandes suivantes :

  • datenum() ou now : indique le numéro du jour actuel, le premier jour étant le 1er janvier de l'an 0 s'il avait existé (en fait, l'an -1, puisque le calendrier grégorien passe de l'an -1 à l'an 1) ; la partie décimale correspond à l'heure ;
  • getdate() : idem, mais avec un format différent : année, mois, numéro de la semaine dans le mois, numéro du jour d'ans l'année, numéro du jour dans la semaine (1 pour dimanche, 2 pour lundi, …), numéro du jour dans le mois, heure, minutes, secondes (entier), millisecondes ;
  • datenum([aaaa, mm, jj]) indique le numéro du jour jj du mois mm de l'année aaaa ; on peut aussi ajouter l'heure pour avoir la partie décimale : datenum([aaaa, mm, jj, hh, mm, ss]) ;
  • clock donne directement la date et l'heure au format [aaaa, mm, jj, hh, mm, ss] ;
  • datevec(numjour) : convertit le numéro du jour en un vecteur contenant l'année, le mois, le jour, l'heure, les minutes et les secondes (opération inverse de la commande précédente) ;
  • weekday(numjour) indique le jour de la semaine du jour numjour : 1 pour dimanche, 2 pour lundi, … Avec la syntaxe [N, S] = weekday(numjour), Scilab indique le numéro du jour de la semaine, dans N, ainsi que le nom du jour selon l'abréviation anglaise, dans S ;
  • eomday(année, mois) : donne le nombre de jours du mois (indiquer sous forme de nombre, de 1 à 12) et de l'année indiqué ;
  • calendar() : donne le calendrier du mois en cours ;
  • calendar(aaaa, mm) : donne le calendrier du mois mm de l'année aaaa.

Par exemple :

aujourdhui = datenum()
datevec(aujourdhui)
clock
calendar()
[N, S] = weekday(aujourdhui)
datenum([1969, 07, 21])
eomday([1969, 07])

Si l'on s'intéresse à l'intervalle entre deux dates :

  • etime(t1, t2) : indique le nombre de secondes écoulées (elapsed time) entre les dates t1 et t2, sous la forme de vecteurs [aaaa, mm, jj, hh, mm, ss] ;
  • getdate("s") : indique le nombre de secondes écoulée entre le 1er janvier 1970 et l'instant actuel ;
  • getdate(numsec) : indique le jour et l'heure, au format getdate(), à partir du nombre de secondes écoulées depuis le 1er janvier 1970 ;
  • tic() déclenche un chronomètre, et toc() renvoie le nombre de secondes écoulées depuis ce déclenchement.

Précision des calculs[modifier | modifier le wikicode]

Dans Scilab, les nombres sont stockés sour la forme de nombre décimaux de double précision (doubles), selon la norme IEEE 754[5]. Ce sont des nombres à virgule flottante (floating point) codés sur 64 bits (8 octets). Il s'ensuit que :

  • la précision d'un résultat est au mieux de 10-16 ; cette précision est indiquée par la variable %eps ;
  • la valeur absolue la plus petite que l'on puisse représenter est environ 10-307 ; un nombre ayant une valeur absolue plus petite sera considéré comme nul, 0 (soupassement, underflow) ;
  • la valeur absolue la plus grande que l'on puisse représenter est environ 10308 ; un nombre ayant une valeur absolue plus grande sera considéré comme infini, inf (dépassement, overflow).

Un certain nombre de fonction internes de Scilab utilisent des méthodes qui minimisent l'erreur de calcul, et font donc mieux qu'une codage « naïf » de ces méthode. Par exemple, pour résoudre une équation du second degré, il vaut mieux utiliser la commande roots() que calculer le résultat à partir du discriminant ; la différence ne sera pas visible pour la plupart des cas, mais peut être flagrante pour des cas extrêmes.

Lire les valeurs dans un fichier[modifier | modifier le wikicode]

Il existe plusieurs manières de lire des données écrites dans un fichier. Nous nous limitons ici à deux cas :

  • fichier CSV, par exemple créé avec un tableur comme LibreOffice Calc ou Microsoft Office Excel : M = csvRead("nom_du_fichier.csv") ;
  • fichier de texte pur (ASCII ou Unicode) dans lequel les données sont sous forme de tableau : M = fscanfMat("nom_du_fichier.txt").

Dans les deux cas, les donées sont placées dans la matrice M.


Pour plus de détails voir : Découvrir Scilab/Gestion des fichiers.

Voir aussi[modifier | modifier le wikicode]

Dans Wikipédia

Notes[modifier | modifier le wikicode]

  1. la documentation utilise le terme de « primitive » ; il y a une ambiguité en français que nous avons décidé de lever dans le présent ouvrage en changeant le terme (en anglais, le problème ne se pose pas puisque la primitive d'une fonction est appelée « antiderivative »)
  2. dans la représentation des nombres à virgule flottante IEEE 754, un NaN est un nombre dont l'exposant ne contient que des 1 et la mantisse n'est pas nulle
  3. En particulier, il faut éviter le test a == [], car la notion de « variable vide » peut évoluer. Par exemple, on peut envisager qu'il existera un jour un tableau de n lignes et 0 colonne, qui est équivalent à la variable vide d'un point de vue mathématiques mais pas à [] d'un point de vue informatique. La commande isempty() est donc plus robuste aux évolutions du logiciel [1].
  4. les commandes code2str() et str2code() permettent aussi de passer de nombres aux caractères et vice versa (il s'agit par contre d'une table de conversion propre à Scilab), mais sont déclarées obsolètes et ont été supprimées dans la version 6
  5. voir Scilab is not naive, Scilab Enterprises, décembre 2010

Interface < > Calcul numérique


Calcul numérique

Table des matièresIndex



3. Calcul numérique


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().

Le problème de la précision peut parfois se manifester sur des calculs « banals ». Par exemple, le vecteur [-5.1:0.2:5.1] devrait contenir les 51 éléments {–5,1 ; –4,9 ; … ; 4,9 ; 5,1}. Or, les cinquante additions successives de la quantité 0,2, du fait des erreurs d'arrondis, font que le dernier élément calculé est légèrement supérieur à 5,1 (il vaut environ 5,1 + 2⋅10-15) et donc la liste s'arrête à 4,9. Ce genre de problème peut se gérer avec la commande nearfloat() :

-5.1:0.2:nearfloat("succ", 5.1)

ou bien, pour la génération de vecteurs, avec la commande linspace().

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) )

Pour effectuer des tests, Scilab dispose de la constante spéciale %eps qui est égale à la plus petite valeur absolue représentable.

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
Sur Scilab.org
Autres
  • (en) QL - Quadratic Programming sur Université de Bayreuth, 2010. Consulté le 2013-01-24
  • [pdf] (en) Alan Edelman (MIT), « Open Source and Traditional Technical Computing » sur Scilabtec10, 2010-06-16. Consulté le 2017-07-01
  • [pdf] anglais David Goldberg, « What every computer scientist should know about floating-point arithmetics », dans ACM Computing Surveys, vol. 23, no 1, mars 1991 [texte intégral (page consultée le 2016-08-30)] 

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


Matrices creuses

Table des matièresIndex



4. Matrices creuses


Qu'est-ce qu'une matrice creuse ?[modifier | modifier le wikicode]

Une matrice creuse est une matrice dont seuls les éléments non nuls sont stockés en mémoire, par opposition à une matrice pleine dont tous les termes sont rangés en mémoire. Dans le cas de matrices diagonales, ou de matrices ne comportant qu'un faible nombre d'éléments non nuls, l'économie en terme de mémoire peut être considérable.

Outre cette économie en terme de mémoire, l'utilisation de matrices creuses peut accélérer fortement certains calculs.

Création d'une matrice creuse[modifier | modifier le wikicode]

La fonction sparse permet de créer des matrices creuses, ou de convertir une matrice pleine en matrice creuse.

Syntaxe Description Exemple
sparse (A)
Transforme la matrice A en une matrice creuse
-->A = [0, 1; 2, 0]
 A  =
 
    0.    1.  
    2.    0.  
 
-->sparse (A)
 ans  =
 
(    2,    2) sparse matrix
 
(    1,    2)        1. 
(    2,    1)        2.
sparse (ij, v [,mn])
Crée une matrice creuse :
  • ij : matrice à deux colonnes donnant la position des éléments non nuls ;
  • vv : vecteur donnant la valeur des éléments non nuls ;
  • nm : vecteur à deux éléments donnant la dimension de la matrice.
-->A = sparse ([1, 2; 4, 3], [1, 1])
 A  =
 
(    4,    3) sparse matrix
 
(    1,    2)        1. 
(    4,    3)        1. 
 
-->full (A)
 ans  =
 
    0.    1.    0.  
    0.    0.    0.  
    0.    0.    0.  
    0.    0.    1.
sparse ([], [], [n, m])
Crée une matrice creuse "vide" de dimension n*m. Équivalent au code matlab :
sparse (n, m)
-->A = sparse ([], [], [2, 3])
 A  =
 
(    2,    3) zero sparse matrix
 
 
-->full (A)
 ans  =
 
    0.    0.    0.  
    0.    0.    0.

Matrice diagonale[modifier | modifier le wikicode]

La fonction matlab spdiags n'existe pas sous scilab. Pour créer une matrice diagonale, il faut passer par la fonction générique sparse de création d'une matrice creuse :

sparse ([1:n; 1:n]', d)
  • n : dimension de la matrice;
  • d : vecteur des valeurs de la diagonale.

Exemple :

d = [1, 5, 2, 4]
 d  =
 
    1.    5.    2.    4.
A = sparse ([1:4; 1:4]', d)
 A  =
 
(    4,    4) sparse matrix
 
(    1,    1)        1. 
(    2,    2)        5. 
(    3,    3)        2. 
(    4,    4)        4.
full (A)
 ans  =
 
    1.    0.    0.    0.  
    0.    5.    0.    0.  
    0.    0.    2.    0.  
    0.    0.    0.    4.

Notes[modifier | modifier le wikicode]


Voir aussi[modifier | modifier le wikicode]

Liens externes[modifier | modifier le wikicode]


Calcul numérique < > Graphiques et sons


Graphiques et sons

Table des matièresIndex



5. Graphiques et sons


Ce chapitre traite de la création de graphique et de sons. Pour la création et la lecture de fichiers graphiques et sonores, voir Gestion des fichiers.

Tracé de fonction[modifier | modifier le wikicode]

Tracé 2D[modifier | modifier le wikicode]

En Coordonnées cartésiennes[modifier | modifier le wikicode]

Le tracé d’une fonction se fait en deux étapes

  1. définir l’étendue de la variable abscisse et le pas, sous la forme d’un vecteur ligne ou colonne, par exemple
    x = [début:pas:fin], ou
    x = linspace(début, fin, nombre_de_points) ;
  2. tracer la fonction avec la commande
    plot(x, f(x))
    si ƒ est la fonction.

On note qu'en fait, ƒ(x) est elle-même un vecteur. On peut de manière générale définir un vecteur des valeurs de x et un vecteur des valeurs de y, et la fonction plot(x, y) tracera le nuage de points (x(i ),y(i )).

Si ƒ est une fonction externe (par exemple définie avec deff ou function, voir Calcul numérique), alors on peut tracer directement la fonction avec

fplot2d(x, f).

On peut aussi définir le vecteur y par

y = feval(x, f),

puis tracer avec

plot(x,y).

Les fonctions plot2di, utilisées à la place de plot, permettent de faire varier l’apparence générale du tracé :

  • plot2d() : trait « normal » ; identique à plot, mais l’utilisation de marqueurs et des couleurs est différente (voir ci-après) ;
  • plot2d2() : trait en escalier ;
  • plot2d3() : tracé en barres ;
  • plot2d4() : tracé en « flèches » : les points adjacents sont reliés par une flèche.

Ces fonctions plot2di acceptent des arguments modifiant le tracé, sous la forme

plot2di(x, y, arguments).

Les arguments sont de la forme mot-clef = valeur. L’argument rect = [xmin, ymin, xmax, ymax] permet de limiter le tracé à la zone comprise dans le rectangle défini par les valeurs dans la matrice. D'autres arguments sont exposés plus loin.

Modes de tracé 2d

Exemple

clear; clf;
x = [0:0.5:2*%pi]';

subplot(2,2,1)
plot2d1(x,sin(x),rect=[0,-1,2*%pi,1])
xtitle("plot2d1")

subplot(2,2,2)
plot2d2(x,sin(x),rect=[0,-1,2*%pi,1])
xtitle("plot2d2")

subplot(2,2,3)
plot2d3(x,sin(x),rect=[0,-1,2*%pi,1])
xtitle("plot2d3")

subplot(2,2,4)
plot2d4(x,sin(x),rect=[0,-1,2*%pi,1])
xtitle("plot2d4")

Nuage de points[modifier | modifier le wikicode]

La fonction scatter(x, y) trace un nuage de points. On peut choisir les marqueurs et leur couleur, par exemple :

scatter(x, y, "blue", "+")

place des croix bleues. On peut également utiliser xrects(), xarcs() et xfarcs() pour placer respectivement des rectangles, des ellipses et des ellipses remplies. Les paramètres d'entrée sont des matrices. Pour xarcs(), c'est une matrice de type :

[x1, x2,, xn
 y1, y2,, yn
 l1, l2,, ln
 h1, h2,, hn]

où le rectangle i est défini par son point en haut à gauche (xi, yi), sa largeur li et sa hauteur hi. Par exemple, si x et y sont des vecteurs de même dimension et a est le côté du carré que l'on veut tracer, on peut utiliser

xrects([x - a/2 ; y + a/2 ; a*ones(x) ; a*ones(x)])

On note que la hauteur est comptée positive vers le bas, raison pour laquelle il faut ajouter a/2 à y (pour remonter le coin du rectangle). On peut ajouter un vecteur de couleurs de même dimension que x, que nous appelons ici (vecteur d'entiers correpondant à la carte des couleurs :

xrects([x - a/2 ; y + a/2 ; a*ones(x) ; a*ones(x)], z)

Les commandes xarcs() et xfarcs() sont similaires (elles tracent des ellipses inscrites dans les rectangles définis comme avec xrects()), mais il faut de plus indiquer l'angle de départ et l'angle de fin de l'arc, en 1/64 de degré (donc 0 et 34×360 pour une ellipse complète). Donc :

xarcs([x - a/2 ; y + a/2 ; a*ones(x) ; a*ones(x) ; zeros(x) ; 34*360*ones(x)], z)

Tracé de plusieurs courbes[modifier | modifier le wikicode]

On peut superposer plusieurs courbes. Il suffit d'avoir :

  • pour les abscisses x, un vecteur colonnes ;
  • pour les ordonnées y, une matrice Y dont chaque colonne contient les ordonnées d'une courbe.

Par exemple, si l'on a deux courbes, la matrice Y a deux colonnes, et les courbes tracées sont les courbes (x(i), Y(i, 1)) et (x(i), Y(i, 2)) :

x = (0:0.05:2*%pi)';
Y = [cos(x), sin(x)];
plot(x, Y)

On peut aussi simplement utiliser plot(x, [cos(x), sin(x)]).

Dans ce cas-là, Scilab utilise une couleur différente pour chaque courbe ; on peut inscrire une légende sous le graphique pour chaque couleur, avec la commande legend() :

legend("cosinus", "sinus");

Les légendes des axes peuvent contenir des expressions mathématiques LaTeX :

legend("$\cos \alpha$", "$\sin \alpha$");

Si l'on a de nombreuses courbes, on peut afficher les légendes sur plusieurs colonnes avec la fonction legends_multicolumn() développée par Samuel Gougeon, et téléchargeable sur le Scilab FileExchange.

Si les listes d'abscisses sont différentes, alors on répète la liste des abscisses et ordonnées :

plot(x1, f1(x1), x2, f2(x2), x3, f3(x3))

On peut définir la manière dont est tracée la courbe en mettant une chaîne de caractère après la matrice des y. Cette chaîne de caractère peut contenir :

  • un signe moins « - » pour indiquer que l'on trace la ligne en continu ;
  • un signe pour tracer les ligne en discontinu : deux moins « -- » (tirets longs), deux-points « : » (pointillets), un moins et un point « -. » (trait mixte) ;
  • un signe pour indiquer que l'on met un marqueur : un plus « + », un astérisque « * », la lettre « o » (rond), un accent circonflexe « ^ » (triangle pointe en haut), inférieur « < » (triangle pointe à gauche), supérieur « > » (triangle pointe à droite), la lettre « v » (triangle pointe en bas), la lettre « d » (losange, diamond), la lettre « p » pour une étoile à cinq branches (pentagram), la lettre « s » pour un carré (square), la lettre « x » pour des croix de saint André ;
  • une lettre indiquant la couleur : « c » pour cyan, « g » pour vert (green), « k » pour noir (black), « m » pour magenta, « r » pour rouge, « w » pour blanc (white), « y » pour jaune (yellow).

Ces paramètres ne fonctionnent qu'avec la fonction plot et pas avec les fonctions plot2d.

Exemple

x=linspace(0,%pi,50);
plot(x, sin(x), "+g", x, cos(x), "-^r")
trace un sinus en croix droites vertes et un cosinus en trait plein et avec des triangles en rouge.

Avec les commandes plot2di, on change le type de tracé (couleur, marqueur) en utilisant l’argument style = n ou n est un entier positif ou négatif (par exemple plot2d(x, y, style = 1)) :

  • un nombre négatif remplace les points par des marqueurs : une étoile pour -10, des ronds pour -9, …, des petits points pour 0 ; la liste s'obtient en tapant la commande getsymbol() ;
  • un nombre positif indique un trait plein mais de couleur déterminée ; on peut voir la correspondance entre le nombre et la couleur (carte des couleurs) grâce à la commande getcolor() (voir ci-après pour la modification de la cartes de couleur).

Avec plusieurs courbes, on place les différents styles dans un vecteur ligne, par exemple

plot2d(x, [cos(x), sin(x)], style = [-1, -2])

Notons enfin que pour Scilab, les courbes sont des objets de type polyline. On peut donc ajuster leur apparence (couleur, marqueur, type de tracé, …) en modifiant les propriétés de ces objets, voir la section Propriétés des objets.

On peut également afficher les courbes dans des sous-fenêtres séparées, voir Utilisation de sous-fenêtres.

En coordonnées polaires[modifier | modifier le wikicode]

Courbe polaire.

La fonction polarplot permet en tracé en coordonnées polaires. La fonction r(theta) se trace par la commande : polarplot(theta, r)r et theta sont deux matrices de même dimension.

Par exemple :

theta = linspace(0, 2*%pi, 20)
r = 1 + 0.2*rand(1, 20, "normal")

polarplot(theta, r)

Tracé de statistiques[modifier | modifier le wikicode]

La commande bar() permet de tracer un histogramme à barres verticales, et barh() un histogramme à barres horizontales. Les barres sont équidistantes et de même largeur. Dans sa forme la plus simple, il suffit de fournir un vecteur de valeurs à tracer ; les barres sont centrées sur les indices des valeurs dans le vecteur. Par exemple

y = [1 -3 5];
bar(y)

Trace un histogramme à trois barres centrées sur x = 1, x = 2 et x = 3, ayant pour altitude respective 1, -3 et 5. On peut préciser la largeur l sous la forme d'un nombre l, et éventuellement la couleur en donnant le nom en anglais entre guillemets :

bar(y, l, "green")

La variable y peut aussi être une matrice de m lignes et n colonnes. Alors, la commande trace m histogrammes, un par ligne. Chaque histogramme est centré sur l'abscisse correspondant au numéro de la ligne. Mais si l'on ajoute l'option stacked" (empilé), alors la commande trace un histogramme de m barres, une par ligne, et les barres correspondant aux valeurs de chaque ligne sont empilées. Par exemple

y = [1 3 5 ; 2 4 6];
bar(y, 0.3, "stacked")

La fonction histplot() analyse un jeu de données et trace l'histogramme correspondant. Si x est un vecteur, la fonction histplot(n, x), n étant un entier, va découper l’intervalle de valeurs prises par les coefficients de x en n tranches d’égale largeur, et tracer l’histogramme de répartition des valeurs selon ces tranches. Si n est un vecteur dont les coefficients sont strictements croissants, les valeurs des coefficients de n servent à déterminer les tranches.

Si x est une matrice, hist3d(x) trace un histogramme 3D tel que le parallélépipède situé en (i, j ) a pour hauteur x(i, j ). Comme pour toutes les fonctions de tracé en trois dimension (voir plus loin), on peut définir l’angle de vue avec θ et α.

Carte d’un champ[modifier | modifier le wikicode]

Considérons un champ de scalaires z défini sur un plan (x, y). On a donc :

  • un vecteur colonne x de m éléments,
  • un vecteur colonne y de n éléments et
  • z une matrice m×n.

Ainsi, au point de coordonnées (x(i), y(j)) est associée la valeur z(i, j).

Cartographie couleur[modifier | modifier le wikicode]

La fonction

grayplot(x, y, z)

associe une couleur à chaque valeur de z et trace une carte de couleurs, chaque point (x(i ), y(j )) ayant la couleur associée au coefficient z(i, j ). Si l'on veut afficher la légende des couleurs, il faut utiliser la commande colobar() avant la commande grayplot(), en lui indiquant les valeurs extrêmes :

zmin = min(z); zmax = max(z);
colorbar(zmin, zmax)
grayplot(x, y, z)

Notez que la barre de couleurs est un objet de type composé (compound) encapsulé dans un objet de type axes (axes). Le type axes possède une propriété title (objet de type étiquette, label), qui permet de mettre un titre à la barre :

colorbar(zmin, zmax)
e = gce();
e.parent.title.text = "Série 1"

(voir plus bas Propriétés des objets).

La commande

Sgrayplot(x, y, z)

fait un lissage (smoothing) des valeurs pour avoir des dégradés de couleurs plutôt qu'un pavage. La fonction

Sfgrayplot(x, y, f)

fait une cartographie couleur lissée de ƒ(x, y).

Les fonctions Sgrayplot() et Sfgrayplot() font en fait appel à la fonction fec(). Cette fonction permet de faire une interpolation des valeurs et donc de représenter un dégradé, on considérant une maillage triangulaire comme ceux utilisés dans la méthode des éléments finis (FE, finite elements) :

  • on a n nœuds, le nœud i étant défini par les coordonnes x(i), y(i), z(i) ; x, y et z sont donc des vecteur de dimension n ;
  • un nombre t de mailles triangulaires, une maille étant définie par son numéro j, les numéros i des trois nœuds la composant, et un entier qui n'est pas utilisé ici (utiliser la valeur 1) ; T est donc une matrice de dimension t×5 ;

l'affichage « vu de dessus » se fait par

fec(x, y, T, z);

La commande

Matplot(A)

fait une carte couleur des valeurs de la matrice A. L'abscisse x est le numéro de colonne, l'ordonnée y est le numéro de ligne, la première ligne étant en haut, respectant en cela l'organisation du tableau. Cette fonction utilise la couleur correspondante à la partie entière de la valeur.

Les niveaux de couleur sont indiqués par la fonction

set(gcf(), "color_map", cmap)

cmap est une matrice de trois colonnes dont chaque ligne contient la couleur associée à un niveau, sous la forme RVB (les éléments de la matrice allant de 0 à 1). La première ligne de la matrice correspond au plus bas niveau, la dernière ligne au plus haut.

Exemple de codage RVB de couleurs
Couleur Ligne de
la matrice
Rouge   [1, 0, 0]
Jaune   [1, 1, 0]
Vert   [0, 1, 0]
Cyan   [0, 1, 1]
Bleu   [0, 0, 1]
Magenta vidéo   [1, 0, 1]
Magenta imprimerie   [0.859, 0, 0.451]
Blanc (100 %)   [1, 1, 1]
Gris clair (75 %)   [0.75, 0.75, 0.75]
Gris moyen (50 %)   [0.5, 0.5, 0.5]
Gris foncé (25 %)   [0.25, 0.25, 0.25]
Noir (0 %)   [0, 0, 0]

Cette matrice peut être générée de manière automatique par les fonctions …colormap(n), où n est le nombre de couleurs (entier) :

get(sdf(), "color_map") : carte de couleurs par défaut ;
coolcolormap(n) : dégradé entre le cyan et le magenta ;
coppercolormap(n) : dégradé entre le noir et le cuivre clair (codage RVB [1, 0.8, 0.5]) ;
graycolormap(n) : niveaux de gris entre le noir et le blanc ;
hotcolormap(n) : dégradé entre le rouge et le jaune ;
hsvcolormap(n) : dégradé entre suivant le cercle chromatique (hue, saturation, value) : rouge, jaune, vert, cyan, bleu, magenta, et retour au rouge ;
jetcolormap(n) : dégradé entre le bleu et le rouge ;
oceancolormap(n) : dégradé de bleus ;
pinkcolormap(n) : sépia ;
rainbowcolormap(n) : dégradé selon les couleurs de l'arc-en-ciel (rouge, orange, jaune, vert, bleu, violet) ;
springcolormap(n) : dégradé entre le magenta et le jaune ;
summercolormap(n) : dégradé entre le vert et le jaune ;
whitecolormap(n) : blanc uniforme (matrice remplie de 1, équivalent à ones(n, 3)) ;
wintercolormap(n) : dégradé entre le bleu et le vert.

On peut par exemple utiliser

set(gcf(), "color_map", graycolormap(32))

pour avoir 32 niveaux de gris. On peut réaliser un dégradé du noir vers le rouge avec

cmap = graycolormap(32); cmap(:, 2:3) = 0;
set(gcf(), "color_map", cmap)

ou avec

r = linspace(0, 1, 32)'; cmap = [r zeros(32, 2)];
set(gcf(), "color_map", cmap)

et un dégradé du bleu vers le blanc avec

cmap = graycolormap(32); cmap(:, 2:3) = 1;
set(gcf(), "color_map", cmap)

ou avec

r = linspace(0, 1, 32)'; cmap = [r ones(32, 2)];
set(gcf(), "color_map", cmap)

Les niveaux de couleur sont également utilisés lorsque l’on trace plusieurs courbes sur le même graphique. Dans ce cas-là, des dégradés ne fournissent pas un contraste permettant de distinguer facilement des courbes voisines. La carte par défaut est bien conçue pour ce cas-là, et on peut la restaurer par

cmap = get(sdf(), "color_map");
set(gcf(), "color_map", cmap)

Courbes de niveau[modifier | modifier le wikicode]

On peut tracer des courbes de niveau avec la fonction

contour2d(x, y, z, n)

n est le nombre de niveaux que l’on veut voir figurer. On peut aussi donner les valeurs des niveaux z1, z2, …, zn par un vecteur

contour2d(x, y, z, [z1,z2,...,zn])

Champ de vecteurs[modifier | modifier le wikicode]

On peut également tracer un champ de vecteurs, sous la forme de flèches. Pour cela, il faut un vecteur colonne vx ayant les composante selon x du champ de vecteur, un vecteur colonne vy ayant les composantes selon y de ce champ, et utiliser la fonction

champ(x, y, vx, vy)

le vecteur de composantes (vx(i,j ),vy(i,j )) étant placé en (x(i ), y(j )).

Avec la fonction champ1, les vecteurs tracés ont tous la même longueur, la norme du champ est indiquée par la couleur du vecteur, suivant le principe exposé pour grayplot.

Tracé 3D[modifier | modifier le wikicode]

Scilab permet également le tracé de surfaces à trois dimensions. Il s'agit en général de nuages de points de la forme(x, y, z).

De même que les courbes sont des objets « polyline », les surfaces sont des objets « fac3d » (facettes à trois dimensions).

Tracé de fonctions de deux variables[modifier | modifier le wikicode]

Nous voulons tracer la surface représentative d’une fonction ƒ, on a

z = ƒ(x, y).

Si x est une matrice colonne de m éléments, y une matrice colonne de n éléments, alors z est une matrice m×n définie par

z(i, j ) = ƒ(x(i ), y(j ))

alors la fonction

plot3d(x, y, z)

va tracer la surface générée par le nuage de points (x(i ), y(j ), z(i, j )).

Mise en œuvre

Partons de deux vecteurs-ligne x et y, de taille possiblement différente. On peut les transformer en matrices de même dimension par la multiplication avec un vecteur de uns. Par exemple
x = -2:0.1:2;
y = -3:0.1:3;

X = ones(y')*x;
Y = y'*ones(x);
Z = X.^2 + Y.^2;

plot3d(X, Y, Z)
Ici, la matrice X est une matrice colonne pour laquelle chaque valeur de x est répétée n fois (n étant la taille du vecteur y), la matrice Y est une matrice colonne pour laquelle la suite des valeurs de y est répétée m fois (m étant la taille du vecteur x), et la matrice Z résulte donc de toutes les combinaisons possibles de couples (xi, yj ).