version="$Id: control.lib,v 1.39 2008-10-07 17:55:52 levandov Exp $"; category="System and Control Theory"; info=" LIBRARY: control.lib Algebraic analysis tools for System and Control Theory AUTHORS: Oleksandr Iena yena@mathematik.uni-kl.de @* Markus Becker mbecker@mathematik.uni-kl.de @* Viktor Levandovskyy levandov@mathematik.uni-kl.de SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis' (Project of Dr. E. Zerz and V. Levandovskyy), Uni Kaiserslautern MAIN PROCEDURES: control(R); analysis of controllability-related properties of R (using Ext modules) controlDim(R); analysis of controllability-related properties of R (using dimension) autonom(R); analysis of autonomy-related properties of R (using Ext modules) autonomDim(R); analysis of autonomy-related properties of R (using dimension) COMPONENT PROCEDURES: leftKernel(R); a left kernel of R rightKernel(R); a right kernel of R leftInverse(R); a left inverse of R rightInverse(R); a right inverse of R colrank(M); a column rank of M as of matrix genericity(M); analysis of the genericity of parameters canonize(L); Groebnerification for modules in the output of control or autonomy procs iostruct(R); computes an IO-structure of behavior given by a module R findTorsion(R, I); generators of the submodule of a module R, annihilated by the ideal I AUXILIARY PROCEDURES: controlExample(s); set up an example from the mini database inside of the library view(); well-formatted output of lists, modules and matrices "; LIB "homolog.lib"; LIB "poly.lib"; LIB "primdec.lib"; LIB "matrix.lib"; //--------------------------------------------------------------- static proc Opt_Our() "USAGE: Opt_Our(); RETURN: intvec, where previous options are stored PURPOSE: save previous options and set customized options " { intvec v; v=option(get); option(redSB); option(redTail); return (v); } //------------------------------------------------------------------------- static proc space(int n) "USAGE:spase(n); n is an integer (number of needed spaces) RETURN: string consisting of n spaces NOTE: the procedure is used in the procedure 'view' to have a better formatted output "{ int i; string s=""; for(i=1;i<=n;i++) { s=s+" "; }; return(s); }; //----------------------------------------------------------------------------- proc view(M) "USAGE: view(M); M is of any type RETURN: no return value PURPOSE: procedure for (well-) formatted output of modules, matrices, lists of modules, matrices; shows everything even if entries are long NOTE: in case of other types( not 'module', 'matrix', 'list') works just as standard 'print' procedure EXAMPLE: example view; shows an example "{ // to be replaced with something more feasible if ( (typeof(M)=="module")||(typeof(M)=="matrix") ) { int @R=nrows(M); int @C=ncols(M); int i; int j; list MaxLength=list(); int Size=0; int max; string s; for(i=1;i<=@C;i++) { max=0; for(j=1;j<=@R;j++) { Size=size( string( M[j,i] ) ); if( Size>max ) { max=Size; }; }; MaxLength[i] = max; }; for(i=1;i<=@R;i++) { s=""; for(j=1;j<@C;j++) { s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) ) +","; }; s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) ); if (i!=@R) { s=s+","; }; print(s); }; return(); }; if(typeof(M)=="list") { int sz=size(M); int i; for(i=1;i<=sz;i++) { view(M[i]); print(""); }; return(); }; print(M); return(); } example {"EXAMPLE:";echo = 2; ring r; list L; matrix M[1][3] = x2+x,y3-y,z5-4z+7; L[1] = "a matrix:"; L[2] = M; L[3] = "an ideal:"; L[4] = ideal(M); view(L); }; //-------------------------------------------------------------------------- proc rightKernel(matrix M) "USAGE: rightKernel(M); M a matrix RETURN: module PURPOSE: computes the right kernel of matrix M (a module of all elements v such that Mv=0) EXAMPLE: example rightKernel; shows an example "{ return(modulo(M,std(0))); } example {"EXAMPLE:";echo = 2; ring r = 0,(x,y,z),dp; matrix M[1][3] = x,y,z; print(M); matrix R = rightKernel(M); print(R); // check: print(M*R); }; //------------------------------------------------------------------------- proc leftKernel(matrix M) "USAGE: leftKernel(M); M a matrix RETURN: module PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0) EXAMPLE: example leftKernel; shows an example " { return( transpose( modulo( transpose(M),std(0) ) ) ); } example {"EXAMPLE:";echo = 2; ring r= 0,(x,y,z),dp; matrix M[3][1] = x,y,z; print(M); matrix L = leftKernel(M); print(L); // check: print(L*M); }; //------------------------------------------------------------------------ proc leftInverse(module M) "USAGE: leftInverse(M); M a module RETURN: module PURPOSE: computes such a matrix L, that LM = Id; EXAMPLE: example leftInverse; shows an example NOTE: exists only in the case when M is free submodule " { // it works also for the NC case; int NCols = ncols(M); module Id = freemodule(NCols); module N = transpose(M); intvec old_opt=Opt_Our(); Id = std(Id); matrix T; // check the correctness (Id \subseteq M) // via dimension: dim (M) = -1! int d = dim_Our(N); if (d != -1) { // No left inverse exists return(matrix(0)); } matrix T2 = lift(N, Id); T2 = transpose(T2); option(set,old_opt); // set the options back return(T2); } example { "EXAMPLE:";echo =2; // a trivial example: ring r = 0,(x,z),dp; matrix M[2][1] = 1,x2z; print(M); print( leftInverse(M) ); kill r; // derived from the example TwoPendula: ring r=(0,m1,m2,M,g,L1,L2),Dt,dp; matrix U[3][1]; U[1,1]=(-L2)*Dt^4+(g)*Dt^2; U[2,1]=(-L1)*Dt^4+(g)*Dt^2; U[3,1]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2); module M = module(U); module L = leftInverse(M); print(L); // check print(L*M); }; //----------------------------------------------------------------------- proc rightInverse(module R) "USAGE: rightInverse(M); M a module RETURN: module PURPOSE: computes such a matrix L, that ML = Id EXAMPLE: example rightInverse; shows an example NOTE: exists only in the case when M is free submodule " { return(transpose(leftInverse(transpose(R)))); } example { "EXAMPLE:";echo =2; // a trivial example: ring r = 0,(x,z),dp; matrix M[1][2] = 1,x2+z; print(M); print( rightInverse(M) ); kill r; // derived from the TwoPendula example: ring r=(0,m1,m2,M,g,L1,L2),Dt,dp; matrix U[1][3]; U[1,1]=(-L2)*Dt^4+(g)*Dt^2; U[1,2]=(-L1)*Dt^4+(g)*Dt^2; U[1,3]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2); module M = module(U); module L = rightInverse(M); print(L); // check print(M*L); }; //----------------------------------------------------------------------- static proc dim_Our(module R) { int d; if (attrib(R,"isSB")<>1) { R = std(R); } d = dim(R); return(d); } //----------------------------------------------------------------------- static proc Ann_Our(module R) { return(Ann(R)); } //----------------------------------------------------------------------- static proc Ext_Our(int i, module R, list #) { // mimicking 'Ext_R' from homolog.lib int ExtraArg = ( size(#)>0 ); if (ExtraArg) { return( Ext_R(i,R,#[1]) ); } else { return( Ext_R(i,R) ); } } //------------------------------------------------------------------------ static proc is_zero_Our { //just a copy of 'is_zero' from "poly.lib" patched with GKdim int d = dim_Our(std(#[1])); int a = ( d==-1 ); if( size(#) >1 ) { list L=a,d; return(L); } return(a); // return( is_zero(R) ) ; }; //------------------------------------------------------------------------ static proc control_output(int i, int NVars, module R, module Ext_1, list Gen) "USAGE: control_output(i, NVars, R, Ext_1), PURPOSE: where @* i is integer (number of first nonzero Ext or a number of variables in a basering + 1 in case that all the Exts are zero), @* NVars: integer, number of variables in a base ring, @* R: module R (cokernel representation), @* Ext_1: module, the first Ext(its cokernel representation) RETURN: list with all the contollability properties of the system which is to be returned in 'control' procedure NOTE: this procedure is used in 'control' procedure "{ // TODO: NVars to be replaced with the global hom. dimension of basering!!! // Is not clear what to do with gl.dim of qrings string DofS = "dimension of the system:"; string Fn = "number of first nonzero Ext:"; string Gen_mes = "Parameter constellations which might lead to a non-controllable system:"; module RK = rightKernel(R); int d=dim_Our(std(transpose(R))); if (i==1) { return( list ( Fn, i, "not controllable , image representation for controllable part:", RK, "kernel representation for controllable part:", leftKernel( RK ), "obstruction to controllability", Ext_1, "annihilator of torsion module (of obstruction to controllability)", Ann_Our(Ext_1), DofS, d ) ); }; if(i>NVars) { return( list( Fn, -1, "strongly controllable(flat), image representation:", RK, "left inverse to image representation:", leftInverse(RK), DofS, d, Gen_mes, Gen) ); }; // //now i<=NVars // if( (i==2) ) { return( list( Fn, i, "controllable, not reflexive, image representation:", RK, DofS, d, Gen_mes, Gen) ); }; if( (i>=3) ) { return( list ( Fn, i, "reflexive, not strongly controllable, image representation:", RK, DofS, d, Gen_mes, Gen) ); }; }; //------------------------------------------------------------------------- proc control(module R) "USAGE: control(R); R a module (R is the matrix of the system of equations to be investigated) RETURN: list PURPOSE: compute the list of all the properties concerning controllability of the system (behavior), represented by the matrix R NOTE: the properties and corresponding data like controllability, flatness, dimension of the system, degree of controllability, kernel and image representations, genericity of parameters, obstructions to controllability, annihilator of torsion submodule and left inverse are investigated EXAMPLE: example control; shows an example " { int i; int NVars=nvars(basering); // TODO: NVars to be replaced with the global hom. dimension of basering!!! int ExtIsZero; intvec v=Opt_Our(); module R_std=std(R); module Ext_1 = std(Ext_Our(1,R_std)); ExtIsZero=is_zero_Our(Ext_1); i=1; while( (ExtIsZero) && (i<=NVars) ) { i++; ExtIsZero = is_zero_Our( Ext_Our(i,R_std) ); }; matrix T=lift(R,R_std); list l=genericity(T); option(set,v); return( control_output( i, NVars, R, Ext_1, l ) ); } example {"EXAMPLE:";echo = 2; // a WindTunnel example ring A = (0,a, omega, zeta, k),(D1, delta),dp; module R; R = [D1+a, -k*a*delta, 0, 0], [0, D1, -1, 0], [0, omega^2, D1+2*zeta*omega, -omega^2]; R=transpose(R); view(R); view(control(R)); }; //-------------------------------------------------------------------------- proc controlDim(module R) "USAGE: controlDim(R); R a module (R is the matrix of the system of equations to be investigated) RETURN: list PURPOSE: computes list of all the properties concerning controllability of the system (behavior), represented by the matrix R EXAMPLE: example controlDim; shows an example NOTE: the properties and corresponding data like controllability, flatness, dimension of the system, degree of controllability, kernel and image representations, genericity of parameters, obstructions to controllability, annihilator of torsion submodule and left inverse are investigated. @* This procedure is analogous to 'control' but uses dimension calculations. @* The implemented approach works for full row rank matrices only (the check is done automatically). " { if( nrows(R) != colrank(transpose(R)) ) { return ("controlDim cannot be applied, since R does not have full row rank"); } intvec v = Opt_Our(); module R_std = std(R); int d = dim_Our(R_std); int NVars = nvars(basering); int i = NVars-d; module Ext_1 = std(Ext_Our(1,R_std)); matrix T = lift(R,R_std); list l = genericity(T); option(set, v); return( control_output( i, NVars, R, Ext_1, l)); } example {"EXAMPLE:";echo = 2; //a WindTunnel example ring A = (0,a, omega, zeta, k),(D1, delta),dp; module R; R = [D1+a, -k*a*delta, 0, 0], [0, D1, -1, 0], [0, omega^2, D1+2*zeta*omega, -omega^2]; R=transpose(R); view(R); view(controlDim(R)); }; //------------------------------------------------------------------------ proc colrank(module M) "USAGE: colrank(M); M a matrix/module RETURN: int PURPOSE: compute the column rank of M as of matrix NOTE: this procedure uses Bareiss algorithm "{ // NOte continued: // which might not terminate in some cases module M_red = bareiss(M)[1]; int NCols_red = ncols(M_red); return (NCols_red); } example {"EXAMPLE: ";echo = 2; // de Rham complex ring r=0,(D(1..3)),dp; module R; R=[0,-D(3),D(2)], [D(3),0,-D(1)], [-D(2),D(1),0]; R=transpose(R); colrank(R); }; //------------------------------------------------------------------------ static proc autonom_output( int i, int NVars, module RC, int R_rank ) "USAGE: proc autonom_output(i, NVars, RC, R_rank) i: integer, number of first nonzero Ext or just number of variables in a base ring + 1 in case that all the Exts are zero NVars: integer, number of variables in a base ring RC: module, kernel-representation of controllable part of the system R_rank: integer, column rank of the representation matrix PURPOSE: compute all the autonomy properties of the system which is to be returned in 'autonom' procedure RETURN: list NOTE: this procedure is used in 'autonom' procedure " { int d=NVars-i;//that is the dimension of the system string DofS="dimension of the system:"; string Fn = "number of first nonzero Ext:"; if(i==0) { return( list( Fn, i, "not autonomous", "kernel representation for controllable part", RC, "column rank of the matrix", R_rank, DofS, d ) ); }; if( i>NVars ) { return( list( Fn, -1, "trivial", DofS, d ) ); }; // //now i<=NVars // if( i==1 ) // in case that NVars==1 there is no sense to consider the notion // of strongly autonomous behavior, because it does not imply // that system is overdetermined in this case { return( list ( Fn, i, "autonomous, not overdetermined", DofS, d ) ); }; if( i==NVars ) { return( list( Fn, i, "strongly autonomous(fin. dimensional),in particular overdetermined", DofS, d) ); }; if( iv[i+1]) { temp=v[i]; v[i]=v[i+1]; v[i+1]=temp; b=1; } } } P=R[v]; //generate P for(i=1;i<=NCols;i++) //generate Q { if(elementof(i,v)==1) { i++; continue; } Q=Q,R[i]; } Q=simplify(Q,2); string s="The following components have been chosen as outputs: "; return (list(s,v,P,Q)); } example {"EXAMPLE:";echo = 2; //Example Antenna ring r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), (c,dp); module RR; RR = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0], [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta], [0, 0, Dt, -K1, 0, 0, 0, 0, 0], [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta], [0, 0, 0, 0, Dt, -K1, 0, 0, 0], [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta]; module R = transpose(RR); view(iostruct(R)); }; //--------------------------------------------------------------- static proc smdeg(matrix N) "USAGE: smdeg( N ); N a matrix RETURN: intvec PURPOSE: returns an intvec of length 2 with the index of an element of N with smallest degree " { int n = nrows(N); int m = ncols(N); int d,d_temp; intvec v; int i,j; // counter if (N==0) { v = 1,1; return(v); } for (i=1; i<=n; i++) // hier wird ein Element ausgewaehlt(!=0) und mit dessen Grad gestartet { for (j=1; j<=m; j++) { if( deg(N[i,j])!=-1 ) { d=deg(N[i,j]); break; } } if (d != -1) { break; } } for(i=1; i<=n; i++) { for(j=1; j<=m; j++) { d_temp = deg(N[i,j]); if ( (d_temp < d) && (N[i,j]!=0) ) { d=d_temp; } } } for (i=1; i<=n; i++) { for (j=1; j<=m;j++) { if ( (deg(N[i,j]) == d) && (N[i,j]!=0) ) { v = i,j; return(v); } } } } //--------------------------------------------------------------- static proc NoNon0Pol(vector v) "USAGE: NoNon0Pol(v), v a vector RETURN: int PURPOSE: returns 1, if there is only one non-zero element in v and 0 else "{ int i,j; int n = nrows(v); for( j=1; j<=n; j++) { if (v[j] != 0) { i++; } } if ( i!=1 ) { i=0; } return(i); } //--------------------------------------------------------------- static proc extgcd_Our(poly p, poly q) { ideal J; //for extgcd-computations matrix T; //----------"------------ list L; // the extgcd-command has a bug in versions before 2-0-7 if ( system("version")<=2006 ) { J = p,q; // J = N[k-1,k-1],N[k,k]; //J is of type ideal L[1] = liftstd(J,T); //T is of type matrix if(J[1]==p) //this is just for the case the SINGULAR swaps the // two elements due to ordering { L[2] = T[1,1]; L[3] = T[2,1]; } else { L[2] = T[2,1]; L[3] = T[1,1]; } } else { L=extgcd(p,q); // L=extgcd(N[k-1,k-1],N[k,k]); //one can use this line if extgcd-bug is fixed } return(L); } static proc normalize_Our(matrix N, matrix Q) "USAGE: normalize_Our(N,Q), N, Q are two matrices PURPOSE: normalizes N and divides the columns of Q through the leading coefficients of the columns of N RETURN: normalized matrix N and altered Q(according to the scheme mentioned in purpose). If number of columns of N and Q do not coincide, N and Q are returned unchanged NOTE: number of columns of N and Q must coincide. " { if(ncols(N) != ncols(Q)) { return (N,Q); } module M = module(N); module S = module(Q); int NCols = ncols(N); number n; for(int i=1;i<=NCols;i++) { n = leadcoef(M[i]); if( n != 0 ) { M[i]=M[i]/n; S[i]=S[i]/n; } } N = matrix(M); Q = matrix(S); return (N,Q); } //--------------------------------------------------------------- proc oldsmith( module M ) "USAGE: oldsmith(M); M a module/matrix PURPOSE: computes the Smith normal form of a matrix RETURN: a list of length 4 with the following entries: @* [1]: the Smith normal form S of M, @* [2]: the rank of M, @* [3]: a unimodular matrix U, @* [4]: a unimodular matrix V, such that U*M*V=S. An warning is returned when no Smith form exists. NOTE: Older experimental implementation. The Smith form only exists over PIDs (principal ideal domains). Use global ordering for computations! " { if (nvars(basering)>1) //if more than one variable, return empty list { string s="The Smith-Form only exists for principal ideal domains"; return (s); } matrix N = matrix(M); //Typecasting int n = nrows(N); int m = ncols(N); matrix P = unitmat(n); //left transformation matrix matrix Q = unitmat(m); //right transformation matrix int k, i, j, deg_temp; poly tmp; vector v; list L; //for extgcd-computation intmat f[1][n]; //to save degrees matrix lambda[1][n]; //to save leadcoefficients intmat g[1][m]; //to save degrees matrix mu[1][m]; //to save leadcoefficients int ii; //counter while ((k!=n) && (k!=m) ) { k++; while ((k<=n) && (k<=m)) //outer while-loop for column-operations { while(k<=m ) //inner while-loop for row-operations { if( (n>m) && (k < n) && (kn) { break; } if(NoNon0Pol(transpose(N)[k])==1) { break; } tmp=leadcoef(N[k,k]); deg_temp=ord(N[k,k]); //ord outputs the leading degree of N[k][k] for(ii=k+1;ii<=m;ii++) { mu[1,ii]=leadcoef(N[k,ii])/tmp; g[1,ii]=deg(N[k,ii])-deg_temp; } for(ii=k+1;ii<=m;ii++) { N=addcol(N,k,-mu[1,ii]*var(1)^g[1,ii],ii); Q=addcol(Q,k,-mu[1,ii]*var(1)^g[1,ii],ii); N,Q=normalize_Our(N,Q); } } if( (k!=1) && (knur wenn string in liste dann verbatim { t=L[i]; if(nr_loop==1) { write(l,"\\begin\{center\}"); write(l,"\\begin\{verbatim\}"); } write(l,t); if(nr_loop==0) { write(l,"\\par"); } if(nr_loop==1) { write(l,"\\end\{verbatim\}"); write(l,"\\end\{center\}"); } break; } if(typeof(L[i])=="module") { texobj(name,matrix(L[i])); break; } if(typeof(L[i])=="list") { list_tex(L[i],name,l,1); break; } write(l,"\\begin\{center\}"); texobj(name,L[i]); write(l,"\\end\{center\}"); write(l,"\\par"); break; } } } } example { "EXAMPLE:";echo = 2; } //--------------------------------------------------------------- proc verbatim_tex(string s, link l) "USAGE: verbatim_tex(s,l), where s is a string and l a link PURPOSE: writes the content of s in verbatim-environment in the file specified by link RETURN: nothing " { write(l,"\\begin{verbatim}"); write(l,s); write(l,"\\end{verbatim}"); write(l,"\\par"); } example { "EXAMPLE:";echo = 2; } //--------------------------------------------------------------- proc findTorsion(module R, ideal TAnn) "USAGE: findTorsion(R, I); R an ideal/matrix/module, I an ideal RETURN: module PURPOSE: computes the Groebner basis of the submodule of R, annihilated by I NOTE: especially helpful, when I is the annihilator of the t(R) - the torsion submodule of R. In this case, the result is the explicit presentation of t(R) as the submodule of R EXAMPLE: example findTorsion; shows an example " { // motivation: let R be a module, // TAnn is the annihilator of t(R)\subset R // compute the generators of t(R) explicitly ideal AS = TAnn; module S = R; if (attrib(S,"isSB")<>1) { S = std(S); } if (attrib(AS,"isSB")<>1) { AS = std(AS); } int nc = ncols(S); module To = quotient(S,AS); To = std(NF(To,S)); return(To); } example { "EXAMPLE:";echo = 2; // Flexible Rod ring A = 0,(D1, D2), (c,dp); module R= [D1, -D1*D2, -1], [2*D1*D2, -D1-D1*D2^2, 0]; module RR = transpose(R); list L = control(RR); // here, we have the annihilator: ideal LAnn = D1; // = L[10] module Tr = findTorsion(RR,LAnn); print(RR); // the module itself print(Tr); // generators of the torsion submodule } proc controlExample(string s) "USAGE: controlExample(s); s a string RETURN: ring PURPOSE: set up an example from the mini database by initalizing a ring and a module in a ring NOTE: in order to see the list of available examples, execute @code{controlExample(\"show\");} @* To use ab example, one has to do the following. Suppose one calls the ring, where the example will be activated, A. Then, by executing @* @code{def A = controlExample(\"Antenna\");} and @code{setring A;}, @* A will become a basering from the example \"Antenna\" with the predefined system module R (transposed). After that one can just execute @code{control(R);} respectively @code{autonom(R);} to perform the control resp. autonomy analysis of R. EXAMPLE: example controlExample; shows an example "{ list E, S, D; // E=official name, S=synonym, D=description E[1] = "Cauchy1"; S[1] = "cauchy1"; D[1] = "1-dimensional Cauchy equation"; E[2] = "Cauchy2"; S[2] = "cauchy2"; D[2] = "2-dimensional Cauchy equation"; E[3] = "Control1"; S[3] = "control1"; D[3] = "example of a simple noncontrollable system"; E[4] = "Control2"; S[4] = "control2"; D[4] = "example of a simple controllable system"; E[5] = "Antenna"; S[5] = "antenna"; D[5] = "antenna"; E[6] = "Einstein"; S[6] = "einstein"; D[6] = "Einstein equations in vacuum"; E[7] = "FlexibleRod"; S[7] = "flexible rod"; D[7] = "flexible rod"; E[8] = "TwoPendula"; S[8] = "two pendula"; D[8] = "two pendula mounted on a cart"; E[9] = "WindTunnel"; S[9] = "wind tunnel";D[9] = "wind tunnel"; E[10] = "Zerz1"; S[10] = "zerz1"; D[10] = "example from the lecture of Eva Zerz"; // all the examples so far int i; if ( (s=="show") || (s=="Show") ) { print("The list of examples:"); for (i=1; i<=size(E); i++) { printf("name: %s, desc: %s", E[i],D[i]); } return(); } string t; for (i=1; i<=size(E); i++) { if ( (s==E[i]) || (s==S[i]) ) { t = "def @A = ex"+E[i]+"();"; execute(t); return(@A); } } "No example found"; return(); } example { "EXAMPLE:";echo = 2; controlExample("show"); // let us see all available examples: def B = controlExample("TwoPendula"); // let us set up a particular example setring B; print(R); } //---------------------------------------------------------- // //Some example rings with defined systems //---------------------------------------------------------- //autonomy: // //---------------------------------------------------------- static proc exCauchy1() { ring @r=0,(s1,s2),dp; module R= [s1,-s2], [s2, s1]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- static proc exCauchy2() { ring @r=0,(s1,s2,s3,s4),dp; module R= [s1,-s2], [s2, s1], [s3,-s4], [s4, s3]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- static proc exZerz1() { ring @r=0,(d1,d2),dp; module R=[d1^2-d2], [d2^2-1]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- //control //---------------------------------------------------------- static proc exControl1() { ring @r=0,(s1,s2,s3),dp; module R=[0,-s3,s2], [s3,0,-s1]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- static proc exControl2() { ring @r=0,(s1,s2,s3),dp; module R=[0,-s3,s2], [s3,0,-s1], [-s2,s1,0]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- static proc exAntenna() { ring @r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), dp; module R; R = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0], [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta], [0, 0, Dt, -K1, 0, 0, 0, 0, 0], [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta], [0, 0, 0, 0, Dt, -K1, 0, 0, 0], [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- static proc exEinstein() { ring @r = 0,(D(1..4)),dp; module R = [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)], [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], [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], [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)], [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)], [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], [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)], [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)], [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)], [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]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- static proc exFlexibleRod() { ring @r = 0,(D1, delta), dp; module R; R = [D1, -D1*delta, -1], [2*D1*delta, -D1-D1*delta^2, 0]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- static proc exTwoPendula() { ring @r=(0,m1,m2,M,g,L1,L2),Dt,dp; module R = [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2], [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2], [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2]; R=transpose(R); export R; return(@r); } //---------------------------------------------------------- static proc exWindTunnel() { ring @r = (0,a, omega, zeta, k),(D1, delta),dp; module R = [D1+a, -k*a*delta, 0, 0], [0, D1, -1, 0], [0, omega^2, D1+2*zeta*omega, -omega^2]; R=transpose(R); export R; return(@r); }