Changeset 0ceca2 in git


Ignore:
Timestamp:
Dec 22, 2004, 7:03:16 PM (20 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
Children:
40d02c5bf05d362fee584c360803b8c93c23f239
Parents:
6fe0f6a88ab2e438b1590f5e383517ed59728be1
Message:
*levandov: ncdetection and involution moved to involution-lib, better examples and fixes for LeftInverse, RightKernel etc


git-svn-id: file:///usr/local/Singular/svn/trunk@7638 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/control.lib

    r6fe0f6 r0ceca2  
    1 version="$Id: control.lib,v 1.18 2004-12-16 16:28:53 levandov Exp $";
     1version="$Id: control.lib,v 1.19 2004-12-22 18:03:16 levandov Exp $";
    22category="Applications";
    33info="
     
    2323
    2424AUXILIARY PROCEDURES:
    25 ncdetection(ring r);          computes an ideal, presenting an involution map on non-comm algebra r;
    26 involution(m, map theta);    applies the involution, presented by theta to m of type poly, vector, ideal, module;
    2725declare(string NameOfRing, Variables[,string  Parameters, Ordering]);     defines the ring, optional parametes are strings of parameters and ordering,
    2826view();                      Well-formatted output of lists, modules and matrixes
     
    5553// NOTE: static things should not be shown for end-user
    5654// static Ext_Our(...)                  Copy of Ext_R from 'homolog.lib' in commutative case;
    57 // static is_zero_Our(module R)         Copy of is_zero from 'OBpoly.lib';
     55// static is_zero_Our(module R)         Copy of is_zero from 'poly.lib';
    5856//  static space(int n)           Procedure used inside the procedure 'Print' to have a better formatted output
    5957// static control_output();      Generating the output for the procedure 'control'
     
    139137//  print(" To compute controllability please enter: control(R)");
    140138//  print(" To compute autonomy please enter: autonom(R)");
    141 //
    142 //
    143 //
    144139//-------------------------------------------------------------------------
     140
    145141static proc space(int n)
    146142"USAGE:spase(n);
     
    239235//--------------------------------------------------------------------------
    240236proc RightKernel(matrix M)
    241 "USAGE:  RightKernel(M);
    242            M:  matrix
     237"USAGE:  RightKernel(M);   M a matrix
    243238RETURN:  right kernel of matrix M, i.e., the module of all elements v such that Mv=0
    244 NOTE:  in commutative case it is a left  module, in noncommutative (will be implemented later) it is a right module   
     239NOTE:    in the noncommutative case RightKernel is a right module
    245240EXAMPLE:  example RightKernel; shows an example
    246241"
     
    250245example
    251246{"EXAMPLE:";echo = 2;
    252   ring r;
     247  ring r = 0,(x,y,z),dp;
    253248  matrix M[1][3] = x,y,z;
    254249  print(M);
    255   print( RightKernel(M) );
     250  matrix R = RightKernel(M);
     251  print(R);
     252  // check:
     253  print(M*R);
    256254};
    257255//-------------------------------------------------------------------------
    258256proc LeftKernel(matrix M)
    259 "USAGE:  LeftKernel(M);
    260            M:  matrix
    261 RETURN:  left kernel of matrix M, i.e., the matrix whose rows are generators of left module
    262          (elements of this module are to be rows) of all elements v such that vM=0   
     257"USAGE:  LeftKernel(M);   M a matrix
     258RETURN:  left kernel of matrix M, i.e. the module of all elements v such that vM=0
    263259EXAMPLE:  example LeftKernel; shows an example
    264260"
     
    268264example
    269265{"EXAMPLE:";echo = 2;
    270   ring r;
     266  ring r= 0,(x,y,z),dp;
    271267  matrix M[3][1] = x,y,z;
    272268  print(M);
    273   print( LeftKernel(M) );
     269  matrix L = LeftKernel(M);
     270  print(L);
     271  // check:
     272  print(L*M);
    274273};
    275274//------------------------------------------------------------------------
     
    277276"USAGE:  LeftInverse(M);  M a module
    278277RETURN:  matrix L such that LM == Id;
    279 EXAMPLE:  example LeftInverse; shows an example
     278EXAMPLE: example LeftInverse; shows an example
    280279NOTE: exists only in the case when Id \subseteq M!
    281280"
    282281{
    283282  // now it works also for the NC case;
    284   int NCols=ncols(M);
     283  // two approaches are possible;
     284  // the older one is commented out
     285  int NCols = ncols(M);
    285286  module Id = freemodule(NCols);
    286   M = transpose(M);
     287  module N = transpose(M);
    287288  Id = std(Id);
    288289  matrix T;
    289290  option(redSB);
    290291  option(redTail);
    291   // check the correctness:
    292   // Id \subseteq M
    293   module N = liftstd(M, T);
    294   ideal TT = ideal(NF(N,Id));
    295   TT = simplify(TT,2);
    296 
    297   if (TT[1]!=0)
    298   {
    299     "No left inverse exists";
     292  // check the correctness (Id \subseteq M)
     293  // via dimension: {dim/GKdim} (M) = -1!
     294  int d = dim_Our(N);
     295  if (d != -1)
     296  {
     297    //    "No left inverse exists";
    300298    return(matrix(0));
    301299  }
     300  //  module N = liftstd(M, T);
     301  //  ideal TT = ideal(NF(N,Id));
     302  //  TT = simplify(TT,2);
     303  //  if (TT[1]!=0)
     304  //  {
     305  //    "No left inverse exists";
     306  //    return(matrix(0));
     307  //  }
    302308  matrix T2 = lift(N, Id);
    303   T2 = transpose(T2)*transpose(T);
     309  //  T2 = transpose(T2)*transpose(T);
    304310  T2 =  transpose(T2);
     311  // set options back
     312  option(noredTail);
    305313  return(T2);
    306314}
    307315example
    308 {"EXAMPLE:";echo =2;
    309   ring r;
     316{ // trivial example:
     317"EXAMPLE:";echo =2;
     318  ring r = 0,(x,z),dp;
    310319  matrix M[2][1] = 1,x2z;
    311320  print(M);
    312321  print( LeftInverse(M) ); 
     322  kill r;
     323  // derived from the exTwoPendula
     324  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
     325  matrix U[3][1];
     326  U[1,1]=(-L2)*Dt^4+(g)*Dt^2;
     327  U[2,1]=(-L1)*Dt^4+(g)*Dt^2;
     328  U[3,1]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2);
     329  module M = module(U);
     330  module L = LeftInverse(M);
     331  print(L);
    313332};
    314333//-----------------------------------------------------------------------
     
    316335{
    317336  int d;
     337  if (attrib(R,"isSB")<>1)
     338  {
     339    R = std(R);
     340  }
    318341  if (Is_NC())
    319342  {
     
    328351  if (Is_NC())
    329352  {
     353    // TODO: use preAnn, that is an annihilator of the generators
     354    // what is a left ideal!
    330355    return(-1);
    331356  }
     
    338363//-----------------------------------------------------------------------
    339364static proc Ext_Our(int i, module R, list #)
    340 // just a copy of 'Ext_R' from "homolog.lib" in commutative case
     365// mimicking 'Ext_R' from "homolog.lib" in commutative case
    341366{
    342367  int ISNC = Is_NC();
     
    369394//------------------------------------------------------------------------
    370395static proc control_output(int i, int NVars, module R, module Ext_1)
    371 "USAGE:  proc control_output(i, NVars, R, Ext_1)
     396"USAGE:  control_output(i, NVars, R, Ext_1), where
    372397           i:  integer, number of first nonzero Ext or
    373                just number of variables in a base ring + 1 in case that all the Exts are zero
    374            NVars:  integer, number of variables in a base ring 
    375            R:  module R (cokernel representation)
     398               just number of variables in a base ring + 1 in case that all the Exts are zero,
     399           NVars:  integer, number of variables in a base ring,
     400           R:  module R (cokernel representation),
    376401           Ext_1:  module, the first Ext(its cokernel representation)     
    377402RETURN:  list with all the contollability properties of the system which is to be returned in 'control' procedure
     
    379404"
    380405{
    381   int d = dim_Our( std( transpose(R) ) ) ;; //this is the dimension of the system
    382   string DofS= "dimension of the system:";
    383   string Fn= "number of first nonzero Ext:";
    384   if(i==1)
     406  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
     407  int d = dim_Our( std( transpose(R) ) ); // this is the dimension of the system
     408  string DofS = "dimension of the system:";
     409  string Fn   = "number of first nonzero Ext:";
     410  if (i==1)
    385411  {
    386     module RK = RightKernel(R);
     412    module RK = RightKernel(R); 
    387413    return(
    388414            list ( Fn,
     
    403429 
    404430  if(i>NVars)
    405   { module RK =RightKernel(R);
     431  { module RK = RightKernel(R);
    406432    return( list(  Fn,
    407433                   -1,
     
    439465                   d)
    440466          );
    441   };
    442              
    443  
     467  }; 
    444468}; 
    445469//-------------------------------------------------------------------------
    446470
    447471proc control(module R)
    448 "USAGE:  control(R);
    449            R:  module (R is a matrix of the system of equations which is to be investigated)
    450 RETURN:  list of all the properties concerning controllability of the system(behavior) represented by the  matrix R
     472"USAGE:  control(R); R a module (R is the matrix of the system of equations which is to be investigated)
     473RETURN:  list of all the properties concerning controllability of the system (behavior) represented by the  matrix R
    451474EXAMPLE:  example control; shows an example
    452475"
     
    454477  int i;
    455478  int NVars=nvars(basering);
     479  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
    456480  int ExtIsZero;
    457  
    458        
     481
    459482  module Ext_1 = std(Ext_Our(1,R));
    460483   
     
    480503  view(R);
    481504  view(control(R));
    482 
    483505};
    484506//------------------------------------------------------------------------
     
    521543     
    522544  if( i==1 )
    523   //in case that NVars==1 there is no sence to consider the notion
    524   //of strongly autonomous behavior, because it does not imply
    525   //that system is overdetermined in this case
     545  // in case that NVars==1 there is no sense to consider the notion
     546  // of strongly autonomous behavior, because it does not imply
     547  // that system is overdetermined in this case
    526548  {
    527549    return( list ( Fn,
     
    552574          );
    553575  };
    554    
    555576}; 
    556577//--------------------------------------------------------------------------
    557578proc autonom2(module R)
    558 "USAGE:  autonom2(R);
    559            R:  module (R is a matrix of the system of equations which is to be investigated)
    560 RETURN:  list of all the properties concerning autonomy of the system(behavior) represented by the  matrix R
     579"USAGE:  autonom2(R);   R a module (R is a matrix of the system of equations which is to be investigated)
     580RETURN:  list of all the properties concerning autonomy of the system (behavior) represented by the matrix R
    561581NOTE:  this procedure is an analogue to 'autonom' using dimension calculations
    562582EXAMPLE:  example autonom2; shows an example
     
    581601  view( autonom2(R) );
    582602}; 
    583 //---------------------------------------------------------------------------
    584 
     603//----------------------------------------------------------
    585604proc autonom(module R)
    586 "USAGE:  autonom(R);
    587            R:  module (R is a matrix of the system of equations which is to be investigated)
    588 RETURN:  list of all the properties concerning autonomy of the system(behavior) represented by the  matrix R
     605"USAGE:  autonom(R);   R a module (R is a matrix of the system of equations which is to be investigated)
     606RETURN:  list of all the properties concerning autonomy of the system (behavior) represented by the  matrix R
    589607EXAMPLE:  example autonom; shows an example
    590608"
     
    600618    ExtIsZero = is_zero_Our(Ext_Our(i,R));
    601619  };
    602 
    603620  return(autonom_output(i,NVars));     
    604621}
     
    616633}; 
    617634
    618 //--------------------------------------------------------------------------
     635//----------------------------------------------------------
    619636//
    620637//Some example rings with defined systems
    621 //----------------------------------------------------------------------------
     638//----------------------------------------------------------
    622639//autonomy:
    623640//
    624 //----------------------------------------------------------------------------
     641//----------------------------------------------------------
    625642proc exCauchy()
    626643{
     
    632649  return(@r);       
    633650};
    634 //----------------------------------------------------------------------------
     651//----------------------------------------------------------
    635652proc exCauchy2()
    636653{
     
    644661  return(@r);       
    645662};
    646 //----------------------------------------------------------------------------
     663//----------------------------------------------------------
    647664proc exZerz()
    648665{
     
    654671  return(@r);       
    655672}; 
    656 //----------------------------------------------------------------------------
     673//----------------------------------------------------------
    657674//control
    658675//
     
    666683  return(@r);
    667684};
    668 //----------------------------------------------------------------------------
     685//----------------------------------------------------------
    669686proc ex2()
    670687{
     
    677694  return(@r);
    678695};
    679 //----------------------------------------------------------------------------
     696//----------------------------------------------------------
    680697proc exAntenna()
    681698{
     
    694711};
    695712
    696 //----------------------------------------------------------------------------
     713//----------------------------------------------------------
    697714
    698715proc exEinstein()
     
    717734
    718735
    719 //---------------------------------------------------------------------------------------------
     736//----------------------------------------------------------
    720737
    721738proc exFlexibleRod()
     
    730747};
    731748
    732 //---------------------------------------------------------------------------------------------
     749//----------------------------------------------------------
    733750proc exTwoPendula()
    734751
     
    742759  return(@r); 
    743760};
    744 //---------------------------------------------------------------------------------------------
     761//----------------------------------------------------------
    745762proc exWindTunnel()
    746763{
     
    754771  return(@r); 
    755772};
    756 //--------------------------------------------------------------------------------------------
    757 
    758 
    759 //---------------------------------------------------------------------------
    760 //---------------------------------------------------------------------------
    761 //---------------------------------------------------------------------------
    762 //---------------------------------------------------------------------------
    763 
    764 static proc invo_poly(poly m, map theta)
    765 //applies the involution map theta to m, where m=polynomial
    766 {
    767   int i,j;
    768   intvec v;
    769   poly p,z;
    770   poly n = 0;
    771   i = 1;
    772   while(m[i]!=0)
    773   {
    774     v = leadexp(m[i]);
    775     z =1;
    776     for(j=nvars(basering); j>=1; j--)
    777     {
    778       if (v[j]!=0)
    779       {
    780         p = var(j);
    781         p = theta(p);
    782         z = z*(p^v[j]);
    783       }
    784     }   
    785     n = n + (leadcoef(m[i])*z);
    786     i++; 
    787   }
    788   return(n);
    789 }
    790 
    791 proc involution(m, map theta)
    792 //applies the involution map theta to m, where m=vector, polynomial,
    793 //module,ideal
    794 {
    795   int i,j;
    796   intvec v;
    797   poly p,z;
    798   if (typeof(m)=="poly")
    799   {
    800     return (invo_poly(m,theta)); 
    801   }
    802   if ( typeof(m)=="ideal" )
    803   {
    804     ideal n;
    805     for (i=1; i<=size(m); i++)
    806     {
    807       n[i] = invo_poly(m[i],theta);
    808     }
    809     return(n);
    810   }
    811   if (typeof(m)=="vector")
    812   {
    813     for(i=1;i<=size(m);i++)
    814     {
    815       m[i] = invo_poly(m[i],theta);
    816     }
    817     return (m); 
    818   }
    819  
    820   if ( (typeof(m)=="matrix") || (typeof(m)=="module"))
    821   { 
    822 //    m=transpose(m);
    823     matrix n = matrix(m);
    824     int @R=nrows(n);
    825     int @C=ncols(n);
    826     for(i=1; i<=@R; i++)
    827     {
    828       for(j=1; j<=@C; j++)
    829       {
    830         n[i,j] = invo_poly( m[i,j], theta);
    831       }
    832     }
    833    }
    834   if (typeof(m)=="module")
    835   {
    836     return (module(n));
    837   }
    838   return(n);
    839 }
    840 example
    841 {
    842   "EXAMPLE:";echo = 2;
    843   ring r = 0,(x,d),dp;
    844   ncalgebra(1,1); // Weyl-Algebra
    845   map F = r,x,-d;
    846   poly f =  x*d^2+d;
    847   poly If = involution(f,F);
    848   f-If;
    849   poly g = x^2*d+2*x*d+3*x+7*d;
    850   poly tg = -d*x^2-2*d*x+3*x-7*d;
    851   poly Ig = involution(g,F);
    852   tg-Ig;
    853   ideal I = f,g;
    854   ideal II = involution(I,F);
    855   II;
    856   I - involution(II,F);
    857   module M  = [f,g,0],[g,0,x^2*d];
    858   module IM = involution(M,F);
    859   print(IM);
    860   print(M - involution(IM,F)); 
    861 }
    862 
    863 proc ncdetection(def r)
    864 // in this procedure an involution map is generated from the NCRelations
    865 // that will be used in the function involution
    866 // in dieser proc. wird eine matrix erzeugt, die in der i-ten zeile die indices
    867 // der differential-,shift- oder advance-operatoren enthaelt mit denen die i-te
    868 // variable nicht kommutiert.
    869 {
    870   int i,j,k,LExp;
    871   int NVars=nvars(r);
    872   matrix rel = NCRelations(r)[2];
    873   intmat M[NVars][3];
    874   int NRows = nrows(rel);
    875   intvec v,w;
    876   poly d,d_lead;
    877   ideal I;
    878   map theta;
    879  
    880   for( j=NRows; j>=2; j-- )
    881   {
    882    if( rel[j] == w )       //the whole column is zero
    883     {
    884       j--;
    885       continue;
    886     }
    887    
    888     for( i=1; i<j; i++ )         
    889     {
    890       if( rel[i,j]==1 )        //relation of type var(j)*var(i) = var(i)*var(j) +1
    891       {
    892          M[i,1]=j;
    893       }
    894       if( rel[i,j] == -1 )    //relation of type var(i)*var(j) = var(j)*var(i) -1
    895       {
    896         M[j,1]=i;
    897       }
    898       d = rel[i,j];
    899       d_lead = lead(d);
    900       v = leadexp(d_lead); //in the next lines we check wether we have a  relation of differential or shift type
    901       LExp=0;
    902       for(k=1; k<=NVars; k++)
    903       {
    904         LExp = LExp + v[k];
    905       }
    906       //      if( (d-d_lead != 0) || (LExp > 1) )
    907 if ( ( (d-d_lead) != 0) || (LExp > 1) || ( (LExp==0) && ((d_lead>1) || (d_lead<-1)) ) )
    908       {
    909         return(theta);
    910       }
    911 
    912       if( v[j] == 1)                   //relation of type var(j)*var(i) = var(i)*var(j) -lambda*var(j)
    913       {
    914         if (leadcoef(d) < 0)
    915         {
    916           M[i,2] = j;
    917         }
    918         else
    919         {
    920           M[i,3] = j;
    921         }
    922       }
    923       if( v[i]==1 )                    //relation of type var(j)*var(i) = var(i)*var(j) -lambda*var(i)
    924       {
    925         if (leadcoef(d) > 0)
    926         {
    927           M[j,2] = i;
    928         }
    929         else
    930         {
    931           M[j,3] = i;
    932         }
    933       }
    934     }
    935   }
    936   // from here on, the map is computed
    937   for(i=1;i<=NVars;i++)
    938   {
    939     I=I+var(i);
    940   }
    941 
    942   for(i=1;i<=NVars;i++)
    943   {
    944     if( M[i,1..3]==(0,0,0) )
    945     {
    946       i++;
    947       continue;
    948     }
    949     if( M[i,1]!=0 )
    950     {
    951       if( (M[i,2]!=0) && (M[i,3]!=0) )
    952       {
    953         I[M[i,1]] = -var(M[i,1]);
    954         I[M[i,2]] = var(M[i,3]);
    955         I[M[i,3]] = var(M[i,2]);
    956       }
    957       if( (M[i,2]==0) && (M[i,3]==0) )
    958       {
    959         I[M[i,1]] = -var(M[i,1]);
    960       }                 
    961       if( ( (M[i,2]!=0) && (M[i,3]==0) )|| ( (M[i,2]!=0) && (M[i,3]==0) )
    962 )
    963       {
    964         I[i] = -var(i);
    965       }
    966     }
    967     else
    968     {
    969       if( (M[i,2]!=0) && (M[i,3]!=0) )
    970       {
    971         I[i] = -var(i);
    972         I[M[i,2]] = var(M[i,3]);
    973         I[M[i,3]] = var(M[i,2]);
    974       }
    975       else
    976       {
    977         I[i] = -var(i);
    978       }
    979     }
    980   }
    981   return(I);
    982 }
    983 example
    984 {
    985   "EXAMPLE:"; echo = 2;
    986   ring r=0,(x,y,z,D(1..3)),dp;
    987   matrix D[6][6];
    988   D[1,4]=1;
    989   D[2,5]=1;
    990   D[3,6]=1;
    991   ncalgebra(1,D);
    992   ncdetection(r);
    993   kill r;
    994   //----------------------------------------
    995   ring r=0,(x,S),dp;
    996   ncalgebra(1,-S);
    997   ncdetection(r);
    998   kill r;
    999   //----------------------------------------
    1000   ring r=0,(x,D(1),S),dp;
    1001   matrix D[3][3];
    1002   D[1,2]=1;
    1003   D[1,3]=-S;
    1004   ncalgebra(1,D);
    1005   ncdetection(r);
    1006 }
     773//----------------------------------------------------------
    1007774
    1008775proc genericity(matrix M)
    1009776"USAGE:  genericity(M), M is a matrix
    1010 RETURN:  list of strings with expressions, used as non-zero in SB
    1011 NOTE:  effective with the liftstd procedure
     777RETURN:  list of strings with expressions, which have been assumed to be non-zero in the process of computing the Groebner basis
     778NOTE:  effective with the liftstd procedure for modules with parameters
    1012779"
    1013780{
     
    1085852  module SR = liftstd(R,T);
    1086853  genericity(T);
     854  //-- The result is different when computing reduced bases:
     855  matrix T2;
     856  option(redSB);
     857  option(redTail);
     858  module SR2 = liftstd(R,T2);
     859  genericity(T2);
    1087860}
    1088861
    1089862proc canonize(list L)
    1090863"USAGE:  canonize(L), L a list
    1091 ASSUME: L is the output of control/autonomy procs
    1092 RETURN:  canonized list
     864ASSUME:  L is the output of control/autonomy procs
     865RETURN:  a canonized list, where modules are canonized by computing
     866their reduced minimal Groebner bases
    1093867"
    1094868{
     
    1103877    {
    1104878      M[i] = std(M[i]);
    1105       //      M[i] = prune(M[i]); // mimimal embedding
     879      //      M[i] = prune(M[i]); // mimimal embedding: no need yet
    1106880      //      M[i] = std(M[i]);
    1107881    }
    1108882  }
     883  option(noredTail); // set the initial option back
    1109884  return(M);
    1110885}
     
    13281103  option(redTail);
    13291104  ring r=0,x,dp;
    1330   // This is example is trivial, but we want to see,
     1105  // This example is trivial, but we want to see,
    13311106  // that nothing is changed, when the matrix is already in Smith-Form!
    13321107  module M = [x,0,0],[0,x2,0],[0,0,x3];
Note: See TracChangeset for help on using the changeset viewer.