- Timestamp:
- May 23, 2005, 5:36:46 PM (19 years ago)
- Branches:
- (u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd08f5f0bb3329b8ca19f23b74cb1473686415c3a')
- Children:
- 22674766072105828e687803f9aa5476d378511f
- Parents:
- d2ef5ca4ff890d10ffea61a0fda0ae01e2d7ef7f
- File:
-
- 1 edited
-
Singular/LIB/control.lib (modified) (33 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/control.lib
rd2ef5c ra29de75 1 version="$Id: control.lib,v 1.3 1 2005-05-06 14:38:12 hannesExp $";1 version="$Id: control.lib,v 1.32 2005-05-23 15:36:46 levandov Exp $"; 2 2 category="System and Control Theory"; 3 3 info=" … … 13 13 MAIN PROCEDURES: 14 14 control(R); analysis of controllability-related properties of R (using Ext modules) 15 control 2(R); analysis of controllability-related properties of R (using dimension)15 controlDim(R); analysis of controllability-related properties of R (using dimension) 16 16 autonom(R); analysis of autonomy-related properties of R (using Ext modules) 17 autonom 2(R); analysis of autonomy-related properties of R (using dimension)17 autonomDim(R); analysis of autonomy-related properties of R (using dimension) 18 18 19 19 COMPONENT PROCEDURES: 20 LeftKernel(R); a left kernel of R21 RightKernel(R); a right kernel of R22 LeftInverse(R); a left inverse of R23 RightInverse(R); a right inverse of R20 leftKernel(R); a left kernel of R 21 rightKernel(R); a right kernel of R 22 leftInverse(R); a left inverse of R 23 rightInverse(R); a right inverse of R 24 24 smith(M); a Smith form of a module M 25 25 colrank(M); a column rank of M as of matrix … … 27 27 canonize(L); Groebnerification for modules in the output of control or autonomy procs 28 28 iostruct(R); computes an IO-structure of behavior given by a module R 29 FindTorsion(R, I); generators of the submodule of a module R, annihilated by the ideal I29 findTorsion(R, I); generators of the submodule of a module R, annihilated by the ideal I 30 30 31 31 AUXILIARY PROCEDURES: 32 ControlExample(s); set up an example from the mini database inside of the library 33 declare(N,V,P,O); defines the ring easily 32 controlExample(s); set up an example from the mini database inside of the library 34 33 view(); well-formatted output of lists, modules and matrices 35 34 "; … … 142 141 {"EXAMPLE:";echo = 2; 143 142 ring r; 144 matrix M[1][3] = x,y,z; 145 print(M); 146 view(M); 143 list L; 144 matrix M[1][3] = x2+x,y3-y,z5-4z+7; 145 L[1] = "a matrix:"; 146 L[2] = M; 147 L[3] = "an ideal:"; 148 L[4] = ideal(M); 149 view(L); 147 150 }; 148 151 //-------------------------------------------------------------------------- 149 proc RightKernel(matrix M)150 "USAGE: RightKernel(M); M a matrix152 proc rightKernel(matrix M) 153 "USAGE: rightKernel(M); M a matrix 151 154 PURPOSE: computes the right kernel of matrix M (a module of all elements v such that Mv=0) 152 155 RETURN: module 153 EXAMPLE: example RightKernel; shows an example156 EXAMPLE: example rightKernel; shows an example 154 157 "{ 155 158 return(modulo(M,std(0))); … … 160 163 matrix M[1][3] = x,y,z; 161 164 print(M); 162 matrix R = RightKernel(M);165 matrix R = rightKernel(M); 163 166 print(R); 164 167 // check: … … 166 169 }; 167 170 //------------------------------------------------------------------------- 168 proc LeftKernel(matrix M)169 "USAGE: LeftKernel(M); M a matrix171 proc leftKernel(matrix M) 172 "USAGE: leftKernel(M); M a matrix 170 173 PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0) 171 174 RETURN: module 172 EXAMPLE: example LeftKernel; shows an example175 EXAMPLE: example leftKernel; shows an example 173 176 " 174 177 { … … 180 183 matrix M[3][1] = x,y,z; 181 184 print(M); 182 matrix L = LeftKernel(M);185 matrix L = leftKernel(M); 183 186 print(L); 184 187 // check: … … 186 189 }; 187 190 //------------------------------------------------------------------------ 188 proc LeftInverse(module M)189 "USAGE: LeftInverse(M); M a module190 PURPOSE: computes such a matrix L, that LM = =Id;191 proc leftInverse(module M) 192 "USAGE: leftInverse(M); M a module 193 PURPOSE: computes such a matrix L, that LM = Id; 191 194 RETURN: module 192 EXAMPLE: example LeftInverse; shows an example193 NOTE: exists only in the case when Id belongs to M!195 EXAMPLE: example leftInverse; shows an example 196 NOTE: exists only in the case when M is free submodule 194 197 " 195 198 { … … 221 224 matrix M[2][1] = 1,x2z; 222 225 print(M); 223 print( LeftInverse(M) );226 print( leftInverse(M) ); 224 227 kill r; 225 228 // derived from the example TwoPendula: … … 230 233 U[3,1]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2); 231 234 module M = module(U); 232 module L = LeftInverse(M);235 module L = leftInverse(M); 233 236 print(L); 234 237 // check … … 236 239 }; 237 240 //----------------------------------------------------------------------- 238 proc RightInverse(module R)239 "USAGE: RightInverse(M); M a module240 PURPOSE: computes such a matrix L, that ML = = Id;241 proc rightInverse(module R) 242 "USAGE: rightInverse(M); M a module 243 PURPOSE: computes such a matrix L, that ML = Id 241 244 RETURN: module 242 EXAMPLE: example RightInverse; shows an example243 NOTE: exists only in the case when Id belongs to M!244 " 245 { 246 return(transpose( LeftInverse(transpose(R))));245 EXAMPLE: example rightInverse; shows an example 246 NOTE: exists only in the case when M is free submodule 247 " 248 { 249 return(transpose(leftInverse(transpose(R)))); 247 250 } 248 251 example … … 252 255 matrix M[1][2] = 1,x2+z; 253 256 print(M); 254 print( RightInverse(M) );257 print( rightInverse(M) ); 255 258 kill r; 256 259 // derived from the TwoPendula example: … … 261 264 U[1,3]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2); 262 265 module M = module(U); 263 module L = RightInverse(M);266 module L = rightInverse(M); 264 267 print(L); 265 268 // check … … 323 326 string Gen_mes = "Parameter constellations which might lead to a non-controllable system:"; 324 327 325 module RK = RightKernel(R);328 module RK = rightKernel(R); 326 329 int d=dim_Our(std(transpose(R))); 327 330 … … 334 337 RK, 335 338 "kernel representation for controllable part:", 336 LeftKernel( RK ),339 leftKernel( RK ), 337 340 "obstruction to controllability", 338 341 Ext_1, … … 352 355 RK, 353 356 "left inverse to image representation:", 354 LeftInverse(RK),357 leftInverse(RK), 355 358 DofS, 356 359 d, … … 433 436 }; 434 437 //-------------------------------------------------------------------------- 435 proc control 2(module R)436 "USAGE: control 2(R); R a module (R is the matrix of the system of equations to be investigated)438 proc controlDim(module R) 439 "USAGE: controlDim(R); R a module (R is the matrix of the system of equations to be investigated) 437 440 PURPOSE: computes list of all the properties concerning controllability of the system (behavior), represented by the matrix R 438 441 RETURN: list 439 EXAMPLE: example control 2; shows an example442 EXAMPLE: example controlDim; shows an example 440 443 NOTE: this procedure is analogous to 'control' but uses dimension calculations.This approach works for full row rank matrices only. 441 444 " … … 443 446 if( nrows(R) != colrank(transpose(R)) ) 444 447 { 445 return ("control 2cannot be applied, since R does not have full row rank");448 return ("controlDim cannot be applied, since R does not have full row rank"); 446 449 } 447 450 intvec v=Opt_Our(); … … 466 469 R=transpose(R); 467 470 view(R); 468 view(control 2(R));471 view(controlDim(R)); 469 472 }; 470 473 //------------------------------------------------------------------------ 471 474 proc colrank(module M) 472 "USAGE: proc colrank(M),M a matrix/module475 "USAGE: colrank(M); M a matrix/module 473 476 PURPOSE: compute the column rank of M as of matrix 474 477 RETURN: int 475 NOTE: this procedure uses bareiss-algorithm which might not terminate in some cases 476 " 477 { 478 NOTE: this procedure uses Bareiss algorithm 479 480 " 481 { 482 // NOte continued: 483 // which might not terminate in some cases 478 484 module M_red = bareiss(M)[1]; 479 485 int NCols_red = ncols(M_red); … … 571 577 }; 572 578 //-------------------------------------------------------------------------- 573 proc autonom 2(module R)574 "USAGE: autonom 2(R); R a module (R is a matrix of the system of equations which is to be investigated)579 proc autonomDim(module R) 580 "USAGE: autonomDim(R); R a module (R is a matrix of the system of equations which is to be investigated) 575 581 PURPOSE: computes the list of all the properties concerning autonomy of the system (behavior), represented by the matrix R 576 582 RETURN: list 577 583 NOTE: this procedure is analogous to 'autonom' but uses dimension calculations 578 EXAMPLE: example autonom 2; shows an example584 EXAMPLE: example autonomDim; shows an example 579 585 " 580 586 { … … 588 594 if( d==0 ) 589 595 { 590 RC= LeftKernel(RightKernel(R));596 RC=leftKernel(rightKernel(R)); 591 597 R_rank=colrank(R); 592 598 } … … 603 609 R=transpose(R); 604 610 view( R ); 605 view( autonom 2(R) );611 view( autonomDim(R) ); 606 612 }; 607 613 //---------------------------------------------------------- … … 627 633 if (i==0) 628 634 { 629 RC= LeftKernel(RightKernel(R));635 RC=leftKernel(rightKernel(R)); 630 636 R_rank=colrank(R); 631 637 } … … 648 654 //---------------------------------------------------------- 649 655 proc genericity(matrix M) 650 "USAGE: genericity(M) ,M is a matrix/module656 "USAGE: genericity(M); M is a matrix/module 651 657 PURPOSE: determine parametric expressions which have been assumed to be non-zero in the process of computing the Groebner basis 652 658 RETURN: list (of strings) … … 869 875 //--------------------------------------------------------------- 870 876 proc canonize(list L) 871 "USAGE: canonize(L) ,L a list877 "USAGE: canonize(L); L a list 872 878 PURPOSE: modules in the list are canonized by computing their reduced minimal (= unique up to constant factor w.r.t. the given ordering) Groebner bases 873 879 RETURN: list … … 926 932 RETURN: list L with entries: string s, intvec v, module P and module Q 927 933 PURPOSE: if R is the kernel-representation-matrix of some system, then we output a input-ouput representation Py=Qu of the system, the components that have been chosen as outputs(intvec v) and a comment s 928 NOTE: the procedure uses Bareiss algorithm which might not terminate in some cases934 NOTE: the procedure uses Bareiss algorithm 929 935 EXAMPLE: example iostruct; shows an example 930 936 " 931 937 { 938 // NOTE cont'd 939 //which might not terminate in some cases 932 940 list L = bareiss(R); 933 941 int R_rank = ncols(L[1]); … … 1127 1135 //--------------------------------------------------------------- 1128 1136 proc smith( module M ) 1129 "USAGE: smith(M) , M a module or a matrix,1130 PURPOSE: computes the Smith form of a matrix1137 "USAGE: smith(M); M a module/matrix 1138 PURPOSE: computes the Smith normal form of a matrix 1131 1139 RETURN: a list of length 4 with the following entries: 1132 @* [1]: The Smith-Form S of M,1140 @* [1]: the Smith normal form S of M, 1133 1141 @* [2]: the rank of M, 1134 1142 @* [3]: a unimodular matrix U, 1135 1143 @* [4]: a unimodular matrix V, 1136 such that U*M*V=S. An warning is returned when no Smith Form exists.1144 such that U*M*V=S. An warning is returned when no Smith form exists. 1137 1145 NOTE: The Smith form only exists over PIDs (principal ideal domains). Use global ordering for computations! 1138 1146 " … … 1256 1264 option(redSB); 1257 1265 option(redTail); 1258 ring r =0,x,dp;1259 // see what happens when the matrix is already in Smith-Form1260 module M = [x,0,0],[0,x2,0],[0,0,x3];1261 list L= smith(M);1262 print( L[1]);1263 matrix N =matrix(M);1264 matrix B =L[3]*N*L[4];1266 ring r = 0,x,dp; 1267 module M = [x2,x,3x3-4], [2x2-1,4x,5x2], [2x5,3x,4x]; 1268 print(M); 1269 list P = smith(M); 1270 print(P[1]); 1271 matrix N = matrix(M); 1272 matrix B = P[3]*N*P[4]; 1265 1273 print(B); 1266 //------- and yet another example -------------- 1267 module M2=[x2,x,3x3-4],[2x2-1,4x,5x2],[2x5,3x,4x]; 1268 print(M2); 1269 list P=smith(M2); 1270 print(P[1]); 1271 matrix N2=matrix(M2); 1272 matrix B2=P[3]*N2*P[4]; 1273 print(B2); 1274 } 1274 } 1275 // see what happens when the matrix is already in Smith-Form 1276 // module M = [x,0,0],[0,x2,0],[0,0,x3]; 1277 // list L = smith(M); 1278 // print(L[1]); 1279 //matrix N=matrix(M); 1280 //matrix B=L[3]*N*L[4]; 1281 //print(B); 1275 1282 //--------------------------------------------------------------- 1276 1283 static proc list_tex(L, string name,link l,int nr_loop) … … 1355 1362 } 1356 1363 //--------------------------------------------------------------- 1357 proc FindTorsion(module R, ideal TAnn)1358 "USAGE: FindTorsion(R, I); R an ideal/matrix/module, I an ideal1364 proc findTorsion(module R, ideal TAnn) 1365 "USAGE: findTorsion(R, I); R an ideal/matrix/module, I an ideal 1359 1366 PURPOSE: computes the Groebner basis of the submodule of R, annihilated by I 1360 1367 ETURN: module 1361 1368 NOTE: especially helpful, when I is the annihilator of the t(R) - the torsion submodule of R. In this case, the result is the explicit presentation of t(R) as 1362 1369 the submodule of R 1363 EXAMPLE: example FindTorsion; shows an example1370 EXAMPLE: example findTorsion; shows an example 1364 1371 " 1365 1372 { … … 1392 1399 // here, we have the annihilator: 1393 1400 ideal LAnn = D1; // = L[10] 1394 module Tr = FindTorsion(RR,LAnn);1401 module Tr = findTorsion(RR,LAnn); 1395 1402 print(RR); // the module itself 1396 1403 print(Tr); // generators of the torsion submodule … … 1398 1405 1399 1406 1400 proc ControlExample(string s)1401 "USAGE: ControlExample(s); s a string1407 proc controlExample(string s) 1408 "USAGE: controlExample(s); s a string 1402 1409 PURPOSE: set up an example from the mini database by initalizing a ring and a module in a ring 1403 1410 RETURN: ring 1404 NOTE: in order to see the list of available examples, execute @code{ ControlExample(\"show\");}1411 NOTE: in order to see the list of available examples, execute @code{controlExample(\"show\");} 1405 1412 @* To use ab example, one has to do the following. Suppose one calls the ring, where the example will be activated, A. Then, by executing 1406 @* @code{def A = ControlExample(\"Antenna\");} and @code{setring A;},1413 @* @code{def A = controlExample(\"Antenna\");} and @code{setring A;}, 1407 1414 @* A will become a basering from the example \"Antenna\" with 1408 1415 the predefined system module R (transposed). 1409 1416 After that one can just execute @code{control(R);} respectively 1410 1417 @code{autonom(R);} to perform the control resp. autonomy analysis of R. 1411 EXAMPLE: example ControlExample; shows an example1418 EXAMPLE: example controlExample; shows an example 1412 1419 "{ 1413 1420 list E, S, D; // E=official name, S=synonym, D=description … … 1449 1456 { 1450 1457 "EXAMPLE:";echo = 2; 1451 ControlExample("show"); // let us see all available examples:1452 def B = ControlExample("TwoPendula"); // let us set up a particular example1458 controlExample("show"); // let us see all available examples: 1459 def B = controlExample("TwoPendula"); // let us set up a particular example 1453 1460 setring B; 1454 1461 print(R); … … 1591 1598 return(@r); 1592 1599 }; 1593 1594 1595 //---------------------------------------------------------------1596 proc declare(string NameOfRing, string Variables, list #)1597 "USAGE: declare(NameOfRing, Variables,[Parameters, Ordering]); where1598 @* NameOfRing is string with name of ring,1599 @* Variables is a string with names of variables separated by commas (e.g. \"x,y,z\"),1600 @* Parameters is string of parameters in the ring separated by commas (e.g. \"a,b,c\"),1601 @* Ordering is string with name of ordering (by default, the ordering (dp,C) is used).1602 PURPOSE: define the ring easily1603 RETURN: no return value1604 EXAMPLE: example declare; shows an example1605 "1606 {1607 if(size(#)==0)1608 {1609 execute("ring "+NameOfRing+"=0,("+Variables+"),dp;");1610 }1611 else1612 {1613 if(size(#)==1)1614 {1615 execute( "ring " + NameOfRing + "=(0," + #[1] + "),(" + Variables + "),dp;" );1616 }1617 else1618 {1619 if( (size(#[1])!=0)&&(#[1]!=" ") )1620 {1621 execute( "ring " + NameOfRing + "=(0," + #[1] + "),(" + Variables + "),("+#[2]+");" );1622 }1623 else1624 {1625 execute( "ring " + NameOfRing + "=0,("+Variables+"),("+#[2]+");" );1626 };1627 };1628 };1629 keepring(basering);1630 }1631 example1632 {"EXAMPLE:";echo = 2;1633 string v="x,y,z";1634 string p="q,p";1635 string Ord ="c,lp";1636 //----------------------------------1637 declare("Ring_1",v);1638 print(nameof(basering));1639 print(basering);1640 //----------------------------------1641 declare("Ring_2",v,p);1642 print(basering);1643 print(nameof(basering));1644 //----------------------------------1645 declare("Ring_3",v,p,Ord);1646 print(basering);1647 print(nameof(basering));1648 //----------------------------------1649 declare("Ring_4",v,"",Ord);1650 print(basering);1651 print(nameof(basering));1652 //----------------------------------1653 declare("Ring_5",v," ",Ord);1654 print(basering);1655 print(nameof(basering));1656 }1657 //1658 // maybe reasonable to add this in declare1659 //1660 // print("Please enter your representation matrix in the following form:1661 // module R=[1st row],[2nd row],...");1662 // print("Type the command: R=transpose(R)");1663 // print(" To compute controllability please enter: control(R)");1664 // print(" To compute autonomy please enter: autonom(R)");
Note: See TracChangeset
for help on using the changeset viewer.
