Changeset 6045c3 in git


Ignore:
Timestamp:
Dec 22, 2008, 10:45:09 PM (14 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')
Children:
565e866c54dc28cc8085355e58a8d9ba2b7e0d2c
Parents:
3e0c94f2a56b2afdc4e51bebb82b8a047d4487ed
Message:
*levandov: docu and minor bugfixes, more examples, formatting


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/jacobson.lib

    r3e0c94f r6045c3  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="$Id: jacobson.lib,v 1.5 2008-12-09 16:50:21 levandov Exp $";
     2version="$Id: jacobson.lib,v 1.6 2008-12-22 21:45:09 levandov Exp $";
    33category="System and Control Theory";
    44info="
    … …  
    77
    88THEORY: We work over a ring R, that is a principal ideal domain.
    9 @*  If R is commutative, we suppose R to be a polynomial ring in one variable.
    10 @*  If R is non-commutative, we suppose R to be in two variables, say x and d. We treat
    11 @* then the Ore localization of R with respect to the mult.closed set S = K[x] without 0.
    12 @*  Given a rectangular matrix M over R, one can compute unimodular matrices U, V
    13 @*  such that U*M*V=D is a diagonal matrix.
    14 @*  Depending on the ring, the diagonal entries of D have certain properties.
     9@*   If R is commutative, we suppose R to be a polynomial ring in one variable.
     10@* If R is non-commutative, we suppose R to be in two variables, say x and d.
     11@* We treat then the basering as principal ideal ring with d a polynomial
     12@* variable and x an invertible one. That is, we work in the Ore localization of R
     13@* with respect to the mult. closed set S = K[x] without 0.
     14@* Note, that in computations no division by x will actually happen.
     15@*   Given a rectangular matrix M over R, one can compute unimodular (that is invertible)
     16@* square matrices U and V, such that U*M*V=D is a diagonal matrix.
     17@*   Depending on the ring, the diagonal entries of D have certain properties.
     18
     19REFERENCES:
     20@* N. Jacobson, 'The theory of rings', AMS, 1943.
     21@* Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner: Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'. PhD thesis, Universidad de Santiago de Compostela, 2005.
     22
    1523
    1624PROCEDURES:
    17 smith(R,M[,eng1,eng2]);  compute the Smith Normal Form of M over commutative ring
    18 jacobson(R,M[,eng]);     compute a weak Jacobson Normal Form of M over non-commutative ring, i.e. a diagonal matrix
     25smith(M[,eng1,eng2]);  compute the Smith Normal Form of M over commutative ring
     26jacobson(M[,eng]);       compute a weak Jacobson Normal Form of M over non-commutative ring
    1927
    2028SEE ALSO: control_lib
    2129";
     30
    2231  LIB "poly.lib";
    2332  LIB "involut.lib";
    2433
    25 proc smith(R, matrix MA, list #)
    26 "USAGE:  smith(R, M[, eng1, eng2]);  R ring, M matrix, eng1 and eng2 are optional int
    27 RETURN: list
    28 NOTE: By default, the diagonal matrix, the Smith normal form, is returned.
    29 @* R stays for the ring, where computations will take place.
    30 @* If optional integer eng1 is nonzero, smith returns
    31 @* a list of matrices L such that L[1]*M*L[3]=L[2] with L[2]_11|L[2]_22|...|L[2]_nn and
    32 @* L[1], L[3] unimodular matrices.
    33 @* eng2 determines the engine, that computes the Groebner basis. By default eng2 equals zero.
    34 @* If optional integer eng2 = 0 than std is used to caculate a Groebner basis ,
    35 @* If eng2 = 1 than groebner is used to caculate a Groebner basis
    36 @* If eng2 = 2 than slimgb is used to caculate a Groebner basis
    37 @* If @code{printlevel}=1, progress debug messages will be printed,
     34proc smith(matrix MA, list #)
     35"USAGE: smith(M[, eng1, eng2]);  M matrix, eng1 and eng2 are optional integers
     36RETURN: matrix or list
     37ASSUME: The current ring is assumed to be the commutative polynomial ring in
     38        one variable
     39PURPOSE: compute the Smith Normal Form of M with transformation matrices (optionally)
     40NOTE: If the optional integer eng1 is non-zero, returns the list {U,D,V},
     41where U*M*V = D and the diagonal field entries of D are not normalized.
     42@*    Otherwise only the matrix D, that is the Smith Normal Form of M, is returned.
     43@*    Here, U and V are square unimodular (invertible) matrices. The procedure works for rectangular matrix M.
     44@* The optional integer eng2 determines the engine, that computes the Groebner basis:
     45@* 0 means 'std' (default), 1 means 'groebner' and 2 means 'slimgb'.
     46DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
    3847@* if @code{printlevel}>=2, all the debug messages will be printed.
    3948EXAMPLE: example smith; shows examples
    4049"
    4150{
     51  def R = basering;
     52  // check assume: R must be a polynomial ring in 1 variable
     53  setring R;
     54  if (nvars(R) >1 )
     55  {
     56    ERROR("Basering must have exactly one variable");
     57  }
     58
    4259  int eng = 0;
    4360  int BASIS;
    … …  
    102119  ring r = 0,x,Dp;
    103120  matrix m[3][2]=x, x^4+x^2+21, x^4+x^2+x, x^3+x, 4*x^2+x, x;
    104   list s=smith(r,m,1);
    105   print(s[2]);
    106   print(s[1]*m*s[3]);
     121  list s=smith(m,1);
     122  print(s[2]);  // Smith form of m
     123  print(s[1]*m*s[3] - s[2]); // check U*M*V = D
    107124}
    108125
    … …  
    554571}
    555572
    556 proc jacobson(R, matrix MA, list #)
    557  "USAGE:  jacobson(R, matrix MA, eng);  R ring, M matrix, eng an optional int
     573proc jacobson(matrix MA, list #)
     574 "USAGE:  jacobson(M, eng); M matrix, eng an optional int
    558575RETURN: list
    559 NOTE: A list of matrices L such that L[1]*M*L[3]=L[2] such that L[2] is a diagonal matrix and
    560 @* L[1], L[3] unimodular matrices.
    561 @* R stays for the ring, where computations will take place.
    562 @* eng determines the engine, that computes the Groebner basis. By default eng equals zero.
    563 @* If eng = 0 than std is used to caculate a Groebner basis
    564 @* If eng = 1 than groebner is used to caculate a Groebner basis
    565 @* If eng = 2 than slimgb is used to caculate a Groebner basis
    566 @* If @code{printlevel}=1, progress debug messages will be printed,
     576ASSUME: Basering is a noncommutative ring in two variables.
     577PURPOSE: compute a weak Jacobson Normal Form of M over a noncommutative ring
     578NOTE: A list L of matrices {U,D,V} is returned. That is L[1]*M*L[3]=L[2], where
     579@*      L[2] is a diagonal matrix and L[1], L[3] square invertible (unimodular) matrices.
     580@*      Note, that M can be rectangular.
     581@* The optional integer @code{eng} determines the engine, that computes the Groebner basis:
     582@* 0 means 'std' (default), 1 means 'groebner' and 2 means 'slimgb'.
     583DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
    567584@* if @code{printlevel}>=2, all the debug messages will be printed.
    568585EXAMPLE: example jacobson; shows examples
    569586"
    570587{
     588  def R = basering;
     589  // check assume: R must be a polynomial ring in 2 variables
     590  setring R;
     591  if ( (nvars(R) !=2 ) )
     592  {
     593    ERROR("Basering must have exactly two variables");
     594  }
     595
     596  // check if MA is zero matrix and return corr. U,V
     597  if ( (size(module(MA))==0) and (size(transpose(module(MA)))==0) )
     598  {
     599    int nr = nrows(MA);
     600    int nc = ncols(MA);
     601    matrix U = matrix(freemodule(nr));
     602    matrix V = matrix(freemodule(nc));
     603    list L; L[1]=U; L[2] = MA; L[3] = V;
     604    return(L);
     605  }
     606
    571607  int B=0;
    572608  if ( size(#)>0 )
    … …  
    615651}
    616652example
    617 { "EXAMPLE:"; echo = 2;
     653{
     654  "EXAMPLE:"; echo = 2;
    618655  ring r = 0,(x,d),Dp;
    619   def R=nc_algebra(1,1); // the Weyl algebra
    620   setring R;
    621   matrix m[2][2]=d,x,0,d; print(m);
    622   list J=jacobson(R,m); // returns a list with 3 entries
    623   print(J[2]); // a Jacobson Form D
    624   print(J[1]*m*J[3]); // check that U*M*V = D
    625 }
    626 
    627 
     656  def R = nc_algebra(1,1);   setring R; // the 1st Weyl algebra
     657  matrix m[2][2] = d,x,0,d; print(m);
     658  list J = jacobson(m); // returns a list with 3 entries
     659  print(J[2]); // a Jacobson Form D for m
     660  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D
     661  /*   now, let us do the same for the shift algebra  */
     662  ring r2 = 0,(x,s),Dp;
     663  def R2 = nc_algebra(1,s);   setring R2; // the 1st shift algebra
     664  matrix m[2][2] = s,x,0,s; print(m); // matrix of the same for as above
     665  list J = jacobson(m);
     666  print(J[2]); // a Jacobson Form D, quite different from above
     667  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D
     668}
    628669
    629670static proc triangle( matrix MA, int B)
    … …  
    844885/////interesting examples for smith////////////////
    845886
    846 static proc Ex_One_wheeled_bicycle()
     887/*
     888
     889//static proc Ex_One_wheeled_bicycle()
    847890{
    848891ring RA=(0,m), D, lp;
    849892matrix bicycle[2][3]=(1+m)*D^2, D^2, 1, D^2, D^2-1, 0;
    850 list s=smith(RA,bicycle, 1,0);
     893list s=smith(bicycle, 1,0);
    851894print(s[2]);
    852895print(s[1]*bicycle*s[3]-s[2]);
    … …  
    854897
    855898
    856 static proc Ex_RLC-circuit()
     899//static proc Ex_RLC-circuit()
    857900{
    858901ring RA=(0,m,R1,R2,L,C), D, lp;
    859902matrix  circuit[2][3]=D+1/(R1*C), 0, -1/(R1*C), 0, D+R2/L, -1/L;
    860 list s=smith(RA,circuit, 1,0);
     903list s=smith(circuit, 1,0);
    861904print(s[2]);
    862905print(s[1]*circuit*s[3]-s[2]);
    … …  
    864907
    865908
    866 static proc Ex_two_pendula()
     909//static proc Ex_two_pendula()
    867910{
    868911ring RA=(0,m,M,L1,L2,m1,m2,g), D, lp;
    869 
    870912matrix pendula[3][4]=m1*L1*D^2,m2*L2*D^2,(M+m1+m2)*D^2,-1,m1*L1^2*D^2-m1*L1*g,0,m1*L1*D^2,0,0,
    871913m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0;
    872 list s=smith(RA,pendula, 1,0);
     914list s=smith(pendula, 1,0);
    873915print(s[2]);
    874916print(s[1]*pendula*s[3]-s[2]);
    875917}
    876918
    877 
    878 
    879 static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit()
     919//static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit()
    880920{
    881921ring RA=(0,m,omega,r,a,b), D, lp;
    882 
    883922matrix satellite[4][6]=
    884923D,-1,0,0,0,0,
    … …  
    8869250,0,D,-1,0,0,
    8879260,2*omega/r,0,D,0,-b/(m*r);
    888 
    889 list s=smith(RA,satellite, 1,0);
     927list s=smith(satellite, 1,0);
    890928print(s[2]);
    891929print(s[1]*satellite*s[3]-s[2]);
    892930}
    893931
    894 static proc Ex_flexible_one_link_robot()
     932//static proc Ex_flexible_one_link_robot()
    895933{
    896934ring RA=(0,M11,M12,M13,M21,M22,M31,M33,K1,K2), D, lp;
    897 
    898935matrix robot[3][4]=M11*D^2,M12*D^2,M13*D^2,-1,M21*D^2,M22*D^2+K1,0,0,M31*D^2,0,M33*D^2+K2,0;
    899 list s=smith(RA,robot, 1,0);
     936list s=smith(robot, 1,0);
    900937print(s[2]);
    901938print(s[1]*robot*s[3]-s[2]);
    902939}
    903940
     941*/
    904942
    905943
    906944/////interesting examples for jacobson////////////////
    907945
    908 static proc Ex_compare_output_with_maple_package_JanetOre() 
    909 {   ring w = 0,(x,d),Dp;
     946/*
     947
     948//static proc Ex_compare_output_with_maple_package_JanetOre() 
     949{ 
     950    ring w = 0,(x,d),Dp;
    910951    def W=nc_algebra(1,1);
    911952    setring W;
    912953    matrix m[3][3]=[d2,d+1,0],[d+1,0,d3-x2*d],[2d+1, d3+d2, d2];
    913     list J=jacobson(W,m,0);
     954    list J=jacobson(m,0);
    914955    print(J[1]*m*J[3]);
    915956    print(J[2]);
    … …  
    919960}
    920961
    921 
    922 static proc Ex_cyclic()
    923 {   ring w = 0,(x,d),Dp;
    924     def W=nc_algebra(1,1);
     962// modification for shift algebra
     963{
     964    ring w = 0,(x,s),Dp;
     965    def W=nc_algebra(1,s);
    925966    setring W;
    926     matrix m[3][3]=d,0,0,x*d+1,d,0,0,x*d,d;
    927     list J=jacobson(W,m,0);
     967    matrix m[3][3]=[s^2,s+1,0],[s+1,0,s^3-x^2*s],[2*s+1, s^3+s^2, s^2];
     968    list J=jacobson(m,0);
    928969    print(J[1]*m*J[3]);
    929970    print(J[2]);
    … …  
    933974}
    934975
     976//static proc Ex_cyclic()
     977{   
     978    ring w = 0,(x,d),Dp;
     979    def W=nc_algebra(1,1);
     980    setring W;
     981    matrix m[3][3]=d,0,0,x*d+1,d,0,0,x*d,d;
     982    list J=jacobson(m,0);
     983    print(J[1]*m*J[3]);
     984    print(J[2]);
     985    print(J[1]);
     986    print(J[3]);
     987    print(J[1]*m*J[3]-J[2]);
     988}
     989
     990// modification for shift algebra: gives a module with ann = s^2
     991{
     992    ring w = 0,(x,s),Dp;
     993    def W=nc_algebra(1,s);
     994    setring W;
     995    matrix m[3][3]=s,0,0,x*s+1,s,0,0,x*s,s;
     996    list J=jacobson(m,0);
     997    print(J[1]*m*J[3]);
     998    print(J[2]);
     999    print(J[1]);
     1000    print(J[3]);
     1001    print(J[1]*m*J[3]-J[2]);
     1002}
     1003
     1004// yet another modification for shift algebra: gives a module with ann = s
     1005// xs+1 is replaced by x*s+s
     1006{
     1007    ring w = 0,(x,s),Dp;
     1008    def W=nc_algebra(1,s);
     1009    setring W;
     1010    matrix m[3][3]=s,0,0,x*s+s,s,0,0,x*s,s;
     1011    list J=jacobson(m,0);
     1012    print(J[1]*m*J[3]);
     1013    print(J[2]);
     1014    print(J[1]);
     1015    print(J[3]);
     1016    print(J[1]*m*J[3]-J[2]);
     1017}
     1018
     1019// variations for above setup with shift algebra:
     1020
     1021// easy but instructive one: see the difference to Weyl case!
     1022matrix m[2][2]=s,x,0,s; print(m);
     1023
     1024// more interesting:
     1025matrix n[3][3]=s,x,0,0,s,x,s,0,x;
     1026
     1027// things blow up
     1028matrix m[3][4] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3, (x+s)^2, x*s;
     1029
     1030// this one is quite nasty and time consuming
     1031matrix m[3][4] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3, (x+s)^2, x*s,x,x^2,x^3,s;
     1032
     1033*/
Note: See TracChangeset for help on using the changeset viewer.