# Kapitel 7_h:
# Freie viskos gedämpfte Schwingungen
# mit endlich vielen Freiheitsgraden
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> interface(displayprecision =4): 
> restart: with(LinearAlgebra): with(plots): with(plottools):
# Hinweis: Zur Berechnung der Zustandsgrößen sind die Wurzeln M-1/2 bzw. M1/2 einer symmetrischen Matrix ( hier Massenmatrix) zu berechnen (s.h. Buch TM3). Diese Ausdrücke sind als Matrizenfunktionen zu verstehen, zu deren Berechnung die Eigenwerte mj (j = 1,…,n) und Eigenvektoren ej aus dem speziellen Eigenwertproblem (M - m1)e benötigt werden. Wir stellen dazu im Vorfeld folgende Prozedur zur Verfügung:
# 
> Proc_Calc_30:=proc(M::Matrix)
> #------------------------------------------------------
# #  Eingabe:
# #       M : Symmetrische und positiv definite Matrix
> 
> #  Ausgabe:
> #      Mhp: Matrix M 1/2
> #      Mhm: Matrix M-1/2
> #------------------------------------------------------
> local  DET,n,DMpeh,DMmeh,EW,Phi,PhiI,DM,i;
> global Mhp, Mhm;
> description "Berechnung von M 1/2 und M -1/2 einer symmetrischen positiv definiten Matrix M";
> DET:=LinearAlgebra[Determinant](M);
> if abs(DET) < 10.^(-Digits+2) then
>   print(`Die Matrix M ist singulär, Abbruch der Berechnung`);
>   return;
> end if;
> n     :=LinearAlgebra[RowDimension](M);
> DMpeh :=Matrix(n,n);
> DMmeh :=Matrix(n,n);
> EW,Phi:=Eigenvectors(M);
> EW:=Re(EW); Phi:=Re(Phi);
> PhiI:= MatrixInverse(Phi);
> DM:=PhiI.M.Phi;
> for i to n do
>   DMpeh[i,i]:= sqrt(DM[i,i]);
>   DMmeh[i,i]:=  1/DMpeh[i,i];
> end do;
> Mhp:=Phi.DMpeh.PhiI;
> Mhm:=Phi.DMmeh.PhiI;
> print(`M 1/2  = `, Mhp);
> print(`M-1/2  = `, Mhm);
> end proc:
# Wir testen die Prozedur und wählen
> MT:= Matrix(4,{(1, 1)= 2.,(1, 2) =  7.,(1, 3) =  2.,(1, 4) = 1.,
>                           (2, 2) = 15.,(2, 3) =  4.,(2, 4) = 1.,
>                                        (3, 3) =  6.,(3, 4) = 1.,
>                                                     (4, 4) = 16.}, shape = symmetric);
> Proc_Calc_30(MT);
# Das Matrixprodukt  M 1/2 . M -1/2  muss die Einheitsmatrix ergeben. Wir kontrollieren das:
> Mhp.Mhm;
# Beispiel 7-7:
# 
# Für das in der obigen Abbildung skizzierte System sind sämtliche Zustandsgrößen zu berechnen. Die Massen m1 und m2 erhalten zum Zeitpunkt t = 0 die Anfangsauslenkungen x10 = 1 cm und x20 = 5 cm sowie die Anfangsgeschwindigkeiten v10 = 1 m/s und v20 = -1 m/s. 
> M:=Matrix([[m1,0],[0,m2]]);        #Massenmatrix
;
> C:=Matrix([[c1+c2,-c2],[-c2,c2]]); #Dämpfungsmatrix
;
> K:=Matrix([[k1+k2,-k2],[-k2,k2]]); #Steifigkeitsmatrix
;
# Entsprechend der Aufgabenstellung werden folgende Systemwerte gewählt:
> m1:=20.; c1:=10.; k1:=1600.;
> m2:= 5.; c2:= 3.; k2:=1000.;
> M; C; K;
> n:=LinearAlgebra[RowDimension](M):
# Wir gehen in Schritten vor:
# 1. Berechnung von  M -1/2 und M 1/2
> Proc_Calc_30(M);
# 2. Berechnung von ” und ©2
> Delta:=0.5*Mhm.C.Mhm; OMQ:=Mhm.K.Mhm;
# 3. Berechnung der Eigenwerte
> E:=IdentityMatrix(2):
> DSM:=zeta^2*E + 2*zeta*Delta + OMQ;
> Poly:=factor(LinearAlgebra[Determinant](DSM));
> EW:=solve(Poly,zeta);
# 4. Berechnung der Eigenvektoren
> Phi:=Matrix(n,2*n):
> for i to 4 do
>   DSMk:=evalf(subs(zeta=EW[i],DSM)):
>   Ek:=Vector([a,1]):
>   gl:=DSMk.Ek: 
>   a:=solve(gl[1],a):
>   Ek:=LinearAlgebra[Normalize](Ek,Euclidean);
>     for m from 1 to 2 do
>       Phi[m,i]:=Ek[m];
>     end do;
>   a:='a': b:=1:
> end do:
> Phi;
# 5. Transformation der Anfangswerte
> qd0:=Vector([0.01,0.05]); qdp0:=Vector([1.,-1.]);
> pd0:=Mhp.qd0; pdp0:=Mhp.qdp0;
# 6. Berechnung der Matrizen Z0 und Z1
> W1:=Matrix(n): Z1:=Matrix(n):
> for i to 2*n do
>   CV := LinearAlgebra[Column](Phi,i);
>   CVT:= LinearAlgebra[Transpose](CV);
>   W1S:= CV.CVT/(CVT.OMQ.CV-EW[i]*EW[i]*(CVT.CV))*exp(EW[i]*t);
>   W1 := W1 + W1S;
>   Z1 := Z1 - EW[i]*W1S;
> end do:
> Z0:=W1.OMQ:
# 7. Berechnung von pd
> pd:=evalc(Z0.pd0 + Z1.pdp0);
# 8. Rücktransformation in physikalische Koordinaten: x(t) = M -1/2 Åpd
> x:=Mhm.pd;         #Auslenkung
;
> v:=map(diff,x,t);  #Geschwindigkeit
;
> 
;
# Grafische Ausgabe der Ergebnisse
> para:=axes=boxed,gridlines, titlefont = ["ARIAL", 15], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont=["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom]:
> pl1:=plot([x[1],x[2]],t=0..1,para,title = "\nAuslenkungen\n", legend = ["x1(t)", "x2(t)"], labels = ["Zeit t [s]", "Verschiebungen [m]"],color=[black,blue]):
> pl2:=plot([v[1],v[2]],t=0..1,para,title = "\nGeschwindigkeiten\n", legend = ["v1(t)", "v2(t)"], labels = ["Zeit t [s]", "Geschwindigkeiten [m/s]"],color=[black,blue]):
# Darstellung der Schwingungen im Phasenraum
> pl3:=plot([x[1],x[2],t=0..5],para,title = "\nPfad in der (x1,x2)-Ebene\n", labels = ["x1 [m]"  , "x2 [m]"]):
> pl4:=plot([v[1],v[2],t=0..5],para,title = "\nPfad in der (v1,v2)-Ebene\n", labels = ["v1 [m/s]", "v2 [m/s]"]):
> plots[display](Matrix(1,4,[pl1,pl2,pl3,pl4]));
