Introduction à Scilab
Semaine logiciels, année 2000-2001
vendredi 20 octobre 2000

Jean-Philippe Chancelier & Michel Cohen de Lara

\fbox{\textbf{http://clubinfo.enpc.fr/scilab/}}


Table des matières

9h00 $\longrightarrow$ 10h15. Manipulations vectorielles

Les lignes de commande suivantes peuvent être directement tapées sous Scilab pour être exécutées immédiatement.

Elles peuvent également être écrites dans un fichier nom_du_fichier.sce, puis exécutées en tapant sous Scilab exec("nom_du_fichier.sce","c").

Objets vectoriels 

// scalaires
5
ans^2
l=100
y=%e
l*log(sqrt(y))

// vecteurs
v=[3,4]
norm(v)
size(v)

w=[%pi,v]'
sin(w)
sum(w)

t=[4:9]
u=sqrt(2)*[1:2:8]'
size(u)
s=ones(u)

// matrices
A=ones(3,7)
B=eye(6,7)
C=[A;(-6)*B]

Exercice

Construire une matrice (4,9) (à 4 lignes et 9 colonnes) dont la première ligne est formée de 1, et dont tous les autres termes sont nuls.

Construire une matrice (3,5) dont la première colonne est formée de 2, la deuxième colonne des entiers de 1 à 3, et le reste de -1.

Opérations vectorielles 

//addition
A=rand(5,6)
B=sin(A)
A+B

//transposition
A'

// rang
rank(A)

//multiplication
A'*A
A*A'
C=eye(A)
A*C

//déterminant
A=[1,2;3,4]
det(A)

//exponentiation
D=rand(3,3)
expm(D)

//extraction de sous-matrices
A=[11:19;21:29;31:39;41:49;51:59;61:69]
A(1,1)
A(3,4)
A(1,:)
A(:,5)
A(2:3,:)
A(2:3,7:9)
A([1 3 5],[2 4 6 8])

Exercice

Cliquer sur la touche Help

Taper \ dans la case A propos

Choisir des matrices $A$ et $b$

Résoudre numériquement le système linéaire $Ax=b$

Recommencer avec / et le système linéaire $xb=A$

Valeurs propres et vecteurs propres 

//valeurs propres
A=eye(5,5)
spec(A)
A=[1,2;3,4]
spec(A)

//vecteurs propres
[Ab,X]=bdiag(A)
(A - Ab(1,1)*eye(A))*X(:,1)
(A - Ab(2,2)*eye(A))*X(:,2)
inv(X)*A*X

Exercice

Diagonaliser une matrice rand(4,4).

10h45 $\longrightarrow$ 12h00. Fonctions Scilab, graphiques, etc.

Graphiques 

// plot
x=0.1:0.1:10
y=sin(x)./x
plot(x,y)

// plot2d
r=rand(1,100)
plot2d(r)
xbasc()
plot2d(1:100,r,-5)
s=rand(1,100);
xbasc()
plot2d(r,s,-2)

// histogrammes
xbasc()
histplot(10,r)
xbasc()
histplot(0:0.2:1,r)

Exercice

Engendrer un vecteur aléatoire $r$ de taille 200.

Faire help cumsum

Créer, à partir de la suite $r$, la suite $s$ des sommes partielles.

Représenter les couples de points $(r,s)$.

Programmation 

 // 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(0:N,[0,x])

// fonction (d'expression simple)
deff('[y]=densite(x)','y=1/(%pi*sqrt(x*(1-x)))');
densite(0.8)
deff('[y]=Gauss(x,mu,sigma)','y=1/(sigma*sqrt(2*%pi))*...
exp(-((x-mu)/sigma)^2/2)'); 
// mettre ... avant de passer à la ligne

Exercice

Écrire la fonction sinus cardinal ($\sin x /x$).

Fonctions Scilab

Une fonction est une suite d'instructions, éventuellement paramétrées, écrites dans un fichier nom_du_fichier.sci. 

function [x]=permutation_lent(n)
// permutation aléatoire de {1,...,n}
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;

function [x]=permutation_rapide(n)
// permutation aléatoire de {1,...,n}
x=grand(1,'prm',[1:n]')

Une fonction est chargée en tapant sous Scilab getf("nom_du_fichier.sci","c"). Toute fonction est ensuite appelée sous Scilab par son nom. 

// chargement de la fonction permutation_lent
getf("nom_du_fichier.sci","c")
// appel de la fonction permutation_lent
permutation_lent(100)

Exercice

Écrire dans un fichier

- la fonction sinus cardinal ($\sin x /x$)

- une fonction d'entrée 2 réels $a <b$, qui produit son graphe sur l'intervalle $[a,b]$.

Saisie de tableaux de chiffres

Par copier-coller, transférer les lignes suivantes dans un fichier donnees.txt : 

 21.132487    56.084861    30.760907    50.153416    28.06498  
 75.604385    66.235694    93.296162    43.685876    12.800585 
 0.0221135    72.635068    21.460079    26.931248    77.831286 
 33.032709    19.851438    31.2642      63.257449    21.190304 
 66.53811     54.425732    36.16361     40.51954     11.213547 
 62.839179    23.207479    29.222666    91.847078    68.56896  
 84.974524    23.122372    56.642488    4.3733433    15.312167 
 68.573102    21.646326    48.26472     48.185089    69.708506 
 87.821648    88.338878    33.217189    26.39556     84.155184 
 6.8374037    65.251349    59.350947    41.481037    40.620248


// fscanfMat
// que fait la fonction fscanfMat ?
help fscanfMat               
// lecture et stockage dans une matrice M      
M= fscanfMat('donnees.txt');

Un exemple récapitulatif : analyse de données de température 

//saisie de données de température
M= fscanfMat('SCI/contrib/tp_enpc/temperatures_globe.txt')
//on extrait les moyennes annuelles
[m,n]=size(M)
y=M(:,n)
//on programme une régression linéaire
t=[1:m]'; // représente le temps (années)
//formule de régression
x=[ones(M),t]; b=(x'*x)^(-1)*x'*y;
//superposition graphique des données et de la droite de régression
xbasc();plot2d(t',y');      
plot2d(t',(x*b)');
//examen des résidus
r=y-x*b
plot(r)
histplot(r)

Un exemple récapitulatif : analyse de données de débit

Les debits moyen journaliers de la Gervanne à Beaufort sur Gervanne (en l/s) sont stockés dans un fichier.

On peut par exemple lire les données correspondante dans Scilab avec la fonction fscanfMat comme suit : 

// lecture et stockage dans une matrice M
M= fscanfMat('SCI/contrib/tp_enpc/gervanne.txt')

Combien de données avons-nous lu ? 

[m,n]=size(M)    // On a lu une seule colonne de chiffre : on doit avoir n=1
                 // m nous indique le nombre de valeurs lues 
M=M/1000;        // on se ramène a des m3/s

Visualisons les débits journaliers 

 
plot(M);         // visualisation de toutes les données 
xbasc();plot(M(1:100));  // les 100 premières valeurs

Analyse de la série temporelle des débits journaliers 

 
[cor,moy]=corr(M,100)  // les 100 premiers coefficiente de corrélation

Transformée de Fourier (discrète) du signal 

 
y=fft(M,-1);
xbasc();plot(abs(y(1:100)));

À quoi correspond le premier « pic » dans la transformée de Fourier discrète ? 

 
[m,ind]=maxi(abs(y(2:100))); // on cherche l'indice qui réalise le max 
                             // y étant complexe, on prend son module (abs)
y1=y;
y1(ind+1)=0;
y1($-ind+1)=0;               // on élimine la fréquence du signal y 
M1=real(fft(y1,1));          // transformée de Fourier inverse 
per=M-M1;                    // signal périodique (de période l'année)

xbasc();plot(per);
xbasc();plot(per(1:360));

13h00 $\longrightarrow$ 14h15. Probabilités

Générateurs aléatoires 

// rand
r=rand(1,100)
plot2d(r)
xbasc()
plot2d(1:100,r,-5)

// histogrammes
xbasc()
histplot(10,r)
xbasc()
histplot(0:0.2:1,r)

Simulations de lancers de pièces et de dés 

// lancer de pièce
int(2*rand)

// 50 lancers de pièce
int(2*rand(1,50))

// somme des faces de 2 dés
n=100;
s=int(1+6*rand(1,n))+int(1+6*rand(1,n));
xbasc();histplot(1.5:1:12.5,s,1,"011"," ",[1.5,0,12.5,1])

Exercice

Écrire dans un fichier une fonction de $n$ entier qui effectue $n$ simulations comme suit.

- Tirer un nombre $x$ au hasard dans $[0,1/37]$.

- Calculer la partie entière $y$ du (grand) scalaire $1/x$.

- Tirer un entier $z$ au hasard dans $\{1,...,y\}$.

- Calculer $r = z$ (modulo $37$).

Tracer un histogramme du résultat.

Quel jeu de hasard a-t-on ainsi approximativement simulé ?

Convergences 

 
// loi forte des grands nombres (avec loi uniforme)
rand('uniform');
N=1000;
x=rand(1:N);
// x est un tirage de N=1000 v.a. uniformes sur [0,1]
for i=2:N,x(i)=(i-1)*x(i-1)/i+x(i)/i;end;
// calcul récursif des moyennes partielles
xbasc();
plot2d(1:N,x,1,"011"," ",[0,0,N,1],[0,N/50,0,10]);

// une programmation plus rapide qui évite une boucle for
N=1000;
x=rand(1:N);
s=cumsum(x); // somme cumulée
x=s./[1:N]; // division terme à terme
xbasc();
plot2d(1:N,x,1,"011"," ",[0,0,N,1],[0,N/50,0,10]);

Autres travaux pratiques

Voir le chapitre [*] Simulations de loi du TD Calcul des probabilités.

14h45 $\longrightarrow$ 16h00. Zéros de fonctions, équations différentielles, optimisation, graphes

Zéros avec fsolve 

// fsolve
help fsolve

Zéro de fonction scalaire 

 
// zéro de polynôme
  deff('[y]=fct(x)','y=2*x^3-30*x^2-3*x+200');
x=[-3:0.1:15];xbasc();plot(x,fct(x));
  x1=fsolve(-1,fct)
fct(x1)
  x2=fsolve(1,fct)
fct(x2)
  x3=fsolve(11,fct)
fct(x3)

// zéro de polynôme, avec gradient
  deff('[y]=fct(x)','y=2*x^3-30*x^2-3*x+200');
  deff('[y]=grad_fct(x)','y=6*x^2-60*x-3');
x=[-3:0.1:15];xbasc();plot(x,fct(x));
  x1=fsolve(-1,fct,grad_fct)
fct(x1)
  x2=fsolve(1,fct)
fct(x2)
  x3=fsolve(11,fct)
fct(x3)

Intersection de quadriques 

// intersection de quadriques
deff('[z]=quadrique1(x,y)','z=2*x^2+ 5*y^2-30*x+20');
deff('[z]=quadrique2(x,y)','z=2*x^2 -y^2-3*y-20');

// dessin 
x=-2:10;
y=-10:10;
// on trace le contour 0 pour avoir la quadrique 
//  quadrique1(x,y)=0 
xbasc();
fcontour2d(x,y,quadrique1,[0,0],[9,9])
// On superpose le deuxieme 
fcontour2d(x,y,quadrique2,[0,0],[12,12],"000")

// recherche d'un point d'intersection 
deff('[Y]=quadriques(X)','Y=[quadrique1(X(1),X(2)),quadrique2(X(1),X(2))]');
rep=fsolve([-1,1],quadriques)

// on rajoute le point sur le dessin 
xpolys(rep(1),rep(2),-1)
 
// on vérifie le calcul 
quadriques(rep)

Autres travaux pratiques

Voir le chapitre [*] La recherche de zéro avec fsolve du TD Optimisation.

Intégration d'équations différentielles avec ode


\begin{displaymath} \frac{dy}{dt}=f(t,y) \, , \quad y(t_0)=y_0 \in {\mathbb{R}}^n \end{displaymath} 

// ode
help ode

Équations différentielles scalaires autonomes


\begin{displaymath} \frac{dy}{dt}=f(y) \, , \quad y(t_0)=y_0 \in {\mathbb{R}} \end{displaymath} 

 
// y dot = sin(y)
deff("[ydot]=f(t,y)","ydot=sin(y)")
//attention ! on écrit f(t,y) même si f ne dépend pas de t
y0=0.2;t0=0;t=0:0.1:15;
y=ode(y0,t0,t,f);
xbasc();  plot(t,y)

// y dot = -y^2
deff("[ydot]=f(t,y)","ydot=-y^2")
y0=0.2;t0=0;t=0:0.1:30;
y=ode(y0,t0,t,f);
xbasc();  plot(t,y)

// y dot = y^2
deff("[ydot]=f(t,y)","ydot=y^2")
y0=0.2;t0=0;t=0:0.1:30;
y=ode(y0,t0,t,f);
xbasc();  plot(t,y)

Équation différentielle scalaire non autonome


\begin{displaymath} \frac{dy}{dt}=f(t,y) \, , \quad y(t_0)=y_0 \in {\mathbb{R}} \end{displaymath} 

 // y dot = sin(t*y)
deff("[ydot]=f(t,y)","ydot=sin(t*y)")
y0=0.2;t0=0;t=0:0.1:15;
y=ode(y0,t0,t,f);
xbasc();  plot(t,y)

Système différentiel


\begin{displaymath} \left\{ \begin{array}{rcl} \displaystyle \frac{dy_1}{dt} &=&... ...(y_1(t_0),y_2(t_0))=y_0 \in {\mathbb{R}}^2 \end{array} \right. \end{displaymath} 

// y dot dot = -sin(y)
deff('[z]=fct1(y1,y2)','z=y2');
deff('[z]=fct2(y1,y2)','z=-sin(y1)');
deff('[Z]=fct(Y)','Z=[fct1(Y(1),Y(2)),fct2(Y(1),Y(2))]');
y0=[0.3,0.2]';t0=0;t=0:0.1:30;
y=ode(y0,t0,t,f);
xbasc();  plot(t,y(1,:))

Autres travaux pratiques

Voir le chapitre [*] Systèmes dynamiques linéaires et stabilité du TD Systèmes dynamiques.

Voir le chapitre [*] Modèles biologiques et chimiques du TD Systèmes dynamiques.

Voir le chapitre [*] Attracteurs de systèmes dynamiques du TD Systèmes dynamiques.

Optimisation avec optim 

 //optim
help optim

Minimum d'une fonction de deux variables, sans contraintes


\begin{displaymath} \min_{x,y} J(x,y) \end{displaymath} 

 
// un exemple simple : optimiser x^2 + y^2 
deff('[f,g,ind]=cost(x,ind)',['f=x(1)^2+x(2)^2';'g=[2*x(1);2*x(2)]']);
// g est le gradient
// ici, ind est un paramètre non utilisé mais qui doit être présent
[f,xopt]=optim(cost,[1;2])
// le coût est quasi nul

Minimum d'une fonction de deux variables, avec contraintes


\begin{displaymath} \min_{x_{min} \leq x \leq x_{max}, y_{min} \leq y \leq y_{max}} J(x,y) \end{displaymath} 

//la même chose sous contrainte x dans [2,10] et y dans [-10,10] 
  //  le min est atteint sur un bord 

deff('[z]=C(x,y)','z=x^2+y^2');
x=2:10;y=-10:10;
z=feval(x,y,C);
xbasc();
plot3d(x,y,z);

[f,xopt,gopt]=optim(cost,'b',[2;-10],[10;10],[5;5]);
// f n'est pas nul 
// le gradient en xopt est perpendiculaire au bord

Autres travaux pratiques

Voir le sous-chapitre [*] Une minimisation de fonction différentiable du TD Optimisation.

Graphes avec Metanet 

 
Help

Metanet

showgraph

shortest_path

À propos de ce document...

Introduction à Scilab
Semaine logiciels, année 2000-2001
vendredi 20 octobre 2000

This document was generated using the LaTeX2HTML translator Version 99.2beta8 (1.43)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split 0 -nonavigation TPscilab_20oct2000

The translation was initiated by The Scilab Group on 2000-11-30

The Scilab Group 2000-11-30
| Français | English |