source: git/Singular/LIB/control.lib @ 359b79

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