# Lektion 7
# Matrizeneigenwertprobleme
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(PolynomialTools):
# In der Grundvorlesung zur Technischen Mechanik treten Matrizeneigenwertprobleme u.a. bei der Berechnung der Hauptträgheitsmomente eines ebenen Querschnitts und der Ermittlung der Hauptspannungen eines ebenen oder auch räumlichen Spannungszustandes auf. Diese Hauptwerte zeichnen sich dadurch aus, dass sie Extremwerte annehmen. Als einführendes Beispiel betrachten wir die Spiegelung eines Punktes P nach P` an einer Ebene mit dem Stellungsvektor n.

# 
#  Der Vektor n habe den Betrag 1 (Einheitsvektor). Dann entnehmen wir der obigen Grafik: rsp = r - 2 (r · n) n = r ·[ I - 2n—n]  = rÅA(n),  mit A(n) = I - 2n—n =A(-n).
#  Wir interessieren uns nun für alle gespiegelten Vektoren rsp, die ein Vielfaches von r darstellen. Diese Vektoren besitzten dann dieselbe Orientierung oder Gegenorientierung wie r. Mit rsp =  » r folgt nach Zusammenfassung, unter Beachtung von  A = I - 2 n—n , das homogene lineare Gleichungssystem (A - » I)År = 0, das nur dann eine nichttriviale Lösung besitzt, wenn die Determinante der charakteristischen Matrix  A - » I  verschwindet, also det (A - » I) = 0 erfüllt ist. Der noch freie Parameter  »  heißt Eigenwert der Matrix A. 
> N:=Vector(3,symbol=n);
> A:=IdentityMatrix(3)-2*N.Transpose(N);
# Die Summe der Hauptdiagonalelemente von  A  heißt Sp(A) (engl. trace)
> SPA:= Trace(A);
# Die charakteristische Matrix A - » I  und das charakteristische Polynom ermitteln wir mit den Befehlen  CharacteristicMatrix  bzw.  CharacteristicPolynomial  aus der Bibliothek LinearAlgebra.
> CA:=CharacteristicMatrix(A,lambda);
> CP:=CharacteristicPolynomial(A,lambda);
# Die Koeffizienten bei  -»2,  »1  und  -»0 bezeichnet man als erste, zweite und dritte Hauptinvariante der Matrix A.
> CL:=CoefficientList(CP,lambda,'termorder' = 'reverse');
> I1:=-CL[2];I2:=CL[3];I3:=-CL[4];
# Für die Invarianten von A  gilt:
# I[1](A) = Sp(A)*I[2](A) and Sp(A)*I[2](A) = (1/2)[Sp(A)^2-Sp(A^2)]*I[3](A) and (1/2)[Sp(A)^2-Sp(A^2)]*I[3](A) = det*A;
> I1:=Trace(A);
> I2:=simplify(1/2*(SPA^2 - Trace(A.A)));
> I3:=Determinant(A);
# Unter Einbeziehung der Invarianten kann das charakteristische Polynom  in der Form  lambda^3-I[1]*lambda^2+Typesetting[delayDotProduct](I[2], lambda, true)-I[3] = 0; geschrieben werden.
> sort(lambda^3-I1*lambda^2+I2*lambda-I3,lambda,descending);
# Sind  »1, »2 und »3 die Wurzeln (Nullstellen) des charakteristischen Polynoms, dann gilt nach dem Vietaschen Wurzelsatz
> CPV:=(lambda-lambda1)*(lambda-lambda2)*(lambda-lambda3);
> CLV:=CoefficientList(CPV,lambda,'termorder' = 'reverse');
# und für die Hauptinvarianten folgt:
> I1V:=-CLV[2];I2V:=expand(CLV[3]);I3V:=-CLV[4];
# Faktorisieren wir das charakteristische Polynom CP, denn sehen wir, dass » = 1  doppelte Nullstelle ist. 
> factor(CP);
# Die dritte Nullstelle beschaffen wir uns aus dem rechts stehenden Klammerausdruck
> solve(op(%)[2],lambda);
# Sind die Eigenwerte bekannt, dann lassen sich durch  manuelle Rechnungen  zu jedem Eigenwert die Eigenvektoren berechnen. Diese Aufgabe übernimmt Maple mit den  Befehlen Eigenvalues  und  Eigenvectors 
> Eigenvalues(A);
> EW,Phi:=Eigenvectors(A);
# Die Eigenwerte stehen zeilenweise im Vektor EW. Die den Eigenwerten zugeordneten  Eigenvektoren finden sich spaltenweise in der Eigenvektor-Matrix ¦. Wir setzen die berechneten Eigenwerte und Eigenvektoren in die linke Seite (A - » I)År ein und erhalten erwartungsgemäß
> seq(simplify((A-EW[i]*IdentityMatrix(3)).Column(Phi,i)),i=1..3);
# Wir konkretisieren unser Problem und wählen als Spiegelungsebene die (1,3)-Ebene.  Der Normaleneinheitsvektor n hat dann nur eine Komponente in 2-Richtung.
> n[1]:=0: n[2]:=1: n[3]:=0:
# Für die Matrix A und das charakteristische Polynom erhalten wir
> A;
# Die Hauptinvarianten der Matrix A  sind
> I1; I2; I3;
# und für das charakteristische Polynom folgt:
> CP:=factor(CP);
# Neben dem einfachen Eigenwerte » = -1  tritt der doppelte Eigenwert  » = 1  auf. Wir stellen das charakteristische Polynom grafisch dar:
> plot(CP,lambda=-1.5..2,-1..2,axes=boxed,gridlines=true);
# Wir berechnen die Eigenwerte und Eigenvektoren 
> EW,EV:=Eigenvectors(A);Determinant(EV);
# Die drei Eigenvektoren sind linear unabhängig , da die Matrix EV wegen det (EV) = 1 regulär ist.
# Die Lösungen können wie folgt gedeutet werden: Der zum Eigenwert  » = -1 gehörende Eigenvektor e (» = -1) = [0,1,0]T  stellt alle  Punkte P mit den Ortsvektoren r = [0, ±, 0]T  auf der 2-Achse (mit  Ausnahme von ± = 0) dar, die bei Spiegelung an der (1,3)-Ebene in den Gegenvektor rsp = [0, -±, 0]T  übergehen. Der  Eigenwert  » = 1 tritt doppelt auf. Dazu gehören die Eigenvektoren e (» = 1) = [0,0,1]T  und e (» = 1) = [1,0,0]T  . Der Eigenvektor  e = [0,0,1]T   bezeichnet alle diejenigen Punkte P mit den Ortsvektoren r = [0, 0, ²]T  auf der 3-Achse (mit  Ausnahme von ² = 0),  die bei Spiegelung an der (1,3)-Ebene in sich selbst übergehen.  Die analoge Argumentation trifft für den Vektor e = [1,0,0]T zu. 
# Für die Berechnung von Spur und Determinante  gelten folgende Rechenregeln (Beweis durch Ausrechnen)
# 1.  Sp(A+alpha B)=Sp(A)+alpha Sp(B) 2.  Sp(A)=Sp(A^(T));
# 3.  Sp(I*alpha) = 3*alpha*im*Dreidimensionalen;
# 4.  det(A) = det(A^T);
# 5.  det(alpha A)=alpha^(3)det(A)        (im Dreidimensionalen) 6.  det(A*B)=det(A) det(B);
# 7. 1/[det(A)] = det(1/A);
# 8. det(I) = 1;
# Jede quadratische Matrix kann eindeutig additiv in ihre Deviatormatrix A'  und ihre  Kugelmatrix A[K];aufgespalten werden: A(x) = diff(A(x), x)+A(x)[K]; . Für die Deviatormatrix gilt: diff(A(x), x) = A(x)-(1/3)*Sp(A(x))*I;mit  Sp( A' ) = 0. 
# Wir wählen folgende (3×3)-Matrix A, von der wir ihre Deviatormatrix und Kugelmatrix berechnen:
> A:=Matrix(3,symbol=a);
> As:=A-1/3*Trace(A)*IdentityMatrix(3);
> Trace(As);
# Die Kugelmatrix A[K];  ist ein skalares Vielfaches der Einheitsmatrix I
> Ak:=A-As;
# Die Eigenwerte einer  Diagonalmatrix  sind die Diagonalelemente selbst, und die Eigenvektoren sind die kartesischen Einheitsvektoren. Wir zeigen das an folgendem Beispiel
> v:=Vector(3,symbol=d): DM:=DiagonalMatrix(v,3);
> EW,Phi:=Eigenvectors(DM);
# Ist  A eine symmetrische (n×n)-Matrix, dann existiert eine Orthonormalbasis aus reellen Eigenvektoren e[i]; (i=1...n),  ;die wir spaltenweise angeordnet in der Eigenvektormatrix  ¦  vorfinden. Für die Eigenvektormatrix gilt: Phi^T = und*damit*Phi^T/Phi*Phi and und*damit*Phi^T/Phi*Phi = I;. Die Matrix ¦ ist somit eine orthogonale Matrix. Ersetzen wir in dem Eigenwertproblem (A-Typesetting[delayDotProduct](lambda, I, true))*a = 0; den Vektor a durch den neuen Vektor v mit a=Phi*v, ; dann liefert eine Linksmultiplikation der Eigenwertgleichung mit Phi^T; die transformierte Gleichung (Phi^T*A*Phi-Typesetting[delayDotProduct](lambda, I, true))*v = 0.;Die Spektralmatrix  ›=Phi^(T)*A*Phi ; ist aufgrund der Orthogonalitätseigenschaften der Eigenvektoren eine Diagonalmatrix, deren reelle Eigenwerte auf der Hauptdiagonale stehen. Wir zeigen das an einem Zahlenbeispiel für n = 3:
> A:=Matrix(3,{(1,2)=3.,(1,3)=4.,(2,2)=1.},fill=2.,shape=symmetric);
> EW,Phi:=Eigenvectors(A);
# Wir zeigen, dass die Eigenvektormatrix ¦ orthogonal ist:
> Transpose(Phi).Phi;
# Wir geben die Eigenvektoren grafisch aus:
> Student[LinearAlgebra][EigenPlot](A,showsphere=false,showimagevectors=false,eigenoptions=[shape=cylindrical_arrow,width = 0.3,head_length=0.3],axes=boxed,labels=[x,y,z],scaling = constrained);
# Wir berechnen die Spektralmatrix
> Lambda:=Transpose(Phi).A.Phi;
# Für die Matrix A gilt dann umgekehrt
> A:=Phi.Lambda.Transpose(Phi);
> restart: with(LinearAlgebra):
# Wir betrachten abschließend die Kurven  2. Ordnung 
# F(x[1], x[2];) = a2*b*x[1]*x[2]+c*x[2]^2+2*d*x[1]+2*e*x[2]+x[1]^2+f = 0, oder in Matrizenschreibweise;  F(x[1],x[2])=x^(T)*A*x+2 a^(T)*x+f=0,    A=[[[a,b],[b,c]]],   a=[[[d],[e]]],   x=[[[x[1]],[x[2]]]].  Kurven, die aus Geraden bestehen, werden singuläre Kurven 2. Ordnung genannt, die anderen heißen reguläre Kurven oder Kegelschnitte. ;
# Um eine Klassifizierung vornehmen zu können, transformieren wir die Kurve zweiter Ordnung auf die  Normalform, was einer Transformation in ein Hauptzentralachsensystem (HZA) gleichkommt. Das passiert in zwei Schritten:
# 1. Translation des Koordinatensystems mit dem Vektor  y = x-u;.  Der Vektor u wird aus der Gleichung  A*u=-a  ermittelt. ;Damit erscheint die Kurve ohne die linearen Glieder in der Form: 
#  F(y[1], y[2];) = y^T*A*y+g = 0, mit*g = u^T*A*u+2*a^T*u+Typesetting[delayDotProduct](f . 2., Drehung, true)*des*Koordinatensystems and u^T*A*u+2*a^T*u+Typesetting[delayDotProduct](f . 2., Drehung, true)*des*Koordinatensystems in Typesetting[delayDotProduct](Typesetting[delayDotProduct](Richtung*der*Eigenvektoren*von*A . Mit, der, true)*Eigenvektormatrix, Phi, true)*erfolgt*die*Drehung*y and Typesetting[delayDotProduct](Typesetting[delayDotProduct](Richtung*der*Eigenvektoren*von*A . Mit, der, true)*Eigenvektormatrix, Phi, true)*erfolgt*die*Drehung*y = Typesetting[delayDotProduct](Phi*z . Mit, R, true) and Typesetting[delayDotProduct](Phi*z . Mit, R, true) = Phi^T*A*Phi*erhalten*wir*die*Normalform*der*Kurve*zweiter*Ordnung;
# F(z[1], z[2];)  = z^T*R*z+g = ((0.*Die)*Einklassifizierung*Name*der*Kurve*erfolgt*über*die*Vorzeichen*der)*Eigenwerte*und*der*Konstanten*g*mittels*einer*Tabelle;(s.h. Bronstein-Semendjajew, Taschenbuch der Mathematik). ;
# Ist*die*Matrix*A*singulär, dann*ist*das*Gleichungssystem*A*u = -Typesetting[delayDotProduct](a*nicht*lösbar . Es, existiert, true)*ein*Nulleigenwert;, und in diesem Fall müssen programmtechnisch Sonderüberlegungen angestellt  werden.
# Wir stellen zur Transformation von Kurven 2. Ordnung auf Normalform die folgende Prozedur zur Verfügung:
> NORMALFORM:=proc(A::Matrix,a::Vector,f)
> local k,m,n,eps,AT,x,y,z,DET,Phi1,Phi2,R,DA,h,par,p1,p2,p3;
> global FU,F1,FG,lambda,g,Phi,u,phi;
> description "Berechnung der Normalform einer Kurve 2.Ordnung";
> k,m:= Dimension(A);
> n  := Dimension(a);
> eps:=abs(10.^(-Digits+1));
> if k <> 2 or m <> 2 or n <> 2 then 
>   print(`Matrix- oder Vektordimensionen überprüfen`); return;
> end if;
> AT:=LinearAlgebra[Transpose](A);
> if LinearAlgebra[MatrixNorm](A-AT) > eps then
>   print(`Matrix ist nicht symmetrisch`); return; 
> end if;
> x:=Vector([x1,x2]); y:=Vector([y1,y2]); z:=Vector([z1,z2]); u:=Vector(2);
> FU:=LinearAlgebra[Transpose](x).A.x + 2*LinearAlgebra[Transpose](a).x + f;
> lambda,Phi:= LinearAlgebra[Eigenvectors](A);
> phi := convert(arctan(Phi[2,1],Phi[1,1]),degrees);
> DET := LinearAlgebra[Determinant](Phi);
> Phi1:= LinearAlgebra[Normalize](LinearAlgebra[Column](Phi,1),Euclidean,conjugate=false);
> Phi2:= LinearAlgebra[Normalize](LinearAlgebra[Column](Phi,2),Euclidean,conjugate=false);
> if DET > 0 then
>  Phi:= Matrix([ Phi1,Phi2]):
> else
>  Phi:= Matrix([-Phi1,Phi2]):
> end if:
> R :=LinearAlgebra[Transpose](Phi).A.Phi:
> DA:=LinearAlgebra[Determinant](A);
> if abs(DA) > eps then
>   u :=LinearAlgebra[LinearSolve](A,-a);
>   g :=LinearAlgebra[Transpose](u).A.u + 2*LinearAlgebra[Transpose](a).u + f;
>   F1:=LinearAlgebra[Transpose](y).A.y + g;
>   FG:=LinearAlgebra[Transpose](z).R.z + g;
> else
>    h :=evalf(2*LinearAlgebra[Transpose](a).Phi);
> F1:=LinearAlgebra[Transpose](y).R.y + h.y + f;
> if   abs(lambda[1]) < eps  and  abs(h[1]) > eps then
>   u[1]:=-f/h[1]+h[2]^2/(4*h[1]*lambda[2]); u[2]:=-h[2]/2/lambda[2];
>   g:=0:
>   FG:=lambda[2]*z2^2 + h[1]*z1; 
> elif abs(lambda[1]) < eps  and  abs(h[1]) < eps then
>   u[1]:=0; u[2]:=-h[2]/2/lambda[2];
>   g:= f - h[2]^2/4/lambda[2];
>   FG:=lambda[2]*z2^2  + g;
> elif abs(lambda[2]) < eps  and  abs(h[2]) > eps then
>   u[1]:=-h[1]/2/lambda[1]; u[2]:=-f/h[2]+h[1]^2/(4*h[2]*lambda[1]);
>   g:=0:
>   FG:=lambda[1]*z1^2 + h[2]*z2;
> elif abs(lambda[2]) < eps  and  abs(h[2]) < eps then
>   u[1]:=-h[1]/2/lambda[1]; u[2]:=0;
>   g:= f - h[1]^2/4/lambda[1]:
>   FG:=lambda[1]*z1^2  + g;
> end if:
> end if:
> #Ausgabe der Ergebnisse
> par:=axes=normal,gridrefine=4,scaling=constrained,gridlines=true,thickness=2:
> p1:=plots[implicitplot](FU,x1=-10..10,x2=-10..10,par,color=blue):
> p2:=plots[implicitplot](F1,y1=-10..10,y2=-10..10,par,color= red):
> p3:=plots[implicitplot](FG,z1=-10..10,z2=-10..10,par,color=gold):
> print(plots[display](Matrix([p1,p2,p3]),scaling=constrained));
> print(`Eigenwerte lambda = `,lambda, ` g =`,g);
> print(`Gleichung der Kurve 2.Ordnung in Normalform F(z1,z2) = 0 = `,FG);
> print(`Translation des Koordinatensystems mit u = `,u);
> print(`Drehung des Koordinatensystems mit phi = `, phi);
> end proc:
# 1. Beispiel:  Hyperbel
> AM:=Matrix(2,2,[a,b,b,c]): av:=Vector([d,e]):
> a:=1: b:=2: c:=1: d:=-1: e:=1: f:=-4:
> print(AM,av);
> NORMALFORM(AM,av,f);
# 2. Beispiel: Ellipse
> a:=2: b:=-1: c:=2: d:=1: e:=1: f:=1:
> print(AM,av);
> NORMALFORM(AM,av,f);
# 3. Beispiel: Reelles, sich schneidendes Geradenpaar
> a:=1: b:=0: c:=-1: d:=1: e:=-1: f:=0:
> print(AM,av);
> NORMALFORM(AM,av,f);
> 
;
