Lösung Beispiel 4, Helmut Hörner Mat.Nr. 8850092

Translatorische kinetische Energie der Masse 1

TT1=(1/2)m1 ((x1'[t])^2);

Translatorische kinetische Energie der Masse 2

TT2=(1/2)m2 ((v2[t])^(2 ))/. {((v2[t])^(2 ))->  ((D[x1[t] + l Sin[β[t]],t])^2)+((D[- l Cos[β[t]],t])^2)}

Rotatorische kinetische Energie der Masse 2

TR2=(1/2)II ((β'[t])^2) /. {II -> m2 (l^2)}

Potentielle Federenergie der Masse 1

U1= (1/2)k ((x1[t])^2)

Potentielle Gravitationsenergie der Masse 2

U2=m2 g h /. {h->-l Cos[β[t]]}

Gesamte kinetische Energie

T=TT1+TT2+TR2

Gesamte potentielle Energie

U = U1 + U2

Lagrange Funktion:

L = T - U

Definiere konkrete Werte (u.A.  m1=0,5 und m2=1)

Werte1={g->10, l->5, k->1, m1-> 0.5, m2->1};

Anfangsbedingungen

AB={x1[0]==0, β[0]==(π/4), x1'[0]==0, β'[0]==0};

Definiere Lagrange-DGL

DGL1=D[D[L,x1'[t]],t]-D[L,x1[t]]==0;
DGL2=D[D[L,β'[t]],t]-D[L,β[t]]==0;

Löse DGL numerisch

Solution1=NDSolve[{DGL1,DGL2,AB} //. Werte1,{x1[t],β[t]},{t,0,200}]

Animation mit m1=0,5 und m2=1

Animate[
Graphics[
{Line[{{-7,0},{10,0}}],
Line[{{x1[t]/.First[Solution1] /. t->tt,0},
{(x1[t]/.First[Solution1] /. t->tt)+l*Sin[β[t]/.First[Solution1] /. t->tt],
-l*Cos[β[t]/.First[Solution1] /. t->tt]}}] ,
Rectangle[{x1[t]-0.5/.First[Solution1] /. t->tt,0.7},{x1[t]+0.5/.First[Solution1] /. t->tt,0}],
Disk[{(x1[t]/.First[Solution1] /. t->tt)+l*Sin[β[t]/.First[Solution1] /. t->tt],-l*Cos[β[t]/.First[Solution1] /. t->tt]},0.5`]},

PlotRange->{{-7,10},{-7,7}}],
{tt,0,200},AnimationRate->3,AnimationRunning->False]/. Werte1

Definiere neue Werte (u.A.  m1=5 und m2=0.1)

Werte2={g->10, l->5, k->1, m1-> 5, m2->0.1};

Löse DGL mit neuen Werten numerisch

Solution2=NDSolve[{DGL1,DGL2,AB} //. Werte2,{x1[t],β[t]},{t,0,200}]

Animation mit m1=5 und m2=0,1

Animate[
Graphics[
{Line[{{-7,0},{10,0}}],
Line[{{x1[t]/.First[Solution2] /. t->tt,0},
{(x1[t]/.First[Solution2] /. t->tt)+l*Sin[β[t]/.First[Solution2] /. t->tt],
-l*Cos[β[t]/.First[Solution2] /. t->tt]}}] ,
Rectangle[{x1[t]-0.5/.First[Solution2] /. t->tt,0.7},{x1[t]+0.5/.First[Solution2] /. t->tt,0}],
Disk[{(x1[t]/.First[Solution2] /. t->tt)+l*Sin[β[t]/.First[Solution2] /. t->tt],-l*Cos[β[t]/.First[Solution2] /. t->tt]},0.5`]},

PlotRange->{{-7,10},{-7,7}}],
{tt,0,200},AnimationRate->3,AnimationRunning->False]/. Werte2

Linearisieren der DGL:

DGL1LIN=FullSimplify[Expand[DGL1 /. {Sin[β[t]]-> β[t], Cos[β[t]]-> 1} ]/.{((β[t])^2)-> 0,((β'[t])^2)-> 0} /. {β[t]-> (((x2[t]-x1[t]))/l), β''[t]-> D[(((x2[t]-x1[t]))/l),{t,2}]}]
DGL2LIN=FullSimplify[Expand[DGL2/. {Sin[β[t]]-> β[t], Cos[β[t]]-> 1} ]/.{((β[t])^2)-> 0, ((β'[t])^2)-> 0}/. {β[t]-> (((x2[t]-x1[t]))/l), β''[t]-> D[(((x2[t]-x1[t]))/l),{t,2}]},Assumptions->m2>0]

Definiere neue Werte für linearisierte DGL

Werte3={g->9.8, l->10, k->0.2, m1-> 1, m2->2};

Exakte Lösung der linearisierten DGLs

Solution3=DSolve[{DGL1LIN,DGL2LIN,AB /. {β[0]-> (((x2[t]-x1[t]))/l)/. t->0, β'[0]-> D[(((x2[t]-x1[t]))/l),t]/. t->0}}//. Werte3,{x1[t],x2[t]},t] // FullSimplify;

Animation mit linearisierten DGLs

Animate[
Graphics[
{Line[{{-11,0},{11,0}}],
Line[{{x1[t]/.First[Solution3] /. t->tt,0},{x2[t]/.First[Solution3] /. t->tt,-l*Cos[(((x2[t]/.First[Solution3])-(x1[t]/.First[Solution3]))/l) /. t->tt]}}],
Rectangle[{x1[t]-0.5/.First[Solution3] /. t->tt,0.7},{x1[t]+0.5/.First[Solution3] /. t->tt,0}],
Disk[{x2[t]/.First[Solution3] /. t->tt,-l*Cos[(((x2[t]/.First[Solution3])-(x1[t]/.First[Solution3]))/l) /. t->tt]},0.5`]},

PlotRange->{{-11,10},{-11,7}}],
{tt,0,200},AnimationRate->3,AnimationRunning->False]/. Werte3

Eigenschwingungen: 
Normalmodenansatz x1[t]=a Sin[ω t]; x2[t]=b Sin[ω t]

gl1=DGL1LIN /.{x1[t]->a Sin[ω t], x1''[t]->-a (ω^2) Sin[ω t],x2[t]->b Sin[ω t], x2''[t]->-b (ω^2) Sin[ω t]}
gl2=DGL2LIN/.{x1[t]->a Sin[ω t], x1''[t]->-a (ω^2) Sin[ω t],x2[t]->b Sin[ω t], x2''[t]->-b (ω^2) Sin[ω t]}

Bestimme Koeffizienten

gl1=gl1[[1]]/Sin[ω t] //FullSimplify //ExpandAll
gl2=gl2[[1]]/Sin[ω t] -gl2[[2]]/Sin[ω t]//FullSimplify//ExpandAll

Bestimme Eigenmoden

det=Solve[Det[(({{Coefficient[gl1,a], Coefficient[gl1,b]}, {Coefficient[gl2,a], Coefficient[gl2,b]}}))]==0,ω]

(Nur positive Frequenzen berücksichtigen)

(Subscript[ω, 1])=det[[2,1,2]]
(Subscript[ω, 2])=det[[4,1,2]]

(Subscript[a, 1])=Solve[{gl1==0,gl2==0}/.{ω-> (Subscript[ω, 1]),b->1},a][[1,1,2]]
(Subscript[a, 2])=Solve[{gl1==0,gl2==0}/.{ω-> (Subscript[ω, 2]),b->1},a][[1,1,2]]

Grafische Darstellung

Plot[{(Subscript[a, 1]) Sin[(Subscript[ω, 1]) t] ,Sin[(Subscript[ω, 1]) t]}/.Werte3,{t,0,100}]

Plot[{(Subscript[a, 2]) Sin[(Subscript[ω, 2]) t] /.Werte3,Sin[(Subscript[ω, 2]) t]/.Werte3},{t,0,100}]