# source:git/Singular/LIB/control.lib@a09c2a0

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