[a09c2a0] | 1 | version="$Id: control.lib,v 1.28 2005-04-29 14:53:50 levandov Exp $"; |
---|
| 2 | category="System and Control Theory"; |
---|
[4713c6] | 3 | info=" |
---|
| 4 | LIBRARY: control.lib Procedures for System and Control Theory |
---|
[f52c2f] | 5 | AUTHORS: Oleksandr Iena yena@mathematik.uni-kl.de |
---|
| 6 | @* Markus Becker mbecker@mathematik.uni-kl.de |
---|
| 7 | @* Viktor Levandovskyy levandov@mathematik.uni-kl.de |
---|
[4713c6] | 8 | |
---|
[5831fa9] | 9 | SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis' (Project of Dr. E. Zerz |
---|
| 10 | and V. Levandovskyy), Uni Kaiserslautern |
---|
[4713c6] | 11 | |
---|
[5831fa9] | 12 | NOTE: This library provides algebraic analysis tools for System and Control Theory |
---|
[4713c6] | 13 | |
---|
| 14 | PROCEDURES: |
---|
[6e9d48a] | 15 | control (R); analysis of controllability-related properties of R (using Ext modules), |
---|
| 16 | control2 (R); analysis of controllability-related properties of R (using dimension), |
---|
| 17 | autonom (R); analysis of autonomy-related properties of R (using Ext modules), |
---|
| 18 | autonom2 (R); analysis of autonomy-related properties of R (using dimension), |
---|
| 19 | LeftKernel (R); a left kernel of R, |
---|
| 20 | RightKernel (R); a right kernel of R, |
---|
| 21 | LeftInverse (R); a left inverse of R, |
---|
| 22 | RightInverse (R); a right inverse of R, |
---|
| 23 | smith (M); a Smith form of a module M, |
---|
| 24 | colrank (M); a column rank of M as of matrix, |
---|
| 25 | genericity (M); analysis of the genericity of parameters, |
---|
| 26 | canonize (L); Groebnerification for modules in the output of control/autonomy procs, |
---|
| 27 | iostruct (R); computes an I/O-structure of behavior given by a module R |
---|
| 28 | FindTorsion (R, I); generators of the submodule of a module R, annihilated by the ideal I. |
---|
[0c824aa] | 29 | |
---|
[4713c6] | 30 | |
---|
[5831fa9] | 31 | AUXILIARY PROCEDURES: |
---|
[d5f9d98] | 32 | declare(NameOfRing,Variables [,Parameters, Ordering]); define the ring easily |
---|
| 33 | view(); well-formatted output of lists, modules and matrices |
---|
[f691ff] | 34 | |
---|
| 35 | NOTE (EXAMPLES): In order to use examples below, execute the commands |
---|
| 36 | @* def A = exAntenna(); setring A; |
---|
[d5f9d98] | 37 | @* Thus A will become a basering from the example with the predefined module R (transposed), corresponding to the system. After that you can just type in |
---|
[a09c2a0] | 38 | @* control(R); // respectively autonom(R); |
---|
[f691ff] | 39 | and check the result. |
---|
| 40 | |
---|
| 41 | |
---|
| 42 | EXAMPLES (AUTONOMY): |
---|
| 43 | exCauchy(); example of 1-dimensional Cauchy equation, |
---|
| 44 | exCauchy2(); example of 2-dimensional Cauchy equation, |
---|
| 45 | exZerz(); example from the lecture of Eva Zerz, |
---|
| 46 | |
---|
| 47 | EXAMPLES (CONTROLLABILITY): |
---|
| 48 | ex1(); example of noncontrollable system, |
---|
| 49 | ex2(); example of controllable system , |
---|
| 50 | exAntenna(); Antenna, |
---|
| 51 | exEinstein(); Einstein equations, |
---|
| 52 | exFlexibleRod(); Flexible Rod, |
---|
[d5f9d98] | 53 | exTwoPendula(); Two Pendula mounted on a cart, |
---|
[f691ff] | 54 | exWindTunnel(); Wind Tunnel. |
---|
[4713c6] | 55 | "; |
---|
| 56 | |
---|
[5831fa9] | 57 | // NOTE: static things should not be shown for end-user |
---|
| 58 | // static Ext_Our(...) Copy of Ext_R from 'homolog.lib' in commutative case; |
---|
[0ceca2] | 59 | // static is_zero_Our(module R) Copy of is_zero from 'poly.lib'; |
---|
[5831fa9] | 60 | // static space(int n) Procedure used inside the procedure 'Print' to have a better formatted output |
---|
| 61 | // static control_output(); Generating the output for the procedure 'control' |
---|
| 62 | // static autonom_output(); Generating the output for the procedure 'autonom' and 'autonom2' |
---|
[6e9d48a] | 63 | // static extgcd_Our(poly p, poly q) Computes extgcd of p and q. for versions ealier than 2006 extgcd has a bug and is therefore not used |
---|
| 64 | // static normalize_Our(matrix N, matrix Q) normalizes the columns of N and divides the columns of Q through the leading coefficients of the columns of N |
---|
[5831fa9] | 65 | |
---|
[4713c6] | 66 | LIB "homolog.lib"; |
---|
| 67 | LIB "poly.lib"; |
---|
| 68 | LIB "primdec.lib"; |
---|
[501f7c5] | 69 | LIB "matrix.lib"; |
---|
[0e200a] | 70 | |
---|
[1deec7] | 71 | //--------------------------------------------------------------- |
---|
[0e200a] | 72 | static proc Opt_Our() |
---|
| 73 | "USAGE: Opt_Our(); |
---|
| 74 | RETURN: intvec, where previous options are stored |
---|
| 75 | PURPOSE: save previous options and set customized options |
---|
| 76 | " |
---|
| 77 | { |
---|
| 78 | intvec v; |
---|
| 79 | v=option(get); |
---|
| 80 | option(redSB); |
---|
| 81 | option(redTail); |
---|
| 82 | return (v); |
---|
| 83 | } |
---|
[0c824aa] | 84 | //--------------------------------------------------------------- |
---|
[4713c6] | 85 | proc declare(string NameOfRing, string Variables, list #) |
---|
[2e4e6ec] | 86 | "USAGE: declare(NameOfRing, Variables,[Parameters, Ordering]); |
---|
[d5f9d98] | 87 | @* NameOfRing: string with name of ring, |
---|
| 88 | @* Variables: string with names of variables separated by commas(e.g. \"x,y,z\"), |
---|
| 89 | @* Parameters: string of parameters in the ring separated by commas(e.g. \"a,b,c\"), |
---|
[a09c2a0] | 90 | @* Ordering: string with name of ordering (by default, the ordering (dp,C) is used) |
---|
[d5f9d98] | 91 | PURPOSE: define the ring easily |
---|
[2e4e6ec] | 92 | RETURN: no return value |
---|
| 93 | EXAMPLE: example declare; shows an example |
---|
| 94 | " |
---|
[4713c6] | 95 | { |
---|
| 96 | if(size(#)==0) |
---|
| 97 | { |
---|
| 98 | execute("ring "+NameOfRing+"=0,("+Variables+"),dp;"); |
---|
| 99 | } |
---|
| 100 | else |
---|
| 101 | { |
---|
| 102 | if(size(#)==1) |
---|
| 103 | { |
---|
| 104 | execute( "ring " + NameOfRing + "=(0," + #[1] + "),(" + Variables + "),dp;" ); |
---|
| 105 | } |
---|
| 106 | else |
---|
| 107 | { |
---|
| 108 | if( (size(#[1])!=0)&&(#[1]!=" ") ) |
---|
| 109 | { |
---|
| 110 | execute( "ring " + NameOfRing + "=(0," + #[1] + "),(" + Variables + "),("+#[2]+");" ); |
---|
| 111 | } |
---|
| 112 | else |
---|
| 113 | { |
---|
| 114 | execute( "ring " + NameOfRing + "=0,("+Variables+"),("+#[2]+");" ); |
---|
| 115 | }; |
---|
| 116 | }; |
---|
| 117 | }; |
---|
| 118 | keepring(basering); |
---|
[5831fa9] | 119 | } |
---|
[4713c6] | 120 | example |
---|
| 121 | {"EXAMPLE:";echo = 2; |
---|
| 122 | string v="x,y,z"; |
---|
| 123 | string p="q,p"; |
---|
| 124 | string Ord ="c,lp"; |
---|
[a09c2a0] | 125 | //---------------------------------- |
---|
[4713c6] | 126 | declare("Ring_1",v); |
---|
| 127 | print(nameof(basering)); |
---|
| 128 | print(basering); |
---|
[a09c2a0] | 129 | //---------------------------------- |
---|
[4713c6] | 130 | declare("Ring_2",v,p); |
---|
| 131 | print(basering); |
---|
| 132 | print(nameof(basering)); |
---|
[a09c2a0] | 133 | //---------------------------------- |
---|
[4713c6] | 134 | declare("Ring_3",v,p,Ord); |
---|
| 135 | print(basering); |
---|
| 136 | print(nameof(basering)); |
---|
[a09c2a0] | 137 | //---------------------------------- |
---|
[4713c6] | 138 | declare("Ring_4",v,"",Ord); |
---|
| 139 | print(basering); |
---|
| 140 | print(nameof(basering)); |
---|
[a09c2a0] | 141 | //---------------------------------- |
---|
[4713c6] | 142 | declare("Ring_5",v," ",Ord); |
---|
| 143 | print(basering); |
---|
| 144 | print(nameof(basering)); |
---|
| 145 | }; |
---|
| 146 | // |
---|
[a09c2a0] | 147 | // maybe reasonable to add this in declare |
---|
[4713c6] | 148 | // |
---|
| 149 | // print("Please enter your representation matrix in the following form: |
---|
| 150 | // module R=[1st row],[2nd row],..."); |
---|
| 151 | // print("Type the command: R=transpose(R)"); |
---|
| 152 | // print(" To compute controllability please enter: control(R)"); |
---|
| 153 | // print(" To compute autonomy please enter: autonom(R)"); |
---|
| 154 | //------------------------------------------------------------------------- |
---|
[0ceca2] | 155 | |
---|
[4713c6] | 156 | static proc space(int n) |
---|
[2e4e6ec] | 157 | "USAGE:spase(n); |
---|
| 158 | n: integer, number of needed spaces |
---|
| 159 | RETURN: string consisting of n spaces |
---|
| 160 | NOTE: the procedure is used in the procedure 'view' to have a better formatted output |
---|
| 161 | " |
---|
[4713c6] | 162 | { |
---|
| 163 | int i; |
---|
| 164 | string s=""; |
---|
| 165 | for(i=1;i<=n;i++) |
---|
| 166 | { |
---|
| 167 | s=s+" "; |
---|
| 168 | }; |
---|
| 169 | return(s); |
---|
| 170 | }; |
---|
| 171 | //----------------------------------------------------------------------------- |
---|
[2e4e6ec] | 172 | proc view(M) |
---|
| 173 | "USAGE: view(M); |
---|
| 174 | M: any type |
---|
| 175 | RETURN: no return value |
---|
| 176 | PURPOSE: procedure for ( well-) formatted output of modules, matrices, lists of modules, matrices; |
---|
| 177 | shows everything even if entries are long |
---|
| 178 | NOTE: in case of other types( not 'module', 'matrix', 'list') works just as standard 'print' procedure |
---|
| 179 | EXAMPLE: example view; shows an example |
---|
| 180 | " |
---|
[4713c6] | 181 | { |
---|
[d5f9d98] | 182 | // to be replaced with something more feasible |
---|
[4713c6] | 183 | if ( (typeof(M)=="module")||(typeof(M)=="matrix") ) |
---|
| 184 | { |
---|
| 185 | int @R=nrows(M); |
---|
| 186 | int @C=ncols(M); |
---|
| 187 | int i; |
---|
| 188 | int j; |
---|
| 189 | list MaxLength=list(); |
---|
| 190 | int Size=0; |
---|
| 191 | int max; |
---|
| 192 | string s; |
---|
| 193 | |
---|
| 194 | for(i=1;i<=@C;i++) |
---|
| 195 | { |
---|
| 196 | max=0; |
---|
| 197 | |
---|
| 198 | for(j=1;j<=@R;j++) |
---|
| 199 | { |
---|
| 200 | Size=size( string( M[j,i] ) ); |
---|
| 201 | if( Size>max ) |
---|
| 202 | { |
---|
| 203 | max=Size; |
---|
| 204 | }; |
---|
| 205 | }; |
---|
| 206 | MaxLength[i] = max; |
---|
| 207 | }; |
---|
| 208 | |
---|
| 209 | for(i=1;i<=@R;i++) |
---|
| 210 | { |
---|
| 211 | s=""; |
---|
| 212 | for(j=1;j<@C;j++) |
---|
| 213 | { |
---|
| 214 | s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) ) +","; |
---|
| 215 | }; |
---|
| 216 | |
---|
| 217 | s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) ); |
---|
| 218 | |
---|
| 219 | if (i!=@R) |
---|
| 220 | { |
---|
| 221 | s=s+","; |
---|
| 222 | }; |
---|
| 223 | print(s); |
---|
| 224 | }; |
---|
| 225 | |
---|
| 226 | return(); |
---|
| 227 | }; |
---|
| 228 | |
---|
| 229 | if(typeof(M)=="list") |
---|
| 230 | { |
---|
| 231 | int sz=size(M); |
---|
| 232 | int i; |
---|
| 233 | for(i=1;i<=sz;i++) |
---|
| 234 | { |
---|
[2e4e6ec] | 235 | view(M[i]); |
---|
[4713c6] | 236 | print(""); |
---|
| 237 | }; |
---|
| 238 | |
---|
| 239 | return(); |
---|
| 240 | }; |
---|
| 241 | print(M); |
---|
| 242 | return(); |
---|
[5831fa9] | 243 | } |
---|
[2e4e6ec] | 244 | example |
---|
| 245 | {"EXAMPLE:";echo = 2; |
---|
| 246 | ring r; |
---|
| 247 | matrix M[1][3] = x,y,z; |
---|
| 248 | print(M); |
---|
| 249 | view(M); |
---|
| 250 | }; |
---|
[4713c6] | 251 | //-------------------------------------------------------------------------- |
---|
| 252 | proc RightKernel(matrix M) |
---|
[0ceca2] | 253 | "USAGE: RightKernel(M); M a matrix |
---|
[d5f9d98] | 254 | PURPOSE: computes the right kernel of matrix M (a module of all elements v such that Mv=0) |
---|
| 255 | RETURN: module |
---|
| 256 | EXAMPLE: example RightKernel; shows an example |
---|
[2e4e6ec] | 257 | " |
---|
[4713c6] | 258 | { |
---|
[596ab3] | 259 | return(modulo(M,std(0))); |
---|
[5831fa9] | 260 | } |
---|
[2e4e6ec] | 261 | example |
---|
| 262 | {"EXAMPLE:";echo = 2; |
---|
[0ceca2] | 263 | ring r = 0,(x,y,z),dp; |
---|
[2e4e6ec] | 264 | matrix M[1][3] = x,y,z; |
---|
| 265 | print(M); |
---|
[0ceca2] | 266 | matrix R = RightKernel(M); |
---|
| 267 | print(R); |
---|
| 268 | // check: |
---|
| 269 | print(M*R); |
---|
[2e4e6ec] | 270 | }; |
---|
[4713c6] | 271 | //------------------------------------------------------------------------- |
---|
| 272 | proc LeftKernel(matrix M) |
---|
[0ceca2] | 273 | "USAGE: LeftKernel(M); M a matrix |
---|
[d5f9d98] | 274 | PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0) |
---|
| 275 | RETURN: module |
---|
[2e4e6ec] | 276 | EXAMPLE: example LeftKernel; shows an example |
---|
| 277 | " |
---|
[4713c6] | 278 | { |
---|
[596ab3] | 279 | return( transpose( modulo( transpose(M),std(0) ) ) ); |
---|
[5831fa9] | 280 | } |
---|
[2e4e6ec] | 281 | example |
---|
| 282 | {"EXAMPLE:";echo = 2; |
---|
[0ceca2] | 283 | ring r= 0,(x,y,z),dp; |
---|
[2e4e6ec] | 284 | matrix M[3][1] = x,y,z; |
---|
| 285 | print(M); |
---|
[0ceca2] | 286 | matrix L = LeftKernel(M); |
---|
| 287 | print(L); |
---|
| 288 | // check: |
---|
| 289 | print(L*M); |
---|
[2e4e6ec] | 290 | }; |
---|
[4713c6] | 291 | //------------------------------------------------------------------------ |
---|
[501f7c5] | 292 | proc LeftInverse(module M) |
---|
| 293 | "USAGE: LeftInverse(M); M a module |
---|
[d5f9d98] | 294 | PURPOSE: computes such a matrix L, that LM == Id; |
---|
| 295 | RETURN: module |
---|
[0ceca2] | 296 | EXAMPLE: example LeftInverse; shows an example |
---|
[501f7c5] | 297 | NOTE: exists only in the case when Id \subseteq M! |
---|
[2e4e6ec] | 298 | " |
---|
[e50c50] | 299 | { |
---|
[0c824aa] | 300 | // it works also for the NC case; |
---|
[0ceca2] | 301 | int NCols = ncols(M); |
---|
[0ae788] | 302 | module Id = freemodule(NCols); |
---|
[0ceca2] | 303 | module N = transpose(M); |
---|
[0e200a] | 304 | intvec old_opt=Opt_Our(); |
---|
[501f7c5] | 305 | Id = std(Id); |
---|
| 306 | matrix T; |
---|
[0ceca2] | 307 | // check the correctness (Id \subseteq M) |
---|
[0c824aa] | 308 | // via dimension: dim (M) = -1! |
---|
[0ceca2] | 309 | int d = dim_Our(N); |
---|
| 310 | if (d != -1) |
---|
[501f7c5] | 311 | { |
---|
[0c824aa] | 312 | // No left inverse exists |
---|
[501f7c5] | 313 | return(matrix(0)); |
---|
| 314 | } |
---|
| 315 | matrix T2 = lift(N, Id); |
---|
| 316 | T2 = transpose(T2); |
---|
[0c824aa] | 317 | option(set,old_opt); // set the options back |
---|
[501f7c5] | 318 | return(T2); |
---|
[5831fa9] | 319 | } |
---|
[e50c50] | 320 | example |
---|
[0ceca2] | 321 | { // trivial example: |
---|
| 322 | "EXAMPLE:";echo =2; |
---|
| 323 | ring r = 0,(x,z),dp; |
---|
[2e4e6ec] | 324 | matrix M[2][1] = 1,x2z; |
---|
| 325 | print(M); |
---|
| 326 | print( LeftInverse(M) ); |
---|
[0ceca2] | 327 | kill r; |
---|
| 328 | // derived from the exTwoPendula |
---|
| 329 | ring r=(0,m1,m2,M,g,L1,L2),Dt,dp; |
---|
| 330 | matrix U[3][1]; |
---|
| 331 | U[1,1]=(-L2)*Dt^4+(g)*Dt^2; |
---|
| 332 | U[2,1]=(-L1)*Dt^4+(g)*Dt^2; |
---|
| 333 | U[3,1]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2); |
---|
| 334 | module M = module(U); |
---|
| 335 | module L = LeftInverse(M); |
---|
| 336 | print(L); |
---|
[2e4e6ec] | 337 | }; |
---|
[999541] | 338 | //----------------------------------------------------------------------- |
---|
[0e200a] | 339 | proc RightInverse(module R) |
---|
[d5f9d98] | 340 | "USAGE: RightInverse(M); M a module |
---|
| 341 | PURPOSE: computes such a matrix L, that ML == Id; |
---|
| 342 | RETURN: module |
---|
| 343 | EXAMPLE: example RightInverse; shows an example |
---|
| 344 | NOTE: exists only in the case when Id \subseteq M! |
---|
| 345 | " |
---|
[0e200a] | 346 | { |
---|
| 347 | return(transpose(LeftInverse(transpose(R)))); |
---|
| 348 | } |
---|
[d5f9d98] | 349 | example |
---|
[0c824aa] | 350 | { "EXAMPLE:";echo =2; |
---|
| 351 | // trivial example: |
---|
[d5f9d98] | 352 | ring r = 0,(x,z),dp; |
---|
| 353 | matrix M[1][2] = 1,x2+z; |
---|
| 354 | print(M); |
---|
| 355 | print( RightInverse(M) ); |
---|
| 356 | kill r; |
---|
| 357 | // derived from the exTwoPendula |
---|
| 358 | ring r=(0,m1,m2,M,g,L1,L2),Dt,dp; |
---|
| 359 | matrix U[1][3]; |
---|
| 360 | U[1,1]=(-L2)*Dt^4+(g)*Dt^2; |
---|
| 361 | U[1,2]=(-L1)*Dt^4+(g)*Dt^2; |
---|
| 362 | U[1,3]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2); |
---|
| 363 | module M = module(U); |
---|
| 364 | module L = RightInverse(M); |
---|
| 365 | print(L); |
---|
| 366 | }; |
---|
[0e200a] | 367 | //----------------------------------------------------------------------- |
---|
[f52c2f] | 368 | static proc dim_Our(module R) |
---|
| 369 | { |
---|
[501f7c5] | 370 | int d; |
---|
[0ceca2] | 371 | if (attrib(R,"isSB")<>1) |
---|
| 372 | { |
---|
| 373 | R = std(R); |
---|
| 374 | } |
---|
[0e200a] | 375 | d = dim(R); |
---|
[501f7c5] | 376 | return(d); |
---|
[f52c2f] | 377 | } |
---|
| 378 | //----------------------------------------------------------------------- |
---|
| 379 | static proc Ann_Our(module R) |
---|
| 380 | { |
---|
[0e200a] | 381 | return(Ann(R)); |
---|
[f52c2f] | 382 | } |
---|
| 383 | //----------------------------------------------------------------------- |
---|
[501f7c5] | 384 | static proc Ext_Our(int i, module R, list #) |
---|
[999541] | 385 | { |
---|
[0c824aa] | 386 | // mimicking 'Ext_R' from homolog.lib |
---|
[501f7c5] | 387 | int ExtraArg = ( size(#)>0 ); |
---|
| 388 | if (ExtraArg) |
---|
[999541] | 389 | { |
---|
| 390 | return( Ext_R(i,R,#[1]) ); |
---|
[501f7c5] | 391 | } |
---|
| 392 | else |
---|
| 393 | { |
---|
| 394 | return( Ext_R(i,R) ); |
---|
| 395 | } |
---|
[5831fa9] | 396 | } |
---|
[999541] | 397 | //------------------------------------------------------------------------ |
---|
[f52c2f] | 398 | static proc is_zero_Our |
---|
[999541] | 399 | { |
---|
[0c824aa] | 400 | //just a copy of 'is_zero' from "poly.lib" patched with GKdim |
---|
[f52c2f] | 401 | int d = dim_Our(std(#[1])); |
---|
| 402 | int a = ( d==-1 ); |
---|
| 403 | if( size(#) >1 ) { list L=a,d; return(L); } |
---|
| 404 | return(a); |
---|
| 405 | // return( is_zero(R) ) ; |
---|
[999541] | 406 | }; |
---|
[e50c50] | 407 | //------------------------------------------------------------------------ |
---|
[0e200a] | 408 | static proc control_output(int i, int NVars, module R, module Ext_1, list Gen) |
---|
| 409 | //static proc control_output(int i, int NVars, module R, module Ext_1) |
---|
[0ceca2] | 410 | "USAGE: control_output(i, NVars, R, Ext_1), where |
---|
[2e4e6ec] | 411 | i: integer, number of first nonzero Ext or |
---|
[d5f9d98] | 412 | just a number of variables in a base ring + 1 in case that all the Exts are zero, |
---|
[0ceca2] | 413 | NVars: integer, number of variables in a base ring, |
---|
| 414 | R: module R (cokernel representation), |
---|
[2e4e6ec] | 415 | Ext_1: module, the first Ext(its cokernel representation) |
---|
| 416 | RETURN: list with all the contollability properties of the system which is to be returned in 'control' procedure |
---|
| 417 | NOTE: this procedure is used in 'control' procedure |
---|
| 418 | " |
---|
[4713c6] | 419 | { |
---|
[0ceca2] | 420 | // TODO: NVars to be replaced with the global hom. dimension of basering!!! |
---|
[0c824aa] | 421 | // Is not clear what to do with gl.dim of qrings |
---|
[0ceca2] | 422 | string DofS = "dimension of the system:"; |
---|
| 423 | string Fn = "number of first nonzero Ext:"; |
---|
[0c824aa] | 424 | string Gen_mes = "Parameter constellations which might lead to a non-controllable system:"; |
---|
[0e200a] | 425 | |
---|
| 426 | module RK = RightKernel(R); |
---|
| 427 | int d=dim_Our(std(transpose(R))); |
---|
| 428 | |
---|
[0ceca2] | 429 | if (i==1) |
---|
[4713c6] | 430 | { |
---|
| 431 | return( |
---|
[999541] | 432 | list ( Fn, |
---|
| 433 | i, |
---|
[4713c6] | 434 | "not controllable , image representation for controllable part:", |
---|
| 435 | RK, |
---|
| 436 | "kernel representation for controllable part:", |
---|
| 437 | LeftKernel( RK ), |
---|
| 438 | "obstruction to controllability", |
---|
| 439 | Ext_1, |
---|
[501f7c5] | 440 | "annihilator of torsion module (of obstruction to controllability)", |
---|
[f52c2f] | 441 | Ann_Our(Ext_1), |
---|
[4713c6] | 442 | DofS, |
---|
| 443 | d |
---|
| 444 | ) |
---|
| 445 | ); |
---|
| 446 | }; |
---|
| 447 | |
---|
| 448 | if(i>NVars) |
---|
[0e200a] | 449 | { |
---|
[999541] | 450 | return( list( Fn, |
---|
| 451 | -1, |
---|
[0e200a] | 452 | "strongly controllable(flat), image representation:", |
---|
[e50c50] | 453 | RK, |
---|
| 454 | "left inverse to image representation:", |
---|
| 455 | LeftInverse(RK), |
---|
[0e200a] | 456 | DofS, |
---|
| 457 | d, |
---|
| 458 | Gen_mes, |
---|
| 459 | Gen) |
---|
[4713c6] | 460 | ); |
---|
| 461 | }; |
---|
| 462 | |
---|
| 463 | // |
---|
| 464 | //now i<=NVars |
---|
| 465 | // |
---|
| 466 | |
---|
| 467 | if( (i==2) ) |
---|
| 468 | { |
---|
[999541] | 469 | return( list( Fn, |
---|
| 470 | i, |
---|
| 471 | "controllable, not reflexive, image representation:", |
---|
[0e200a] | 472 | RK, |
---|
[4713c6] | 473 | DofS, |
---|
[0e200a] | 474 | d, |
---|
| 475 | Gen_mes, |
---|
| 476 | Gen) |
---|
[4713c6] | 477 | ); |
---|
| 478 | }; |
---|
| 479 | |
---|
| 480 | if( (i>=3) ) |
---|
| 481 | { |
---|
[999541] | 482 | return( list ( Fn, |
---|
| 483 | i, |
---|
[4713c6] | 484 | "reflexive, not strongly controllable, image representation:", |
---|
[0e200a] | 485 | RK, |
---|
[4713c6] | 486 | DofS, |
---|
[0e200a] | 487 | d, |
---|
| 488 | Gen_mes, |
---|
| 489 | Gen) |
---|
[4713c6] | 490 | ); |
---|
[0ceca2] | 491 | }; |
---|
[4713c6] | 492 | }; |
---|
| 493 | //------------------------------------------------------------------------- |
---|
| 494 | |
---|
| 495 | proc control(module R) |
---|
[d5f9d98] | 496 | "USAGE: control(R); R a module (R is the matrix of the system of equations to be investigated) |
---|
| 497 | PURPOSE: compute the list of all the properties concerning controllability of the system (behavior), represented by the matrix R |
---|
| 498 | RETURN: list |
---|
[2e4e6ec] | 499 | EXAMPLE: example control; shows an example |
---|
| 500 | " |
---|
[4713c6] | 501 | { |
---|
| 502 | int i; |
---|
| 503 | int NVars=nvars(basering); |
---|
[0ceca2] | 504 | // TODO: NVars to be replaced with the global hom. dimension of basering!!! |
---|
[4713c6] | 505 | int ExtIsZero; |
---|
[0e200a] | 506 | intvec v=Opt_Our(); |
---|
| 507 | module R_std=std(R); |
---|
| 508 | module Ext_1 = std(Ext_Our(1,R_std)); |
---|
[4713c6] | 509 | |
---|
[999541] | 510 | ExtIsZero=is_zero_Our(Ext_1); |
---|
[4713c6] | 511 | i=1; |
---|
| 512 | while( (ExtIsZero) && (i<=NVars) ) |
---|
| 513 | { |
---|
| 514 | i++; |
---|
[0e200a] | 515 | ExtIsZero = is_zero_Our( Ext_Our(i,R_std) ); |
---|
[4713c6] | 516 | }; |
---|
[0e200a] | 517 | matrix T=lift(R,R_std); |
---|
| 518 | list l=genericity(T); |
---|
| 519 | option(set,v); |
---|
| 520 | |
---|
| 521 | return( control_output( i, NVars, R, Ext_1, l ) ); |
---|
[5831fa9] | 522 | } |
---|
[4713c6] | 523 | example |
---|
| 524 | {"EXAMPLE:";echo = 2; |
---|
[82d3aec] | 525 | //Wind Tunnel |
---|
| 526 | ring A = (0,a, omega, zeta, k),(D1, delta),dp; |
---|
| 527 | module R; |
---|
| 528 | R = [D1+a, -k*a*delta, 0, 0], |
---|
| 529 | [0, D1, -1, 0], |
---|
| 530 | [0, omega^2, D1+2*zeta*omega, -omega^2]; |
---|
| 531 | R=transpose(R); |
---|
| 532 | view(R); |
---|
| 533 | view(control(R)); |
---|
[4713c6] | 534 | }; |
---|
[0c824aa] | 535 | //-------------------------------------------------------------------------- |
---|
[0e200a] | 536 | proc control2(module R) |
---|
[d5f9d98] | 537 | "USAGE: control2(R); R a module (R is the matrix of the system of equations to be investigated) |
---|
| 538 | PURPOSE: computes list of all the properties concerning controllability of the system (behavior), represented by the matrix R |
---|
| 539 | RETURN: list |
---|
[eccbb1] | 540 | EXAMPLE: example control2; shows an example |
---|
[0c824aa] | 541 | NOTE: same as control(R); but using dimensions, this procedure only works for full row rank matrices |
---|
[0e200a] | 542 | " |
---|
| 543 | { |
---|
[6e9d48a] | 544 | if( nrows(R) != colrank(transpose(R)) ) |
---|
[2b9b32] | 545 | { |
---|
| 546 | return ("control2 cannot be applied, since R does not have full row rank"); |
---|
| 547 | } |
---|
[6e9d48a] | 548 | intvec v=Opt_Our(); |
---|
[0e200a] | 549 | module R_std=std(R); |
---|
| 550 | int d=dim_Our(R_std); |
---|
| 551 | int NVars=nvars(basering); |
---|
| 552 | int i=NVars-d; |
---|
[6e9d48a] | 553 | module Ext_1=std(Ext_Our(1,R_std)); |
---|
[0e200a] | 554 | matrix T=lift(R,R_std); |
---|
| 555 | list l=genericity(T); |
---|
| 556 | option(set, v); |
---|
| 557 | return( control_output( i, NVars, R, Ext_1, l)); |
---|
| 558 | } |
---|
| 559 | example |
---|
| 560 | {"EXAMPLE:";echo = 2; |
---|
| 561 | //Wind Tunnel |
---|
| 562 | ring A = (0,a, omega, zeta, k),(D1, delta),dp; |
---|
| 563 | module R; |
---|
| 564 | R = [D1+a, -k*a*delta, 0, 0], |
---|
| 565 | [0, D1, -1, 0], |
---|
| 566 | [0, omega^2, D1+2*zeta*omega, -omega^2]; |
---|
| 567 | R=transpose(R); |
---|
| 568 | view(R); |
---|
| 569 | view(control2(R)); |
---|
| 570 | }; |
---|
[4713c6] | 571 | //------------------------------------------------------------------------ |
---|
[6e9d48a] | 572 | proc colrank(module M) |
---|
| 573 | "USAGE: proc colrank(M), M a matrix/module |
---|
[2b9b32] | 574 | PURPOSE: compute the column rank of M as of matrix |
---|
| 575 | RETURN: int |
---|
[6e9d48a] | 576 | NOTE: this procedure uses bareiss-algorithm which might not terminate in some cases |
---|
[2b9b32] | 577 | " |
---|
[eccbb1] | 578 | { |
---|
[2b9b32] | 579 | module M_red = bareiss(M)[1]; |
---|
| 580 | int NCols_red = ncols(M_red); |
---|
[eccbb1] | 581 | return (NCols_red); |
---|
| 582 | } |
---|
[2b9b32] | 583 | example |
---|
| 584 | {"EXAMPLE: ";echo = 2; |
---|
[0c824aa] | 585 | // de Rham complex |
---|
[2b9b32] | 586 | ring r=0,(D(1..3)),dp; |
---|
| 587 | module R; |
---|
| 588 | R=[0,-D(3),D(2)], |
---|
| 589 | [D(3),0,-D(1)], |
---|
| 590 | [-D(2),D(1),0]; |
---|
| 591 | R=transpose(R); |
---|
[6e9d48a] | 592 | colrank(R); |
---|
[2b9b32] | 593 | }; |
---|
| 594 | |
---|
[eccbb1] | 595 | //------------------------------------------------------------------------ |
---|
| 596 | static proc autonom_output( int i, int NVars, module RC, int R_rank ) |
---|
[2b9b32] | 597 | "USAGE: proc autonom_output(i, NVars, RC, R_rank) |
---|
[2e4e6ec] | 598 | i: integer, number of first nonzero Ext or |
---|
| 599 | just number of variables in a base ring + 1 in case that all the Exts are zero |
---|
| 600 | NVars: integer, number of variables in a base ring |
---|
[2b9b32] | 601 | RC: module, kernel-representation of controllable part of the system |
---|
[6e9d48a] | 602 | R_rank: integer, column rank of the representation matrix |
---|
[2b9b32] | 603 | PURPOSE: compute all the autonomy properties of the system which is to be returned in 'autonom' procedure |
---|
| 604 | RETURN: list |
---|
[2e4e6ec] | 605 | NOTE: this procedure is used in 'autonom' procedure |
---|
| 606 | " |
---|
[4713c6] | 607 | { |
---|
| 608 | int d=NVars-i;//that is the dimension of the system |
---|
| 609 | string DofS="dimension of the system:"; |
---|
[999541] | 610 | string Fn = "number of first nonzero Ext:"; |
---|
[4713c6] | 611 | if(i==0) |
---|
| 612 | { |
---|
[999541] | 613 | return( list( Fn, |
---|
| 614 | i, |
---|
[4713c6] | 615 | "not autonomous", |
---|
[0e200a] | 616 | "kernel representation for controllable part", |
---|
| 617 | RC, |
---|
[6e9d48a] | 618 | "column rank of the matrix", |
---|
[eccbb1] | 619 | R_rank, |
---|
[4713c6] | 620 | DofS, |
---|
| 621 | d ) |
---|
| 622 | ); |
---|
| 623 | }; |
---|
| 624 | |
---|
| 625 | if( i>NVars ) |
---|
| 626 | { |
---|
[999541] | 627 | return( list( Fn, |
---|
| 628 | -1, |
---|
[4713c6] | 629 | "trivial", |
---|
| 630 | DofS, |
---|
| 631 | d ) |
---|
| 632 | ); |
---|
| 633 | }; |
---|
| 634 | |
---|
| 635 | // |
---|
| 636 | //now i<=NVars |
---|
| 637 | // |
---|
| 638 | |
---|
| 639 | |
---|
| 640 | if( i==1 ) |
---|
[0ceca2] | 641 | // in case that NVars==1 there is no sense to consider the notion |
---|
| 642 | // of strongly autonomous behavior, because it does not imply |
---|
| 643 | // that system is overdetermined in this case |
---|
[4713c6] | 644 | { |
---|
[999541] | 645 | return( list ( Fn, |
---|
| 646 | i, |
---|
[4713c6] | 647 | "autonomous, not overdetermined", |
---|
| 648 | DofS, |
---|
| 649 | d ) |
---|
| 650 | ); |
---|
| 651 | }; |
---|
| 652 | |
---|
| 653 | if( i==NVars ) |
---|
| 654 | { |
---|
[999541] | 655 | return( list( Fn, |
---|
| 656 | i, |
---|
[0e200a] | 657 | "strongly autonomous(fin. dimensional),in particular overdetermined", |
---|
[4713c6] | 658 | DofS, |
---|
| 659 | d) |
---|
| 660 | ); |
---|
| 661 | }; |
---|
| 662 | |
---|
| 663 | if( i<NVars ) |
---|
| 664 | { |
---|
[999541] | 665 | return( list ( Fn, |
---|
| 666 | i, |
---|
[4713c6] | 667 | "overdetermined, not strongly autonomous", |
---|
| 668 | DofS, |
---|
| 669 | d) |
---|
| 670 | ); |
---|
| 671 | }; |
---|
| 672 | }; |
---|
| 673 | //-------------------------------------------------------------------------- |
---|
| 674 | proc autonom2(module R) |
---|
[0ceca2] | 675 | "USAGE: autonom2(R); R a module (R is a matrix of the system of equations which is to be investigated) |
---|
[d5f9d98] | 676 | PURPOSE: computes the list of all the properties concerning autonomy of the system (behavior), represented by the matrix R |
---|
| 677 | RETURN: list |
---|
| 678 | NOTE: this procedure is analogous to 'autonom' but uses dimension calculations |
---|
[2e4e6ec] | 679 | EXAMPLE: example autonom2; shows an example |
---|
| 680 | " |
---|
[4713c6] | 681 | { |
---|
| 682 | int d; |
---|
| 683 | int NVars = nvars(basering); |
---|
[0e200a] | 684 | module RT = transpose(R); |
---|
| 685 | module RC; //for computation of controllable part if if exists |
---|
[0c824aa] | 686 | int R_rank = ncols(R); |
---|
[2b9b32] | 687 | d = dim_Our( std(RT) ); //this is the dimension of the system |
---|
| 688 | int i = NVars-d; //First non-zero Ext |
---|
| 689 | if( d==0 ) |
---|
[0e200a] | 690 | { |
---|
| 691 | RC=LeftKernel(RightKernel(R)); |
---|
[6e9d48a] | 692 | R_rank=colrank(R); |
---|
[0e200a] | 693 | } |
---|
[0c824aa] | 694 | return( autonom_output(i,NVars,RC,R_rank) ); |
---|
[5831fa9] | 695 | } |
---|
[4713c6] | 696 | example |
---|
| 697 | {"EXAMPLE:"; echo = 2; |
---|
[6a02de] | 698 | //Cauchy |
---|
[4713c6] | 699 | ring r=0,(s1,s2,s3,s4),dp; |
---|
| 700 | module R= [s1,-s2], |
---|
| 701 | [s2, s1], |
---|
| 702 | [s3,-s4], |
---|
| 703 | [s4, s3]; |
---|
| 704 | R=transpose(R); |
---|
[2e4e6ec] | 705 | view( R ); |
---|
| 706 | view( autonom2(R) ); |
---|
[4713c6] | 707 | }; |
---|
[0ceca2] | 708 | //---------------------------------------------------------- |
---|
[4713c6] | 709 | proc autonom(module R) |
---|
[0ceca2] | 710 | "USAGE: autonom(R); R a module (R is a matrix of the system of equations which is to be investigated) |
---|
[2b9b32] | 711 | PURPOSE: find all the properties concerning autonomy of the system (behavior) represented by the matrix R |
---|
| 712 | RETURN: list |
---|
| 713 | EXAMPLE: example autonom; shows an example |
---|
[2e4e6ec] | 714 | " |
---|
[4713c6] | 715 | { |
---|
| 716 | int NVars=nvars(basering); |
---|
[0e200a] | 717 | int ExtIsZero; |
---|
| 718 | module RT=transpose(R); |
---|
| 719 | module RC; |
---|
[eccbb1] | 720 | int R_rank=ncols(R); |
---|
[0e200a] | 721 | ExtIsZero=is_zero_Our(Ext_Our(0,RT)); |
---|
[5831fa9] | 722 | int i=0; |
---|
[4713c6] | 723 | while( (ExtIsZero)&&(i<=NVars) ) |
---|
| 724 | { |
---|
| 725 | i++; |
---|
[0e200a] | 726 | ExtIsZero = is_zero_Our(Ext_Our(i,RT)); |
---|
[4713c6] | 727 | }; |
---|
[0c824aa] | 728 | if (i==0) |
---|
| 729 | { |
---|
| 730 | RC=LeftKernel(RightKernel(R)); |
---|
[6e9d48a] | 731 | R_rank=colrank(R); |
---|
[0c824aa] | 732 | } |
---|
[eccbb1] | 733 | return(autonom_output(i,NVars,RC,R_rank)); |
---|
[5831fa9] | 734 | } |
---|
[4713c6] | 735 | example |
---|
| 736 | {"EXAMPLE:"; echo = 2; |
---|
[0c824aa] | 737 | // Cauchy |
---|
[82d3aec] | 738 | ring r=0,(s1,s2,s3,s4),dp; |
---|
| 739 | module R= [s1,-s2], |
---|
| 740 | [s2, s1], |
---|
| 741 | [s3,-s4], |
---|
| 742 | [s4, s3]; |
---|
| 743 | R=transpose(R); |
---|
| 744 | view( R ); |
---|
| 745 | view( autonom(R) ); |
---|
| 746 | }; |
---|
| 747 | |
---|
[0ceca2] | 748 | //---------------------------------------------------------- |
---|
[82d3aec] | 749 | // |
---|
| 750 | //Some example rings with defined systems |
---|
[0ceca2] | 751 | //---------------------------------------------------------- |
---|
[82d3aec] | 752 | //autonomy: |
---|
| 753 | // |
---|
[0ceca2] | 754 | //---------------------------------------------------------- |
---|
[82d3aec] | 755 | proc exCauchy() |
---|
| 756 | { |
---|
| 757 | ring @r=0,(s1,s2),dp; |
---|
[4713c6] | 758 | module R= [s1,-s2], |
---|
| 759 | [s2, s1]; |
---|
| 760 | R=transpose(R); |
---|
[82d3aec] | 761 | export R; |
---|
| 762 | return(@r); |
---|
| 763 | }; |
---|
[0ceca2] | 764 | //---------------------------------------------------------- |
---|
[82d3aec] | 765 | proc exCauchy2() |
---|
| 766 | { |
---|
| 767 | ring @r=0,(s1,s2,s3,s4),dp; |
---|
[4713c6] | 768 | module R= [s1,-s2], |
---|
| 769 | [s2, s1], |
---|
| 770 | [s3,-s4], |
---|
| 771 | [s4, s3]; |
---|
| 772 | R=transpose(R); |
---|
[82d3aec] | 773 | export R; |
---|
| 774 | return(@r); |
---|
| 775 | }; |
---|
[0ceca2] | 776 | //---------------------------------------------------------- |
---|
[f691ff] | 777 | proc exZerz() |
---|
[82d3aec] | 778 | { |
---|
| 779 | ring @r=0,(d1,d2),dp; |
---|
[4713c6] | 780 | module R=[d1^2-d2], |
---|
| 781 | [d2^2-1]; |
---|
| 782 | R=transpose(R); |
---|
[82d3aec] | 783 | export R; |
---|
| 784 | return(@r); |
---|
[4713c6] | 785 | }; |
---|
[0ceca2] | 786 | //---------------------------------------------------------- |
---|
[82d3aec] | 787 | //control |
---|
[0c824aa] | 788 | //---------------------------------------------------------- |
---|
[6a02de] | 789 | proc ex1() |
---|
| 790 | { |
---|
| 791 | ring @r=0,(s1,s2,s3),dp; |
---|
| 792 | module R=[0,-s3,s2], |
---|
| 793 | [s3,0,-s1]; |
---|
| 794 | R=transpose(R); |
---|
| 795 | export R; |
---|
| 796 | return(@r); |
---|
| 797 | }; |
---|
[0ceca2] | 798 | //---------------------------------------------------------- |
---|
[6a02de] | 799 | proc ex2() |
---|
| 800 | { |
---|
| 801 | ring @r=0,(s1,s2,s3),dp; |
---|
| 802 | module R=[0,-s3,s2], |
---|
[0c824aa] | 803 | [s3,0,-s1], |
---|
| 804 | [-s2,s1,0]; |
---|
[6a02de] | 805 | R=transpose(R); |
---|
| 806 | export R; |
---|
| 807 | return(@r); |
---|
| 808 | }; |
---|
[0ceca2] | 809 | //---------------------------------------------------------- |
---|
[394dad9] | 810 | proc exAntenna() |
---|
| 811 | { |
---|
| 812 | ring @r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), dp; |
---|
| 813 | module R; |
---|
| 814 | R = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0], |
---|
| 815 | [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta], |
---|
| 816 | [0, 0, Dt, -K1, 0, 0, 0, 0, 0], |
---|
| 817 | [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta], |
---|
| 818 | [0, 0, 0, 0, Dt, -K1, 0, 0, 0], |
---|
| 819 | [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta]; |
---|
| 820 | |
---|
| 821 | R=transpose(R); |
---|
| 822 | export R; |
---|
| 823 | return(@r); |
---|
| 824 | }; |
---|
| 825 | |
---|
[0ceca2] | 826 | //---------------------------------------------------------- |
---|
[394dad9] | 827 | |
---|
| 828 | proc exEinstein() |
---|
| 829 | { |
---|
| 830 | ring @r = 0,(D(1..4)),dp; |
---|
| 831 | module R = |
---|
| 832 | [D(2)^2+D(3)^2-D(4)^2, D(1)^2, D(1)^2, -D(1)^2, -2*D(1)*D(2), 0, 0, -2*D(1)*D(3), 0, 2*D(1)*D(4)], |
---|
| 833 | [D(2)^2, D(1)^2+D(3)^2-D(4)^2, D(2)^2, -D(2)^2, -2*D(1)*D(2), -2*D(2)*D(3), 0, 0, 2*D(2)*D(4), 0], |
---|
| 834 | [D(3)^2, D(3)^2, D(1)^2+D(2)^2-D(4)^2, -D(3)^2, 0, -2*D(2)*D(3), 2*D(3)*D(4), -2*D(1)*D(3), 0, 0], |
---|
| 835 | [D(4)^2, D(4)^2, D(4)^2, D(1)^2+D(2)^2+D(3)^2, 0, 0, -2*D(3)*D(4), 0, -2*D(2)*D(4), -2*D(1)*D(4)], |
---|
| 836 | [0, 0, D(1)*D(2), -D(1)*D(2), D(3)^2-D(4)^2, -D(1)*D(3), 0, -D(2)*D(3), D(1)*D(4), D(2)*D(4)], |
---|
| 837 | [D(2)*D(3), 0, 0, -D(2)*D(3),-D(1)*D(3), D(1)^2-D(4)^2, D(2)*D(4), -D(1)*D(2), D(3)*D(4), 0], |
---|
| 838 | [D(3)*D(4), D(3)*D(4), 0, 0, 0, -D(2)*D(4), D(1)^2+D(2)^2, -D(1)*D(4), -D(2)*D(3), -D(1)*D(3)], |
---|
| 839 | [0, D(1)*D(3), 0, -D(1)*D(3), -D(2)*D(3), -D(1)*D(2), D(1)*D(4), D(2)^2-D(4)^2, 0, D(3)*D(4)], |
---|
| 840 | [D(2)*D(4), 0, D(2)*D(4), 0, -D(1)*D(4), -D(3)*D(4), -D(2)*D(3), 0, D(1)^2+D(3)^2, -D(1)*D(2)], |
---|
| 841 | [0, D(1)*D(4), D(1)*D(4), 0, -D(2)*D(4), 0, -D(1)*D(3), -D(3)*D(4), -D(1)*D(2), D(2)^2+D(3)^2]; |
---|
| 842 | |
---|
| 843 | R=transpose(R); |
---|
| 844 | export R; |
---|
| 845 | return(@r); |
---|
| 846 | }; |
---|
| 847 | |
---|
[0ceca2] | 848 | //---------------------------------------------------------- |
---|
[394dad9] | 849 | proc exFlexibleRod() |
---|
| 850 | { |
---|
| 851 | ring @r = 0,(D1, delta), dp; |
---|
| 852 | module R; |
---|
| 853 | R = [D1, -D1*delta, -1], [2*D1*delta, -D1-D1*delta^2, 0]; |
---|
| 854 | |
---|
| 855 | R=transpose(R); |
---|
| 856 | export R; |
---|
| 857 | return(@r); |
---|
| 858 | }; |
---|
| 859 | |
---|
[0ceca2] | 860 | //---------------------------------------------------------- |
---|
[394dad9] | 861 | proc exTwoPendula() |
---|
| 862 | { |
---|
| 863 | ring @r=(0,m1,m2,M,g,L1,L2),Dt,dp; |
---|
| 864 | module R = [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2], |
---|
| 865 | [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2], |
---|
| 866 | [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2]; |
---|
| 867 | |
---|
| 868 | R=transpose(R); |
---|
| 869 | export R; |
---|
| 870 | return(@r); |
---|
| 871 | }; |
---|
[0ceca2] | 872 | //---------------------------------------------------------- |
---|
[394dad9] | 873 | proc exWindTunnel() |
---|
| 874 | { |
---|
| 875 | ring @r = (0,a, omega, zeta, k),(D1, delta),dp; |
---|
| 876 | module R = [D1+a, -k*a*delta, 0, 0], |
---|
| 877 | [0, D1, -1, 0], |
---|
| 878 | [0, omega^2, D1+2*zeta*omega, -omega^2]; |
---|
| 879 | |
---|
| 880 | R=transpose(R); |
---|
| 881 | export R; |
---|
| 882 | return(@r); |
---|
| 883 | }; |
---|
[0ceca2] | 884 | //---------------------------------------------------------- |
---|
[0ae788] | 885 | proc genericity(matrix M) |
---|
[0c824aa] | 886 | "USAGE: genericity(M), M is a matrix/module |
---|
[d5f9d98] | 887 | PURPOSE: determine parametric expressions which have been assumed to be non-zero in the process of computing the Groebner basis |
---|
[eccbb1] | 888 | RETURN: list (of strings) |
---|
[a4f5ce] | 889 | NOTE: we strongly recommend to switch on the redSB and redTail options; |
---|
[d5f9d98] | 890 | @* the procedure is effective with the lift procedure for modules with parameters |
---|
[a4f5ce] | 891 | EXAMPLE: example genericity; shows an example |
---|
[0ae788] | 892 | " |
---|
[2b9b32] | 893 | { |
---|
| 894 | // returns "-", if there are no parameters! |
---|
| 895 | if (npars(basering)==0) |
---|
| 896 | { |
---|
| 897 | return("-"); |
---|
| 898 | } |
---|
[0c824aa] | 899 | list RT = evas_genericity(M); // list of strings |
---|
| 900 | if ((size(RT)==1) && (RT[1] == "")) |
---|
| 901 | { |
---|
| 902 | return("-"); |
---|
| 903 | } |
---|
[2b9b32] | 904 | return(RT); |
---|
| 905 | } |
---|
| 906 | example |
---|
| 907 | { // TwoPendula |
---|
| 908 | "EXAMPLE:"; echo = 2; |
---|
| 909 | ring r=(0,m1,m2,M,g,L1,L2),Dt,dp; |
---|
| 910 | module RR = |
---|
| 911 | [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2], |
---|
| 912 | [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2], |
---|
| 913 | [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2]; |
---|
| 914 | module R = transpose(RR); |
---|
| 915 | module SR = std(R); |
---|
| 916 | matrix T = lift(R,SR); |
---|
| 917 | genericity(T); |
---|
| 918 | //-- The result might be different when computing reduced bases: |
---|
| 919 | matrix T2; |
---|
| 920 | option(redSB); |
---|
| 921 | option(redTail); |
---|
| 922 | module SR2 = std(R); |
---|
| 923 | T2 = lift(R,SR2); |
---|
| 924 | genericity(T2); |
---|
| 925 | } |
---|
| 926 | //--------------------------------------------------------------- |
---|
[0c824aa] | 927 | static proc victors_genericity(matrix M) |
---|
[0ae788] | 928 | { |
---|
[eccbb1] | 929 | // returns "-", if there are no parameters! |
---|
[d5f9d98] | 930 | if (npars(basering)==0) |
---|
| 931 | { |
---|
[eccbb1] | 932 | return("-"); |
---|
[d5f9d98] | 933 | } |
---|
| 934 | int plevel = printlevel-voice+2; |
---|
[0ae788] | 935 | // M is a matrix over a ring with params and vars; |
---|
| 936 | ideal I = ideal(M); // a list of entries |
---|
[a4f5ce] | 937 | I = simplify(I,2); // delete 0's |
---|
| 938 | // decompose every coeff in every poly |
---|
| 939 | int i; |
---|
[0ae788] | 940 | int s = size(I); |
---|
[a4f5ce] | 941 | ideal NM; |
---|
[0ae788] | 942 | poly p; |
---|
[a4f5ce] | 943 | number num; |
---|
| 944 | int cl=1; |
---|
| 945 | intvec ZeroVec; ZeroVec[nvars(basering)] = 0; |
---|
| 946 | intvec W; |
---|
[d5f9d98] | 947 | ideal Numero, Denomiro; |
---|
| 948 | int cNu=0; int cDe=0; |
---|
[0ae788] | 949 | for (i=1; i<=s; i++) |
---|
| 950 | { |
---|
[a4f5ce] | 951 | // remove contents and add them as polys |
---|
| 952 | p = I[i]; |
---|
| 953 | W = leadexp(p); |
---|
| 954 | if (W == ZeroVec) // i.e. just a coef |
---|
| 955 | { |
---|
| 956 | num = denominator(leadcoef(p)); // from poly.lib |
---|
[d5f9d98] | 957 | NM[cl] = numerator(leadcoef(p)); |
---|
| 958 | dbprint(p,"numerator:"); |
---|
| 959 | dbprint(p, string(NM[cl])); |
---|
| 960 | cNu++; Numero[cNu]= NM[cl]; |
---|
[a4f5ce] | 961 | cl++; |
---|
[d5f9d98] | 962 | NM[cl] = num; // denominator |
---|
| 963 | dbprint(p,"denominator:"); |
---|
| 964 | dbprint(p, string(NM[cl])); |
---|
| 965 | cDe++; Denomiro[cDe]= NM[cl]; |
---|
[a4f5ce] | 966 | cl++; |
---|
| 967 | p = p - lead(p); // for the next cycle |
---|
| 968 | } |
---|
| 969 | if ( p!= 0) |
---|
| 970 | { |
---|
| 971 | num = content(p); |
---|
| 972 | p = p/num; |
---|
| 973 | NM[cl] = denominator(num); |
---|
[d5f9d98] | 974 | dbprint(p,"content denominator:"); |
---|
| 975 | dbprint(p, string(NM[cl])); |
---|
| 976 | cNu++; Numero[cNu]= NM[cl]; |
---|
[a4f5ce] | 977 | cl++; |
---|
| 978 | NM[cl] = numerator(num); |
---|
[d5f9d98] | 979 | dbprint(p,"content numerator:"); |
---|
| 980 | dbprint(p, string(NM[cl])); |
---|
| 981 | cDe++; Denomiro[cDe]= NM[cl]; |
---|
[a4f5ce] | 982 | cl++; |
---|
| 983 | } |
---|
[d5f9d98] | 984 | // it seems that the next elements will not have real influence |
---|
[0ae788] | 985 | while( p != 0) |
---|
| 986 | { |
---|
[d5f9d98] | 987 | NM[cl] = leadcoef(p); // should be all integer, i.e. non-rational |
---|
| 988 | dbprint(p,"coef:"); |
---|
| 989 | dbprint(p, string(NM[cl])); |
---|
[0ae788] | 990 | cl++; |
---|
| 991 | p = p - lead(p); |
---|
[a4f5ce] | 992 | } |
---|
| 993 | } |
---|
| 994 | NM = simplify(NM,4); // delete identical |
---|
[0ae788] | 995 | string newvars = parstr(basering); |
---|
| 996 | def save = basering; |
---|
| 997 | string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;"; |
---|
| 998 | execute(NewRing); |
---|
| 999 | // get params as variables |
---|
[f52c2f] | 1000 | // create a list of non-monomials |
---|
[a4f5ce] | 1001 | ideal @L; |
---|
[0ae788] | 1002 | ideal F; |
---|
[a4f5ce] | 1003 | ideal NM = imap(save,NM); |
---|
| 1004 | NM = simplify(NM,8); //delete multiples |
---|
[0ae788] | 1005 | poly p,q; |
---|
| 1006 | cl = 1; |
---|
| 1007 | int j, cf; |
---|
| 1008 | for(i=1; i<=size(NM);i++) |
---|
| 1009 | { |
---|
| 1010 | p = NM[i] - lead(NM[i]); |
---|
| 1011 | if (p!=0) |
---|
| 1012 | { |
---|
| 1013 | // L[cl] = p; |
---|
| 1014 | F = factorize(NM[i],1); //non-constant factors only |
---|
| 1015 | cf = 1; |
---|
| 1016 | // factorize every polynomial |
---|
[d5f9d98] | 1017 | // throw away every monomial from factorization (also constants from above ring) |
---|
[0ae788] | 1018 | for (j=1; j<=size(F);j++) |
---|
| 1019 | { |
---|
| 1020 | q = F[j]-lead(F[j]); |
---|
| 1021 | if (q!=0) |
---|
| 1022 | { |
---|
[a4f5ce] | 1023 | @L[cl] = F[j]; |
---|
[0ae788] | 1024 | cl++; |
---|
| 1025 | } |
---|
| 1026 | } |
---|
| 1027 | } |
---|
| 1028 | } |
---|
[501f7c5] | 1029 | // return the result [in string-format] |
---|
[d5f9d98] | 1030 | @L = simplify(@L,2+4+8); // skip zeroes, doubled and entries, diff. by a constant |
---|
[0ae788] | 1031 | list SL; |
---|
[a4f5ce] | 1032 | for (j=1; j<=size(@L);j++) |
---|
[0ae788] | 1033 | { |
---|
[a4f5ce] | 1034 | SL[j] = string(@L[j]); |
---|
[0ae788] | 1035 | } |
---|
| 1036 | setring save; |
---|
| 1037 | return(SL); |
---|
| 1038 | } |
---|
[0c824aa] | 1039 | //--------------------------------------------------------------- |
---|
| 1040 | static proc evas_genericity(matrix M) |
---|
| 1041 | { |
---|
| 1042 | // called from the main genericity proc |
---|
| 1043 | ideal I = ideal(M); |
---|
| 1044 | I = simplify(I,2+4); |
---|
| 1045 | int s = size(I); |
---|
| 1046 | ideal Den; |
---|
| 1047 | poly p; |
---|
| 1048 | int i; |
---|
| 1049 | for (i=1; i<=s; i++) |
---|
| 1050 | { |
---|
| 1051 | p = I[i]; |
---|
| 1052 | while (p !=0) |
---|
| 1053 | { |
---|
| 1054 | Den = Den, denominator(leadcoef(p)); |
---|
| 1055 | p = p-lead(p); |
---|
| 1056 | } |
---|
| 1057 | } |
---|
| 1058 | Den = simplify(Den,2+4); |
---|
| 1059 | string newvars = parstr(basering); |
---|
| 1060 | def save = basering; |
---|
| 1061 | string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;"; |
---|
| 1062 | execute(NewRing); |
---|
| 1063 | ideal F; |
---|
| 1064 | ideal Den = imap(save,Den); |
---|
| 1065 | Den = simplify(Den,2); |
---|
| 1066 | int s1 = size(Den); |
---|
| 1067 | for (i=1; i<=s1; i++) |
---|
| 1068 | { |
---|
| 1069 | if (Den[i] !=1) |
---|
| 1070 | { |
---|
| 1071 | F= F, factorize(Den[i],1); |
---|
| 1072 | } |
---|
| 1073 | } |
---|
| 1074 | F = simplify(F, 2+4+8); |
---|
| 1075 | ideal @L = F; |
---|
| 1076 | list SL; |
---|
| 1077 | int c,j; |
---|
| 1078 | string Mono; |
---|
| 1079 | c = 1; |
---|
| 1080 | for (j=1; j<=size(@L);j++) |
---|
| 1081 | { |
---|
| 1082 | if (leadcoef(@L[j]) <0) |
---|
| 1083 | { |
---|
| 1084 | @L[j] = -1*@L[j]; |
---|
| 1085 | } |
---|
| 1086 | if ( (@L[j] - lead(@L[j]))==0 ) //@L[j] is a monomial |
---|
| 1087 | { |
---|
| 1088 | Mono = Mono + string(@L[j])+ ","; // concatenation |
---|
| 1089 | } |
---|
| 1090 | else |
---|
| 1091 | { |
---|
| 1092 | c++; |
---|
| 1093 | SL[c] = string(@L[j]); |
---|
| 1094 | } |
---|
| 1095 | } |
---|
| 1096 | if (Mono!="") |
---|
| 1097 | { |
---|
| 1098 | Mono = Mono[1..size(Mono)-1]; // delete the last semicolon |
---|
| 1099 | } |
---|
| 1100 | SL[1] = Mono; |
---|
| 1101 | setring save; |
---|
| 1102 | return(SL); |
---|
[0ae788] | 1103 | } |
---|
[0c824aa] | 1104 | |
---|
[eccbb1] | 1105 | //--------------------------------------------------------------- |
---|
[0ae788] | 1106 | proc canonize(list L) |
---|
[f52c2f] | 1107 | "USAGE: canonize(L), L a list |
---|
[eccbb1] | 1108 | 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 |
---|
[d5f9d98] | 1109 | RETURN: list |
---|
[a4f5ce] | 1110 | ASSUME: L is the output of control/autonomy procedures |
---|
| 1111 | EXAMPLE: example canonize; shows an example |
---|
[0ae788] | 1112 | " |
---|
| 1113 | { |
---|
| 1114 | list M = L; |
---|
[0e200a] | 1115 | intvec v=Opt_Our(); |
---|
[0ae788] | 1116 | int s = size(L); |
---|
| 1117 | int i; |
---|
| 1118 | for (i=2; i<=s; i=i+2) |
---|
| 1119 | { |
---|
| 1120 | if (typeof(M[i])=="module") |
---|
| 1121 | { |
---|
| 1122 | M[i] = std(M[i]); |
---|
[0ceca2] | 1123 | // M[i] = prune(M[i]); // mimimal embedding: no need yet |
---|
[501f7c5] | 1124 | // M[i] = std(M[i]); |
---|
[0ae788] | 1125 | } |
---|
| 1126 | } |
---|
[0e200a] | 1127 | option(set, v); //set old values back |
---|
[0ae788] | 1128 | return(M); |
---|
| 1129 | } |
---|
| 1130 | example |
---|
| 1131 | { // TwoPendula with L1=L2=L |
---|
| 1132 | "EXAMPLE:"; echo = 2; |
---|
| 1133 | ring r=(0,m1,m2,M,g,L),Dt,dp; |
---|
| 1134 | module RR = |
---|
| 1135 | [m1*L*Dt^2, m2*L*Dt^2, -1, (M+m1+m2)*Dt^2], |
---|
| 1136 | [m1*L^2*Dt^2-m1*L*g, 0, 0, m1*L*Dt^2], |
---|
| 1137 | [0, m2*L^2*Dt^2-m2*L*g, 0, m2*L*Dt^2]; |
---|
| 1138 | module R = transpose(RR); |
---|
| 1139 | list C = control(R); |
---|
| 1140 | list CC = canonize(C); |
---|
| 1141 | view(CC); |
---|
| 1142 | } |
---|
[6e9d48a] | 1143 | |
---|
| 1144 | //---------------------------------------------------------------- |
---|
| 1145 | |
---|
| 1146 | static proc elementof (int i, intvec v) |
---|
| 1147 | { |
---|
| 1148 | int b=0; |
---|
| 1149 | for(int j=1;j<=nrows(v);j++) |
---|
| 1150 | { |
---|
| 1151 | if(v[j]==i) |
---|
| 1152 | { |
---|
| 1153 | b=1; |
---|
| 1154 | return (b); |
---|
| 1155 | } |
---|
| 1156 | } |
---|
| 1157 | return (b); |
---|
| 1158 | } |
---|
| 1159 | //----------------------------------------------------------------- |
---|
| 1160 | proc iostruct(module R) |
---|
| 1161 | "USAGE: iostruct( R ); R a module |
---|
| 1162 | RETURN: list L with entries: string s, intvec v, module P and module Q |
---|
| 1163 | 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 |
---|
| 1164 | NOTE: the procedure uses Bareiss algorithm which might not terminate in some cases |
---|
| 1165 | EXAMPLE: example iostruct; shows an example |
---|
| 1166 | " |
---|
| 1167 | { |
---|
| 1168 | list L = bareiss(R); |
---|
| 1169 | int R_rank = ncols(L[1]); |
---|
| 1170 | int NCols=ncols(R); |
---|
| 1171 | intvec v=L[2]; |
---|
| 1172 | int temp; |
---|
| 1173 | int NRows=nrows(v); |
---|
| 1174 | int i,j; |
---|
| 1175 | int b=1; |
---|
| 1176 | module P; |
---|
| 1177 | module Q; |
---|
| 1178 | int n=0; |
---|
| 1179 | |
---|
| 1180 | while(b==1) //sort v through bubblesort |
---|
| 1181 | { |
---|
| 1182 | b=0; |
---|
| 1183 | for(i=1;i<NRows;i++) |
---|
| 1184 | { |
---|
| 1185 | if(v[i]>v[i+1]) |
---|
| 1186 | { |
---|
| 1187 | temp=v[i]; |
---|
| 1188 | v[i]=v[i+1]; |
---|
| 1189 | v[i+1]=temp; |
---|
| 1190 | b=1; |
---|
| 1191 | } |
---|
| 1192 | } |
---|
| 1193 | } |
---|
| 1194 | P=R[v]; //generate P |
---|
| 1195 | for(i=1;i<=NCols;i++) //generate Q |
---|
| 1196 | { |
---|
| 1197 | if(elementof(i,v)==1) |
---|
| 1198 | { |
---|
| 1199 | i++; |
---|
| 1200 | continue; |
---|
| 1201 | } |
---|
| 1202 | Q=Q,R[i]; |
---|
| 1203 | } |
---|
| 1204 | Q=simplify(Q,2); |
---|
| 1205 | string s="The following components have been chosen as outputs: "; |
---|
| 1206 | return (list(s,v,P,Q)); |
---|
| 1207 | } |
---|
| 1208 | example |
---|
| 1209 | {"EXAMPLE:";echo = 2; |
---|
| 1210 | //Example Antenna |
---|
| 1211 | ring r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), (c,dp); |
---|
| 1212 | |
---|
| 1213 | module RR; |
---|
| 1214 | RR = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0], |
---|
| 1215 | [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta], |
---|
| 1216 | [0, 0, Dt, -K1, 0, 0, 0, 0, 0], |
---|
| 1217 | [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta], |
---|
| 1218 | [0, 0, 0, 0, Dt, -K1, 0, 0, 0], |
---|
| 1219 | [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta]; |
---|
| 1220 | module R = transpose(RR); |
---|
| 1221 | view(R); |
---|
| 1222 | view(iostruct(R)); |
---|
| 1223 | }; |
---|
| 1224 | |
---|
[501f7c5] | 1225 | //--------------------------------------------------------------- |
---|
| 1226 | static proc smdeg(matrix N) |
---|
| 1227 | // returns an intvec of length 2 with the index of an element of N with smallest degree |
---|
| 1228 | { |
---|
| 1229 | int n = nrows(N); |
---|
| 1230 | int m = ncols(N); |
---|
| 1231 | int d,d_temp; |
---|
| 1232 | intvec v; |
---|
| 1233 | int i,j; // counter |
---|
[0e200a] | 1234 | |
---|
[0c824aa] | 1235 | if (N==0) |
---|
[0e200a] | 1236 | { |
---|
[0c824aa] | 1237 | v = 1,1; |
---|
[0e200a] | 1238 | return(v); |
---|
| 1239 | } |
---|
[501f7c5] | 1240 | |
---|
[0c824aa] | 1241 | for (i=1; i<=n; i++) |
---|
| 1242 | // hier wird ein Element ausgewaehlt(!=0) und mit dessen Grad gestartet |
---|
[501f7c5] | 1243 | { |
---|
[0c824aa] | 1244 | for (j=1; j<=m; j++) |
---|
[501f7c5] | 1245 | { |
---|
| 1246 | if( deg(N[i,j])!=-1 ) |
---|
| 1247 | { |
---|
[0c824aa] | 1248 | d=deg(N[i,j]); |
---|
| 1249 | break; |
---|
[501f7c5] | 1250 | } |
---|
| 1251 | } |
---|
[0c824aa] | 1252 | if (d != -1) |
---|
[501f7c5] | 1253 | { |
---|
| 1254 | break; |
---|
| 1255 | } |
---|
| 1256 | } |
---|
[0c824aa] | 1257 | for(i=1; i<=n; i++) |
---|
[501f7c5] | 1258 | { |
---|
[0c824aa] | 1259 | for(j=1; j<=m; j++) |
---|
[501f7c5] | 1260 | { |
---|
| 1261 | d_temp = deg(N[i,j]); |
---|
[0c824aa] | 1262 | if ( (d_temp < d) && (N[i,j]!=0) ) |
---|
[501f7c5] | 1263 | { |
---|
[0c824aa] | 1264 | d=d_temp; |
---|
[501f7c5] | 1265 | } |
---|
| 1266 | } |
---|
| 1267 | } |
---|
| 1268 | for (i=1; i<=n; i++) |
---|
| 1269 | { |
---|
[0c824aa] | 1270 | for (j=1; j<=m;j++) |
---|
[501f7c5] | 1271 | { |
---|
| 1272 | if ( (deg(N[i,j]) == d) && (N[i,j]!=0) ) |
---|
| 1273 | { |
---|
[0c824aa] | 1274 | v = i,j; |
---|
| 1275 | return(v); |
---|
[501f7c5] | 1276 | } |
---|
| 1277 | } |
---|
| 1278 | } |
---|
| 1279 | } |
---|
[d5f9d98] | 1280 | //--------------------------------------------------------------- |
---|
[501f7c5] | 1281 | static proc NoNon0Pol(vector v) |
---|
| 1282 | // returns 1, if there is only one non-zero element in v and 0 else |
---|
| 1283 | { |
---|
| 1284 | int i,j; |
---|
| 1285 | int n = nrows(v); |
---|
[0c824aa] | 1286 | for( j=1; j<=n; j++) |
---|
[501f7c5] | 1287 | { |
---|
[0c824aa] | 1288 | if (v[j] != 0) |
---|
[501f7c5] | 1289 | { |
---|
| 1290 | i++; |
---|
| 1291 | } |
---|
| 1292 | } |
---|
| 1293 | if ( i!=1 ) |
---|
| 1294 | { |
---|
| 1295 | i=0; |
---|
| 1296 | } |
---|
| 1297 | return(i); |
---|
| 1298 | } |
---|
[d5f9d98] | 1299 | //--------------------------------------------------------------- |
---|
[6e9d48a] | 1300 | static proc extgcd_Our(poly p, poly q) |
---|
[0c824aa] | 1301 | { |
---|
| 1302 | ideal J; //for extgcd-computations |
---|
| 1303 | matrix T; //----------"------------ |
---|
| 1304 | list L; |
---|
| 1305 | // the extgcd-command has a bug in versions before 2-0-7 |
---|
| 1306 | if ( system("version")<=2006 ) |
---|
| 1307 | { |
---|
| 1308 | J = p,q; // J = N[k-1,k-1],N[k,k]; //J is of type ideal |
---|
| 1309 | L[1] = liftstd(J,T); //T is of type matrix |
---|
| 1310 | if(J[1]==p) //this is just for the case the SINGULAR swaps the |
---|
| 1311 | // two elements due to ordering |
---|
| 1312 | { |
---|
| 1313 | L[2] = T[1,1]; |
---|
| 1314 | L[3] = T[2,1]; |
---|
| 1315 | } |
---|
| 1316 | else |
---|
| 1317 | { |
---|
| 1318 | L[2] = T[2,1]; |
---|
| 1319 | L[3] = T[1,1]; |
---|
| 1320 | } |
---|
| 1321 | } |
---|
| 1322 | else |
---|
| 1323 | { |
---|
| 1324 | L=extgcd(p,q); |
---|
| 1325 | // L=extgcd(N[k-1,k-1],N[k,k]); |
---|
| 1326 | //one can use this line if extgcd-bug is fixed |
---|
| 1327 | } |
---|
| 1328 | return(L); |
---|
| 1329 | } |
---|
[6e9d48a] | 1330 | static proc normalize_Our(matrix N, matrix Q) |
---|
| 1331 | "USAGE: normalize_Our(N,Q), N, Q are two matrices |
---|
| 1332 | PURPOSE: normalizes N and divides the columns of Q through the leading coefficients of the columns of N |
---|
| 1333 | RETURN: normalized matrix N and altered Q(according to the scheme mentioned in purpose). If number of columns of N and Q do not coincide, N and Q are returned unchanged |
---|
| 1334 | NOTE: number of columns of N and Q must coincide. |
---|
| 1335 | " |
---|
| 1336 | { |
---|
| 1337 | if(ncols(N) != ncols(Q)) |
---|
| 1338 | { |
---|
| 1339 | return (N,Q); |
---|
| 1340 | } |
---|
| 1341 | module M = module(N); |
---|
| 1342 | module S = module(Q); |
---|
| 1343 | int NCols = ncols(N); |
---|
| 1344 | number n; |
---|
| 1345 | for(int i=1;i<=NCols;i++) |
---|
| 1346 | { |
---|
| 1347 | n = leadcoef(M[i]); |
---|
| 1348 | if( n != 0 ) |
---|
| 1349 | { |
---|
| 1350 | M[i]=M[i]/n; |
---|
| 1351 | S[i]=S[i]/n; |
---|
| 1352 | } |
---|
| 1353 | } |
---|
| 1354 | N = matrix(M); |
---|
| 1355 | Q = matrix(S); |
---|
| 1356 | return (N,Q); |
---|
| 1357 | } |
---|
| 1358 | |
---|
[0c824aa] | 1359 | //--------------------------------------------------------------- |
---|
[501f7c5] | 1360 | proc smith( module M ) |
---|
| 1361 | "USAGE: smith(M), M a module or a matrix, |
---|
[d5f9d98] | 1362 | PURPOSE: computes the Smith form of a matrix |
---|
[501f7c5] | 1363 | RETURN: a list of length 4 with the following entries: |
---|
[d5f9d98] | 1364 | @* [1]: The Smith-Form S of M, |
---|
| 1365 | @* [2]: the rank of M, |
---|
| 1366 | @* [3]: a unimodular matrix U, |
---|
| 1367 | @* [4]: a unimodular matrix V, |
---|
[6e9d48a] | 1368 | such that U*M*V=S. An warning is returned when no Smith Form exists. |
---|
| 1369 | NOTE: The Smith form only exists over PIDs (principal ideal domains). Use global ordering for computations! |
---|
[501f7c5] | 1370 | " |
---|
| 1371 | { |
---|
[0e200a] | 1372 | if (nvars(basering)>1) //if more than one variable, return empty list |
---|
[0c824aa] | 1373 | { |
---|
[6e9d48a] | 1374 | string s="The Smith-Form only exists for principal ideal domains"; |
---|
| 1375 | return (s); |
---|
[0c824aa] | 1376 | } |
---|
[501f7c5] | 1377 | matrix N = matrix(M); //Typecasting |
---|
| 1378 | int n = nrows(N); |
---|
| 1379 | int m = ncols(N); |
---|
| 1380 | matrix P = unitmat(n); //left transformation matrix |
---|
| 1381 | matrix Q = unitmat(m); //right transformation matrix |
---|
| 1382 | int k, i, j, deg_temp; |
---|
| 1383 | poly tmp; |
---|
| 1384 | vector v; |
---|
| 1385 | list L; //for extgcd-computation |
---|
| 1386 | intmat f[1][n]; //to save degrees |
---|
| 1387 | matrix lambda[1][n]; //to save leadcoefficients |
---|
| 1388 | intmat g[1][m]; //to save degrees |
---|
| 1389 | matrix mu[1][m]; //to save leadcoefficients |
---|
| 1390 | int ii; //counter |
---|
| 1391 | |
---|
| 1392 | while ((k!=n) && (k!=m) ) |
---|
| 1393 | { |
---|
| 1394 | k++; |
---|
| 1395 | while ((k<=n) && (k<=m)) //outer while-loop for column-operations |
---|
| 1396 | { |
---|
| 1397 | while(k<=m ) //inner while-loop for row-operations |
---|
| 1398 | { |
---|
[0c824aa] | 1399 | if( (n>m) && (k < n) && (k<m)) |
---|
| 1400 | { |
---|
| 1401 | if( simplify((ideal(submat(N,k+1..n,k+1..m))),2)== 0) |
---|
| 1402 | { |
---|
| 1403 | return(N,k-1,P,Q); |
---|
| 1404 | } |
---|
| 1405 | } |
---|
[501f7c5] | 1406 | i,j = smdeg(submat(N,k..n,k..m)); //choose smallest degree in the remaining submatrix |
---|
| 1407 | i=i+(k-1); //indices adjusted to the whole matrix |
---|
| 1408 | j=j+(k-1); |
---|
| 1409 | if(i!=k) //take the element with smallest degree in the first position |
---|
| 1410 | { |
---|
| 1411 | N=permrow(N,i,k); |
---|
| 1412 | P=permrow(P,i,k); |
---|
| 1413 | } |
---|
| 1414 | if(j!=k) |
---|
| 1415 | { |
---|
| 1416 | N=permcol(N,j,k); |
---|
| 1417 | Q=permcol(Q,j,k); |
---|
| 1418 | } |
---|
| 1419 | if(NoNon0Pol(N[k])==1) |
---|
| 1420 | { |
---|
| 1421 | break; |
---|
| 1422 | } |
---|
| 1423 | tmp=leadcoef(N[k,k]); |
---|
| 1424 | deg_temp=ord(N[k,k]); //ord outputs the leading degree of N[k,k] |
---|
| 1425 | for(ii=k+1;ii<=n;ii++) |
---|
| 1426 | { |
---|
| 1427 | lambda[1,ii]=leadcoef(N[ii,k])/tmp; |
---|
| 1428 | f[1,ii]=deg(N[ii,k])-deg_temp; |
---|
| 1429 | } |
---|
| 1430 | for(ii=k+1;ii<=n;ii++) |
---|
| 1431 | { |
---|
[6e9d48a] | 1432 | N = addrow(N,k,-lambda[1,ii]*var(1)^f[1,ii],ii); |
---|
[501f7c5] | 1433 | P = addrow(P,k,-lambda[1,ii]*var(1)^f[1,ii],ii); |
---|
[6e9d48a] | 1434 | N,Q=normalize_Our(N,Q); |
---|
[0c824aa] | 1435 | } |
---|
[501f7c5] | 1436 | } |
---|
| 1437 | if (k>n) |
---|
| 1438 | { |
---|
[0c824aa] | 1439 | break; |
---|
[501f7c5] | 1440 | } |
---|
| 1441 | if(NoNon0Pol(transpose(N)[k])==1) |
---|
| 1442 | { |
---|
| 1443 | break; |
---|
| 1444 | } |
---|
| 1445 | tmp=leadcoef(N[k,k]); |
---|
| 1446 | deg_temp=ord(N[k,k]); //ord outputs the leading degree of N[k][k] |
---|
| 1447 | |
---|
| 1448 | for(ii=k+1;ii<=m;ii++) |
---|
| 1449 | { |
---|
| 1450 | mu[1,ii]=leadcoef(N[k,ii])/tmp; |
---|
| 1451 | g[1,ii]=deg(N[k,ii])-deg_temp; |
---|
| 1452 | } |
---|
| 1453 | for(ii=k+1;ii<=m;ii++) |
---|
| 1454 | { |
---|
| 1455 | N=addcol(N,k,-mu[1,ii]*var(1)^g[1,ii],ii); |
---|
| 1456 | Q=addcol(Q,k,-mu[1,ii]*var(1)^g[1,ii],ii); |
---|
[6e9d48a] | 1457 | N,Q=normalize_Our(N,Q); |
---|
[501f7c5] | 1458 | } |
---|
| 1459 | } |
---|
| 1460 | if( (k!=1) && (k<n) && (k<m) ) |
---|
| 1461 | { |
---|
[6e9d48a] | 1462 | L = extgcd_Our(N[k-1,k-1],N[k,k]); |
---|
[0c824aa] | 1463 | if ( N[k-1,k-1]!=L[1] ) //means that N[k-1,k-1] is not a divisor of N[k,k] |
---|
| 1464 | { |
---|
| 1465 | N=addrow(N,k-1,L[2],k); |
---|
| 1466 | P=addrow(P,k-1,L[2],k); |
---|
[6e9d48a] | 1467 | N,Q=normalize_Our(N,Q); |
---|
[0c824aa] | 1468 | |
---|
[6e9d48a] | 1469 | N=addcol(N,k,-L[3],k-1); |
---|
[0c824aa] | 1470 | Q=addcol(Q,k,-L[3],k-1); |
---|
[6e9d48a] | 1471 | N,Q=normalize_Our(N,Q); |
---|
[0c824aa] | 1472 | k=k-2; |
---|
| 1473 | } |
---|
[501f7c5] | 1474 | } |
---|
| 1475 | } |
---|
| 1476 | if( (k<=n) && (k<=m) ) |
---|
[0c824aa] | 1477 | { |
---|
| 1478 | if( N[k,k]==0) |
---|
[501f7c5] | 1479 | { |
---|
[0c824aa] | 1480 | return(N,k-1,P,Q); |
---|
[501f7c5] | 1481 | } |
---|
[0c824aa] | 1482 | } |
---|
[501f7c5] | 1483 | return(N,k,P,Q); |
---|
| 1484 | } |
---|
| 1485 | example |
---|
| 1486 | { |
---|
| 1487 | "EXAMPLE:";echo = 2; |
---|
| 1488 | option(redSB); |
---|
| 1489 | option(redTail); |
---|
| 1490 | ring r=0,x,dp; |
---|
[0c824aa] | 1491 | // see what happens when the matrix is already in Smith-Form |
---|
[501f7c5] | 1492 | module M = [x,0,0],[0,x2,0],[0,0,x3]; |
---|
| 1493 | print(M); |
---|
| 1494 | list L = smith(M); |
---|
| 1495 | print(L[1]); |
---|
| 1496 | matrix N=matrix(M); |
---|
| 1497 | matrix B=L[3]*N*L[4]; |
---|
| 1498 | print(B); |
---|
| 1499 | //------- and yet another example -------------- |
---|
| 1500 | module M2=[x2,x,3x3-4],[2x2-1,4x,5x2],[2x5,3x,4x]; |
---|
| 1501 | print(M2); |
---|
| 1502 | list P=smith(M2); |
---|
| 1503 | print(P[1]); |
---|
| 1504 | matrix N2=matrix(M2); |
---|
| 1505 | matrix B2=P[3]*N2*P[4]; |
---|
| 1506 | print(B2); |
---|
| 1507 | } |
---|
[d5f9d98] | 1508 | //--------------------------------------------------------------- |
---|
[eccbb1] | 1509 | proc list_tex(L, string name,link l,int nr_loop) |
---|
| 1510 | "USAGE: list_tex(L,name,l), where L is a list, name a string, l a link |
---|
| 1511 | writes the content of list L in a tex-file 'name' |
---|
[d5f9d98] | 1512 | RETURN: nothing |
---|
| 1513 | " |
---|
| 1514 | { |
---|
[eccbb1] | 1515 | if(typeof(L)!="list") //in case L is not a list |
---|
| 1516 | { |
---|
| 1517 | texobj(name,L); |
---|
| 1518 | } |
---|
| 1519 | if(size(L)==0) |
---|
| 1520 | { |
---|
| 1521 | } |
---|
[d5f9d98] | 1522 | else |
---|
[eccbb1] | 1523 | { |
---|
| 1524 | string t; |
---|
| 1525 | for (int i=1;i<=size(L);i++) |
---|
[d5f9d98] | 1526 | { |
---|
[eccbb1] | 1527 | while(1) |
---|
| 1528 | { |
---|
| 1529 | if(typeof(L[i])=="string") //Fehler hier fuer normalen output->nur wenn string in liste dann verbatim |
---|
| 1530 | { |
---|
| 1531 | t=L[i]; |
---|
| 1532 | if(nr_loop==1) |
---|
| 1533 | { |
---|
| 1534 | write(l,"\\begin\{center\}"); |
---|
| 1535 | write(l,"\\begin\{verbatim\}"); |
---|
| 1536 | } |
---|
| 1537 | write(l,t); |
---|
| 1538 | if(nr_loop==0) |
---|
| 1539 | { |
---|
| 1540 | write(l,"\\par"); |
---|
| 1541 | } |
---|
| 1542 | if(nr_loop==1) |
---|
| 1543 | { |
---|
| 1544 | write(l,"\\end\{verbatim\}"); |
---|
| 1545 | write(l,"\\end\{center\}"); |
---|
| 1546 | } |
---|
| 1547 | break; |
---|
| 1548 | } |
---|
| 1549 | if(typeof(L[i])=="module") |
---|
[d5f9d98] | 1550 | { |
---|
[eccbb1] | 1551 | texobj(name,matrix(L[i])); |
---|
| 1552 | break; |
---|
[d5f9d98] | 1553 | } |
---|
[eccbb1] | 1554 | if(typeof(L[i])=="list") |
---|
| 1555 | { |
---|
| 1556 | list_tex(L[i],name,l,1); |
---|
| 1557 | break; |
---|
| 1558 | } |
---|
| 1559 | write(l,"\\begin\{center\}"); |
---|
| 1560 | texobj(name,L[i]); |
---|
| 1561 | write(l,"\\end\{center\}"); |
---|
| 1562 | write(l,"\\par"); |
---|
| 1563 | break; |
---|
| 1564 | } |
---|
[d5f9d98] | 1565 | } |
---|
[eccbb1] | 1566 | } |
---|
[d5f9d98] | 1567 | } |
---|
| 1568 | //--------------------------------------------------------------- |
---|
| 1569 | proc verbatim_tex(string s, link l) |
---|
| 1570 | "USAGE: verbatim_tex(s,l), where s is a string and l a link |
---|
| 1571 | PURPOSE: writes the content of s in verbatim-environment in the file |
---|
| 1572 | specified by link |
---|
| 1573 | RETURN: nothing |
---|
| 1574 | " |
---|
| 1575 | { |
---|
| 1576 | write(l,"\\begin{verbatim}"); |
---|
| 1577 | write(l,s); |
---|
| 1578 | write(l,"\\end{verbatim}"); |
---|
| 1579 | write(l,"\\par"); |
---|
| 1580 | } |
---|
[eccbb1] | 1581 | //--------------------------------------------------------------- |
---|
| 1582 | proc FindTorsion(module R, ideal TAnn) |
---|
| 1583 | "USAGE: FindTorsion(R, I); R an ideal/matrix/module, I an ideal |
---|
| 1584 | PURPOSE: computes the Groebner basis of the submodule of R, annihilated by I |
---|
[6e9d48a] | 1585 | ETURN: module |
---|
[eccbb1] | 1586 | 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 |
---|
| 1587 | the submodule of R |
---|
| 1588 | EXAMPLE: example FindTorsion; shows an example |
---|
| 1589 | " |
---|
| 1590 | { |
---|
| 1591 | // motivation: let R be a module, |
---|
| 1592 | // TAnn is the annihilator of t(R)\subset R |
---|
| 1593 | // compute the generators of t(R) explicitly |
---|
| 1594 | ideal AS = TAnn; |
---|
| 1595 | module S = R; |
---|
| 1596 | if (attrib(S,"isSB")<>1) |
---|
| 1597 | { |
---|
| 1598 | S = std(S); |
---|
| 1599 | } |
---|
| 1600 | if (attrib(AS,"isSB")<>1) |
---|
| 1601 | { |
---|
| 1602 | AS = std(AS); |
---|
| 1603 | } |
---|
| 1604 | int nc = ncols(S); |
---|
[0c824aa] | 1605 | module To = quotient(S,AS); |
---|
| 1606 | To = std(NF(To,S)); |
---|
[eccbb1] | 1607 | return(To); |
---|
| 1608 | } |
---|
| 1609 | example |
---|
| 1610 | { |
---|
| 1611 | "EXAMPLE:";echo = 2; |
---|
| 1612 | // Flexible Rod |
---|
| 1613 | ring A = 0,(D1, D2), (c,dp); |
---|
| 1614 | module R= [D1, -D1*D2, -1], [2*D1*D2, -D1-D1*D2^2, 0]; |
---|
| 1615 | module RR = transpose(R); |
---|
| 1616 | list L = control(RR); |
---|
| 1617 | // here, we have the annihilator: |
---|
| 1618 | ideal LAnn = D1; // = L[10] |
---|
| 1619 | module Tr = FindTorsion(RR,LAnn); |
---|
| 1620 | print(RR); // the module itself |
---|
| 1621 | print(Tr); // generators of the torsion submodule |
---|
| 1622 | } |
---|