Next: 8. System interconnexion
Up: No Title
Previous: 6.6.2.0.5 Scilab functions
We give here the Scilab session corresponding to the first demo.
-->// SCILAB OBJECTS
-->// 1. SCALARS
-->a=1 //real constant
a =
1.
-->1==1 //boolean
ans =
T
-->'string' //character string
ans =
string
-->z=poly(0,'z') // polynomial with variable 'z' and with one root at zero
z =
z
-->p=1+3*z+4.5*z^2 //polynomial
p =
2
1 + 3z + 4.5z
-->r=z/p //rational
r =
z
-------------
2
1 + 3z + 4.5z
-->// 2. MATRICES
-->A=[a+1 2 3
--> 0 0 atan(1)
--> 5 9 -1] //3 x 3 constant matrix
A =
! 2. 2. 3. !
! 0. 0. 0.7853982 !
! 5. 9. - 1. !
-->b=[%t,%f] //1 x 2 boolean matrix
b =
! T F !
-->Mc=['this','is';
--> 'a' ,'matrix'] //2 x 2 matrix of strings
Mc =
!this is !
! !
!a matrix !
-->Mp=[p,1-z;
--> 1,z*p] //2 x 2 polynomial matrix
Mp =
! 2 !
! 1 + 3z + 4.5z 1 - z !
! !
! 2 3 !
! 1 z + 3z + 4.5z !
-->F=Mp/poly([1+%i 1-%i 1],'z') //rational matrix
F =
! 2 !
! 1 + 3z + 4.5z - 1 !
! --------------- --------- !
! 2 3 2 !
! - 2 + 4z - 3z + z 2 - 2z + z !
! !
! 2 3 !
! 1 z + 3z + 4.5z !
! --------------- --------------- !
! 2 3 2 3 !
! - 2 + 4z - 3z + z - 2 + 4z - 3z + z !
-->Sp=sparse([1,2;4,5;3,10],[1,2,3]) //sparse matrix
Sp =
( 4, 10) sparse matrix
( 1, 2) 1.
( 3, 10) 3.
( 4, 5) 2.
-->Sp(1,10)==Sp(1,1) //boolean sparse matrix
ans =
( 1, 1) sparse matrix
( 1, 1) T
-->// 3. LISTS
-->L=list(a,-(1:5), Mp,['this','is';'a','list']) //list
L =
L(1)
1.
L(2)
! - 1. - 2. - 3. - 4. - 5. !
L(3)
! 2 !
! 1 + 3z + 4.5z 1 - z !
! !
! 2 3 !
! 1 z + 3z + 4.5z !
L(4)
!this is !
! !
!a list !
-->L(2)(3) //sub-entry in list
ans =
- 3.
-->Lt=tlist(['mylist','color','position','weight'],'blue',[0,1],10) //typed-list
Lt =
Lt(1)
!mylist color position weight !
Lt(2)
blue
Lt(3)
! 0. 1. !
Lt(4)
10.
-->Lt('color') //extracting
ans =
blue
-->Lt('weight') //extracting
ans =
10.
-->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 =
Sl(1) (state-space system:)
!lss A B C D X0 dt !
Sl(2) = A matrix =
! 2. 0. 0. !
! 0. 3. 0. !
! 0. 0. 4. !
Sl(3) = B matrix =
! 1. 0. !
! 0. 1. !
! 0. 0. !
Sl(4) = C matrix =
! 1. - 1. 0. !
Sl(5) = D matrix =
! 0. 0. !
Sl(6) = X0 (initial state) =
! 0. !
! 0. !
! 0. !
Sl(7) = Time domain =
c
-->Sl("A"), Sl("C") //Retrieving elements of a typed list
ans =
! 2. 0. 0. !
! 0. 3. 0. !
! 0. 0. 4. !
ans =
! 1. - 1. 0. !
-->Slt=ss2tf(Sl) // Transfer matrix
Slt =
! 1 - 1 !
! ----- ----- !
! - 2 + s - 3 + s !
-->Slt('num'), Slt('den')
ans =
! 1 - 1 !
ans =
! - 2 + s - 3 + s !
-->// OPERATIONS
-->v=1:5;W=v'*v //constant matrix multiplication
W =
! 1. 2. 3. 4. 5. !
! 2. 4. 6. 8. 10. !
! 3. 6. 9. 12. 15. !
! 4. 8. 12. 16. 20. !
! 5. 10. 15. 20. 25. !
-->W(1,:) //extracting first row
ans =
! 1. 2. 3. 4. 5. !
-->W(:,$) //extracting last column
ans =
! 5. !
! 10. !
! 15. !
! 20. !
! 25. !
-->Mp'*Mp+eye() //polynomial matrix
ans =
column 1
! 2 3 4 !
! 3 + 6z + 18z + 27z + 20.25z !
! !
! 2 !
! 1 + 3z + 4.5z !
column 2
! 2 !
! 1 + 3z + 4.5z !
! !
! 2 3 4 5 6 !
! 2 - 2z + 2z + 6z + 18z + 27z + 20.25z !
-->Mp1=Mp(1,1)+4.5*%i //complex
Mp1 =
real part
2
1 + 3z + 4.5z
imaginary part
4.5
-->Fi=C*(z*eye()-A)^(-1)*B; //transfer function evaluation
-->F(:,1)*Fi //operations with rationals
ans =
! 2 2 !
! 1 + 3z + 4.5z - 1 - 3z - 4.5z !
! --------------------- --------------------- !
! 2 3 4 2 3 4 !
! 4 - 10z + 10z - 5z + z 6 - 14z + 13z - 6z + z !
! !
! 1 - 1 !
! --------------------- --------------------- !
! 2 3 4 2 3 4 !
! 4 - 10z + 10z - 5z + z 6 - 14z + 13z - 6z + z !
-->M=[Mp -Mp; Mp' Mp+eye()] //concatenation of polynomial matrices
M =
column 1 to 3
! 2 2 !
! 1 + 3z + 4.5z 1 - z - 1 - 3z - 4.5z !
! !
! 2 3 !
! 1 z + 3z + 4.5z - 1 !
! !
! 2 2 !
! 1 + 3z + 4.5z 1 2 + 3z + 4.5z !
! !
! 2 3 !
! 1 - z z + 3z + 4.5z 1 !
column 4
! - 1 + z !
! !
! 2 3 !
! - z - 3z - 4.5z !
! !
! 1 - z !
! !
! 2 3 !
! 1 + z + 3z + 4.5z !
-->[Fi, Fi(:,1)] // ... or rationals
ans =
! 1 - 1 1 !
! ----- ----- ----- !
! - 2 + z - 3 + z - 2 + z !
-->F=syslin('c',F);
-->Num=F('num');Den=F('den'); //operation on transfer matrix
-->// SOME NUMERICAL PRIMITIVES
-->inv(A) //Inverse
ans =
! 0.5 0. 0. !
! 0. 0.3333333 0. !
! 0. 0. 0.25 !
-->inv(Mp) //Inverse
ans =
column 1
! 2 3 !
! z + 3z + 4.5z !
! ------------------------------- !
! 2 3 4 5 !
! - 1 + 2z + 6z + 18z + 27z + 20.25z !
! !
! - 1 !
! ------------------------------- !
! 2 3 4 5 !
! - 1 + 2z + 6z + 18z + 27z + 20.25z !
column 2
! - 1 + z !
! ------------------------------- !
! 2 3 4 5 !
! - 1 + 2z + 6z + 18z + 27z + 20.25z !
! !
! 2 !
! 1 + 3z + 4.5z !
! ------------------------------- !
! 2 3 4 5 !
! - 1 + 2z + 6z + 18z + 27z + 20.25z !
-->inv(Sl*Sl') //Product of two linear systems and inverse
ans =
ans(1) (state-space system:)
!lss A B C D X0 dt !
ans(2) = A matrix =
! 2.8641369 - 0.9304438 0. 0. !
! 0.4111970 2.1358631 0. 0. !
! 0. - 9.339D-16 4. 0. !
! 0. 0. 0. 4. !
ans(3) = B matrix =
! 0.7027925 !
! 0.3620534 !
! 1.665D-16 !
! 0. !
ans(4) = C matrix =
! - 0.3238304 0.6285968 1.890D-15 0. !
ans(5) = D matrix =
2
2.75 - 2.5s + 0.5s
ans(6) = X0 (initial state) =
! 0. !
! 0. !
! 0. !
! 0. !
ans(7) = Time domain =
c
-->w=ss2tf(ans) //Transfer function representation
w =
2 3 4
18 - 30s + 18.5s - 5s + 0.5s
----------------------------
2
6.5 - 5s + s
-->w1=inv(ss2tf(Sl)*ss2tf(Sl')) //Product of two transfer functions and inverse
w1 =
2 3 4
36 - 60s + 37s - 10s + s
------------------------
2
13 - 10s + 2s
-->clean(w-w1)
ans =
1.730D-09 - 6.605D-10s
----------------------
2
6.5 - 5s + s
-->A=rand(3,3);;B=rand(3,1);n=contr(A,B) //Controllability
n =
3.
-->K=ppol(A,B,[-1-%i -1+%i -1]) //Pole placement
K =
! 7.1638394 7.2295307 0.3176982 !
-->poly(A-B*K,'z')-poly([-1-%i -1+%i -1],'z') //Check...
ans =
2
- 8.882D-16 + 1.776D-15z - 1.332D-15z
-->s=sin(0:0.1:5*%pi);
-->ss=fft(s(1:128),-1); //FFT
-->xbasc();
-->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')
Warning: obsolete use of = instead of ==
if n=0 then x=1,else x=n*fact(n-1),end
!
at line 2 of function fact called by :
deff('[x]=fact(n)','if n=0 then x=1,else x=n*fact(n-1),end')
-->10+fact(5)
ans =
130.
-->// 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')
g =
! 0. !
! 0. !
x =
! 1. !
! 1. !
f =
0.
-->// SIMULATION
-->a=rand(3,3)
a =
! 0.7263507 0.2320748 0.8833888 !
! 0.1985144 0.2312237 0.6525135 !
! 0.5442573 0.2164633 0.3076091 !
-->e=expm(a)
e =
! 2.6034702 0.5788017 1.7895052 !
! 0.6508072 1.4242213 1.1108378 !
! 1.0616116 0.4238856 1.8909166 !
-->deff('[ydot]=f(t,y)','ydot=a*y');
-->e(:,1)-ode([1;0;0],0,1,f)
ans =
! - 1.425D-07 !
! - 7.368D-08 !
! - 8.683D-08 !
-->// SYSTEM DEFINITION
-->s=poly(0,'s');
-->h=[1/s,1/(s+1);1/s/(s+1),1/(s+2)/(s+2)]
h =
! 1 1 !
! - ----- !
! s 1 + s !
! !
! 1 1 !
! ----- --------- !
! 2 2 !
! s + s 4 + 4s + s !
-->w=tf2ss(h);
-->ss2tf(w)
ans =
! 1 1 !
! ------------- ----- !
! - 4.710D-16 + s 1 + s !
! !
! 1 + 6.935D-16s 1 + 2.448D-16s !
! ---------------- -------------- !
! 2 2 !
! - 1.610D-15 + s + s 4 + 4s + s !
-->h1=clean(ans)
h1 =
! 1 1 !
! - ----- !
! s 1 + s !
! !
! 1 1 !
! ----- --------- !
! 2 2 !
! s + s 4 + 4s + s !
-->// EXAMPLE: SECOND ORDER SYSTEM ANALYSIS
-->sl=syslin('c',1/(s*s+0.2*s+1))
sl =
1
-----------
2
1 + 0.2s + s
-->instants=0:0.05:20;
-->// step response:
-->y=csim('step',instants,sl);
-->xbasc();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');
-->// Impulse response;
-->yi=csim('imp',instants,sl);xbasc();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);
Warning :redefining function: u
-->yyy=flts(u,sld);
-->xbasc();plot(instants,yyy)
-->// Impulse response
-->u=0*ones(instants);u(1)=1/dt;
-->yy=flts(u,sld);
-->xbasc();plot(instants,yy)
-->// system interconnexion
-->w1=[w,w];
-->clean(ss2tf(w1))
ans =
! 1 1 1 1 !
! - ----- - ----- !
! s 1 + s s 1 + s !
! !
! 1 1 1 1 !
! ----- --------- ----- --------- !
! 2 2 2 2 !
! s + s 4 + 4s + s s + s 4 + 4s + s !
-->w2=[w;w];
-->clean(ss2tf(w2))
ans =
! 1 1 !
! - ----- !
! s 1 + s !
! !
! 1 1 !
! ----- --------- !
! 2 2 !
! s + s 4 + 4s + s !
! !
! 1 1 !
! - ----- !
! s 1 + s !
! !
! 1 1 !
! ----- --------- !
! 2 2 !
! s + s 4 + 4s + s !
-->// change of variable
-->z=poly(0,'z');
-->horner(h,(1-z)/(1+z)) //bilinear transform
ans =
! 1 + z 1 + z !
! ----- ----- !
! 1 - z 2 !
! !
! 2 2 !
! 1 + 2z + z 1 + 2z + z !
! ---------- ---------- !
! 2 !
! 2 - 2z 9 + 6z + z !
-->// PRIMITIVES
-->H=[1. 1. 1. 0.;
--> 2. -1. 0. 1;
--> 1. 0. 1. 1.;
--> 0. 1. 2. -1];
-->ww=spec(H)
ww =
! 2.7320508 !
! - 2.7320508 !
! 0.7320508 !
! - 0.7320508 !
-->// STABLE SUBSPACES
-->[X,d]=schur(H,'cont');
-->X'*H*X
ans =
! - 2.7320508 - 1.110D-15 0. 1. !
! 0. - 0.7320508 - 1. - 7.772D-16 !
! 7.216D-16 0. 2.7320508 0. !
! 0. - 6.106D-16 0. 0.7320508 !
-->[X,d]=schur(H,'disc');
-->X'*H*X
ans =
! 0.7320508 0. 7.772D-16 1. !
! 0. - 0.7320508 - 1. 8.604D-16 !
! 0. 0. 2.7320508 - 1.166D-15 !
! 7.772D-16 1.110D-15 - 1.277D-15 - 2.7320508 !
-->//Selection of user-defined eigenvalues (# 3 and 4 here);
-->deff('[flg]=sel(x)',...
-->'flg=0,ev=x(2)/x(3),...
--> if abs(ev-ww(3))<0.0001|abs(ev-ww(4))<0.00001 then flg=1,end')
-->[X,d]=schur(H,sel)
d =
2.
X =
! - 0.5705632 - 0.2430494 - 0.6640233 - 0.4176813 !
! - 0.4176813 0.6640233 - 0.2430494 0.5705632 !
! 0.5705632 - 0.2430494 - 0.6640233 0.4176813 !
! 0.4176813 0.6640233 - 0.2430494 - 0.5705632 !
-->X'*H*X
ans =
! 0.7320508 0. 7.772D-16 1. !
! 0. - 0.7320508 - 1. 8.604D-16 !
! 0. 0. 2.7320508 - 1.166D-15 !
! 7.772D-16 1.110D-15 - 1.277D-15 - 2.7320508 !
-->// With matrix pencil
-->[X,d]=gschur(H,eye(H),sel)
d =
2.
X =
! 0.5705632 0.2430494 0.6640233 0.4176813 !
! 0.4176813 - 0.6640233 0.2430494 - 0.5705632 !
! - 0.5705632 0.2430494 0.6640233 - 0.4176813 !
! - 0.4176813 - 0.6640233 0.2430494 0.5705632 !
-->X'*H*X
ans =
! 0.7320508 0. 9.576D-16 1. !
! 0. - 0.7320508 - 1. 0. !
! 8.882D-16 0. 2.7320508 0. !
! 0. 0. 0. - 2.7320508 !
-->// block diagonalization
-->[ab,x,bs]=bdiag(H);
-->inv(x)*H*x
ans =
! 2.7320508 1.610D-15 0. 0. !
! - 3.664D-15 - 2.7320508 0. 6.661D-16 !
! 0. 0. 0.7320508 - 7.910D-16 !
! 0. 0. 0. - 0.7320508 !
-->// Matrix pencils
-->E=rand(3,2)*rand(2,3);
-->A=rand(3,2)*rand(2,3);
-->s=poly(0,'s');
-->w=det(s*D-A) //determinant
w =
2
- 0.0149837s + 0.0004193s
-->[al,be]=gspec(A,E);
-->al./(be+%eps*ones(be))
ans =
! 9.202D+14 !
! 35.734043 !
! - 1.170D-16 !
-->roots(w)
ans =
! 0 !
! 35.734043 !
-->[Ns,d]=coffg(s*D-A); //inverse of polynomial matrix;
-->clean(Ns/d*(s*D-A))
ans =
! 1 0 0 !
! - - - !
! 1 1 1 !
! !
! 0 1 0 !
! - - - !
! 1 1 1 !
! !
! 0 0 1 !
! - - - !
! 1 1 1 !
-->[Q,M,i1]=pencan(E,A); // Canonical form;
rank A^k rcond
2. 0.169D-01
rank A^k rcond
2. 0.847D+00
-->clean(M*E*Q)
ans =
! 1. 0. 0. !
! 0. 1. 0. !
! 0. 0. 0. !
-->clean(M*A*Q)
ans =
! 35.774235 - 3.0560234 0. !
! 0.4704929 - 0.0401920 0. !
! 0. 0. 1. !
-->// PAUSD-RESUME
-->write(%io(2),'pause command...');
pause command...
-->write(%io(2),'TO CONTINUE...');
TO CONTINUE...
-->write(%io(2),'ENTER ''resume (or return) or click on resume!!''');
ENTER 'resume (or return) or click on resume!!'
-->//pause;
-->// CALLING EXTERNAL ROUTINE
-->foo=['void foo(a,b,c)';
--> 'double *a,*b,*c;';
--> '{ *c = *a + *b; }' ];
-->unix_s('rm -f foo.c')
-->write('foo.c',foo);
-->unix_s('make foo.o') //Compiling...(needs a compiler)
-->link('foo.o','foo','C') //Linking to Scilab
ans =
0.
-->//5+7 by C function
-->fort('foo',5,1,'d',7,2,'d','out',[1,1],3,'d')
ans =
12.
Next: 8. System interconnexion
Up: No Title
Previous: 6.6.2.0.5 Scilab functions
Scilab Group