source: git/Singular/LIB/control.lib @ 058269

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