Fermer X

Introduction à Scilab
Graphiques, fonctions Scilab, programmation, saisie de données

Jean-Philippe Chancelier & Michel De Lara
cermics, École des Ponts ParisTech
(last modification date: April 12, 2018)
Version pdf de ce document
Version sans bandeaux

Contents

Ouvrir une fenêtre Scilab

Pour ces travaux pratiques d'introduction à Scilab, il vous faut lancer le logiciel Scilab et disposer ainsi d'une fenêtre permettant de saisir et d'exécuter des instructions.

Taper des instructions Scilab

Dans ces premiers travaux pratiques, vous trouverez une série de lignes de commandes Scilab précédées du signe -->. Pour commencer, il vous suffit de les recopier ou de les saisir par copier-coller (sans -->) pour les exécuter immédiatement dans la fenêtre Scilab.

Charger des instructions Scilab

Pour des successions de lignes de commandes, il est préférable de disposer d'un éditeur (Emacs), d'ouvrir un fichier (par exemple nom_du_fichier.sce), d'y écrire les lignes de commandes puis de les exécuter en tapant sous Scilab l'instruction :

 exec("nom_du_fichier.sce")

On peut également utiliser la rubrique File Operations du menu File.

Commentaires

Toute ligne débutant par // est une ligne de commentaires.

1 Graphiques 2D

Pour effectuer des graphiques en deux dimensions, on utilise la commande plot2d.

 
x=[1,3,4,9];  
y=[2,7,-6,-3];  
plot2d(x,y)  
// relie les points de coordonnées (x(1),y(1)), (x(2),y(2)), etc.  

PIC

Tracé d'une unique courbe
x=0.1:0.1:10;  
y=sin(x)./x;  
xbasc(); // efface le contenu de la fenêtre graphique  
plot2d(x,y)  

PIC

Tracé de deux courbes (par superposition)
 
x=0.1:0.1:10;  
y=sin(x)./x;  
z=cos(x)./x;  
xset("window",1);  
// ouvre une nouvelle fenêtre 1  
// (la fenêtre par défaut est 0)  
plot2d(x,y);plot2d(x,z) // les deux graphiques sont superposés  

Tracé de deux courbes (simultanément)
 
x=x';y=y';z=z';  
// x, y et z sont maintenant des vecteurs colonnes  
// ce qui est obligatoire quand on veut dessiner plus de deux courbes.  
xbasc();plot2d(x,[y z])  
// idem, mais les deux graphiques sont tracés chacun avec sa couleur  

PIC

Tracé de plusieurs courbes (simultanément)
 
