# Lektion 5
# Vektoren und Matrizen
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): unprotect(D,I): 
# Für Vektoren sind folgende Rechenregeln definiert:
# a+b = a*kommutativ+b;
# a+b+c = a+(b+c)*assoziativ;
# a+0 = a*Nullvektor;
# a-a = 0*Negativelement;
# alpha(a+b) = alpha*a+Typesetting[delayDotProduct](alpha, b, true)*distributiv;
# a = a*Einselement;
# Vektorräume, die ein Skalarprodukt oder inneres Produkt besitzen, werden Euklidische Vektorräume genannt, in denen folgende Rechenregeln gelten:
# a*b = b*a*kommutativ;
# alpha*a*b = alpha(a*b)*assoziativ;
# (a+b)*c = b*c*distributiv+a*c;
# a*a > (0*für)*a and (0*für)*a <> 0*(positiv-definit);
# Zur Eingabe von Vektoren  und  Matrizen  stehen mehrere Möglichkeiten zur Verfügung, u.a. mittels der  Palette .   
# Hinweis: Falls nichts anderes vereinbart, werden in Maple Vektoren immer als Spaltenvektoren interpretiert.
> A:=Vector(3,symbol=a): B:=Vector(3,symbol=b):
> C:=Vector(3,symbol=c): D:=Vector(3,symbol=d):   
> print(A,B,C,D); 
# Addition und Subtraktion a - b + c
> A-B+C; 
# Multiplikation eines Vektors mit einem Skalar ±
> alpha*B;
# Das geht auch mit dem Befehl ScalarMultiply
> ScalarMultiply(B,alpha);
# Zahlenmäßig kann man einen Vektor definieren durch
> ZV:=Vector([4,7,9]);
# Den zum Spaltenvektor  a   transponierten Vektor   aT (Zeilenvektor) erhalten wir mit dem Befehl 
> AT:=Transpose(A);
# oder auch kürzer
> AT:=A^+;  AT:=A^%T;
# Das Skalarprodukt a*b;  der  beiden Vektoren  a  und  b   berechnet sich als Produkt des Zeilenvektors a  mit dem Spaltenvektor  b. Das Ergebnis ist ein Skalar
> Transpose(A).B;
# Dasselbe Ergebnis liefert der Befehl DotProduct mit dem Zusatz conjugate = false
> DotProduct(A,B,conjugate = false);
# Für das Skalarprodukt gilt das kommutative Gesetz a*b = a*b;
> DotProduct(A,B,conjugate=false)-DotProduct(B,A,conjugate=false);
# Hinweis: Um im Folgenden Schreibarbeit zu sparen, wird  der Ausdruck conjugate = false der Variablen cf zugewiesen.
> cf:=conjugate = false:
# Der folgende Befehl, ohne den Zusatz  cf , unterstellt, dass die zum Produkt kommenden Vektoren komplex sind. Der übergesetzte Querstrich bezeichnet die konjugiert komplexen Koordianten von a  
> DotProduct(A,B);
# Das Skalarprodukt des Vektors  a  mit sich selbst ist
> DotProduct(A,A);DotProduct(A,A,cf);
# Die Wurzel aus dem obigen Ausdruck ist der Betrag des Vektors A (oder auch seine Länge). Er wird auch als 2-Norm bezeichnet, die in Maple mit der Euklidischen- oder auch Frobenius-Norm identisch ist
> Norm(A,Euclidean); Norm(A,Frobenius);
# Dagegen bedeutet  abs(A)  in Maple etwas anderes:
> abs(A);
# Das dyadische Produkt a—b zweier Vektoren a und b  ist ein Spezialfall eines Tensors zweiter Stufe. Die Komponentenmatrix ergibt sich als Matrizenprodukt der Spaltenmatrix von  a  mit der Zeilenmatrix von  b, was oft auch in der Form ab^T;geschrieben wird.
> DYAB:=A.Transpose(B); 
> Rank(DYAB),Determinant(DYAB);
# Weil der Rang eines dyadischen Produkts gleich 1 und damit kleiner als 3  ist, verschwindet dessen  Determinante.
# Das Vektorprodukt zweier Vektoren  a und b  ist ein Vektor, dessen Berechnung veranlasst wird durch den  Befehl:
> V:=CrossProduct(A,B);
# oder auch kürzer
> V:=A &x B;
# Dieser Vektor v steht senkrecht auf der Ebene, die durch a und b aufgespannt wird, was man wie folgt prüfen kann:
> DotProduct(A,V,cf);  simplify(%);
> DotProduct(B,V,cf);  simplify(%);
# Hinweis: Das Vektorprodukt ist nur im 3  definiert.
# Das Vektorprodukt ist nicht kommutativ sondern alterniert, und es gilt a × b = - b × a 
> CrossProduct(A,B) + CrossProduct(B,A);
# Wir bilden das zweifache Vektorprodukt a×(b×b) und  beweisen damit die Rechenregel a×(b×c) = b (a·c) - c (a·b)
> AKBKC:=CrossProduct(A,CrossProduct(B,C));
> AC:=DotProduct(A,C,cf); AB:=DotProduct(A,B,cf);
> AKBKC-(AC*B-AB*C); simplify(%);
# Die Division des Vektors a durch seinen Betrag liefert seinen Einheitsvektor mit dem Betrag 1. Das leistet der Befehl Normalize
> Normalize(A,Euclidean,cf);
# oder auch
> Normalize(A,Frobenius,cf);
# Des Weitern enthält Maple die sog. Unendlichkeitsnorm. Dabei werden die Koordinaten durch das betragsgrößte Element dividiert. Achtung: Der so entstandene Vektor hat nicht den Betrag 1. 
> Normalize(A,infinity,cf);
# Der Flächeninhalt FL des durch die beiden Vektoren a und b aufgespannten Paralleleogramms wird durch den Betrag des Vektorproduktes beschrieben, also FL =|a×b| und die Orientierung der Fläche im Raum ist im Sinne der Rechtsschraubenregel durch den Einheitsvektor von a×b festgelegt. 
> FL:=Norm(V,Euclidean);
> OR:=Normalize(V,Euclidean,cf);
# Nun bilden wir das Spatprodukt [a, b, c] := (`&x`(a, b);)·c der drei Vektoren  a, b, c  als Skalarprodukt der Vektoren `&x`(a, b);  und  c:
> sp1:=DotProduct(V,C,cf); expand(%);
# Das Spatprodukt entspricht dem Volumen des von den drei Vektoren a,b,c aufgespannten Parallelepipeds (Spat).
# Wir definieren nun eine quadratische Matrix, die spaltenweise die Komponenten der Vektoren  a, b  und c enthält:
> Dsp1:=Matrix([A,B,C]);
# Ein Vergleich zeigt: Der Wert der  Determinante  von  Dsp1  ist identisch mit dem Spatprodukt:
> sp2:=Determinant(Dsp1);  
> simplify(sp2-sp1);
# Den Winkel ± , den die Vektoren a und b  miteinander einschließen, erhalten wir mit dem Befehl 
> alpha:=VectorAngle(A,B,cf);
# Maple gibt diesen Wert im Bogenmaß aus. Soll der Winkel in Gradmaß umgerechnet werden, dann gelingt das mit dem Befehl convert.
> convert(alpha,degrees);
# Im Folgenden werden einige Zahlenbeispiele vorgestellt. Da wir mit reellen Komponenten rechnen, kann bis auf Weiteres auf den Parameter conjugate = false verzichtet werden.
> A:=Vector([ 3, 5,-1]):  B:=Vector([1, 2, 2]):
> C:=Vector([ 1, 1, 6]):  D:=Vector([4, 3, 6]):
> print(A,B,C,D);
> A - B + C;
# Der Vektor a wird mit dem Skalar  3 multipliziert 
> 3*A;
# Das Skalarprodukt der beiden Vektoren a und b ist 
> AB:=DotProduct(A,B);
# Hier gilt das kommutative Gesetz a·b = b·a
> BA:=DotProduct(B,A);
# Die Beträge (Längen) der Vektoren sind
> NA:=Norm(A,Euclidean);NB:=Norm(B,Euclidean);
> NC:=Norm(C,Euclidean);ND:=Norm(D,Euclidean);
# Dagegen ist
> abs(A);
# Zur grafischen Darstellung dreidimensionaler Vektoren als Vektorpfeile wird der Befehl  arrow  aus der  Bibliothek plots benötigt.
> par:=shape=cylindrical_arrow,width=0.2, head_length=0.8:
> pa:=plots[arrow](A,par,color =   red):
> pb:=plots[arrow](B,par,color = green):
> pc:=plots[arrow](C,par,color =  blue):
> plots[display]([pa,pb,pc],axes=boxed,orientation=[-50,60,0],scaling=constrained);
# Das dyadische Produkt der Vektoren a und b errechnet sich zu
> DYAB:=A.Transpose(B);
# Das Vektorprodukt der beiden Vektoren a und b ist
> AKB:=CrossProduct(A,B);
# Wir berechnen den Einheitsvektor von a 
> evalf(Normalize(A,Euclidean));
# Wir prüfen seine Länge
> Norm(%,Euclidean);
# Der Flächeninhalt des von den beiden Vektoren a und b aufgespannten Parallelogramms ermitteln wir wie folgt:
> FL:=evalf(Norm(AKB,Euclidean));
# Der im Sinne der Rechtsschraubenregel senkrecht auf der von a und b aufgespannten Ebene stehende Einheitsvektor ist
> OR:=evalf(Normalize(AKB,Euclidean,cf));
# Der Vektor or steht senkrecht auf der durch die Vektoren  a und b aufgespannten Ebene, wie die folgenden Rechnungen zeigen:
> DotProduct(OR,A,cf); DotProduct(OR,B,cf);
# Wir berechnen das Spatprodukt [a, b, c];, das dem Volumen des von a, b und c aufgespannten Parallelepipeds (Spat) entspricht.
> sp1:=Determinant(Matrix([A,B,C]));
# Wir berechnen den Winkel zwischen den Vektoren  a und b  im Bogen- sowie im Gradmaß:
> alpha:=evalf(VectorAngle(A,B));
> evalf(convert(alpha,degrees));
# Die Defintion des Skalarproduktes der beiden Vektoren a und b ist a*b = abs(a || b)*cos(alpha);. Wir überprüfen das 
> evalf(NA*NB*cos(alpha)-AB);
# Für den  Betrag des Vektors  a×b  gilt |a×b| = |a||b| sin±. 
> simplify(NA*NB*sin(alpha)-FL);
# Einfache geometrische Anwendungen linearer Gebilde sind die Gerade und die Ebene. Die Parameterdarstellung der Geraden, die durch die Punkte r0  und  r1 verläuft, ist r = r0 + t (r1 - r0)  und mit dem Richtungsvektor  g = r1 - r0  folgt die Punkt-Richtungs-Form r = r0 + t g  der Geradengleichung. Jeder t-Wert bedeutet genau einen Punkt auf der Geraden. Um beispielsweise denjenigen Punkt r[n] = g*t[n]+r[0];  zu ermitteln, der dem Nullpunkt am nächsten ist, muss  r[n]; der Bedingung r[n] . g = 0;  genügen. Das ist gleichbedeutend mit  t[n]=-(r[0]*g)/(g^(2)). ;
> r0:=Vector([1,1,1]): r1:=Vector([2,3,4]): g:=r1-r0:
> print(r0,r1,g);
> r:=unapply(r0 + t*g,t);
> tn:=-DotProduct(r0,g)/DotProduct(g,g);
> rn:=r(tn);
# Der Abstand vom Nullpunkt ist dann
> evalf(Norm(rn,Euclidean));
# Die Komponente eines Vektors b in Richtung eines zweiten Vektors a ist der Vektor b[a]=(a*b)/(a^(2))a. ;Wir rechnen dazu ein Zahlenbeispiel:
> A :=Vector([1,1,1]):  B :=Vector([1,2,3]):
> AB:=DotProduct(A,B):  AQ:=DotProduct(A,A):
> KBA:=AB/AQ*A;
# Wir stellen das Ergebnis grafisch dar und sehen, dass der Vektor kba (blau)  in die Richtung des Vektors a (rot)  fällt:
> r0:=Vector([0,0,0]):
> par:= shape=cylindrical_arrow,width=0.1, head_length=0.8:
> pa:=plots[arrow](r0,A  ,par,color =   red):
> pb:=plots[arrow](r0,B  ,par,color = green):
> pc:=plots[arrow](r0,KBA,par,color =  blue):
> plots[display]([pa,pb,pc],axes=boxed,orientation=[-50,60,0],scaling=constrained);
# Die Berechnung des  Abstandes  d  eines Punktes Q von einer Geraden  r = g*t+r[0];erfolgt in zwei Schritten. Im ersten Schritt wird derjenige Punkt Q' mit dem Ortsvektor  r[diff(Q(x), x)] = r[0]+(diff(t(x), x))*g;  auf der Geraden ermittelt, für den der Vektor  r[diff(Q(x), x)]-r[Q(x)];  = r0  - rQ  +  t' g  senkrecht auf der Geraden mit dem Richtungsvektor g steht. Das erfordert (r0  - rQ  +  t' g) · g = 0   und damit:   t'  = [(rQ  - r0 )·g ]/g2.
# Der gesuchte Abstand d zwischen den  Punkten Q und Q'  ergibt sich dann aus der Länge des Verbindungsvektors r[Q']- r[Q]. ;Wir rechnen dazu ein Zahlenbeispiel. Die Gerade  g  wird  durch die  Vektoren r[0] und r[1]  in  der Form   g=r[1] - r[0]    beschrieben. ;
> r0:=Vector([1,1,1]): r1:=Vector([2,3,4]): g:=r1-r0:
> r:=unapply(r0+t*g,t):
# Der Punkt Q hat den Ortsvektor
> rQ:=Vector([2,1,3]):
> tst:=DotProduct(rQ-r0,g)/DotProduct(g,g): rQs:=r(tst):
> print(tst,rQs);
# Der gesuchte Abstand d ergibt sich dann aus der Beziehung
> d:=Norm(rQs-rQ,Euclidean);
# Eine Ebene, die den Punkt r0  enthält und von den Vektoren p und q aufgespannt wird, lässt sich mittels zweier Parameter s und t durch die Gleichung r = r0  + s p + t q darstellen (Ebenengleichung). Wir rechnen dazu ein Zahlenbeispiel:
> r0:=Vector([6,3,6]): p:=Vector([6,-1,1]):  q:=Vector([6, 5,2]):
> print(r0,p,q);
> r:=unapply(r0 + s*p + t*q,s,t):
# Wir stellen die durch die Vektoren r0, p und q festgelegte Ebene grafisch dar:
> par:=shape=cylindrical_arrow,width=0.4, head_length=0.8:
> r0t:=op(convert(r0,list)):r0p:=op(convert(r0+p,list)): r0q:=op(convert(r0+q,list)):
> p0:=plots[arrow](r0  ,par,color=black):
> pp:=plots[arrow](r0,p,par,color=  red):
> pq:=plots[arrow](r0,q,par,color=green):
> pl:=plot3d(r(s,t),t=-1.1..1.1,s=-1.1..1.1,style=patchnogrid,color=blue,transparency=0.5):
> plt:=plots[textplot3d]({[r0t,"r0",'align'='above'],[r0p,"p",'align'='above'],[r0q,"q",'align'='above']},'font'=["times",16]):
> plots[display]([p0,pp,pq,pl,plt],axes=normal,scaling=constrained,labels=[x, y, z],orientation=[-55,75,50]);
# Der Vektor p × q steht senkrecht auf der Ebene, die durch  p und q aufgespannt wird, und sein im Sinne der Rechtsschraubregel definierter Einheitsvektor n  legt die Orientierung der Ebene im Raum fest:
> v:=CrossProduct(p,q);
# Wir bilden daraus den Einheitsvektor
> n:=evalf(Normalize(v,Euclidean));
# Wir wählen nun im Raum einen Punkt W mit dem Ortsvektor w.  Wir wollen den Abstand dieses Punktes von der Ebene berechnen. Der Punkt  W  kann dabei oberhalb,unterhalb oder auch in der  Ebene selbst liegen. Was oben bzw.  unten ist, wird durch die Orientierung des Normalenvektors  n festgelegt. 
> w:=Vector([10,0,15]);
# Das Lot  vom Punkt  W  auf die Ebene trifft diese in einem Punkt P mit dem Ortsvektor r. Der Verbindungsvektor  von P  nach W  muss parallel zu n sein, also  w - r  =  » n. Weil n senkrecht auf p und q steht, erhalten wir zwei Gleichungen:
# 1.)  r*p; =w*p;,
# 2.)  r*q; = w*q.;
> gl1:=DotProduct(r(sp,tp),p,cf)= DotProduct(w,p);
> gl2:=DotProduct(r(sp,tp),q,cf)= DotProduct(w,q);
> gln:=gl1,gl2;
> loe:=evalf(solve({gln},{sp,tp}));assign(loe);
# Wir haben mit dem Befehl assign die berechneten Werte den Parametern sp und tp zugeordnet und bestimmen den Ortsvektor  des Punktes  P
> rp:=r(sp,tp);
# Der Lotvektor vom Punkt P zum Punkt W ist dann:
> lot:=w-rp;
# Der Abstand des Punktes  W  von der Ebene ist 
> d:=evalf(Norm(lot,Euclidean));
# Um die Lage des Punktes W zur Ebene zu bestimmen, berechnen wir den Faktor  »  aus der Beziehung » = (w - r )Ån =  lotÅn
# » < 0:  W liegt unterhalb der Ebene,
# » > 0:  W liegt oberhalb der Ebene,
# » = 0:  W liegt in der Ebene.
> lambda:=DotProduct(lot,n);
# Der Punkt W liegt also oberhalb der Ebene. Wir stellen die Ebene und den Vektor lot grafisch dar. Dazu benötigen wir den Befehle plot3d  sowie den Befehl  display3d  aus der Bibliothek plots
> rpt:=op(convert(rp,list)): wt:=op(convert(w,list)):
> pl:=plot3d(r(s,t),t=-1.1..1.1,s=-1.1..2,style=patchnogrid,color=blue,transparency=0.5):
> plt:=plots[textplot3d]({[rpt,"P",'align'='left'],[wt,"W",'align'='above']},'font'=["times",16]):
> pw :=plots[arrow](rp,lot,shape=cylindrical_arrow,width=0.6, head_length=1.,color = green):
> plots[display3d]([pl,plt,pw],axes=framed,scaling=constrained,labels=[x, y, z],orientation=[-55,75,50]);
# 
# Wir wenden uns nun der Verarbeitung von Matrizen zu. 
> restart:with(LinearAlgebra):with(plots): unprotect(D,I):
# Zur Erzeugung einer Matrix stehen in Maple neben der Option Matrix-Palette die Befehle Matrix  und   MVshortcut  aus der Bibliothek LinearAlgebra  zur Verfügung. Spezielle Matrizen sind
# 1.) Die Nullmatrix (hier mit 2 Zeilen und 3 Spalten)
> C:=Matrix(2,3);
# 2.) Die quadratische Matrix (hier mit 3 Zeilen und 3 Spalten)
> A:=Matrix(3,symbol=a): B:=Matrix(3,symbol=b):
> print(A,B);
# 3.) Die Diagonalmatrix
> NULL;D:=Vector(3,symbol=d): DM:=DiagonalMatrix(D,3):
> print(D,DM);
# 4.)  Die Einheitsmatrix  ist eine quadratische Matrix, die auf der Hauptdiagonale jeweils mit einer 1 besetzt ist.
> ID:=IdentityMatrix(3);
# Die Einheitsmatrix bildet jede Matrix auf sich selbst ab:
> A.ID;
> alpha*ID.Vector(3,symbol=a);
# 5.) Die obere Rechtsdreiecksmatrix. Nach Fortran-Konvention gilt (i: Zeilenidex, j: Spaltenindex) . Wir konstruieren uns die folgende obere Rechtsdreiecksmatrix
> g:=(i,j)->i^(i+j-1);
> G:=Matrix(3,g);
> OR:=Matrix(3,g,shape=triangular[upper]);
# 6.) Die untere Linksdreiecksmatrix
> UL:=Matrix(3,g,shape=triangular[lower]);
# Zu jeder Matrix A existiert eine transponierte Matrix AT,  in der gegenüber der Ausgangsmatrix A die  Zeilen und Spalten vertauscht sind.
> AT:=Transpose(A);
# Es gelten die folgenden Rechenregeln für Matrizen (I: Einheitsmatrix)
# I*A = I*A and I*A = A;
# A*0 = 0*A and 0*A = 0*(A^T)^T and 0*(A^T)^T = A;
# (A+B)^T = A^T+B^T;
# (A*B)^T = B^T*A^T;
> A.ID;ID.A;
> Transpose(AT);
> Transpose(A+B);Transpose(A)+Transpose(B);
> Transpose(A.B);Transpose(B).Transpose(A);
# Eine quadratische Matrix S heißt symmetrisch, wenn gilt  ST = S  
> S:=Matrix(3,3,symbol=s,shape=symmetric);
> Transpose(S)-S;
# Eine quadratische Matrix S heißt schiefsymmetrisch oder antisymmetrisch (antimetrisch), wenn gilt  ST = -S 
> S:=Matrix(3,3,symbol=s,shape=antisymmetric);
> Transpose(S)+S;
# Jede quadratische Matrix A kann eindeutig in einen symmetrischen und einen antimetrischen Anteil zerlegt werden:  
# A=A[S]+A[A]  mit  A[S]=1/(2)(A+A^(T))  und  A[A]=1/(2)(A-A^(T)).; Dazu ein Beispiel:
> A:=Matrix([[5,1,-4],[3,7,8],[-2,0,3]]): AT:=Transpose(A):
> print(A,AT);
> AS:=1/2*(A + AT): AA:=1/2*(A - AT):
> print(AS,AA);
# Kontrolle:
> AS + AA - A;
# Zwei Matrizen A und B von je m Zeilen und n Spalten werden addiert (subtrahiert), indem ihre entsprechenden Elemente addiert (subtrahiert) werden. Das Ergebnis ist eine (m×n)-Matrix:
> A:=Matrix(4,3,symbol=a): B:=Matrix(4,3,symbol=b):
> print(A,B);
> print(A+B,A-B);
# Eine Matrix A wird mit einem Skalar » multipliziert, indem jedes Element der Matrix mit » multipliziert wird
> lambda*A;
# oder
> ScalarMultiply(A,lambda);
# Es sei A eine (m×n)-Matrix und B eine (n×p)-Matrix. Dann hat die Produktmatrix C = A*B;  die Dimension (m×p) und ihre Elemente cij errechnen sich als Skalarprodukt aus dem i-ten Zeilenvektor von A und dem j-ten Spaltenvektor von B:
> A:=Matrix(2,3,symbol=a):  B:=Matrix(3,4,symbol=b):
> print(A,B);
> C:=A.B;
> C[2,3];
# Für die Matrizenmultiplikation gilt nicht das Kommutativgesetz, es ist also i. Allg.  A*B <> A*B;
> A:=Matrix([[1,2,3],[2,1,3],[2,3,1]]):
> B:=Matrix([[1,0,0],[0,1,2],[1,0,1]]):
> print(A,B);
> print(A.B,B.A,A.B-B.A);
# Eine quadratische Matrix A heißt orthogonal, wenn gilt  A*A^T = I;. Die Spalten- und Zeilenvektoren einer orthogonalen Matrix A sind untereinander orthonormal, stehen als senkrecht aufeinander, und die Skalarprodukte der Vektoren mit sich selbst sind 1. Die folgenden drei Vektoren erfüllen diese Bedingung. 
> a1:=Vector([1/2*sqrt(2),1/2*sqrt(2),0]):a2:=Vector([0,0,1]):
> a3:=Vector([-1/2*sqrt(2),1/2*sqrt(2),0]):
> print(a1,a2,a3);
# Wir bilden damit die orthogonale Matrix ORTH
> ORTH:=Matrix([a1,a2,a3]):
> print(ORTH,Transpose(ORTH),ORTH.Transpose(ORTH));
# Die Determinante einer quadratischen Matrix A ist eine Zahl. Für eine (3×3)-Matrix erfolgt beispielsweise die Berechnung nach folgender Vorschrift
> A:=Matrix(3,symbol=a);Determinant(A);
# Die Determinante einer oberen Rechtsdreiecksmatrix (untere Linksdreiecksmatrix) ist identisch mit dem Produkt der Hauptdiagonalglieder.
> A:=Matrix(3,symbol=a,shape=triangular[upper]);Determinant(A);
> A:=Matrix(3,symbol=a,shape=triangular[lower]);Determinant(A);
# Hinweis: Eine Matrix A mit det(A)  ` 0 wird reguläre (sonst singuläre) Matrix genannt.
# Unter dem Rang einer Matrix versteht man die höchste Ordnung, die deren nichtverschwindende Unterdeterminanten haben können, wobei für eine nichtquadratische (m×n)-Matrix A gilt: Rang(A) d min{m,n}.
> A:=Matrix(4,5,[[2,-4,3,1,0],[1,-2,1,-4,2],[0,1,-1,3,1],[4,-7,4,-4,5]]);
> Rank(A);
# Der Rang dieser Matrix ist 3 < 4. Zur weiteren Untersuchung wird eine L,U-Zerlegung der Matrix A durchgeführt. Damit kann die Matrix A in der Form  A = P*L*U*notiert*werden;
# P:  Permutationsmatrix, 
# L:  untere (Lower) quadratische Linksdreiecksmatrix 
# U: obere (Upper) Rechtsdreiecksmatrix, hat dieselbe Dimension wie A
> P,L,U:=LUDecomposition(A);
> print(A,P.A,P.P);
# Kontrolle:
> P.L.U - A;
# Der Permutationsmatrix P entnehmen wir folgende Informationen: Die 1. und 4. Zeile bleiben unverändert. Die ursprünglich 2. Zeile wird zur 3. Zeile und die 3. Zeile wird zur 2. Zeile.
# Zur Ausgabe der mit dem Gauß-Algorithmus reduzierten Matrix A ist der folgende Befehl zu verwenden :
> LUDecomposition(A,output = 'U');
# Die letzte Zeile ist eine Null-Zeile, damit reduziert sich der Rang der Matrix um 1.
# Wir wählen im Folgenden eine quadratische Matrix und führen  wieder eine L,U-Zerlegung durch.
> A:=Matrix([[2,3,2],[-1,-1,-3],[3,5,5]]);
> P,L,U:=LUDecomposition(A);
# Die Gauß-Elimination der Matrix A liefert 
> LUDecomposition(A,output='U');Rank(A);
# Die quadratische (m×m)-Matrix X mit der Eigenschaft X*A = X*A and X*A = I;  heißt  inverse  Matrix von A. Dabei ist A eine reguläre (m×m)-Matrix und  I  ist die (m×m)-Einheitsmatrix.  Wir berechnen die Inverse der quadratischen Matrix A
> AI:=MatrixInverse(A);
# oder auch kürzer
> AI:=A^(-1);
# Kontrolle
> A.AI;
# Besonders einfach ist die Invertierung einer Diagonalmatrix
> MatrixInverse(DiagonalMatrix(D));
# Für die Inverse eines Matrizenproduktes gilt (A*B)^(-1)=B^(-1)*A^(-1). ;Wir zeigen das an einem Beispiel
> A:=Matrix([[2,5],[1,7]]); B:=Matrix([[0,-4],[4,3]]); AB:=A.B;
> MatrixInverse(AB);MatrixInverse(B).MatrixInverse(A);
# Im Folgenden werden Zeilen- und Spaltenmanipulationen an einer Matrix vorgenommen. Dazu verwenden wir folgende (4×5)-Matrix
> A:=Matrix([[2,-4,3,1,0],[1,-2,1,-4,2],[0,1,-1,3,1],[4,-7,4,-4,5]]);
# Die Anzahl der Zeilen und Spalten von A erhalten wir mit den Befehlen
> RowDimension(A); ColumnDimension(A);
# Wir bilden aus A die Submatrix B, die die Zeilen 2..3 und die Spalten 3..4 von A enthält
> B:=SubMatrix(A,[2..3],[3..4]);
# Herausgreifen der 1. Zeile von A
> Z1_A := Row(A,1); whattype(%);
# Herausgreifen der 3. Spalte von A :
> S3_A := Column(A,3);whattype(%);
# Vertauschung der 2. und  3.  Zeile, wobei die alte Matrix A erhalten bleibt:
> RowOperation(A,[2,3]):
> print(%,A);
# Soll die alte Matrix A überschrieben werden, dann ist der Zusatz  inplace  zu benutzen:
> RowOperation(A,[2,3],inplace):
> print(%,A);
# Die 3. Zeile wird mit dem Skalar  ´ multipliziert
> RowOperation(A,3,delta,inplace);
# Die neue zweite Zeile setzt sich zusammen aus der alten 2. Zeile plus  5 mal der 1. Zeile
> RowOperation(A,[2,1],5,inplace);
# Von der 4. Zeile wird das Doppelte der 1. Zeile subtrahiert
> RowOperation(A,[4,1],-2,inplace);
# Um Matrizen, Vektoren und Arrays miteinander zu verketten, wird der Befehl Concatenate aus der Biblothek ArrayTools benötigt.
> print(B,Transpose(B));
> ArrayTools[Concatenate](1,B,Transpose(B));
> ArrayTools[Concatenate](2,B,Transpose(B));
# Hinweis: Das Vektorprodukt a × b  zweier Vektoren a und b , die in einem orthonormalen Basissystem vorliegen, kann formal immer als reines Matrizenprodukt  a × b = SAÅb  mit einer schiefsymmetrischen Matrix SA  berechnet werden.
> A:=Vector(3,symbol=a):B:=Vector(3,symbol=b):
# Setzen wir 
> SA:=Matrix([[0,-a[3],a[2]],[a[3],0,-a[1]],[-a[2],a[1],0]]);
# dann ist, wie behauptet,
> CrossProduct(A,B)=SA.B;
# Wir berechnen die Determinante und den Rang der schiefsymmetrischen Matrix SA
> Determinant(SA);Rank(SA);
# Das zweifache Vektorprodukt  wird dann a × (a × b ) =  SAÅ(SAÅb) = SA2Åb = HAÅb  mit der durch den Vektor a eindeutig festgelegten symmetrischen Matrix HA
> HA:=SA.SA;
> simplify(CrossProduct(A,CrossProduct(A,B))) = simplify(HA.B);
# Wir berechnen noch die Determinante und den Rang der symmetrischen Matrix HA
> Determinant(HA);Rank(HA);
# In der Worksheet-Version von Maple werden nur Matrizen angezeigt, deren Reihen- und Spaltendimensionen  im Bereich 1..10 liegen. Matrizen, deren Dimensionen größer sind, werden durch einen Platzhalter angedeutet. Größere Matrizen können durch den Befehl  interface (rtablesize = value) angezeigt werden, wobei value die obere Grenze der Reihen- und Spaltendimension angibt. Mit einem Doppel-Klick auf den Platzhalter kann die Matrix temporär angesehen werden. 
> E:=Matrix(10,11,symbol = e);
> interface(rtablesize=12);E;
> interface(rtablesize=10);
# Abbildungen
> A:=Matrix(3,symbol=a);
# Der k-te Spaltenvektor der Matrix A ist der Bildvektor y der linearen Abbildung y=A*e[k]=a[k].; Wir zeigen das am Beispiel einer (3×3)-Matrix. 
> e1:=UnitVector(1,3): e2:=UnitVector(2,3): e3:=UnitVector(3,3):
> print(e1,e2,e3);
> print(A.e1,A.e2,A.e3);
# Bei einem allgemeinen Vektor x ist die Abbildung y = A*x;eine Linearkombination der Spaltenvektoren a[k];  der Abbildungsmatrix A .
> X:=Vector(3,symbol=x):
> y:=A.X; whattype(%);
# Mit der Transponierten von A kann die sog.  Linksmultiplikation  y = x^T*A^T;dargestellt werden. Das Ergebnis ist ein Zeilenvektor 
> y:=Transpose(X).Transpose(A);whattype(%);
# Ist die Abbildung speziell eine ebene Drehung um die 3-Achse mit dem Winkel Õ, dann gehen die Einheitsvektoren e[1], e[2]; über in die drei um den Winkel Õ gedrehten Einheitsvektoren, wobei der Einheitsvektor e[3]; unverändert bleibt: 
> a1:=Vector([ cos(phi),sin(phi),0]):
> a2:=Vector([-sin(phi),cos(phi),0]):
> a3:=Vector([0        ,0       ,1]):
> R:=Matrix([a1,a2,a3]);
# Die Koordinaten eines Punktes P mit den Koordinaten x[1],;x[2];,x[3];  im Ausgangssystem sind dann im gedrehten Koordinatensystem
> y:=R.X;
# Ist die Inverse von R identisch mit  R^(T), also  R^(-1)=R^(T), ;dann handelt es sich bei R um eine orthogonale Matrix
# Hinweis: Gilt det(R) =  +1, dann handelt es sich bei R um eine Drehmatrix
> simplify(MatrixInverse(R))-Transpose(R);
# Der Rang und die Determinante von R sind:
> Rank(R);simplify(Determinant(R));
# Der Ausdruck  x^T*A*x;  mit einer symmetrischen Matrix A wird  quadratische Form genannt, und   q(x)=;x^T*A*x;+b^(T)*x+c;   heißt quadratisches Polynom. Dazu ein Zahlenbeispiel:
> X:=Vector(2,symbol=x): A:=Matrix(2,[1,-1,-1,3]): b:=Vector([1,1]): c:=-1:
> print(X,A,b,c);
> f:=expand(Transpose(X).A.X + Transpose(b).X + c);
> plots[implicitplot](f,x[1]=-3..1,x[2]=-2..1,thickness=2,coloring=[blue,red], filledregions=true);
# Im Dreidimensionalen sieht das so aus:
> X:=Vector(3,symbol=x): A:=Matrix(3,[-1,0,0,0,2,0,0,0,2]); 
> b:=Vector([1,1,1]): c:=-1:
> f:=expand(Transpose(X).A.X + Transpose(b).X + c);
> plots[implicitplot3d](f,x[1]=-2..2,x[2]=-2.5..2.5,x[3]=-2.5..2.5,axes=framed,scaling=constrained,style=patchnogrid,grid=[60,60,60],lightmodel=light3,shading=zgrayscale,orientation=[-45,55,0]);
# Für das Vektorprodukt einer Matrix  A  mit einem Vektor c,  also  A × c  bzw.  c × A  gelten folgende Zusammenhänge: 
# A × c = - (c × AT)T,    c × A = - (AT × c)T .
# Zur Berechnung des Vektorproduktes A × c stellen wir die folgend Prozedur zur Verfügung. Dabei muss die Matrix A von der Dimension (3×3) sein, und der Vektor C hat die Dimension (3×1). Das Ergebnis ist eine (3×3)-Matrix.
> Physics[Setup](dimension = 3,signature =`+`): 
# Wir sehen uns die aktuelle Rechnerumgebung für das Paket Physics an.
> Physics[Setup](query);
> A_kreuz_C:=proc(A::Matrix,C::Vector)
> #----------------------------------------------------------------------
# #  Eingabe:
# #    A  :         Matrix (3 x 3)
> 
> #    C  :         Vektor (3 x 1)
> #  Ausgabe 
> #     B = A x C:  Matrix (3 x 3)
> #----------------------------------------------------------------------
> local l,i,k,j;
> global B;
> description "Vektorprodukt einer Matrix A mit einem Vektor C: B = A x C";
> B:=Matrix(3,3):
> for l from 1 to 3 do
>   for i from 1 to 3 do
>     for k from 1 to 3 do
>       for j from 1 to 3 do
>         B[i,l]:= B[i,l] + A[i,j]*C[k]*Physics[LeviCivita][j,k,l]:
>       end do:
>     end do:
>   end do:
> end do:
> end proc:
> a:='a':c:='c':
> A:=Matrix(3,3,symbol=a): C:=Vector(3,symbol=c):
> print(A,C);
> A_kreuz_C(A,C):
> print(`A_kreuz_C = `,B);
# Um das Kreuzprodukt c × A = - (AT × c)T zu berechnen, ermitteln wir zunächst  AT × c:
> A_kreuz_C(Transpose(A),C):
> print(B);
# Abschließend multiplizieren wir die transponierte Matrix B mit -1.
> C_kreuz_A:=-Transpose(B);
# Wir zeigen mit A = I  die Beziehung:  I × c  =  c × I  (antisymmetrische Matrix). 
> A:=Matrix(3,shape=identity);
> A_kreuz_C:=A_kreuz_C(A,C): print(B);
> A_kreuz_C(Transpose(A),C);
> C_kreuz_A:=-Transpose(B);
# Ist A eine symmetrische Matrix, dann gilt:  A × c = - (c × A)T,    c × A = - (A × c)T
> A:=Matrix(3,symbol=a,shape=symmetric);
> A_kreuz_C(A,C):print(B);
> C_kreuz_A:=-Transpose(B);
> 
;
# 
# 
