source: git/Singular/LIB/control.lib @ 82d3aec

spielwiese
Last change on this file since 82d3aec was 82d3aec, checked in by Viktor Levandovskyy <levandov@…>, 20 years ago
*yena: new example rings(systems) for autonomy git-svn-id: file:///usr/local/Singular/svn/trunk@7351 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 24.0 KB
Line 
1version="$Id: control.lib,v 1.12 2004-08-10 14:47:02 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 typ poly, vector, ideal, module;
24declare(string NameOfRing, string Variables[,string  Parameters, string Ordering]);     defines the ring, optional parametes are a string of parameters and a sting of 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  //-------------------------------
389  ring r=0,(s1,s2,s3),dp;
390  module R=[0,-s3,s2],
391           [s3,0,-s1];
392  R=transpose(R);         
393  view( R ); 
394  view( control(R) );
395
396  module Rc=[0,-s3,s2],
397            [s3,0,-s1],
398            [-s2,s1,0];
399  Rc=transpose(Rc);         
400  view( Rc );
401  view( control(Rc) );
402  kill r;
403  //----------------------------------
404  //reflector antenna
405  ring A = (0, K1, K2, Te, Kp, Kc),(Dt, delta), dp;
406  module R;
407  R = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
408       [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
409       [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
410       [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
411       [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
412       [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];
413  R=transpose(R);
414  view(R);
415  view(control(R));
416  kill A;
417  //----------------------------------
418  //Flexible Rod
419  ring A = 0,(D1, delta), dp;
420  module R;
421  R = [D1, -D1*delta, -1],
422      [2*D1*delta, -D1-D1*delta^2, 0];
423  R=transpose(R);
424  view(R);
425  view(control(R));
426  kill A;
427  //-------------------------------------
428  //TwoPendula
429  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
430  module R;
431  R =  [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
432       [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
433       [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
434  R=transpose(R);
435  view(R);
436  view(control(R));
437  kill r;
438  //----------------------------------------
439  //Wind Tunnel
440  ring A = (0,a, omega, zeta, k),(D1, delta),dp;
441  module R;
442  R = [D1+a, -k*a*delta, 0, 0],
443      [0, D1, -1, 0],
444      [0, omega^2, D1+2*zeta*omega, -omega^2];
445  R=transpose(R);
446  view(R);
447  view(control(R));
448  kill A;
449  //-------------------------------------------
450
451};
452//------------------------------------------------------------------------
453static proc autonom_output( int i, int NVars )
454"USAGE:  proc autonom_output(i, NVars)
455           i:  integer, number of first nonzero Ext or
456               just number of variables in a base ring + 1 in case that all the Exts are zero
457           NVars:  integer, number of variables in a base ring 
458RETURN:  list with all the autonomy properties of the system which is to be returned in 'autonom' procedure
459NOTE:  this procedure is used in 'autonom' procedure
460"
461{
462  int d=NVars-i;//that is the dimension of the system
463  string DofS="dimension of the system:";
464  string Fn = "number of first nonzero Ext:";
465  if(i==0)
466  {
467    return( list(  Fn,
468                   i,
469                  "not autonomous",
470                   DofS,
471                   d )
472          );
473  };
474 
475  if( i>NVars )
476  { 
477    return( list( Fn,
478                  -1,
479                  "trivial",
480                  DofS,
481                  d )
482          );
483  };
484 
485  //
486  //now i<=NVars
487  //
488   
489     
490  if( i==1 )
491  //in case that NVars==1 there is no sence to consider the notion
492  //of strongly autonomous behavior, because it does not imply
493  //that system is overdetermined in this case
494  {
495    return( list ( Fn,
496                   i,
497                  "autonomous, not overdetermined",
498                   DofS,
499                   d )
500          );
501  };
502   
503  if( i==NVars )
504  { 
505    return( list(  Fn,
506                   i,
507                  "strongly autonomous,in particular overdetermined",
508                   DofS,
509                   d)
510          );
511  };
512 
513  if( i<NVars )
514  {
515    return( list ( Fn,
516                   i,
517                  "overdetermined, not strongly autonomous",
518                   DofS,
519                   d)
520          );
521  };
522   
523}; 
524//--------------------------------------------------------------------------
525proc autonom2(module R)
526"USAGE:  autonom2(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
529NOTE:  this procedure is an analogue to 'autonom' using dimension calculations
530EXAMPLE:  example autonom2; shows an example
531"
532{
533  int d;
534  int NVars = nvars(basering);
535  R=transpose(R);
536  d=dim( std( Ann(R) ) );
537  return( autonom_output(NVars-d,NVars) );
538}
539example
540{"EXAMPLE:"; echo = 2;
541  ring r=0,(s1,s2,s3,s4),dp;
542  module R= [s1,-s2],
543            [s2, s1],
544            [s3,-s4],
545            [s4, s3];       
546  R=transpose(R);
547  view( R );
548  view( autonom2(R) );
549}; 
550//---------------------------------------------------------------------------
551
552proc autonom(module R)
553"USAGE:  autonom(R);
554           R:  module (R is a matrix of the system of equations which is to be investigated)
555RETURN:  list of all the properties concerning autonomy of the system(behavior) represented by the  matrix R
556EXAMPLE:  example autonom; shows an example
557"
558{
559  int NVars=nvars(basering);
560  int ExtIsZero;   
561  R=transpose(R);
562  ExtIsZero=is_zero_Our(Ext_Our(0,R));
563  int i=0;
564  while( (ExtIsZero)&&(i<=NVars) )
565  {
566    i++;
567    ExtIsZero = is_zero_Our(Ext_Our(i,R));
568  };
569
570  return(autonom_output(i,NVars));     
571}
572example
573{"EXAMPLE:"; echo = 2;
574  ring r=0,(s1,s2,s3,s4),dp;
575  module R= [s1,-s2],
576            [s2, s1],
577            [s3,-s4],
578            [s4, s3];       
579  R=transpose(R);
580  view( R );
581  view( autonom(R) );
582}; 
583
584//--------------------------------------------------------------------------
585//
586//Some example rings with defined systems
587//--------------------------------------------------------------------------
588//autonomy:
589//
590proc exAut1()
591{ 
592  ring @r=0,(s1,s2,s3),dp;
593  module R=[0,-s3,s2],
594           [s3,0,-s1];
595  R=transpose(R);         
596  export R;
597  return(@r);
598};
599//-------------------------------------------------------
600proc exAut2()
601{
602  ring @r=0,(s1,s2,s3),dp;
603  module R = [0,-s3,s2],
604             [s3,0,-s1],
605             [-s2,s1,0];
606  R=transpose(R);
607  export R;
608  return(@r);       
609}; 
610
611//-----------------------------------
612proc exCauchy()
613{
614  ring @r=0,(s1,s2),dp;
615  module R= [s1,-s2],
616            [s2, s1];
617  R=transpose(R);
618  export R;
619  return(@r);       
620};
621//-----------------------------------
622proc exCauchy2()
623{
624  ring @r=0,(s1,s2,s3,s4),dp;
625  module R= [s1,-s2],
626            [s2, s1],
627            [s3,-s4],
628            [s4, s3];       
629  R=transpose(R);
630  export R;
631  return(@r);       
632};
633//----------------------------------------
634proc exEvasLecture()
635{
636  ring @r=0,(d1,d2),dp;
637  module R=[d1^2-d2],
638           [d2^2-1];
639  R=transpose(R);
640  export R;
641  return(@r);       
642}; 
643//---------------------------------------------------------------------------
644//control
645//
646proc exAntenna()
647{
648  ring @r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), dp;
649  module R;
650  R = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
651      [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
652      [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
653      [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
654      [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
655      [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];
656
657  R=transpose(R);
658  export R;
659  return(@r); 
660};
661
662//---------------------------------------------------------------------------------------------
663
664proc exEinstein()
665{
666  ring @r = 0,(D(1..4)),dp;
667  module R =
668  [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)],
669  [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],
670  [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],
671  [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)],
672  [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)],
673  [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],
674  [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)],
675  [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)],
676  [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)],
677  [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];
678
679  R=transpose(R);
680  export R;
681  return(@r); 
682};
683
684
685//---------------------------------------------------------------------------------------------
686
687proc exFlexibleRod()
688{
689  ring @r = 0,(D1, delta), dp;
690  module R;
691  R = [D1, -D1*delta, -1], [2*D1*delta, -D1-D1*delta^2, 0];
692 
693  R=transpose(R);
694  export R;
695  return(@r); 
696};
697
698//---------------------------------------------------------------------------------------------
699proc exTwoPendula()
700{ 
701  ring @r=(0,m1,m2,M,g,L1,L2),Dt,dp;
702  module R = [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
703             [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
704             [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
705
706  R=transpose(R);
707  export R;
708  return(@r); 
709};
710//---------------------------------------------------------------------------------------------
711proc exWindTunnel()
712{
713  ring @r = (0,a, omega, zeta, k),(D1, delta),dp;
714  module R = [D1+a, -k*a*delta, 0, 0],
715             [0, D1, -1, 0],
716             [0, omega^2, D1+2*zeta*omega, -omega^2];
717
718  R=transpose(R);
719  export R;
720  return(@r); 
721};
722//--------------------------------------------------------------------------------------------
723
724proc ex()
725{
726};
727
728//---------------------------------------------------------------------------
729//---------------------------------------------------------------------------
730//---------------------------------------------------------------------------
731//---------------------------------------------------------------------------
732
733static proc invo_poly(poly m, map theta)
734//applies the involution map theta to m, where m=polynomial
735{
736  int i,j;
737  intvec v;
738  poly p,z;
739  poly n = 0;
740  i = 1;
741  while(m[i]!=0)
742  {
743    v = leadexp(m[i]);
744    z =1;
745    for(j=nvars(basering); j>=1; j--)
746    {
747      if (v[j]!=0)
748      {
749        p = var(j);
750        p = theta(p);
751        z = z*(p^v[j]);
752      }
753    }   
754    n = n + (leadcoef(m[i])*z);
755    i++; 
756  }
757  return(n);
758}
759
760proc involution(m, map theta)
761//applies the involution map theta to m, where m=vector, polynomial,
762//module,ideal
763{
764  int i,j;
765  intvec v;
766  poly p,z;
767  if (typeof(m)=="poly")
768  {
769    return (invo_poly(m,theta)); 
770  }
771  if ( typeof(m)=="ideal" )
772  {
773    ideal n;
774    for (i=1; i<=size(m); i++)
775    {
776      n[i] = invo_poly(m[i],theta);
777    }
778    return(n);
779  }
780  if (typeof(m)=="vector")
781  {
782    for(i=1;i<=size(m);i++)
783    {
784      m[i] = invo_poly(m[i],theta);
785    }
786    return (m); 
787  }
788 
789  if ( (typeof(m)=="matrix") || (typeof(m)=="module"))
790  { 
791//    m=transpose(m);
792    matrix n = matrix(m);
793    int @R=nrows(n);
794    int @C=ncols(n);
795    for(i=1; i<=@R; i++)
796    {
797      for(j=1; j<=@C; j++)
798      {
799        n[i,j] = invo_poly( m[i,j], theta);
800      }
801    }
802   }
803  if (typeof(m)=="module")
804  {
805    return (module(n));
806  }
807  return(n);
808}
809example
810{
811  "EXAMPLE:";echo = 2;
812  ring r = 0,(x,d),dp;
813  ncalgebra(1,1); // Weyl-Algebra
814  map F = r,x,-d;
815  poly f =  x*d^2+d;
816  poly If = involution(f,F);
817  f-If;
818  poly g = x^2*d+2*x*d+3*x+7*d;
819  poly tg = -d*x^2-2*d*x+3*x-7*d;
820  poly Ig = involution(g,F);
821  tg-Ig;
822  ideal I = f,g;
823  ideal II = involution(I,F);
824  II;
825  I - involution(II,F);
826  module M  = [f,g,0],[g,0,x^2*d];
827  module IM = involution(M,F);
828  print(IM);
829  print(M - involution(IM,F)); 
830}
831
832proc ncdetection( r)
833//in dieser proc. wird eine matrix erzeugt, die in der i-ten zeile die indices
834//der differential-,shift- oder advance-operatoren enthaelt mit denen die i-te
835//variable nicht kommutiert.
836{
837  int i,j,k,LExp;
838  int NVars=nvars(r);
839  matrix rel = NCRelations(r)[2];
840  intmat M[NVars][3];
841  int NRows = nrows(rel);
842  intvec v,w;
843  poly d,d_lead;
844  ideal I;
845  map theta;
846 
847  for( j=NRows;j>=2;j-- )
848  {
849   if( rel[j] == w )       //the whole column is zero
850    {
851      j--;
852      continue;
853    }
854   
855    for( i=1;i<j;i++ )         
856    {
857      if( rel[i,j]==1 )        //relation of type var(j)*var(i) = var(i)*var(j) +1
858      {
859         M[i,1]=j;
860      }
861      if( rel[i,j] == -1 )    //relation of type var(i)*var(j) = var(j)*var(i) -1
862      {
863        M[j,1]=i;
864      }
865      d = rel[i,j];
866      d_lead = lead(d);
867      v=leadexp(d_lead); //in the next lines we check wether we have a  relation of differential or shift type
868      LExp=0;
869      for( k=1;k<=NVars;k++)
870      {
871        LExp = LExp + v[k];
872      }
873      if( (d-d_lead != 0) || (LExp > 1) )
874      {
875        return( "wrong input" );
876      }
877      if( v[j] == 1)                   //relation of type var(j)*var(i) = var(i)*var(j) -lambda*var(j)
878      {
879        if (leadcoef(d) < 0)
880        {
881          M[i,2] = j;
882        }
883        else
884        {
885          M[i,3] = j;
886        }
887      }
888      if( v[i]==1 )                    //relation of type var(j)*var(i) = var(i)*var(j) -lambda*var(i)
889      {
890        if (leadcoef(d) > 0)
891        {
892          M[j,2] = i;
893        }
894        else
895        {
896          M[j,3] = i;
897        }
898      }
899    }
900  }
901  //ab hier wird die map ausgerechnet
902  for(i=1;i<=NVars;i++)
903  {
904    I=I+var(i);
905  }
906
907  for(i=1;i<=NVars;i++)
908  {
909    if( M[i,1..3]==(0,0,0) )
910    {
911      i++;
912      continue;
913    }
914    if( M[i,1]!=0 )
915    {
916      if( (M[i,2]!=0) && (M[i,3]!=0) )
917      {
918        I[M[i,1]] = -var(M[i,1]);
919        I[M[i,2]] = var(M[i,3]);
920        I[M[i,3]] = var(M[i,2]);
921      }
922      if( (M[i,2]==0) && (M[i,3]==0) )
923      {
924        I[M[i,1]] = -var(M[i,1]);
925      }                 
926      if( ( (M[i,2]!=0) && (M[i,3]==0) )|| ( (M[i,2]!=0) && (M[i,3]==0) )
927)
928      {
929        I[i] = -var(i);
930      }
931    }
932    else
933    {
934      if( (M[i,2]!=0) && (M[i,3]!=0) )
935      {
936        I[i] = -var(i);
937        I[M[i,2]] = var(M[i,3]);
938        I[M[i,3]] = var(M[i,2]);
939      }
940      else
941      {
942        I[i] = -var(i);
943      }
944    }
945  }
946  return(I);
947
948}
949example
950{
951  "EXAMPLE:"; echo = 2;
952  ring r=0,(x,y,z,D(1..3)),dp;
953  matrix D[6][6];
954  D[1,4]=1;
955  D[2,5]=1;
956  D[3,6]=1;
957  ncalgebra(1,D);
958  ncdetection(r);
959  kill r;
960  //----------------------------------------
961  ring r=0,(x,S),dp;
962  ncalgebra(1,-S);
963  ncdetection(r);
964  kill r;
965  //----------------------------------------
966  ring r=0,(x,D(1),S),dp;
967  matrix D[3][3];
968  D[1,2]=1;
969  D[1,3]=-S;
970  ncalgebra(1,D);
971  ncdetection(r);
972}
Note: See TracBrowser for help on using the repository browser.