M=[(0:0.1:10)' (5:0.1:15)' (10:0.1:20)'];  
N=[0.1*M(:,1).^2 log(M(:,2)) 0.5*M(:,3)];  
xbasc();plot2d(M,N);  

Créneaux plot2d2, barres plot2d3
 
x=1:10;  
xbasc();plot2d2(x,x);  

PIC

 
xbasc();plot2d3(x,x);  

PIC

Titre xtitle, légende legends
 
 annees=[1985 1986 1987 1988 1989 1990 1991 1992];  
 annee_ref=2004 ;  
// les prix seront tous exprimés en euros de cette année  
 prix_annee_ref=[1.27 1.05 1.03 1.00 1.04 1.04 1.01 0.96];  
 prix_courant=[0.85 0.72 0.73 0.73 0.78 0.81 0.81 0.79];  
 xbasc();plot2d2(annees',[prix_annee_ref ; prix_courant]');  
 xtitle('Litre de supercarburant ...  
en euros de l''annee ' +string(annee_ref),...  
'(annees)','(prix)')  
// Noter le '' pour signifier l'apostrophe dans l''annee.  
// Dans une instruction trop longue pour tenir dans une ligne,  
// mettre ... avant de passer à la ligne  
 legends(['prix en euros '+string(annee_ref);'prix courant'],[1,2],'ur')  

PIC

Exporter un graphique

Dans le menu de la fenêtre graphique, sélectionner File, puis Export pour exporter le graphique (au format Postscript par exemple, ou au format Postscript-Latex qui permet une insertion facile dans un fichier LATEX).

2 Graphiques 2D : options

La fonction Scilab plot2d a de nombreux arguments optionnels permettant de spécifier des attributs d'un graphique : couleur, style et épaisseur des traits, échelle, dimensions du cadre, etc. Pour découvrir ces arguments optionnels, taper help plot2d.

Pour découvrir les attributs associés aux objets graphiques, utiliser le menu Edit de la fenêtre graphique. À travers ce menu, on peut également modifier en temps réel des attributs.

À titre indicatif, noter que les manipulations des attributs d'un graphique peuvent également être effectuées au moyen de commandes Scilab (taper help graphics_entities).

Nous allons illustrer des manières de changer les attributs d'un graphique 2D.

Option style de plot2d

L'option style permet notamment de changer la couleur.

x=(1:100)';  
y=rand(100,1);  
xbasc();plot2d(x,y)  
// Par défaut, les points sont reliés par des droites de couleur noire  
xbasc();plot2d(x,y,style=5)  
// Avec l'option style=5, on obtient du rouge.  
// Le code par défaut des couleurs se découvre par la commande getcolor()  

PIC

Ce changement de couleur peut également être obtenu en utilisant le menu Edit de la fenêtre graphique.

Question 1 Essayez, à l'aide du menu Edit, de modifier la couleur et de rajouter des losanges aux points de la courbe. Pour cela, on cherchera Polyline dans le menu déroulant des entités graphiques.

Option rect de plot2d

L'option rect permet de fixer les échelles.

 
y=[6:-1:0,1:5] ; x =[y,y,y,y,y];  
xbasc();plot2d(1:60,x);  
// Un graphique simple dans la fenêtre courante.  

PIC

 
xset('window',2);  
// Ouvre une deuxième fenêtre graphique.  
xbasc();plot2d(1:60,x,rect=[0,0,70,10]);  
// Les échelles sont fixées par rect =[xmin,ymin,xmax,ymax]  

PIC

Question 2 Effectuer la même manipulation à l'aide du menu Edit.

3 Saisie de tableaux de chiffres

Lecture de données : le pouvoir d'achat du franc

Le tableau de chiffres suivant représente l'indicateur de pouvoir d'achat du franc de 1985 à 2001, exprimé en euros 2005 (source Insee). Il permet de traduire en euros 2005 des valeurs exprimées en francs du passé. Sur une ligne, on trouve

1.
l'année ;
2.
la valeur équivalente en euros, en pouvoir d'achat de 2005, d'1 FF de l'année.

Ainsi, 1 000 F de 1985 équivalent, en pouvoir d'achat, à 231,36 euros de 2005.

 
2001 0.16492  
2000 0.16766  
1999 0.17050  
1998 0.17135  
1997 0.17254  
1996 0.17466  
1995 0.17811  
1994 0.18119  
1993 0.18421  
1992 0.18804  
1991 0.19249  
1990 0.19866  
1989 0.20536  
1988 0.21277  
1987 0.21849  
1986 0.22537  
1985 0.23136  

Le tableau de chiffres suivant représente l'indicateur de pouvoir d'achat de l'euro de 2002 à 2005, exprimé en euros 2005 (source Insee). Il permet de traduire en euros 2005 des valeurs exprimées en euros du passé.

 
2005    1.000  
2004  1.018  
2003  1.040  
2002  1.061  

Par copier-coller, transférer les deux séries de lignes précédentes dans deux fichiers francs.txt et euros.txt. On prendra garde à ce que le fichier finisse par un retour à la ligne, sinon la dernière ligne ne sera pas lue par fscanfMat. On va à présent lire les contenus de ces fichiers et stocker les résultats de la lecture dans deux matrices Mf et Me (à 2 colonnes et n lignes). Pour cela utiliser help fscanfMat :

 
Mf= fscanfMat('francs.txt');  
size(Mf)  
Mf(1,:)  
Me= fscanfMat('euros.txt');  

Question 3 À partir des matrices Mf et Me, construire une matrice N  dont la première colonne soit les années de 2005 à 1985, et la deuxième l'indicateur de pouvoir d'achat de l'euro exprimé en euros 2005. On rappele que 1 EUR= 6,55957 FF.

En utilisant M=N($:(-1):1,:), transformer cette matrice en une dont la première colonne soit les années de 1985 à 2005.

Lecture de données : le prix des carburants
 
1985 85.737913076923 65.298018096154  
1986 72.559868307692 51.697806788462  
1987 73.706167788462 50.176247980769  
1988 73.489055150943 48.769303018868  
1989 78.860116961538 51.618650134615  
1990 81.440023403846 54.186829923077  
1991 81.495726230769 54.59726975  
1992 79.991758057692 52.726837153846  
1993 82.922916660377 55.966047339623  
1994 85.538556557692 58.813071173077  
1995 89.112313596154 58.692871  
1996 94.692211192308 65.334312673077  
1997 98.006277230769 67.640015865385  
1998 95.804121788462 64.239142269231  
1999 100.24807667711 68.998516558797  
2000 116.88560606538 84.682718019723  
2001 111.85613846808 79.601834232065  
2002 109.5216625 77.242352692308  
2003 109.60326923077 79.347115384615  
2004 114.0958490566 88.470943396226  
2005 127.46173076923 102.68038461538  

Question 4 Le tableau précédent représente, en euros courants, les prix des carburants (source direction générale de l'énergie et des matières premières, ministère de l'Industrie) en hectolitres. On touve sur une ligne : année, supercarburant, gazole.

Par copier-coller, transférer la série de lignes précédente dans un fichiers carburants.txt. Saisir le tableau dans une matrice Cc. Convertir ces prix, exprimés en euros courants,en euros de 2005 et les ramener de l'hectolitre au litre. Tracer sur un même graphique les évolutions des prix (en euros de 2005) des deux carburants de 1985 à 2005.

PIC

4 Programmation

Scilab est un langage interprété où la transmission des arguments se fait par valeur, même pour les objets de type matrices. Il faut donc éviter les boucles qui peuvent être très inefficaces en terme de temps de calcul. Pour cela, il faut utiliser si possible des opérateurs vectoriels qui font la même chose.

Nous recommandons ici de disposer d'un éditeur (Emacs), d'ouvrir un fichier (par exemple nom_du_fichier.sce), d'y écrire les lignes de commandes puis de les exécuter en tapant sous Scilab exec("nom_du_fichier.sce").

  • condition if
       p=0.4  
       y=rand(1)  
       if y< p then x=1;  
       elseif y>= p then x=0;  
       end  
       x  

  • boucles while
       p=0.2  
       y=rand(1);  
       compteur=1;  
       while y >=p do compteur=compteur+1; y=rand(1); end  
       compteur  

  • boucles for
       N=500;  
       y=rand(1,N)-0.5*ones(1,N);  
       x=sign(y);  
       for i=2:N,x(i)=x(i-1)+x(i);end;  
       xbasc();plot2d(1:N,x)  
       x1=cumsum(sign(y));  
       norm(x1-x)  

5 Fonctions Scilab

On peut définir une fonction Scilab de façon interactive, mais il est souvent plus naturel d'écrire le code d'une fonction dans un fichier au moyen d'un éditeur de texte.

On peut mélanger des instructions Scilab et des définitions de fonctions dans un même fichier, qu'on fera terminer par l'extension .sce.

L'entête d'une fonction est constituée de la séquence
function [<arguments de retour>]=<nom>(<arguments d'entrée>)
et le corps de la fonction est terminé par le mot clef endfunction.

   function [y]=cube(x), y=x.^3, endfunction;  
   // le nom de la fonction est cube  
   // elle a pour argument un vecteur x, et retourne le vecteur y=x.^3

On notera l'utilisation de l'opérateur .^ (et non pas simplement ^) car l'argument x  est généralement un vecteur (et non pas un scalaire).

   x=0.01:0.01:0.99;y=cube(x);  
   xbasc(); plot2d(x,y);

Voici un autre exemple de fonction.

   function [H] = Heavyside(x)  
     // Fonction de Heavyside  
     H = bool2s(x>=0)  
     // on notera que Heavyside(0)=1  
   endfunction  
 
   x=-1:0.2:1; xbasc(); plot2d2(x,Heavyside(x),rect=[-1,-0.1,1,1.2]);  

Question 5 Écrire la fonction sinus cardinal (sin(x)∕x  ) de telle sorte qu'elle puisse admettre un argument x  vectoriel et retourne un vecteur de même taille contenant le sinus cardinal de chacun de ses éléments. Deux problèmes se posent : comment faire une division vectorielle ? comment éviter une divison par zéro ?1

Une fonction dans Scilab peut avoir plusieurs arguments, scalaires ou vectoriels.

 
   function [y]=Gauss(x,mu,sigma)  
     // mu et sigma sont des scalaires  
     // x peut être un vecteur de taille quelconque  
     y=1 ./(sigma*sqrt(2*%pi))*exp(-((x-mu)/sigma).^2/2)  
   endfunction;  
 
   function [Y]=Normal(X,Mu,Sigma)  
     // Mu est un vecteur ligne de taille [1,dim]  
     // Sigma est une matrice symétrique définie positive de taille [dim,dim]  
     // x peut seulement être un vecteur de taille [1,dim]  
     [lhs,dim]=size(Mu)  
     Y=1/ sqrt( (2*%pi)^dim * det(Sigma) ) * ...  
           exp(- (X-Mu)*inv(Sigma)*(X-Mu)' / 2 )  
   endfunction;

Une fonction dans Scilab peut renvoyer plusieurs valeurs.

 
   function [H,plus] = Heavyside_plus(x)  
     // Fonction de Heavyside  
     H = maxi(sign(x),0);  
     // on notera que Heavyside(0)=0  
     plus=maxi(x,0)  
     // Fonction +  
   endfunction

Si on tape simplement

 Heavyside_plus(-1:0.2:1)

seule la première valeur est retournée. Pour disposer de toutes les valeurs, il faut utiliser la syntaxe suivante qui attribue des noms à ces valeurs

   [a,b] = Heavyside_plus(-1:0.2:1)

Les deux fonctions ci-dessous retournent une permutation de n  éléments tirée au hasard. Interpréter le mécanisme de simulation de la première des deux. Il apparaît que cette méthode est plus lente que l'appel à la fonction Scilab grand.

function [x]=permutation_lent(n)  
// permutation aléatoire de {1,...,n}  
// version avec des boucles  
  x=ones(1,n);  
  for i=2:n,k=0;  
    for j=1:1+int(rand(1)*(n-i+2));k=k+1;  
      while x(k)>1 do  
        k=k+1;  
      end;  
    end;  
    x(k)=i;  
  end;  
endfunction  
 
function [x]=permutation_rapide(n)  
  // permutation aléatoire de {1,...,n}  
  x=grand(1,'prm',[1:n]')  
endfunction

// chargement de la fonction permutation_lent :  
-->exec('nom_du_fichier.sce')  
// la fonction est maintenant connue :  
-->permutation_lent  
// appel de la fonction permutation_lent :  
-->permutation_lent(100)

6 Application : calculs d'emprunts

On considère un emprunt

  • d'un montant K  , capital emprunté ;
  • au taux annuel θ ∈ [0,1]   ;
  • remboursé en T  périodes (T  = 12Na  pour un remboursement mensuel sur Na  années) ;

et on veut connaître

  • le montant R  du remboursement par période ;
  • le montant total I  des intérêts ;
  • pour chaque période t = 1,...,T
    • Kt  le capital restant dû ;
    • At  le capital amorti ;
    • It  les intérêts dus.

Le taux actuariel τ  est (dans le cas de remboursements mensuels)

soit proportionnel  τp =  θ∕12                                     (1)
   soit ˊequivalent   τ =   (1 + θ )1∕12 − 1 (τ >  τ).                (2)
                     e                      e    p

On a les relations suivantes (t = 1,...,N  ) :

     It =   τKt                                    (3)

  Kt+1  =   Kt − At                                (4)
It + At =   R                                      (5)
Avec K1 =  K  , capital emprunté, on en déduit pour t = 1, ...,T   :
Kt  =   R ∕τ − (1 + τ)t− 1(R ∕τ − K )                        (6)
                   t−1
 It =   R −  (1 + τ ) (R −  τK )                           (7)
At  =   (1 + τ)t−1(R − τK ).                               (8)

Introduisons la fonction suivante (telle que ∂Φ
-∂τ ≥ 0  ) :

Φ (τ,T ) = ------τ------.
           1 − (1 + τ )− T
(9)

Avec K1  = K  et KT +1 = 0  (plus de capital restant dû après la dernière échéance), on obtient alors facilement :

          --τ(1 +-τ)T--
 R  =   K (1 + τ)T − 1 =  K Φ(τ,T )                            (10)

                 T
K   =   R (1-+-τ)--−-1-=  --R-----                             (11)
           τ (1 + τ )T      Φ(τ,T )

        ∑T
 I  =       It = T R − K =  K (TΦ (τ,T) − 1).                  (12)
        t=1
 
function [mensualite,restant,amorti,interets]=...  
         emprunt(capital,taux_annuel,annees)  
  // On considère ici un emprunt de capital  
  // à remboursements mensuels sur (12 * annees) périodes  
  // au taux annuel taux_annuel  
  if taux_annuel > 0.99 then  
    error('Attention : mettre un taux de 4,5 pourcent sous la forme 0.045')  
  end  
  T=12*annees;  
  tau=(1+taux_annuel).^(1/12)-1;  
  // taux actuariel équivalent  
  phi=tau/(1-(1+tau).^(-T));  
  // fonction phi  
  mensualite=capital*phi  
  // mensualité  
  periodes=[1:12*annees];  
  // périodes de remboursement  
  amorti=((1+tau).^(periodes-1)).*(mensualite-capital*tau);  
  // capital amorti : vecteur de taille (1,12*annees)  
  interets=mensualite-amorti;  
  // intérêts : vecteur de taille (1,12*annees)  
  restant=interets/tau;  
  // capital restant : vecteur de taille (1,12*annees)  
endfunction  

Question 6 Copier le code Scilab ci-dessus dans un fichier calculs_emprunt.sce et charger la fonction emprunt par la commande exec("calculs_emprunt.sce"). Choisir des paramètres d'emprunt, puis étudier l'effet d'une variation de 1 % du taux annuel. Modifier le code précédent pour qu'il fournisse également en sortie le montant total des intérêts.

xbasc();  
capital=100000;  
taux_annuel=0.05;  
annees=10;  
[mensualite,restant,amorti,interets]=emprunt(capital,taux_annuel,annees);  
periodes=[1:12*annees];  
plot2d(periodes,restant);  
[mensualite_1,restant_1,amorti_1,interets_1]=...  
        emprunt(capital,taux_annuel+0.01,annees);  
plot2d(periodes,restant_1);  
 
sum(interets_1)/sum(interets)-1  
// variation des intérêts versés

Question 7 Programmer une fonction capital_emprunt de la forme [capital]=capital_emprunt(mensualite,taux_annuel,annees), donnant le capital pouvant être emprunté en fonction de la mensualité, du taux et du nombre d'années.

Si l'on veut calculer, pour un plan d'épargne logement (pel), le capital que l'on peut emprunter et les mensualités, il faut procéder comme suit. Un pel terminé se caractérise par

  • des droits à prêt D  (intérêts acquis sur le pel) ;
  • le taux actuariel τ
 el   du pel.

Le capital empruntable K  dépend du nombre Na  d'années d'emprunt : il est tel que les intérêts perçus I  (par le prêteur pour un prêt au taux τ
 el  ) soient égaux à 2,5  fois les intérêts acquis (sur le pel), soit

I = 2,5 × D
(13)

(2,5 doit être remplacé par 1,5 pour un livret ou compte d'épargne). On en déduit

T =  12 × N a  et  K  = 2,5 ×  ------D-------.
                               TΦ (τel,T ) − 1
(14)

Enfin, le capital est plafonné à Kmax  (92 000 euros en 2003) et le prêt est accordé à un taux (actuariel) τ >  τel   :

K   =  min (K    ;2,5 × ------D------- )                    (15)
              max       T Φ (τel,T ) − 1

R   =  K Φ (τ,T ).                                          (16)

Question 8 Programmer une fonction PEL_emprunt de la forme [capital,mensualite]=PEL_emprunt(droits_a_pret,taux_annuel,annees), donnant le capital pouvant être emprunté et le montant des remboursements, en fonction des droits à prêt, du taux annuel du pel et du nombre d'années d'emprunt.

7 Graphiques 3D

Question 9

  • Saisir le tableau suivant sous forme de matrice
    M=[ 10.66   13.75   16.28   18.17   19.34   19.73   19.34   18.17   16.28   13.75   10.66  
        11.31   14.59   17.28   19.28   20.52   20.93   20.52   19.28   17.28   14.59   11.31  
        10.99   14.17   16.79   18.73   19.93   20.34   19.93   18.73   16.79   14.17   10.99  
         9.00   11.61   13.75   15.34   16.33   16.66   16.33   15.34   13.75   11.61    9.00  
         4.66    6.01    7.12    7.95    8.45    8.63    8.45    7.95    7.12    6.01    4.66  
        -2.01   -2.60   -3.08   -3.43   -3.65   -3.73   -3.65   -3.43   -3.08   -2.60   -2.01  
        -8.98  -11.58  -13.72  -15.31  -16.29  -16.63  -16.29  -15.31  -13.72  -11.58   -8.98  
       -11.02  -14.21  -16.83  -18.78  -19.99  -20.39  -19.99  -18.78  -16.83  -14.21  -11.02  
        -2.63   -3.39   -4.01   -4.48   -4.76   -4.86   -4.76   -4.48   -4.01   -3.39   -2.63  
        10.14   13.08   15.49   17.29   18.40   18.77   18.40   17.29   15.49   13.08   10.14  
         4.42    5.69    6.75    7.53    8.01    8.17    8.01    7.53    6.75    5.69    4.42  
       -11.34  -14.63  -17.33  -19.33  -20.57  -20.99  -20.57  -19.33  -17.33  -14.63  -11.34 ]

  • Faire help plot3d.
  • On suppose que le terme M(i,j) représente la cote z du point (i,j,z). Tracer la surface correspondante.
  • Faire help contour et tracer les courbes de niveau.

PIC
L'École des Ponts ParisTech est membre fondateur de

L'École des Ponts ParisTech est certifiée 

ParisTech ParisEst ParisTech

 

Copyright 2014
École des Ponts ParisTech
All rights reserved