# Lektion 6
# Lineare Gleichungssysteme
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra):
# Existenz und Eindeutigkeit von Lösungen:  Als Beispiel betrachten wir die folgenden drei Gleichungen, die wir zu einer Liste zusammenfassen:
> gl1 :=-x+y+z=0; gl2:=x-3*y-2*z=5; gl3:=5*x+y+4*z=3;
> gls1:=[gl1,gl2,gl3];
# Die Unbekannten sind:
> unb :=[x,y,z];
# Die Auflösung des Gleichungssystems nach den drei Unbekannten  x , y  und z  erfolgt mit dem Befehl  solve . Wir erhalten die eindeutige Lösung
> solve(gls1,unb);
# Wir können uns von Maple auch die  Koeffizientenmatrix  A  und die rechte Seite b  des Gleichungssystems A*x = b; ausgeben lassen:
> A,b:=GenerateMatrix(gls1,unb);
# Damit haben wir
> A.Vector(unb)=b;
# Die um die rechte Seite b erweiterte Koeffizientenmatrix erhalten wir wie folgt:
> AE:= GenerateMatrix(gls1,unb,augmented=true);
# Die Lösung des Gleichungssystems  mit dem Befehl  LinearSolve  liefert
> v:=LinearSolve(A,b);
# Wir kontrollieren das Ergebnis
> A.v-b;
# Ist m die Anzahl der Gleichungen und n die Anzahl der Unbekannten, dann ist das lineare Gleichungssystems A*x = b;lösbar, wenn der Rang der Koeffizientenmatrix des erweiterten Gleichungssystems  Rg(A|b)=Rg(A)=r;  ist. Für r = n  ist das System eindeutig lösbar. Für r < n  ist die Lösung nicht eindeutig. Im obigen Fall liegt ein quadratisches Gleichungssystem mit m = n = 3 vor. 
> Rank(A);Rank(AE);
# Die Koeffizientenmatrix  A und die erweiterte Koeffizientenmatrix AE haben denselben Rang r = 3, der der Anzahl der Unbekannten entspricht. Damit ist das System lösbar, und die Lösung ist eindeutig. 
# Das folgende Gleichungssystem besitzt keine eindeutige Lösung:
> gl3:=-3*x+9*y+6*z=-15;
> gls2:=[gl1,gl2,gl3];
> A,b:=GenerateMatrix(gls2,unb);
# Die dritte Zeile entspricht offensichtlich der mit -3 multiplizierten zweiten Zeile. 
> AE:= GenerateMatrix(gls2,unb,augmented=true);
> A.Vector(unb) = b;
# Wir berechnen den Rang von A und AE
> Rank(A);Rank(AE);
# Wegen Rg(A) = Rg(A|b); = 2  < 3 = n  ist die Lösung nicht eindeutig. Es existieren unendlich viele Lösungen mit einem freien Parameter (hier mit `#msub(mi("t"),mo("&InvisibleTimes;"))`; bezeichnet).
> v:=LinearSolve(A,b,free='t');
# Das folgende System ist unlösbar:
> gl1:=x+y+z=1; gl2:=-x-2*y+z=2; gl3:=x-y+5*z=0;
> gls3:=[gl1,gl2,gl3];
> solve(gls3,unb);
# Die Unlösbarkeit des Systems quittiert Maple mit der Ausgabe der leeren Liste [ ].  Wir beschaffen uns den Rang von A und AE
> A,b:=GenerateMatrix(gls2,unb);
> AE:= GenerateMatrix(gls3,unb,augmented=true);
> Rank(A);Rank(AE);
# Mit  Rg(A) = 2  und  Rg(A|b); = 3  ist das System unlösbar. Eine LU-Zerlegung liefert für dieses Beispiel
> P,L,U := LUDecomposition(A);
# Auf der Hauptdiagonalen der Matrix U steht eine 0. Damit ist A singulär, was wir durch Berechnung der Determinante bestätigen
> Determinant(A);
# In nummerischen Anwendungen erfolgt die Lösung des linearen Gleichungssystems A*x = b; in folgenden Schritten:
# 1. Bestimmung der LU-Zerlegung der Matrix A:  A=P*L*U 2. ;Löse P*z = b;   nach z
# 3. Löse L*y = z;   nach y
# 4. Löse U*x = y*nach*x;
# Wir zeigen das an einem Beispiel
> A:=Matrix(4,[0,5,3,-10,3,7,-3,5,12,4,4,4,0,12,0,-8]): 
> b:=Vector([-10,14,8,-8]):
> print(A,b);
# Zu 1.:  LU-Zerlegung von A
> P,L,U := LUDecomposition(A);
# Zu 2.: Unter Beachtung von 1/P = P^T;  (P ist eine orthogonale Matrix) folgt 
> z:=Transpose(P).b;
# Zu 3.: 
> y:=LinearSolve(L, z);
# Zu 4.: 
> x:=LinearSolve(U, y);
# oder mit einem Befehl, wobei uns die obigen Zwischenschritte verborgen bleiben
> x:=LinearSolve(A, b);
> x:='x':y:='y':z:='z':
# Anwendung beim Dreibock
# Ein Dreibock besteht aus drei Stäben, die an einem Punkt, dem Kopfpunkt, zusammengeführt sind. Zur Ermittlung der Stabkräfte benötigen wir eine Knoten- und eine Elementdatei. 
# 1. Knotendatei: 
# Die Knotendatei, sie enthält die kartesischen Koordinaten der Systemknoten, besteht aus einer zweifach geschachtelten Liste, in Maple  listlist genannt, bei der die Listenelemente wieder Listen sind. 
# Beispiel: KNOTEN:=[[0,1,1],[0,0,0],[2,0,0],[2,2,0]];
# Es wird folgendes vereinbart:
# Die erste Liste der Liste enthält die kartesischen Koordinaten des Kopfpunktes (hier [0,1,1]).  Ihr folgen die Koordinaten der drei Fußpunkte 2, 3 und 4. 
# 2. Elementdatei:
# Die Elementdatei entspricht auch einer zweifach geschachtelten Liste mit drei Listeneinträgen für jeden Stab. 
# Beispiel: ELEM:=[[1,4],[1,3],[1,2]];
# Mit der Liste [1,4] wird der 1. Stab (1. Liste in der listlist) mit seinem Anfangspunkt (1)  und seinem Endpunkt (4) festgelegt. Damit liegt auch die Orientierung des 1. Zugstabes fest. Der zweite Stab verläuft vom Knoten 1 zum Knoten 3. Die entsprechenden Knotenkoordinaten können der Knotendatei  entnommen werden. Wird die Stabkraft bei der oben vorgenommenen Elementorientierung positiv berechnet, dann handelt es sich um eine Zugkraft, andernfalls um eine Druckkraft.
# Auf den Knoten 1 wirkt die äußere Kraft F = [Fx, Fy, Fz]. Aus der Stabgeometrie sind die Stabeinheitsvektoren zu berechnen, die spaltenweise die Geometriematrix  G bilden. Ist S = [S1, S2, S3] der Vektor der gesuchten Stabkräfte, dann ist das Gleichungssystem  GÅS = -F zu lösen.  
# Wir schreiben vorab eine Prozedur bockplot, die uns die  grafische Darstellung des Dreibocks liefert.
> bockplot:=proc(knot::listlist,elem::listlist,f::Vector(3))
> local i,ra,re,ri,mst,stlae,stl,maxl,staebe,fn,
>       lambdaf,kraft,linie,ral,rel,kb,eb,textk,texte;
> description "Grafische Darstellung eines Dreibocks";
> if nops(knot)   <> 4 then
> print(`Die Anzahl der Knoten muss 4 betragen`);
> elif nops(elem) <> 3 then
> print(`Die Anzahl der Elemente muss 3 betragen`);
> return; end if;
> linie:={}: stlae:=[]:
> ra:=Vector(knot[1]);
> for i to 3 do
> re    :=Vector(knot[elem[i,2]]);
> stl   :=LinearAlgebra[Norm](re-ra,Euclidean);
> stlae :=[op(stlae),stl];
> ral   :=convert(ra,list);rel:=convert(re,list);
> staebe:=plottools[line](ral,rel,color=blue,linestyle=1,thickness=3);
> linie :=linie union {staebe}:
> end do:
> maxl:= max(stlae):
> fn:= evalf(LinearAlgebra[Norm](f,Euclidean));
> if fn > 0 then 
>  lambdaf:=1/3*maxl/fn:
> end if:
> kraft:=plots[arrow](ra,lambdaf*f,shape=cylindrical_arrow,width=0.1,head_length=0.3,color=gold):
> kb:={}:
> for i to 4 do           #Knotenbeschriftung
> textk:=plots[textplot3d]([op(knot[i]),i],font =[TIMES,BOLD,14],color=red,align=RIGHT):
> kb:= kb union {textk}:
> end do:
> eb:={}:
> for i to 3 do           #Elementbeschriftung
> ri:=1/2*(op(knot[1])+op(knot[i+1]));
> texte:=plots[textplot3d]([ri,i],font=[TIMES,BOLD,14],color=blue,align=RIGHT):
> eb:= eb union {texte}:
> end do:
> plots[display](linie,kraft,kb,eb,axes=boxed,scaling=constrained,orientation=[-60,65,0],labels=[x, y, z]);
> end proc:
# Wir fassen den Lösungsalgorithmus zur Berechnung der Stabkräfte in der Prozedur stabkr zusammen, die dieselbe Parameterliste wie die Prozedur bockplot besitzt. 
> stabkr:=proc(knot::listlist,elem::listlist,f::Vector(3))
> local eps,i,e1,e1n,e2,e2n,e3,e3n,rGdiff;
> global G,GE,det,rG,rGE,S;
> description "Berechnung der Stabkräfte eines Dreibocks";
> eps:=10^(-(Digits-2)):
> if LinearAlgebra[Norm](f) < eps then
>   print(`Sämtliche Stabkräfte sind null`);
> return end if: 
> e1 :=Vector(3): e2:=Vector(3): e3:=Vector(3):
> G  :=Matrix(3,3);
> for i to 3 do
>   e1[i]:=knot[elem[1,2],i]-knot[elem[1,1],i];
>   e2[i]:=knot[elem[2,2],i]-knot[elem[2,1],i];
>   e3[i]:=knot[elem[3,2],i]-knot[elem[3,1],i];
> end do:
> e1 :=LinearAlgebra[Normalize](e1,Euclidean):
> e2 :=LinearAlgebra[Normalize](e2,Euclidean):
> e3 :=LinearAlgebra[Normalize](e3,Euclidean):
> G  :=Matrix([e1,e2,e3]):
> GE :=Matrix([G,-f]):
> det:=LinearAlgebra[Determinant](G):
> rG :=LinearAlgebra[Rank](G); rGE:=LinearAlgebra[Rank](GE);
> rGdiff:=rGE-rG;
> if rGdiff <> 0 then
>  print(`Das System ist unlösbar`);
>  return end if: 
> S:=LinearAlgebra[LinearSolve](G,-f);
> end proc:
# 1. Beispiel:
> KNOTEN:=[[0.,1.,2.],[0.,0.,0.],[2.,2.,0.],[-2.,2.,0.]]; #Knotendatei:
> ELEM:=[[1,2],[1,3],[1,4]];  #Elementdatei:
> F:=Vector([0.,0.,-1.]);     #Kraftvektor:
# Grafische Ausgabe des Dreibocks mit Belastung
> bockplot(KNOTEN,ELEM,F);
# Berechnung der Stabkräfte
> stabkr(KNOTEN,ELEM,F);
> print(rG);print(rGE);
# Wegen Rg(A|b)=Rg(A)=3 ;ist das System eindeutig lösbar. Alle Stabkräfte haben eine negatives Vorzeichen und sind deshalb Druckkräfte. 
# 2. Beispiel:
> KNOTEN:=[[0.,1.,0.],[0.,0.,0.],[2.,2.,0.],[-2.,2.,0.]];  #Knotendatei:
> ELEM:=[[1,2],[1,3],[1,4]];  #Elementdatei:
> F:=Vector([0.,0.,-1.]);     #Kraftvektor:
# Die grafische Ausgabe zeigt, dass alle Stäbe in einer Ebene liegen. 
> bockplot(KNOTEN,ELEM,F);
> stabkr(KNOTEN,ELEM,F);
> print(rG);print(rGE);
# Da der Rang der Geometriematrix und der Rang der erweiterten Geometriematrix nicht gleich sind, ist das System unlösbar. Weil sämtliche Stäbe in einer Ebene liegen, kann die äußere Kraft F, die eine zur Stabebene senkrechte Komponente besitzt, vom Stabsystem nicht aufgenommen werden. Das geht erst dann wieder, wenn beispielsweise der Kopfpunkt 1 eine (kleine) z-Koordinate besitzt.
# 3. Beispiel:
> KNOTEN:=[[0.,1.,0.],[0.,0.,0.],[2.,2.,0.],[-2.,2.,0.]]; #Knotendatei:
> ELEM:=[[1,2],[1,3],[1,4]];   #Elementdatei:
> F:=Vector([1.,1.,0.]);       #Kraftvektor:
# Die grafische Ausgabe zeigt, dass alle Stäbe, und nun auch die äußere Belastung F , in einer Ebene liegen. Damit liegt eine ebenes Problem vor.
> bockplot(KNOTEN,ELEM,F);
> stabkr(KNOTEN,ELEM,F);
> print(rG); print(rGE);
# Wegen Rg(A|b) =  Rg(A)=2<3; ist die Lösung nicht eindeutig. Die mechanische Begründung liegt hier darin, dass  die  Berechnung der drei Stabkräfte aus reinen Gleichgewichtsbedingungen allein nicht möglich ist, da bei einem ebenen Problem nur zwei Kraftgleichgewichtsbedingungen zur Verfügung stehen. Es handelt sich in diesem Fall um eine statisch unbestimmte Aufgabe, die zusätzlich eine Aussage über das Verfomungsverhalten der Stäbe erfordert. Erst dann wird die Lösung wieder eindeutig.  
> 
;
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
