# Kapitel 6_a:
# Die Lagrangeschen Bewegungsgleichungen
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plottools):
# Beispiel 6-1:
# 
# Es sind die Bewegungsgleichungen für den Zweimassenschwinger in der obigen Abbildung mithilfe der Lagrangeschen Bewegungsgleichungen aufzustellen. Für x1 = 0 und  x2 = 0  sind beide Federn entspannt. Beide Massen m1 und m2 sind reibungsfrei gelagert. 
# Die Beispiele orientieren sich an meinem  Buchtitel:
# F. U. Mathiak: Technische Mechanik 3, Kinematik und Kinetik  mit Maple- und MapleSim-Anwendungen, Verlag Walter de Gruyter, 2015
# Wir benötigen zur Aufstellung der Lagrange-Funktion  L = E - U die potenzielle und die kinetische Energie des Systems.
# 1. Die kinetische Energie:
> E:=1/2*m1*(diff(x1(t),t))^2+1/2*m2*(diff(x2(t),t))^2;
# 2. Die potenzielle Energie:
> U:=1/2*k1*x1(t)^2+1/2*k2*(x2(t)-x1(t))^2;
# Damit ergibt sich die Lagrange-Funktion zu
> L:=E-U;
# Zur Lösung des Variationsproblems verwenden wir die Prozedur EulerLagrange aus der Bibliothek  VariationalCalculus. Zusätzlich wird ein erstes Intgral der Lagrange-Funktion ausgegeben. 
> loe:=VariationalCalculus[EulerLagrange](L, t, [x1(t), x2(t)]);
# Wir entnehmen der Lösung  die beiden Differenzialgleichung 2. Ordnung:
> gl1:=loe[1]; gl2:=loe[2]; 
# Auf die Lösung dieses Gleichungssystem wird im Kapitel 7 näher eingegangen.
# Beispiel 6-2:
# 
# Für das in der obigen Abbildung skizzierte physische Doppelpendel sollen unter Zuhilfenahme des Lagrangeschen Formalismus die Bewegungsgleichungen aufgestellt werden. Das System besitzt mit den beiden Drehwinkel Æ1 und Æ2 genau zwei Freiheitsgrade. Zur Bestimmung des Potenzials der Gewichtskräfte G1 = m1 g  und G2 = m2 g  führen wir das Nullniveau (NN) bei x1 = 0 ein. Dann erhalten wir mit den Strecken  AB =    und  BS2 = s2 
> U:=-((G1*s1+G2*l)*cos(phi1(t))+G2*s2*cos(phi2(t)));
# Die kinetische Energie des Gesamtsystems ist
> E:=1/2*(Theta[A]*diff(phi1(t),t)^2 + m2*vS2q + Theta[S2]*diff(phi2(t),t)^2);
# mit
> vS2q:=(l*diff(phi1(t),t))^2+(s2*diff(phi2(t),t))^2+2*l*s2*diff(phi1(t),t)*diff(phi2(t),t)*cos(phi1(t)-phi2(t));
> L:=simplify(E-U);
> loe:=VariationalCalculus[EulerLagrange](L, t, [phi1(t), phi2(t)]);
> gl1:=loe[1]; gl2:=loe[2]; 
> gl:=gl1,gl2:
# Wir rechnen für dieses gekoppelte nichtlineare Differenzialgleichungssystem ein nummerisches Beispiel:
> l:=1; m1:=1; m2:=1; s1:=0.5*l;s2:=0.5*l;Theta[A]:=1/3; Theta[S2]:=1/12; G1:=9.81; G2:=9.81;
> AB:=phi1(0) = phi10, D(phi1)(0) = phi1p0, phi2(0) = phi20, D(phi2)(0) = phi2p0;
# Wir wählen folgende Anfangsbedingungen
> phi10:=Pi/2; phi1p0:=0.; phi20:=0.; phi2p0:=0.;
> GL:=gl1,gl2:
# Wir lösen die Aufgabe nummerisch, wobei aus Gründen der Genauigkeit die Schrittweite möglichst klein zu wählen ist.
> N:=500;  #Anzahl der Zeitschritte
;
> TE:=10.; #Schwingungszeit
;
> val:=Array(1..N+1): for i to N+1 do val[i]:=(i-1)*TE/N end do:
> LN:=dsolve({GL,AB},type=numeric,output=val); #Lösung des nichtlinearen Differenzialgleichungssystems
;
# Wir animieren den Bewegungsvorgang des Doppelpendels:
> Linie1  :=line([0, 0],[0,  -l],color = red,  thickness = 3):
> Linie2  :=line([0,-l],[0,-2*l],color = blue, thickness = 3):
> Scheibe1:=disk([0,-l],0.04*l,color = red):
> Scheibe2:=disk([0,-2*l],0.04*l,color = blue):
> NT:=LinearAlgebra[RowDimension](LN[2,1]): R:=NULL:
> for j to NT do    #Schleife über jeden Zeitschritt
> P1:= [l*sin(LN[2,1][j,2]),-l*cos(LN[2,1][j,2])];
> L1:= rotate(Linie1  ,LN[2,1][j,2],[0,0]):
> S1:= rotate(Scheibe1,LN[2,1][j,2],[0,0]):
> L2:= rotate(translate(Linie2  ,P1[1],P1[2]+l),LN[2,1][j,4],P1):
> S2:= rotate(translate(Scheibe2,P1[1],P1[2]+l),LN[2,1][j,4],P1):
> P := plots[display](L1,L2,S1,S2):
> R := R,P:
> end do:
> AN:=plots[display](R,insequence=true,axes=boxed,title = "Animation des Bewegungsvorganges für das Doppelpendel",titlefont = ["ARIAL", 15],scaling=constrained,gridlines=true):
> AN;
> para:= titlefont = ["ARIAL", 15],labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],gridlines = true,axes=boxed:
> q1:=plots[odeplot](LN,[t,phi1(t)],0..TE, para,title = " Winkel phi1\n",labels=["Zeit t", "phi1"],color=blue):
> q2:=plots[odeplot](LN,[t,diff(phi1(t),t)],0..TE,para,title ="Winkelgeschwindigkeit phi1'\n",labels = ["Zeit t","phi'"],color=blue):
> q3:=plots[odeplot](LN,[t,phi2(t)],0..TE,para,title = " Winkel phi2\n",labels = ["Zeit t", "phi"],color=blue):
> q4:=plots[odeplot](LN,[t,diff(phi2(t),t)],0..TE,para,title ="Winkelgeschwindigkeit phi2'\n",labels = ["Zeit t","phi2'"],color=blue):
>  A:= Array(1 .. 4): A[1]:= q1: A[2]:= q2: A[3]:= q3: A[4]:= q4:
> plots[display](A);
> 
;
