Bibliothek für Berechnung des Potentials

load("vect")$

Lie-Klammer bzw. Lie-Ableitung eines Vektorfeldes

LieBracket([l]):=block([f,g,x,k,Df,Dg,ad],
    if (length(l)<3) or (length(l)>4) then 
       error ("Aufruf mit 3 oder 4 Argumenten."),
    f:l[1],
    g:l[2],
    x:l[3],
    if length(l)=3 then k:1 else k:l[4],
    if not(nonnegintegerp(k)) then 
       error ("Ordnung k muss natürliche Zahl sein."),
    if k=0 then return(g)
    else 
        Df:jacobian(f,x),
        Dg:jacobian(g,x),
        ad:list_matrix_entries(Dg.f-Df.g),
        return(LieBracket(f,ad,x,k-1))
)$

Steuerbarkeitsmatrix

ControllabilityMatrix(f,g,x):=block([i,n,L],
    n:length(x),
    if length(f)#n then
        error("Dimensionen von f und x müssen übereinstimmen."),
    if length(g)#n then
        error("Dimensionen von g und x müssen übereinstimmen."),
    L:makelist(LieBracket(-f,g,x,i),i,0,n-1),
    transpose(apply(matrix,L))
    )$

Beispiel: Inveses Pendel mit Trägheitsrad

f:[x2,(m0*sin(x1))/(J1-I2),-(m0*sin(x1))/(J1-I2)];
g:[0,-1/(J1-I2),J1/(I2*(J1-I2))];
x:[x1,x2,x3];

Berechnung der Steuerbarkeitsmatrix

Qs:ControllabilityMatrix(f,g,x);

Orthogonales Komplement zu den ersten zwei Spalten

orthogonal_complement(col(Qs,1),col(Qs,2))$
Ω:factor(%);

Differentialfrom/Kovektorfeld

ω:list_matrix_entries(first(args(Ω)));

Potential/Ausgangsabbildung

h:factor(potential(ω,x));