Changeset 6045c3 in git
 Timestamp:
 Dec 22, 2008, 10:45:09 PM (14 years ago)
 Branches:
 (u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')
 Children:
 565e866c54dc28cc8085355e58a8d9ba2b7e0d2c
 Parents:
 3e0c94f2a56b2afdc4e51bebb82b8a047d4487ed
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/jacobson.lib
r3e0c94f r6045c3 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: jacobson.lib,v 1. 5 20081209 16:50:21levandov Exp $";2 version="$Id: jacobson.lib,v 1.6 20081222 21:45:09 levandov Exp $"; 3 3 category="System and Control Theory"; 4 4 info=" … … 7 7 8 8 THEORY: 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 noncommutative, 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 noncommutative, 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 19 REFERENCES: 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 15 23 16 24 PROCEDURES: 17 smith( R,M[,eng1,eng2]); compute the Smith Normal Form of M over commutative ring18 jacobson( R,M[,eng]); compute a weak Jacobson Normal Form of M over noncommutative ring, i.e. a diagonal matrix25 smith(M[,eng1,eng2]); compute the Smith Normal Form of M over commutative ring 26 jacobson(M[,eng]); compute a weak Jacobson Normal Form of M over noncommutative ring 19 27 20 28 SEE ALSO: control_lib 21 29 "; 30 22 31 LIB "poly.lib"; 23 32 LIB "involut.lib"; 24 33 25 proc smith( R,matrix MA, list #)26 "USAGE: smith(R, M[, eng1, eng2]); R ring, M matrix, eng1 and eng2 are optional int27 RETURN: list28 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]_11L[2]_22...L[2]_nn and32 @* 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 basis36 @* If eng2 = 2 than slimgb is used to caculate a Groebner basis37 @*If @code{printlevel}=1, progress debug messages will be printed,34 proc smith(matrix MA, list #) 35 "USAGE: smith(M[, eng1, eng2]); M matrix, eng1 and eng2 are optional integers 36 RETURN: matrix or list 37 ASSUME: The current ring is assumed to be the commutative polynomial ring in 38 one variable 39 PURPOSE: compute the Smith Normal Form of M with transformation matrices (optionally) 40 NOTE: If the optional integer eng1 is nonzero, returns the list {U,D,V}, 41 where 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'. 46 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 38 47 @* if @code{printlevel}>=2, all the debug messages will be printed. 39 48 EXAMPLE: example smith; shows examples 40 49 " 41 50 { 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 42 59 int eng = 0; 43 60 int BASIS; … … 102 119 ring r = 0,x,Dp; 103 120 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 107 124 } 108 125 … … 554 571 } 555 572 556 proc jacobson( R,matrix MA, list #)557 "USAGE: jacobson( R, matrix MA, eng); R ring,M matrix, eng an optional int573 proc jacobson(matrix MA, list #) 574 "USAGE: jacobson(M, eng); M matrix, eng an optional int 558 575 RETURN: 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 basis564 @* If eng = 1 than groebner is used to caculate a Groebner basis565 @* If eng = 2 than slimgb is used to caculate a Groebner basis566 @*If @code{printlevel}=1, progress debug messages will be printed,576 ASSUME: Basering is a noncommutative ring in two variables. 577 PURPOSE: compute a weak Jacobson Normal Form of M over a noncommutative ring 578 NOTE: 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'. 583 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 567 584 @* if @code{printlevel}>=2, all the debug messages will be printed. 568 585 EXAMPLE: example jacobson; shows examples 569 586 " 570 587 { 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 571 607 int B=0; 572 608 if ( size(#)>0 ) … … 615 651 } 616 652 example 617 { "EXAMPLE:"; echo = 2; 653 { 654 "EXAMPLE:"; echo = 2; 618 655 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 } 628 669 629 670 static proc triangle( matrix MA, int B) … … 844 885 /////interesting examples for smith//////////////// 845 886 846 static proc Ex_One_wheeled_bicycle() 887 /* 888 889 //static proc Ex_One_wheeled_bicycle() 847 890 { 848 891 ring RA=(0,m), D, lp; 849 892 matrix bicycle[2][3]=(1+m)*D^2, D^2, 1, D^2, D^21, 0; 850 list s=smith( RA,bicycle, 1,0);893 list s=smith(bicycle, 1,0); 851 894 print(s[2]); 852 895 print(s[1]*bicycle*s[3]s[2]); … … 854 897 855 898 856 static proc Ex_RLCcircuit()899 //static proc Ex_RLCcircuit() 857 900 { 858 901 ring RA=(0,m,R1,R2,L,C), D, lp; 859 902 matrix 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);903 list s=smith(circuit, 1,0); 861 904 print(s[2]); 862 905 print(s[1]*circuit*s[3]s[2]); … … 864 907 865 908 866 static proc Ex_two_pendula()909 //static proc Ex_two_pendula() 867 910 { 868 911 ring RA=(0,m,M,L1,L2,m1,m2,g), D, lp; 869 870 912 matrix pendula[3][4]=m1*L1*D^2,m2*L2*D^2,(M+m1+m2)*D^2,1,m1*L1^2*D^2m1*L1*g,0,m1*L1*D^2,0,0, 871 913 m2*L2^2*D^2m2*L2*g,m2*L2*D^2,0; 872 list s=smith( RA,pendula, 1,0);914 list s=smith(pendula, 1,0); 873 915 print(s[2]); 874 916 print(s[1]*pendula*s[3]s[2]); 875 917 } 876 918 877 878 879 static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit() 919 //static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit() 880 920 { 881 921 ring RA=(0,m,omega,r,a,b), D, lp; 882 883 922 matrix satellite[4][6]= 884 923 D,1,0,0,0,0, … … 886 925 0,0,D,1,0,0, 887 926 0,2*omega/r,0,D,0,b/(m*r); 888 889 list s=smith(RA,satellite, 1,0); 927 list s=smith(satellite, 1,0); 890 928 print(s[2]); 891 929 print(s[1]*satellite*s[3]s[2]); 892 930 } 893 931 894 static proc Ex_flexible_one_link_robot()932 //static proc Ex_flexible_one_link_robot() 895 933 { 896 934 ring RA=(0,M11,M12,M13,M21,M22,M31,M33,K1,K2), D, lp; 897 898 935 matrix 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);936 list s=smith(robot, 1,0); 900 937 print(s[2]); 901 938 print(s[1]*robot*s[3]s[2]); 902 939 } 903 940 941 */ 904 942 905 943 906 944 /////interesting examples for jacobson//////////////// 907 945 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; 910 951 def W=nc_algebra(1,1); 911 952 setring W; 912 953 matrix m[3][3]=[d2,d+1,0],[d+1,0,d3x2*d],[2d+1, d3+d2, d2]; 913 list J=jacobson( W,m,0);954 list J=jacobson(m,0); 914 955 print(J[1]*m*J[3]); 915 956 print(J[2]); … … 919 960 } 920 961 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); 925 966 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^3x^2*s],[2*s+1, s^3+s^2, s^2]; 968 list J=jacobson(m,0); 928 969 print(J[1]*m*J[3]); 929 970 print(J[2]); … … 933 974 } 934 975 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! 1022 matrix m[2][2]=s,x,0,s; print(m); 1023 1024 // more interesting: 1025 matrix n[3][3]=s,x,0,0,s,x,s,0,x; 1026 1027 // things blow up 1028 matrix 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 1031 matrix 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.