Introduction to Scilab

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

1 Scalars, vectors, matrices
2 Usual vectorial operations
3 Term by term vectorial operations
4 2D graphics
5 Programming
6 Scilab functions
7 Simulation of dynamic systems in discrete time

Contents

1 Scalars, vectors, matrices
2 Usual vectorial operations
3 Term by term vectorial operations
4 2D graphics
5 Programming
6 Scilab functions
7 Simulation of dynamic systems in discrete time
 7.1 Population dynamics in ecology
 7.2 Carbon cycle
 7.3 An ecosystem
Please launch software Scilab on your computer. To install Scilab, please consult http://www.scilab.org.

On this Internet page, you will find Scilab instructions lines starting with -->. You will copy them (or paste them by copy-paste) in the Scilab window to see how they are executed. Beware not to copy the prompt -->.

Every line starting with // is a line of comments and is not executed by Scilab.

1 Scalars, vectors, matrices

Question 1

2 Usual vectorial operations

3 Term by term vectorial operations

Question 2

4 2D graphics

For planar graphs, you have Scilab macros plot (elementary) and plot2d (with many options).

plot

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

PIC

plot2d

Execute the instruction help plot2d to have information about the many options of plot2d.

 
-->x=0.1:0.1:10;  
-->y=sin(x)./x;  
-->z=cos(x)./x;  
-->xbasc();  
// erase the graphic content of the window  
-->plot2d(x,y);plot2d(x,z)  
// both graphics are superposed  
-->xbasc();plot2d([x ; x]',[y ; z]')  
// the same, but both graphics are drawn each with its own color  
// [x ; x]' is a matrix with 2 colonnes, as is [y ; z]'  
// so that two curves are drawn  
-->xbasc();plot2d([x' x'],[y' z'])  
// gives the same result since [x ; x]' coincides with [x' x']  
// and [y' z'] avec [y ; z]'  
--> xset("window",2)  
// window 2 is now active and future graphic instructions will  
// appear in this window  
--> plot2d(x,y)  
--> xset("window",3) ; plot2d(x,z)  
--> xdel(0)  
// deletes the window 0  
--> xdel(2:3);  
// deletes the windows from 2 to 3  

PIC

-->r=rand(100,1)  
-->xbasc();plot2d(r)  
-->x=(1:100)';  
-->xbasc();plot2d(x,r,style=-5)  
-->s=rand(100,1);  
-->xbasc()  
-->plot2d(r,s,-2)  

PIC

-->y=[6:-1:0,1:5] ; x =[y,y,y,y,y];  
-->xbasc();plot2d(0:60,[x,6]);  
// a simple graphic  
-->xbasc();plot2d(0:60,[x,6],strf="011",rect=[0,0,60,6])  
// the frame is fixed by rect =[xmin,ymin,xmax,ymax]  
-->xbasc();plot2d(0:60,[x,6],strf="011",rect=[0,0,60,6],...  
nax=[0,5,0,6])  
// when an instruction is too long to hold in a single line  
// interrupt it by ... before changing line  
// nax=[0,5,0,6] fixes divisions and subdivisions on the axes  
// execute instruction help plot2d to see all the options  

PIC

plot2d2

-->x=1:10;  
-->xbasc();plot2d2(x,x);  

PIC

histogrammes

-->histplot(10,r)  
// draws the histogram of 100 reals randomly chosen  in [0,1]  
-->xbasc()  
-->histplot(0:0.2:1,r)  
// the same but in choosing the sizes of the classes  

PIC

Question 3 Plot the years in abcissa and the price of gasoline (euros of 2001) in ordinates. Take care that the ordinates graph ranges from 0 to 1 euro.

annees=[1973 1979 1983:2001]  
// only some years are available  
prix_gazole_EC2001=[0.557 0.715 0.879 0.885 0.902 0.699 0.652 0.618 ...  
0.633 0.643 0.627 0.592 0.614 0.638 0.624 0.681 0.697 0.658 0.702 0.848 0.783]  
// price of gasoline (euros of 2001)

PIC

5 Programming

Scilab is an interpreted language, where the transmission of the arguments is made by value, even for the objects of matrix type. It is thus recommended to avoid the loops, as they can prove ineffective in term of computation time. For that purpose, it is better, if possible, to make use of the vectorial operators who perform the same operations.

We recommend to use a text editor (Emacs), to open a file (for instance file_name.sce), to write the code lines and to execute them inside Scilab with the instruction exec("file_name.sce").

6 Scilab functions

One may define a Scilab function interactively. However, it is more practical and flexible to write the Scilab code in a file with a text editor. One can mix Scilabinstructions Scilab and functions definitions in a same file, named with the extension .sce.

