// // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab // Copyright (C) ????-2008 - INRIA // // This file is distributed under the same license as the Scilab package. // mode(7) // SCILAB OBJECTS // 1. SCALARS a=1 //real constant 1==1 //boolean 'string' //character string z=poly(0,'z') // polynomial with variable 'z' and with one root at zero p=1+3*z+4.5*z^2 //polynomial r=z/p //rational // 2. MATRICES A=[a+1 2 3 0 0 atan(1) 5 9 -1] //3 x 3 constant matrix b=[%t,%f] //1 x 2 boolean matrix Mc=['this','is'; 'a' ,'matrix'] //2 x 2 matrix of strings Mp=[p,1-z; 1,z*p] //2 x 2 polynomial matrix F=Mp/poly([1+%i 1-%i 1],'z') //rational matrix Sp=sparse([1,2;4,5;3,10],[1,2,3]) //sparse matrix Sp(1,10)==Sp(1,1) //boolean sparse matrix // 3. LISTS L=list(a,-(1:5), Mp,['this','is';'a','list']) //list L(2)(3) //sub-entry in list Lt=tlist(['mylist','color','position','weight'],'blue',[0,1],10) //typed-list Lt('color') //extracting Lt('weight') //extracting A=diag([2,3,4]);B=[1 0;0 1;0 0];C=[1 -1 0];D=0*C*B;x0=[0;0;0]; Sl=syslin('c',A,B,C,D,x0) //Standard state-space linear system Sl("A"), Sl("C") //Retrieving elements of a typed list Slt=ss2tf(Sl) // Transfer matrix Slt('num'), Slt('den') // OPERATIONS v=1:5;W=v'*v //constant matrix multiplication W(1,:) //extracting first row W(:,\$) //extracting last column Mp'*Mp+eye() //polynomial matrix Mp1=Mp(1,1)+4.5*%i //complex Fi=C*(z*eye()-A)^(-1)*B; //transfer function evaluation F(:,1)*Fi //operations with rationals M=[Mp -Mp; Mp' Mp+eye()] //concatenation of polynomial matrices [Fi, Fi(:,1)] // ... or rationals F=syslin('c',F); Num=F('num');Den=F('den'); //operation on transfer matrix // SOME NUMERICAL PRIMITIVES inv(A) //Inverse inv(Mp) //Inverse inv(Sl*Sl') //Product of two linear systems and inverse w=ss2tf(ans) //Transfer function representation w1=inv(ss2tf(Sl)*ss2tf(Sl')) //Product of two transfer functions and inverse clean(w-w1) A=rand(3,3);;B=rand(3,1);n=contr(A,B) //Controllability K=ppol(A,B,[-1-%i -1+%i -1]) //Pole placement poly(A-B*K,'z')-poly([-1-%i -1+%i -1],'z') //Check... s=sin(0:0.1:5*%pi); ss=fft(s(1:128),-1); //FFT my_handle = scf(100001);clf(my_handle,"reset"); plot2d3("enn",1,abs(ss)'); //simple plot // ON LINE DEFINITION OF FUNCTION deff('[x]=fact(n)','if n==0 then x=1,else x=n*fact(n-1),end') 10+fact(5) // OPTIMIZATION deff('[f,g,ind]=rosenbro(x,ind)', 'a=x(2)-x(1)^2 , b=1-x(2) ,... f=100.*a^2 + b^2 , g(1)=-400.*x(1)*a , g(2)=200.*a -2.*b '); [f,x,g]=optim(rosenbro,[2;2],'qn') // SIMULATION a=rand(3,3) e=expm(a) deff('[ydot]=f(t,y)','ydot=a*y'); e(:,1)-ode([1;0;0],0,1,f) // SYSTEM DEFINITION s=poly(0,'s'); h=[1/s,1/(s+1);1/s/(s+1),1/(s+2)/(s+2)] w=tf2ss(h); ss2tf(w) h1=clean(ans) // EXAMPLE: SECOND ORDER SYSTEM ANALYSIS sl=syslin('c',1/(s*s+0.2*s+1)) instants=0:0.05:20; // step response: y=csim('step',instants,sl); my_handle = scf(100001);clf(my_handle,"reset"); plot2d(instants',y') // Delayed step response deff('[in]=u(t)','if t<3 then in=0;else in=1;end'); y1=csim(u,instants,sl);plot2d(instants',y1'); clear u; // Impulse response; yi=csim('imp',instants,sl);clf();plot2d(instants',yi'); yi1=csim('step',instants,s*sl);plot2d(instants',yi1'); // Discretization dt=0.05; sld=dscr(tf2ss(sl),0.05); // Step response u=ones(instants); yyy=flts(u,sld); my_handle = scf(100001);clf(my_handle,"reset"); plot(instants,yyy); clear u; // Impulse response u=0*ones(instants);u(1)=1/dt; yy=flts(u,sld); my_handle = scf(100001);clf(my_handle,"reset"); plot(instants,yy); clear u; // system interconnexion w1=[w,w]; clean(ss2tf(w1)) w2=[w;w]; clean(ss2tf(w2)) // change of variable z=poly(0,'z'); horner(h,(1-z)/(1+z)) //bilinear transform // PRIMITIVES H=[1. 1. 1. 0.; 2. -1. 0. 1; 1. 0. 1. 1.; 0. 1. 2. -1]; ww=spec(H) // STABLE SUBSPACES [X,d]=schur(H,'cont'); X'*H*X [X,d]=schur(H,'disc'); X'*H*X //Selection of user-defined eigenvalues (# 3 and 4 here); function [flg]=sel(x) flg=%f if abs(x-ww(3))<0.0001|abs(x-ww(4))<0.00001 then flg=%t,end endfunction [X,d]=schur(H,sel) X'*H*X clear sel // With matrix pencil function [flg]=sel(x,y) flg=%f if abs(x/y-ww(3))<0.0001|abs(x/y-ww(4))<0.00001 then flg=%t,end endfunction [X,d]=schur(H,eye(H),sel) X'*H*X // block diagonalization [ab,x,bs]=bdiag(H); inv(x)*H*x // Matrix pencils E=rand(3,2)*rand(2,3); A=rand(3,2)*rand(2,3); s=poly(0,'s'); w=det(s*E-A) //determinant [al,be]=spec(A,E); al./(be+%eps*ones(be)) roots(w) [Ns,d]=coffg(s*E-A); //inverse of polynomial matrix; clean(Ns/d*(s*E-A)) [Q,M,i1]=pencan(E,A); // Canonical form; clean(M*E*Q) clean(M*A*Q) // PAUSE-RESUME write(%io(2),'pause command...'); write(%io(2),'TO CONTINUE...'); write(%io(2),'ENTER ''resume (or return) or click on resume!!'''); //pause; if haveacompiler() then // CALLING EXTERNAL ROUTINE foo=['void foo(double *a,double *b,double *c)'; '{ *c = *a + *b; }' ]; // we use TMPDIR for compilation if ~c_link('foo') then path = pwd(); chdir(TMPDIR); mputl(foo,'foo.c'); ilib_for_link(['foo'],'foo.c',[],"c"); exec loader.sce chdir(path) end //5+7 by C function call('foo',5,1,'d',7,2,'d','out',[1,1],3,'d') end