Changeset 0ceca2 in git
 Timestamp:
 Dec 22, 2004, 7:03:16 PM (20 years ago)
 Branches:
 (u'fiekerDuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
 Children:
 40d02c5bf05d362fee584c360803b8c93c23f239
 Parents:
 6fe0f6a88ab2e438b1590f5e383517ed59728be1
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/control.lib
r6fe0f6 r0ceca2 1 version="$Id: control.lib,v 1.1 8 20041216 16:28:53levandov Exp $";1 version="$Id: control.lib,v 1.19 20041222 18:03:16 levandov Exp $"; 2 2 category="Applications"; 3 3 info=" … … 23 23 24 24 AUXILIARY PROCEDURES: 25 ncdetection(ring r); computes an ideal, presenting an involution map on noncomm algebra r;26 involution(m, map theta); applies the involution, presented by theta to m of type poly, vector, ideal, module;27 25 declare(string NameOfRing, Variables[,string Parameters, Ordering]); defines the ring, optional parametes are strings of parameters and ordering, 28 26 view(); Wellformatted output of lists, modules and matrixes … … 55 53 // NOTE: static things should not be shown for enduser 56 54 // 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'; 58 56 // static space(int n) Procedure used inside the procedure 'Print' to have a better formatted output 59 57 // static control_output(); Generating the output for the procedure 'control' … … 139 137 // print(" To compute controllability please enter: control(R)"); 140 138 // print(" To compute autonomy please enter: autonom(R)"); 141 //142 //143 //144 139 // 140 145 141 static proc space(int n) 146 142 "USAGE:spase(n); … … 239 235 // 240 236 proc RightKernel(matrix M) 241 "USAGE: RightKernel(M); 242 M: matrix 237 "USAGE: RightKernel(M); M a matrix 243 238 RETURN: 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 module239 NOTE: in the noncommutative case RightKernel is a right module 245 240 EXAMPLE: example RightKernel; shows an example 246 241 " … … 250 245 example 251 246 {"EXAMPLE:";echo = 2; 252 ring r ;247 ring r = 0,(x,y,z),dp; 253 248 matrix M[1][3] = x,y,z; 254 249 print(M); 255 print( RightKernel(M) ); 250 matrix R = RightKernel(M); 251 print(R); 252 // check: 253 print(M*R); 256 254 }; 257 255 // 258 256 proc 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 258 RETURN: left kernel of matrix M, i.e. the module of all elements v such that vM=0 263 259 EXAMPLE: example LeftKernel; shows an example 264 260 " … … 268 264 example 269 265 {"EXAMPLE:";echo = 2; 270 ring r ;266 ring r= 0,(x,y,z),dp; 271 267 matrix M[3][1] = x,y,z; 272 268 print(M); 273 print( LeftKernel(M) ); 269 matrix L = LeftKernel(M); 270 print(L); 271 // check: 272 print(L*M); 274 273 }; 275 274 // … … 277 276 "USAGE: LeftInverse(M); M a module 278 277 RETURN: matrix L such that LM == Id; 279 EXAMPLE: 278 EXAMPLE: example LeftInverse; shows an example 280 279 NOTE: exists only in the case when Id \subseteq M! 281 280 " 282 281 { 283 282 // 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); 285 286 module Id = freemodule(NCols); 286 M= transpose(M);287 module N = transpose(M); 287 288 Id = std(Id); 288 289 matrix T; 289 290 option(redSB); 290 291 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"; 300 298 return(matrix(0)); 301 299 } 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 // } 302 308 matrix T2 = lift(N, Id); 303 T2 = transpose(T2)*transpose(T);309 // T2 = transpose(T2)*transpose(T); 304 310 T2 = transpose(T2); 311 // set options back 312 option(noredTail); 305 313 return(T2); 306 314 } 307 315 example 308 {"EXAMPLE:";echo =2; 309 ring r; 316 { // trivial example: 317 "EXAMPLE:";echo =2; 318 ring r = 0,(x,z),dp; 310 319 matrix M[2][1] = 1,x2z; 311 320 print(M); 312 321 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*L1g*L2)*Dt^2+(g^2); 329 module M = module(U); 330 module L = LeftInverse(M); 331 print(L); 313 332 }; 314 333 // … … 316 335 { 317 336 int d; 337 if (attrib(R,"isSB")<>1) 338 { 339 R = std(R); 340 } 318 341 if (Is_NC()) 319 342 { … … 328 351 if (Is_NC()) 329 352 { 353 // TODO: use preAnn, that is an annihilator of the generators 354 // what is a left ideal! 330 355 return(1); 331 356 } … … 338 363 // 339 364 static proc Ext_Our(int i, module R, list #) 340 // just a copy of'Ext_R' from "homolog.lib" in commutative case365 // mimicking 'Ext_R' from "homolog.lib" in commutative case 341 366 { 342 367 int ISNC = Is_NC(); … … 369 394 // 370 395 static 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 372 397 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), 376 401 Ext_1: module, the first Ext(its cokernel representation) 377 402 RETURN: list with all the contollability properties of the system which is to be returned in 'control' procedure … … 379 404 " 380 405 { 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) 385 411 { 386 module RK = RightKernel(R); 412 module RK = RightKernel(R); 387 413 return( 388 414 list ( Fn, … … 403 429 404 430 if(i>NVars) 405 { module RK = RightKernel(R);431 { module RK = RightKernel(R); 406 432 return( list( Fn, 407 433 1, … … 439 465 d) 440 466 ); 441 }; 442 443 467 }; 444 468 }; 445 469 // 446 470 447 471 proc 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) 473 RETURN: list of all the properties concerning controllability of the system (behavior) represented by the matrix R 451 474 EXAMPLE: example control; shows an example 452 475 " … … 454 477 int i; 455 478 int NVars=nvars(basering); 479 // TODO: NVars to be replaced with the global hom. dimension of basering!!! 456 480 int ExtIsZero; 457 458 481 459 482 module Ext_1 = std(Ext_Our(1,R)); 460 483 … … 480 503 view(R); 481 504 view(control(R)); 482 483 505 }; 484 506 // … … 521 543 522 544 if( i==1 ) 523 // in case that NVars==1 there is no sence to consider the notion524 // of strongly autonomous behavior, because it does not imply525 // that system is overdetermined in this case545 // 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 526 548 { 527 549 return( list ( Fn, … … 552 574 ); 553 575 }; 554 555 576 }; 556 577 // 557 578 proc 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) 580 RETURN: list of all the properties concerning autonomy of the system (behavior) represented by the matrix R 561 581 NOTE: this procedure is an analogue to 'autonom' using dimension calculations 562 582 EXAMPLE: example autonom2; shows an example … … 581 601 view( autonom2(R) ); 582 602 }; 583 // 584 603 // 585 604 proc 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) 606 RETURN: list of all the properties concerning autonomy of the system (behavior) represented by the matrix R 589 607 EXAMPLE: example autonom; shows an example 590 608 " … … 600 618 ExtIsZero = is_zero_Our(Ext_Our(i,R)); 601 619 }; 602 603 620 return(autonom_output(i,NVars)); 604 621 } … … 616 633 }; 617 634 618 // 635 // 619 636 // 620 637 //Some example rings with defined systems 621 // 638 // 622 639 //autonomy: 623 640 // 624 // 641 // 625 642 proc exCauchy() 626 643 { … … 632 649 return(@r); 633 650 }; 634 // 651 // 635 652 proc exCauchy2() 636 653 { … … 644 661 return(@r); 645 662 }; 646 // 663 // 647 664 proc exZerz() 648 665 { … … 654 671 return(@r); 655 672 }; 656 // 673 // 657 674 //control 658 675 // … … 666 683 return(@r); 667 684 }; 668 // 685 // 669 686 proc ex2() 670 687 { … … 677 694 return(@r); 678 695 }; 679 // 696 // 680 697 proc exAntenna() 681 698 { … … 694 711 }; 695 712 696 // 713 // 697 714 698 715 proc exEinstein() … … 717 734 718 735 719 // 736 // 720 737 721 738 proc exFlexibleRod() … … 730 747 }; 731 748 732 // 749 // 733 750 proc exTwoPendula() 734 751 { … … 742 759 return(@r); 743 760 }; 744 // 761 // 745 762 proc exWindTunnel() 746 763 { … … 754 771 return(@r); 755 772 }; 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); // WeylAlgebra 845 map F = r,x,d; 846 poly f = x*d^2+d; 847 poly If = involution(f,F); 848 fIf; 849 poly g = x^2*d+2*x*d+3*x+7*d; 850 poly tg = d*x^22*d*x+3*x7*d; 851 poly Ig = involution(g,F); 852 tgIg; 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 iten zeile die indices 867 // der differential,shift oder advanceoperatoren enthaelt mit denen die ite 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( (dd_lead != 0)  (LExp > 1) ) 907 if ( ( (dd_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 // 1007 774 1008 775 proc genericity(matrix M) 1009 776 "USAGE: genericity(M), M is a matrix 1010 RETURN: list of strings with expressions, used as nonzero in SB1011 NOTE: effective with the liftstd procedure 777 RETURN: list of strings with expressions, which have been assumed to be nonzero in the process of computing the Groebner basis 778 NOTE: effective with the liftstd procedure for modules with parameters 1012 779 " 1013 780 { … … 1085 852 module SR = liftstd(R,T); 1086 853 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); 1087 860 } 1088 861 1089 862 proc canonize(list L) 1090 863 "USAGE: canonize(L), L a list 1091 ASSUME: L is the output of control/autonomy procs 1092 RETURN: canonized list 864 ASSUME: L is the output of control/autonomy procs 865 RETURN: a canonized list, where modules are canonized by computing 866 their reduced minimal Groebner bases 1093 867 " 1094 868 { … … 1103 877 { 1104 878 M[i] = std(M[i]); 1105 // M[i] = prune(M[i]); // mimimal embedding 879 // M[i] = prune(M[i]); // mimimal embedding: no need yet 1106 880 // M[i] = std(M[i]); 1107 881 } 1108 882 } 883 option(noredTail); // set the initial option back 1109 884 return(M); 1110 885 } … … 1328 1103 option(redTail); 1329 1104 ring r=0,x,dp; 1330 // This isexample is trivial, but we want to see,1105 // This example is trivial, but we want to see, 1331 1106 // that nothing is changed, when the matrix is already in SmithForm! 1332 1107 module M = [x,0,0],[0,x2,0],[0,0,x3];
Note: See TracChangeset
for help on using the changeset viewer.