The beginning of a Scilab function is
function [<arguments de retour>]=<nom>(<arguments d'entrée>)
and the code ends with the keyword endfunction.

   function [y]=densite(x), y=1 ./(%pi*sqrt(x.*(1-x))), endfunction;  
 
   x=0.01:0.01:0.99;y=densite(x);  
   xbasc(); plot(x,y);  
 
 
   function [H] = Heavyside(x)  
   // Heavyside function  
   H = bool2s(x>=0)  
   // Heavyside(0)=1  
   endfunction  
 
   x=-1:0.2:1; xbasc(); plot2d2(x,Heavyside(x),rect=[-1,-0.1,1,1.2]);  

Question 4 Write the function sinus cardinal (sin x∕x  ).

A Scilab function in Scilab may have several arguments, scalars or vectors.

 
   function [y]=Gauss(x,mu,sigma)  
   // mu and sigma are scalars  
   // x may be a vector of any dimensions  
   y=1 ./(sigma*sqrt(2*%pi))*exp(-((x-mu)/sigma).^2/2)  
   endfunction;  
 
   function [Y]=Normal(X,Mu,Sigma)  
   // Mu is a row vector with dimensions [1,dim]  
   // Sigma is a positive symmetric matrix with dimensions [dim,dim]  
   // x may only be a vector with dimensions [1,dim]  
   [lhs,dim]=size(Mu)  
   Y=1/ sqrt( (2*%pi)^dim * det(Sigma) ) * ...  
           exp(- (X-Mu)*inv(Sigma)*(X-Mu)' / 2 )  
   endfunction;

A Scilab function in Scilab may have several outputs.

 
   function [H,plus] = Heavyside_plus(x)  
   // Heavyside function  
   H = maxi(sign(x),0);  
   // Heavyside(0)=1  
   plus=maxi(x,0)  
   // Function +  
   endfunction  
 
   [H,plus] = Heavyside_plus(-1:0.2:1)

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

// loading of function permutation_lent :  
-->exec("nom_du_fichier.sce")  
// the function is now loaded  
-->permutation_lent  
// call of function permutation_lent  
-->permutation_lent(100)

7 Simulation of dynamic systems in discrete time

We consider the discrete time dynamics system in   n
R  :

x(t + 1) = f (t,x (t)), for t = t0,..,T − 1,
(1)

starting from initial state        n
x0 ∈ R  at time t0

x (t0) = x0.
(2)

Knowing function f  and (t0,x0)  , we are able to compute the whole sequence x (t0),x(t0 + 1),...,x(T)  solution of the problem.

7.1 Population dynamics in ecology

Consider the population dynamics

             Rx (t)
x (t + 1) =----------,t0 = 1,T = 100, x(t0) = x0 = 1,
           1 + Sx (t)

with parameters R = 1.2  and S = 0.02  .

 
t0=1;x0=1;  
x(t0)=x0; T=100; R=1.2; S=0.02;  
 
function y=f(t,x)  
y=R*x./(1+S*x)  
endfunction  
 
for (t=t0:1:T)  
x(t+1)=f(t,x(t));  
end;  
plot2d(t0:T+1,x(t0:T+1))

Question 5 Change the population dynamics with               0.5
f (t,x ) = (ax )   where a = 10  . Compare the behavior of the solutions x(t)  .

7.2 Carbon cycle

Consider the carbon cycle

M  (t + 1) = M (t) + αEBAU  (t)(1 − a) − δ(M (t) − M −∞ ),

starting from initial conditions M   = 354
  0  (ppm) at year t =  1990
 0  with time horizon T  = 100  . The parameters are α =  0.471  , preindustrial concentration M −∞  = 280  , removal rate δ = 0.01  . It is assumed that the ”business as usual” CO2   emissions path is

E    (t) = E    (t )(1 + g)t−t0 (GtC )
 BAU         BAU  0
with the emissions growth set to g = 1%  and initial emissions to EBAU  (1990 ) = 7.2  (GtC).
 
t0=1990; T=2100; M_0=354;  
 alpha=0.471;  M_infty=280; delta=1/120; sigma=0.519;  
E_BAU=7.2;g=0.01  
a=0;  
 
function y=EBAU(t)  
y=E_BAU*(1+g)^(t-t0);  
endfunction  
 
function y=f(t,M)  
y=M+alpha*EBAU(t)*(1-a)-delta*(M-M_infty);  
endfunction  
 
x(t0)=M_0;  
for (t=t0:1:T)  
x(t+1)=f(t,x(t));  
end;  
plot2d(t0:T+1,x(t0:T+1),rect=[t0,0,T+1,1000])

Question 6 Change the mitigation rate a  (a ∈ [0, 1]  ) to compare the behavior of the concentrations profile M (t)  .

7.3 An ecosystem

We consider two populations x1(t),x2(t)  interacting within a trophic web. The dynamics based on a Lotka-Volterra form is characterized by

{
  x1(t + 1)  =  x1 (t)(1 + r1 − q1x2 (t))
  x2(t + 1)  =  x2 (t)(1 − d2 + q2x1(t))
(3)

where r1 > 0  is the intrinsic growth of prey while d2 > 0  is the intrinsic decrease of predator. Parameters q1 > 0,q2 > 0  are related to the catchability and efficiency of trophic relations.

 
t0=1;x0=[1;1];  
 T=100;  
r_1=0.1; d_2=0.1; q_1=0.1; q_2=0.2;  
 
function y=f(t,x)  
x_1=x(1);x_2=x(2);  
y_1=x_1.*(1+r_1-q_1*x_2);  
y_2=x_2.*(1-d_2+q_2*x_1);  
y=[y_1;y_2];  
endfunction  
 
x=zeros(2,T+1);  
x(:,t0)=x0;  
for (t=t0:1:T)  
x(:,t+1)=f(t,x(:,t));  
end;  
plot2d(t0:T+1,x(:,t0:T+1)')

Question 7 Change parameters r,d  or q  to compare the behavior of the populations (x1(t),x2(t))  .