# Kapitel 7_j:
# Erzwungene ungedämpfte Schwingungen
# mit endlich vielen Freiheitsgraden
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> interface(displayprecision =4): restart: with(LinearAlgebra):
# Beispiel 7-8:
# 
# Es soll eine Maple-Prozedur zur Verfügung gestellt werden, mit der die erzwungenen ungedämpften Bewegungen des Zweimassenschwingers in der obigen Abbildung automatisiert berechnet werden können. Sämtliche Zustandsgrößen sind grafisch auszugeben
# Massen- und Steifigkeitsmatrix
> M:=Matrix([[m[1], 0],[ 0,m[2]]]);
> K:=Matrix([[k[1]+k[2],-k[2]],[-k[2],k[2]]]);
> m[1]:=10.; m[2]:=5.; k[1] :=17.; k[2]:=3.;
> M;K;
# Lösung des allgemeinen Eigenwertproblems
> EW,EV:=LinearAlgebra[Eigenvectors](K,M);
> EW:=Re(EW);EV:=Re(EV);
# Normierung der Eigenvektoren auf die Länge 1
> e1 := LinearAlgebra[Normalize]((LinearAlgebra[Column](EV,1),Euclidean));
> e2 := LinearAlgebra[Normalize]((LinearAlgebra[Column](EV,2),Euclidean));
# Berechnung der Eigenkreisfrequenzen
> omega[1]:=sqrt(Transpose(e1).K.e1/(Transpose(e1).M.e1));
> omega[2]:=sqrt(Transpose(e2).K.e2/(Transpose(e2).M.e2));
# Aufstellen des Kraftvektors
> f:=Vector([A1*cos(Omega1*t)+B1*sin(Omega1*t),A2*cos(Omega2*t)+B2*sin(Omega2*t)]);
# Entkopplung der Bewegungsgleichungen:
> fq1:=Transpose(e1).f/(Transpose(e1).M.e1);
> fq2:=Transpose(e2).f/(Transpose(e2).M.e2);
# Krafterregung:
> A1:=0; B1:=3; Omega1:=1;
> A2:=0; B2:=5; Omega2:=2;
> fq1; fq2;
# Aufstellen der entkoppelten Bewegungsgleichungen
> ode1:=diff(h1(t),t,t)+omega[1]^2*h1(t)-fq1;
> ode2:=diff(h2(t),t,t)+omega[2]^2*h2(t)-fq2;
# Anfangsbedingungen
> AB:=h1(0)=0,h2(0)=0,D(h1)(0)=0,D(h2)(0)=0;
> sys:={ode1,ode2,AB};
> dsolve(sys):assign(%);
> h1:=unapply(evalf(h1(t)),t);h2:=unapply(evalf(h2(t)),t);
# Damit folgen die Verschiebungen
> qd:=unapply(h1(t)*e1 + h2(t)*e2,t):
> qd(t)[1];qd(t)[2];
# und durch Ableitung nach der Zeit die Geschwindigkeiten
> vd:=unapply(map(diff,qd(t),t),t): #Geschwindigkeiten
;
> vd(t)[1];vd(t)[2];
> qd(0);vd(0);
# Grafische Ausgabe der Ergebnisse
> tE:=20:
> para:=titlefont = ["ARIAL", 15], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], linestyle = [solid, solid], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom],axes=boxed,color=[black,blue],gridlines=true:
> p1:=plot([qd(t)[1],qd(t)[2]],t=0..tE,para,title = "\nAuslenkungen\n",legend = ["x1", "x2"], labels = ["Zeit t", "Auslenkungen"]):
> p2:=plot([vd(t)[1],vd(t)[2]],t=0..tE ,para,title = "\nGeschwindigkeiten\n", legend = ["v1", "v2"], labels = ["Zeit t", "Geschwindigkeiten"]):
# Darstellung der Bewegung in der Phasenebene
> p3:=plot([qd(t)[1],qd(t)[2] ,t=0..tE],para,title = "\nPfad in der (x1,x2)-Ebene\n", labels = ["x1", "x2"]):
> A:=Array(1..3): A[1]:=p1: A[2]:=p2: A[3]:=p3: 
> plots[display](A);
> 
;
