source: git/Singular/LIB/control.lib @ 6a02de

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