# Kapitel _03_e: Hauptachsentransformation der Matrix des Massenträgheitsmomententensors
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(Student[LinearAlgebra]):
# Beispiel 3-11:
# Es sollen folgende Aufgaben gelöst werden: 
# 1. Aufstellen des charakteristischen Polynoms p(») = det (˜ - »1) = 0,
# 2.	 Berechnung der Grundinvarianten J1, J2 und J3,
# 3. Berechnung der drei reellen Eigenwerte und normierten Eigenvektoren von ˜,
# 4. Berechnung der Diagonalmatrix  T = ¦TÅ˜Å¦ = diag (»1, »2, »3)
# Laut Aufgabenstellung ist
> Theta:=Matrix([[112. , 0 , 0 ],[ 9. ,100. , 0 ],[ 21. , -36. , 28. ]], shape=symmetric);
# Wir ermitteln die Spur der Matrix des Tensors Theta
> SP:=Trace(Theta);SP2:=Trace(Theta^2); SP3:=Trace(Theta^3);
# 1. Berechnung des charakteristischen Polynoms
> p:=CharacteristicPolynomial(Theta, lambda);
> pl:= PolynomialTools[CoefficientList](p, lambda,'termorder' = 'reverse');
# 2. Berechnung der Grundinvarianten
> J1:=-pl[2]; J2:=pl[3]; J3:=-pl[4];
# Die Grundinvarianten lassen sich alle durch die Spur von ˜ ausdrücken (Buch TM3)
> J1:=SP; J2:=(SP^2-SP2)/2; J3:=(SP^3+2*SP3-3*SP2*SP)/6;
# 3. Lösung des Eigenwertproblems
> lambda,Phi := Eigenvectors(Theta);
# Wir berechnen die Eigenvektoren der Matrix ˜ unter der Verwendung der Option output = list.
> EV:=Eigenvectors(Theta, output = 'list');
# Wir geben den zweiten  Eigenvektor und den zugehörigen (1: einfachen) Eigenwert aus.
> EV[2];
# Wir zeigen, dass die Eigenvektoren die Länge 1 haben und senkrecht aufeinander stehen. Dazu bilden wir die Matrix
> C:=Transpose(Phi).Phi;
# Die symmetrische Matrix C enthält an der Stelle  C[i, j] das Skalarprodukt des i-ten  mit dem j-ten Eigenvektor der Matrix ¦.  
# Darstellung der Eigenvektoren im Raum
> EigenPlot(Theta, showsphere = false, showimagevectors = false, eigenoptions = [shape = harpoon]);
# Wir zeigen, dass die Eigenwertgleichung  ˜·x = »·x für jeden Eigenwert » erfüllt ist.
> for i to 3 do  
>   Theta.Phi[1 .. -1, i] = lambda[i].Phi[1 .. -1, i] 
> end do;
# 4. Berechnung der Diagonalmatrix T (Auf der Hauptdiagonalen stehen die Eigenwerte)
> T:=Transpose(Phi).Theta.Phi;
> 
;
