source: git/Singular/LIB/control.lib @ 0ae788

fieker-DuValspielwiese
Last change on this file since 0ae788 was 0ae788, checked in by Viktor Levandovskyy <levandov@…>, 19 years ago
*levandov: genericity and canonize procs added git-svn-id: file:///usr/local/Singular/svn/trunk@7620 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.3 KB
Line 
1version="$Id: control.lib,v 1.16 2004-12-09 13:56:29 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
8SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis' (Project of Dr. E. Zerz
9and V. Levandovskyy), Uni Kaiserslautern
10
11NOTE: This library provides algebraic analysis tools for System and Control Theory
12
13PROCEDURES:
14control(module R);            analysis of controllability-related properties of R,
15autonom(module R);            analysis of autonomy-related properties of R (using Ext modules),
16autonom2(module R);           analysis of autonomy-related properties of R (using dimension),
17LeftKernel(module R);         a left kernel of R,
18RightKernel(module R);        a right kernel of R,
19LeftInverse(module R)         a left inverse of matrix (module).
20
21AUXILIARY PROCEDURES:
22ncdetection(ring r);          computes an ideal, presenting an involution map on non-comm algebra r;
23involution(m, map theta);    applies the involution, presented by theta to m of type poly, vector, ideal, module;
24declare(string NameOfRing, Variables[,string  Parameters, Ordering]);     defines the ring, optional parametes are strings of parameters and ordering,
25view();                      Well-formatted output of lists, modules and matrixes
26
27NOTE (EXAMPLES): In order to use examples below, execute the commands
28@* def A = exAntenna(); setring A;
29@* Thus A will become a basering from the example with the predefined module R (transposed), corresponding to the system.
30After that you can just type in
31@* control(R); //respectively autonom(R);
32and check the result.
33
34
35EXAMPLES (AUTONOMY):
36
37exCauchy();                  example of 1-dimensional Cauchy equation,
38exCauchy2();                 example of 2-dimensional Cauchy equation,
39exZerz();                    example from the lecture of Eva Zerz,
40
41EXAMPLES (CONTROLLABILITY):
42
43ex1();                       example of noncontrollable system,   
44ex2();                       example of controllable system ,
45exAntenna();                 Antenna,
46exEinstein();                Einstein equations,
47exFlexibleRod();             Flexible Rod,
48exTwoPendula();              Two Pendula,
49exWindTunnel();              Wind Tunnel.
50";
51
52// NOTE: static things should not be shown for end-user
53// static Ext_Our(...)                  Copy of Ext_R from 'homolog.lib' in commutative case;
54// static is_zero_Our(module R)         Copy of is_zero from 'OBpoly.lib';
55//  static space(int n)           Procedure used inside the procedure 'Print' to have a better formatted output
56// static control_output();      Generating the output for the procedure 'control'
57// static autonom_output();      Generating the output for the procedure 'autonom' and 'autonom2'
58
59LIB "homolog.lib";
60LIB "poly.lib";
61LIB "primdec.lib";
62LIB "ncalg.lib";
63//---------------------------------------------------------------
64proc declare(string NameOfRing, string Variables, list #)
65"USAGE: declare(NameOfRing, Variables,[Parameters, Ordering]);
66          NameOfRing:  string with name of ring,
67          Variables:   string with names of variables separated by commas(e.g. "x,y,z"),
68          [Parameters, Ordering]: optional, strings:
69            Parameters: string of parameters in the ring separated by commas(e.g. "a,b,c"),
70            Ordering:   string with name of ordering(by default the ordering "dp,C" is used)
71RETURN:  no return value
72EXAMPLE:  example declare; shows an example
73"
74{
75  if(size(#)==0)
76  {
77    execute("ring "+NameOfRing+"=0,("+Variables+"),dp;");
78  }
79  else
80  {
81    if(size(#)==1)
82    {
83      execute( "ring " + NameOfRing + "=(0," + #[1] + "),(" + Variables + "),dp;" );
84    }
85    else
86    {
87      if( (size(#[1])!=0)&&(#[1]!=" ") )
88      {
89        execute( "ring " + NameOfRing + "=(0," + #[1] + "),(" + Variables + "),("+#[2]+");" );
90      }
91      else
92      {
93        execute( "ring " + NameOfRing + "=0,("+Variables+"),("+#[2]+");" );
94      };
95    };
96  };
97  keepring(basering);
98}
99example
100{"EXAMPLE:";echo = 2;
101  string v="x,y,z";
102  string p="q,p";
103  string Ord ="c,lp";
104 
105  declare("Ring_1",v);
106  print(nameof(basering));
107  print(basering);
108 
109  declare("Ring_2",v,p);
110  print(basering);
111  print(nameof(basering));
112 
113  declare("Ring_3",v,p,Ord);
114  print(basering);
115  print(nameof(basering));
116 
117  declare("Ring_4",v,"",Ord);
118  print(basering);
119  print(nameof(basering));
120 
121  declare("Ring_5",v," ",Ord);
122  print(basering);
123  print(nameof(basering));
124};
125//
126//maybe reasonable to add this in declare
127//
128//  print("Please enter your representation matrix in the following form:
129//  module R=[1st row],[2nd row],...");
130//  print("Type the command: R=transpose(R)");
131//  print(" To compute controllability please enter: control(R)");
132//  print(" To compute autonomy please enter: autonom(R)");
133//
134//
135//
136//-------------------------------------------------------------------------
137static proc space(int n)
138"USAGE:spase(n);
139         n: integer, number of needed spaces
140RETURN:  string consisting of n spaces
141NOTE:  the procedure is used in the procedure 'view' to have a better formatted output
142"
143{
144  int i;
145  string s="";
146  for(i=1;i<=n;i++)
147  {
148    s=s+" ";
149  };
150  return(s);
151};
152//-----------------------------------------------------------------------------
153proc view(M)
154"USAGE:  view(M);
155           M:  any type
156RETURN:  no return value
157PURPOSE:  procedure for ( well-) formatted output of modules, matrices, lists of modules, matrices;
158          shows everything even if entries are long
159NOTE:  in case of other types( not 'module', 'matrix', 'list') works just as standard 'print' procedure 
160EXAMPLE:  example view; shows an example
161"
162{
163  if ( (typeof(M)=="module")||(typeof(M)=="matrix") )
164  {
165  int @R=nrows(M);
166  int @C=ncols(M);
167  int i;
168  int j;
169  list MaxLength=list();
170  int Size=0;
171  int max;
172  string s;
173 
174  for(i=1;i<=@C;i++)
175  {
176    max=0;
177   
178    for(j=1;j<=@R;j++)
179    {
180      Size=size( string( M[j,i] ) );
181      if( Size>max )
182      {
183        max=Size;
184      };
185    };
186    MaxLength[i] = max;
187  };
188 
189  for(i=1;i<=@R;i++)
190  {
191    s="";
192    for(j=1;j<@C;j++)
193    {
194      s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) ) +",";
195    };
196   
197    s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) );
198
199    if (i!=@R)
200    {
201      s=s+",";
202    };
203    print(s);
204  };
205
206  return();   
207  };
208 
209  if(typeof(M)=="list")
210  {
211    int sz=size(M);
212    int i;
213    for(i=1;i<=sz;i++)
214    {
215      view(M[i]);
216      print("");
217    };
218
219    return();
220  };
221  print(M);
222  return();
223}
224example
225{"EXAMPLE:";echo = 2;
226  ring r;
227  matrix M[1][3] = x,y,z;
228  print(M);
229  view(M);
230};
231//--------------------------------------------------------------------------
232proc RightKernel(matrix M)
233"USAGE:  RightKernel(M);
234           M:  matrix
235RETURN:  right kernel of matrix M, i.e., the module of all elements v such that Mv=0
236NOTE:  in commutative case it is a left  module, in noncommutative (will be implemented later) it is a right module   
237EXAMPLE:  example RightKernel; shows an example
238"
239{
240  return(modulo(M,std(0)));
241}
242example
243{"EXAMPLE:";echo = 2;
244  ring r;
245  matrix M[1][3] = x,y,z;
246  print(M);
247  print( RightKernel(M) );
248};
249//-------------------------------------------------------------------------
250proc LeftKernel(matrix M)
251"USAGE:  LeftKernel(M);
252           M:  matrix
253RETURN:  left kernel of matrix M, i.e., the matrix whose rows are generators of left module
254         (elements of this module are to be rows) of all elements v such that vM=0   
255EXAMPLE:  example LeftKernel; shows an example
256"
257{
258  return( transpose( modulo( transpose(M),std(0) ) ) );
259}
260example
261{"EXAMPLE:";echo = 2;
262  ring r;
263  matrix M[3][1] = x,y,z;
264  print(M);
265  print( LeftKernel(M) );
266};
267//------------------------------------------------------------------------
268proc LeftInverse(matrix M)
269"USAGE:  LeftInverse(M);
270           M:  matrix
271RETURN:  left inverse of M if exists, i.e., matrix L such that LM == id;
272EXAMPLE:  example LeftInverse; shows an example
273"
274{
275  int NCols=ncols(M);
276  M=transpose(M);
277  //  matrix I[NCols][NCols];
278  //  I=I+1;
279  //  module Id=I;
280  module Id = freemodule(NCols);
281  return( transpose( lift( module(M),Id ) )  );
282}
283example
284{"EXAMPLE:";echo =2;
285  ring r;
286  matrix M[2][1] = 1,x2z;
287  print(M);
288  print( LeftInverse(M) ); 
289};
290//-----------------------------------------------------------------------
291static proc Ext_Our(int i, module R,list #)
292// just a copy of 'Ext_R' from "homolog.lib" in commutative case
293{
294  if ( size(#)==0 )
295  {
296    return( Ext_R(i,R) );
297  }
298  else
299  {
300    return( Ext_R(i,R,#[1]) );
301  };
302}
303//------------------------------------------------------------------------
304static proc is_zero_Our(module R)
305//just a copy of 'is_zero' from "poly.lib"
306{
307  return( is_zero(R) ) ;
308};
309//------------------------------------------------------------------------
310static proc control_output(int i, int NVars, module R, module Ext_1)
311"USAGE:  proc control_output(i, NVars, R, Ext_1)
312           i:  integer, number of first nonzero Ext or
313               just number of variables in a base ring + 1 in case that all the Exts are zero
314           NVars:  integer, number of variables in a base ring 
315           R:  module R (cokernel representation)
316           Ext_1:  module, the first Ext(its cokernel representation)     
317RETURN:  list with all the contollability properties of the system which is to be returned in 'control' procedure
318NOTE:  this procedure is used in 'control' procedure
319"
320{
321  int d=dim( std( Ann( transpose(R) ) ) ) ;; //this is the dimension of the system
322  string DofS= "dimension of the system:";
323  string Fn= "number of first nonzero Ext:";
324  if(i==1)
325  {
326    module RK=RightKernel(R);
327    return(
328            list ( Fn,
329                   i,
330                  "not controllable , image representation for controllable part:",
331                   RK,     
332                  "kernel representation for controllable part:",
333                   LeftKernel( RK ),
334                  "obstruction to controllability",
335                   Ext_1,
336                  "annihilator of torsion module(of obstruction to controllability)",
337                   Ann(Ext_1),
338                   DofS,
339                   d
340                 )
341          );
342  };
343 
344  if(i>NVars)
345  { module RK =RightKernel(R);
346    return( list(  Fn,
347                   -1,
348                  "strongly controllable, image representation:",
349                   RK,
350                  "left inverse to image representation:",
351                   LeftInverse(RK),
352                   DofS,
353                   d)
354          );
355  };
356 
357  //
358  //now i<=NVars
359  //
360       
361  if( (i==2) )
362  {
363    return( list( Fn,
364                  i,
365                 "controllable, not reflexive, image representation:",
366                  RightKernel(R),
367                  DofS,
368                  d )
369          );
370  };
371   
372  if( (i>=3) )
373  {
374    return( list ( Fn,
375                   i,
376                  "reflexive, not strongly controllable, image representation:",
377                   RightKernel(R),
378                   DofS,
379                   d)
380          );
381  };
382             
383 
384}; 
385//-------------------------------------------------------------------------
386
387proc control(module R)
388"USAGE:  control(R);
389           R:  module (R is a matrix of the system of equations which is to be investigated)
390RETURN:  list of all the properties concerning controllability of the system(behavior) represented by the  matrix R
391EXAMPLE:  example control; shows an example
392"
393{
394  int i;
395  int NVars=nvars(basering);
396  int ExtIsZero;
397 
398       
399  module Ext_1 = std(Ext_Our(1,R));
400   
401  ExtIsZero=is_zero_Our(Ext_1);
402  i=1;
403  while( (ExtIsZero) && (i<=NVars) )
404  {
405    i++;
406    ExtIsZero = is_zero_Our( Ext_Our(i,R) );
407  };
408 
409  return( control_output( i, NVars, R, Ext_1 ) );
410}
411example
412{"EXAMPLE:";echo = 2;
413  //Wind Tunnel
414  ring A = (0,a, omega, zeta, k),(D1, delta),dp;
415  module R;
416  R = [D1+a, -k*a*delta, 0, 0],
417      [0, D1, -1, 0],
418      [0, omega^2, D1+2*zeta*omega, -omega^2];
419  R=transpose(R);
420  view(R);
421  view(control(R));
422
423};
424//------------------------------------------------------------------------
425static proc autonom_output( int i, int NVars )
426"USAGE:  proc autonom_output(i, NVars)
427           i:  integer, number of first nonzero Ext or
428               just number of variables in a base ring + 1 in case that all the Exts are zero
429           NVars:  integer, number of variables in a base ring 
430RETURN:  list with all the autonomy properties of the system which is to be returned in 'autonom' procedure
431NOTE:  this procedure is used in 'autonom' procedure
432"
433{
434  int d=NVars-i;//that is the dimension of the system
435  string DofS="dimension of the system:";
436  string Fn = "number of first nonzero Ext:";
437  if(i==0)
438  {
439    return( list(  Fn,
440                   i,
441                  "not autonomous",
442                   DofS,
443                   d )
444          );
445  };
446 
447  if( i>NVars )
448  { 
449    return( list( Fn,
450                  -1,
451                  "trivial",
452                  DofS,
453                  d )
454          );
455  };
456 
457  //
458  //now i<=NVars
459  //
460   
461     
462  if( i==1 )
463  //in case that NVars==1 there is no sence to consider the notion
464  //of strongly autonomous behavior, because it does not imply
465  //that system is overdetermined in this case
466  {
467    return( list ( Fn,
468                   i,
469                  "autonomous, not overdetermined",
470                   DofS,
471                   d )
472          );
473  };
474   
475  if( i==NVars )
476  { 
477    return( list(  Fn,
478                   i,
479                  "strongly autonomous,in particular overdetermined",
480                   DofS,
481                   d)
482          );
483  };
484 
485  if( i<NVars )
486  {
487    return( list ( Fn,
488                   i,
489                  "overdetermined, not strongly autonomous",
490                   DofS,
491                   d)
492          );
493  };
494   
495}; 
496//--------------------------------------------------------------------------
497proc autonom2(module R)
498"USAGE:  autonom2(R);
499           R:  module (R is a matrix of the system of equations which is to be investigated)
500RETURN:  list of all the properties concerning autonomy of the system(behavior) represented by the  matrix R
501NOTE:  this procedure is an analogue to 'autonom' using dimension calculations
502EXAMPLE:  example autonom2; shows an example
503"
504{
505  int d;
506  int NVars = nvars(basering);
507  R=transpose(R);
508  d=dim( std( Ann(R) ) );
509  return( autonom_output(NVars-d,NVars) );
510}
511example
512{"EXAMPLE:"; echo = 2;
513  //Cauchy
514  ring r=0,(s1,s2,s3,s4),dp;
515  module R= [s1,-s2],
516            [s2, s1],
517            [s3,-s4],
518            [s4, s3];       
519  R=transpose(R);
520  view( R );
521  view( autonom2(R) );
522}; 
523//---------------------------------------------------------------------------
524
525proc autonom(module R)
526"USAGE:  autonom(R);
527           R:  module (R is a matrix of the system of equations which is to be investigated)
528RETURN:  list of all the properties concerning autonomy of the system(behavior) represented by the  matrix R
529EXAMPLE:  example autonom; shows an example
530"
531{
532  int NVars=nvars(basering);
533  int ExtIsZero;   
534  R=transpose(R);
535  ExtIsZero=is_zero_Our(Ext_Our(0,R));
536  int i=0;
537  while( (ExtIsZero)&&(i<=NVars) )
538  {
539    i++;
540    ExtIsZero = is_zero_Our(Ext_Our(i,R));
541  };
542
543  return(autonom_output(i,NVars));     
544}
545example
546{"EXAMPLE:"; echo = 2;
547  //Cauchy
548  ring r=0,(s1,s2,s3,s4),dp;
549  module R= [s1,-s2],
550            [s2, s1],
551            [s3,-s4],
552            [s4, s3];       
553  R=transpose(R);
554  view( R );
555  view( autonom(R) );
556}; 
557
558//--------------------------------------------------------------------------
559//
560//Some example rings with defined systems
561//----------------------------------------------------------------------------
562//autonomy:
563//
564//----------------------------------------------------------------------------
565proc exCauchy()
566{
567  ring @r=0,(s1,s2),dp;
568  module R= [s1,-s2],
569            [s2, s1];
570  R=transpose(R);
571  export R;
572  return(@r);       
573};
574//----------------------------------------------------------------------------
575proc exCauchy2()
576{
577  ring @r=0,(s1,s2,s3,s4),dp;
578  module R= [s1,-s2],
579            [s2, s1],
580            [s3,-s4],
581            [s4, s3];       
582  R=transpose(R);
583  export R;
584  return(@r);       
585};
586//----------------------------------------------------------------------------
587proc exZerz()
588{
589  ring @r=0,(d1,d2),dp;
590  module R=[d1^2-d2],
591           [d2^2-1];
592  R=transpose(R);
593  export R;
594  return(@r);       
595}; 
596//----------------------------------------------------------------------------
597//control
598//
599proc ex1()
600{
601  ring @r=0,(s1,s2,s3),dp;
602  module R=[0,-s3,s2],
603           [s3,0,-s1];
604  R=transpose(R);         
605  export R;
606  return(@r);
607};
608//----------------------------------------------------------------------------
609proc ex2()
610{
611  ring @r=0,(s1,s2,s3),dp;
612  module R=[0,-s3,s2],
613            [s3,0,-s1],
614            [-s2,s1,0];
615  R=transpose(R);         
616  export R;
617  return(@r);
618};
619//----------------------------------------------------------------------------
620proc exAntenna()
621{
622  ring @r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), dp;
623  module R;
624  R = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
625      [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
626      [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
627      [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
628      [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
629      [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];
630
631  R=transpose(R);
632  export R;
633  return(@r); 
634};
635
636//----------------------------------------------------------------------------
637
638proc exEinstein()
639{
640  ring @r = 0,(D(1..4)),dp;
641  module R =
642  [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)],
643  [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],
644  [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],
645  [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)],
646  [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)],
647  [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],
648  [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)],
649  [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)],
650  [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)],
651  [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];
652
653  R=transpose(R);
654  export R;
655  return(@r); 
656};
657
658
659//---------------------------------------------------------------------------------------------
660
661proc exFlexibleRod()
662{
663  ring @r = 0,(D1, delta), dp;
664  module R;
665  R = [D1, -D1*delta, -1], [2*D1*delta, -D1-D1*delta^2, 0];
666 
667  R=transpose(R);
668  export R;
669  return(@r); 
670};
671
672//---------------------------------------------------------------------------------------------
673proc exTwoPendula()
674{ 
675  ring @r=(0,m1,m2,M,g,L1,L2),Dt,dp;
676  module R = [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
677             [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
678             [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
679
680  R=transpose(R);
681  export R;
682  return(@r); 
683};
684//---------------------------------------------------------------------------------------------
685proc exWindTunnel()
686{
687  ring @r = (0,a, omega, zeta, k),(D1, delta),dp;
688  module R = [D1+a, -k*a*delta, 0, 0],
689             [0, D1, -1, 0],
690             [0, omega^2, D1+2*zeta*omega, -omega^2];
691
692  R=transpose(R);
693  export R;
694  return(@r); 
695};
696//--------------------------------------------------------------------------------------------
697
698
699//---------------------------------------------------------------------------
700//---------------------------------------------------------------------------
701//---------------------------------------------------------------------------
702//---------------------------------------------------------------------------
703
704static proc invo_poly(poly m, map theta)
705//applies the involution map theta to m, where m=polynomial
706{
707  int i,j;
708  intvec v;
709  poly p,z;
710  poly n = 0;
711  i = 1;
712  while(m[i]!=0)
713  {
714    v = leadexp(m[i]);
715    z =1;
716    for(j=nvars(basering); j>=1; j--)
717    {
718      if (v[j]!=0)
719      {
720        p = var(j);
721        p = theta(p);
722        z = z*(p^v[j]);
723      }
724    }   
725    n = n + (leadcoef(m[i])*z);
726    i++; 
727  }
728  return(n);
729}
730
731proc involution(m, map theta)
732//applies the involution map theta to m, where m=vector, polynomial,
733//module,ideal
734{
735  int i,j;
736  intvec v;
737  poly p,z;
738  if (typeof(m)=="poly")
739  {
740    return (invo_poly(m,theta)); 
741  }
742  if ( typeof(m)=="ideal" )
743  {
744    ideal n;
745    for (i=1; i<=size(m); i++)
746    {
747      n[i] = invo_poly(m[i],theta);
748    }
749    return(n);
750  }
751  if (typeof(m)=="vector")
752  {
753    for(i=1;i<=size(m);i++)
754    {
755      m[i] = invo_poly(m[i],theta);
756    }
757    return (m); 
758  }
759 
760  if ( (typeof(m)=="matrix") || (typeof(m)=="module"))
761  { 
762//    m=transpose(m);
763    matrix n = matrix(m);
764    int @R=nrows(n);
765    int @C=ncols(n);
766    for(i=1; i<=@R; i++)
767    {
768      for(j=1; j<=@C; j++)
769      {
770        n[i,j] = invo_poly( m[i,j], theta);
771      }
772    }
773   }
774  if (typeof(m)=="module")
775  {
776    return (module(n));
777  }
778  return(n);
779}
780example
781{
782  "EXAMPLE:";echo = 2;
783  ring r = 0,(x,d),dp;
784  ncalgebra(1,1); // Weyl-Algebra
785  map F = r,x,-d;
786  poly f =  x*d^2+d;
787  poly If = involution(f,F);
788  f-If;
789  poly g = x^2*d+2*x*d+3*x+7*d;
790  poly tg = -d*x^2-2*d*x+3*x-7*d;
791  poly Ig = involution(g,F);
792  tg-Ig;
793  ideal I = f,g;
794  ideal II = involution(I,F);
795  II;
796  I - involution(II,F);
797  module M  = [f,g,0],[g,0,x^2*d];
798  module IM = involution(M,F);
799  print(IM);
800  print(M - involution(IM,F)); 
801}
802
803proc ncdetection( r)
804// in this procedure an involution map is generated from the NCRelations
805// that will be used in the function involution
806// in dieser proc. wird eine matrix erzeugt, die in der i-ten zeile die indices
807// der differential-,shift- oder advance-operatoren enthaelt mit denen die i-te
808// variable nicht kommutiert.
809{
810  int i,j,k,LExp;
811  int NVars=nvars(r);
812  matrix rel = NCRelations(r)[2];
813  intmat M[NVars][3];
814  int NRows = nrows(rel);
815  intvec v,w;
816  poly d,d_lead;
817  ideal I;
818  map theta;
819 
820  for( j=NRows;j>=2;j-- )
821  {
822   if( rel[j] == w )       //the whole column is zero
823    {
824      j--;
825      continue;
826    }
827   
828    for( i=1;i<j;i++ )         
829    {
830      if( rel[i,j]==1 )        //relation of type var(j)*var(i) = var(i)*var(j) +1
831      {
832         M[i,1]=j;
833      }
834      if( rel[i,j] == -1 )    //relation of type var(i)*var(j) = var(j)*var(i) -1
835      {
836        M[j,1]=i;
837      }
838      d = rel[i,j];
839      d_lead = lead(d);
840      v=leadexp(d_lead); //in the next lines we check wether we have a  relation of differential or shift type
841      LExp=0;
842      for( k=1;k<=NVars;k++)
843      {
844        LExp = LExp + v[k];
845      }
846      //      if( (d-d_lead != 0) || (LExp > 1) )
847      if( ( d-d_lead != 0) || (LExp > 1) || ((LExp==0)&& !((d_lead==1) ||
848                                                           (d_lead==-1))) )
849      {
850        return( "wrong input" );
851      }
852
853      if( v[j] == 1)                   //relation of type var(j)*var(i) = var(i)*var(j) -lambda*var(j)
854      {
855        if (leadcoef(d) < 0)
856        {
857          M[i,2] = j;
858        }
859        else
860        {
861          M[i,3] = j;
862        }
863      }
864      if( v[i]==1 )                    //relation of type var(j)*var(i) = var(i)*var(j) -lambda*var(i)
865      {
866        if (leadcoef(d) > 0)
867        {
868          M[j,2] = i;
869        }
870        else
871        {
872          M[j,3] = i;
873        }
874      }
875    }
876  }
877  //ab hier wird die map ausgerechnet
878  for(i=1;i<=NVars;i++)
879  {
880    I=I+var(i);
881  }
882
883  for(i=1;i<=NVars;i++)
884  {
885    if( M[i,1..3]==(0,0,0) )
886    {
887      i++;
888      continue;
889    }
890    if( M[i,1]!=0 )
891    {
892      if( (M[i,2]!=0) && (M[i,3]!=0) )
893      {
894        I[M[i,1]] = -var(M[i,1]);
895        I[M[i,2]] = var(M[i,3]);
896        I[M[i,3]] = var(M[i,2]);
897      }
898      if( (M[i,2]==0) && (M[i,3]==0) )
899      {
900        I[M[i,1]] = -var(M[i,1]);
901      }                 
902      if( ( (M[i,2]!=0) && (M[i,3]==0) )|| ( (M[i,2]!=0) && (M[i,3]==0) )
903)
904      {
905        I[i] = -var(i);
906      }
907    }
908    else
909    {
910      if( (M[i,2]!=0) && (M[i,3]!=0) )
911      {
912        I[i] = -var(i);
913        I[M[i,2]] = var(M[i,3]);
914        I[M[i,3]] = var(M[i,2]);
915      }
916      else
917      {
918        I[i] = -var(i);
919      }
920    }
921  }
922  return(I);
923
924}
925example
926{
927  "EXAMPLE:"; echo = 2;
928  ring r=0,(x,y,z,D(1..3)),dp;
929  matrix D[6][6];
930  D[1,4]=1;
931  D[2,5]=1;
932  D[3,6]=1;
933  ncalgebra(1,D);
934  ncdetection(r);
935  kill r;
936  //----------------------------------------
937  ring r=0,(x,S),dp;
938  ncalgebra(1,-S);
939  ncdetection(r);
940  kill r;
941  //----------------------------------------
942  ring r=0,(x,D(1),S),dp;
943  matrix D[3][3];
944  D[1,2]=1;
945  D[1,3]=-S;
946  ncalgebra(1,D);
947  ncdetection(r);
948}
949
950proc genericity(matrix M)
951"USAGE:  genericity(M), M is a matrix
952RETURN:  list of strings with
953NOTE:  effective with the liftstd procedure
954"
955{
956  // M is a matrix over a ring with params and vars;
957  ideal I = ideal(M); // a list of entries
958  I = simplify(I,2); // throw 0's away
959  // decompose every coeff
960  int i; int  cl=1;
961  int s = size(I);
962  list NM;
963  poly p;
964  for (i=1; i<=s; i++)
965  {
966    p = I[i];
967    while( p != 0)
968    {
969      NM[cl] = leadcoef(p);
970      cl++;
971      p = p - lead(p);
972    };
973  };
974  string newvars = parstr(basering);
975  def save = basering;
976  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
977  execute(NewRing);
978  // get params as variables
979  // creat a list of non-monomials
980  ideal L;
981  ideal F;
982  list NM = imap(save,NM);
983  poly p,q;
984  cl = 1;
985  int j, cf;
986  for(i=1; i<=size(NM);i++)
987  {
988    p = NM[i] - lead(NM[i]);
989    if (p!=0)
990    {
991      //      L[cl] = p;
992      F  = factorize(NM[i],1); //non-constant factors only
993      cf = 1;
994      // factorize every polynomial
995      // throw away every non-monomial
996      for (j=1; j<=size(F);j++)
997      {
998        q = F[j]-lead(F[j]);
999        if (q!=0)
1000        {
1001          L[cl] = F[j];
1002          cl++;
1003        }
1004      }
1005    }
1006  }
1007  // return the result [in string=format]
1008  L = simplify(L,2+4); // skip zeroes and double entries
1009  list SL;
1010  for (j=1; j<=size(L);j++)
1011  {
1012    SL[j] = string(L[j]);
1013  }
1014  setring save;
1015  return(SL);
1016}
1017example
1018{  // TwoPendula
1019  "EXAMPLE:"; echo = 2;
1020  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
1021  module RR =
1022    [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
1023    [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
1024    [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
1025  module R = transpose(RR);
1026  matrix T;
1027  module SR = liftstd(R,T);
1028  genericity(T);
1029}
1030
1031proc canonize(list L)
1032"USAGE:  canonize(L), L is a list
1033ASSUME: L is the output of control/autonomy procs
1034RETURN:  canonized list
1035"
1036{
1037  list M = L;
1038  option(redSB);
1039  option(redTail);
1040  int s = size(L);
1041  int i;
1042  for (i=2; i<=s; i=i+2)
1043  {
1044    if (typeof(M[i])=="module")
1045    {
1046      M[i] = std(M[i]);
1047      M[i] = prune(M[i]); // mimimal embedding
1048      M[i] = std(M[i]);
1049    }
1050  }
1051  return(M);
1052}
1053example
1054{  // TwoPendula with L1=L2=L
1055  "EXAMPLE:"; echo = 2;
1056  ring r=(0,m1,m2,M,g,L),Dt,dp;
1057  module RR =
1058    [m1*L*Dt^2, m2*L*Dt^2, -1, (M+m1+m2)*Dt^2],
1059    [m1*L^2*Dt^2-m1*L*g, 0, 0, m1*L*Dt^2],
1060    [0, m2*L^2*Dt^2-m2*L*g, 0, m2*L*Dt^2];
1061  module R = transpose(RR);
1062  list C = control(R);
1063  list CC = canonize(C);
1064  view(CC);
1065}
Note: See TracBrowser for help on using the repository browser.