Vlastní výpočet módů

Vlastní problém malých kmitů:
	- výpočet poloh bodů v rovnovážné poloze (klasický “trojúhelníkový” most) ... Body
	- spojení jednotlivých bodů pružinami ... Pruziny
	- funkce potenciální energie pružiny a její druhá derivace ... Ugen, derUgen
	- funkce kinetické energie hmotné(!!!) pružiny a její druhá derivace ... Tgen, derTgen
	- výpočet matice U ... nejsem schopen příliš okomentovat. ;) sestavená po blocích 2x2 ... U
	- výpočet matice T ... sestavená po blocích 2x2 ... T
	- nalezení vlastních hodnot a vektorů ... EigenSystem
	- vykreslení

(* v systému jsou vždycky fixní (nedynamické) body na souřadnicích (0,0) a (1,0) ty budou "mít" indexy Pocet+1 a Pocet+2 (kde Pocet je celkový počet (dynamických) bodů) *)




(*********** Nastavitelné v rozmezí 1 - ∞. *)
pocetN=3; (* pocetN udává počet bodů na mostovce nepočítaje krajní body *)
(***********)




(* body jsou tvaru {rovnovazna_pozice, hmotnost} *)
Body=1.*Join[
Table[{{i(1/(pocetN+1)),0},1},{i,1,pocetN}],
Table[{{((i-1/2)/(pocetN+1)),(1/(pocetN+1))((Sqrt[3])/2)},1},{i,1,pocetN+1}],
{{{0,0},0},{{1,0},0}}
];
Pocet=Length[Body]-2;
(* pruziny jsou tvaru {{index_bodu_1,index_bodu_2},tuhost} - pružina upevněná na nedynamické body je musí mít uvedené v indexu1 !!! *)
Pruziny=Join[
Table[{{i-1,i},1.},{i,1,pocetN+1}],
Table[{{i,i+1},1.},{i,pocetN+1,2pocetN}],
Table[{{i-1,pocetN+i},1.},{i,1,pocetN+1}],
Table[{{i,pocetN+i},1.},{i,1,pocetN+1}]
];
Pruziny[[1]]={{2pocetN+2,1},1.};
Pruziny[[pocetN+1]]={{2pocetN+3,pocetN},1.};
Pruziny[[2pocetN+2]]={{2pocetN+2,pocetN+1},1.};
Pruziny[[4pocetN+3]]={{2pocetN+3,2pocetN+1},1.};

(* obecný tvar potenciální a kinetické energie pro hmotnou(!!!) pružinu mezi body (x[1],x[2]) a (x[3],x[4]) a její derivace *)
Ugen=(1/2)kk((((Sqrt[(((x[1]-x[3]))^2)+(((x[2]-x[4]))^2)])-l0))^2);
derUgen=Table[D[Ugen,x[i],x[j]],{i,1,4},{j,1,4}]//FullSimplify;
Tgen=(1/2)m((1/12)(((((v[4]-v[2])(x[3]-x[1])-(x[4]-x[2])(v[3]-v[1])))^2)/((((x[3]-x[1]))^2)+(((x[4]-x[2]))^2)))+(1/4)((((v[1]+v[3]))^2)+(((v[2]+v[4]))^2))+(1/12)(((((x[3]-x[1])(v[3]-v[1])+(x[4]-x[2])(v[4]-v[2])))^2)/((((x[3]-x[1]))^2)+(((x[4]-x[2]))^2))));
derTgen=Table[D[Tgen,v[i],v[j]],{i,1,4},{j,1,4}]//FullSimplify;

(* vygenerování matice potenciální energie *)
U=Table[0,{2Pocet},{2Pocet}];
For[i=1,i≤Length[Pruziny],i++,
Which[
Pruziny[[i,1,1]]==Pocet+1,
U[[2Pruziny[[i,1,2]]-1;;2Pruziny[[i,1,2]],2Pruziny[[i,1,2]]-1;;2Pruziny[[i,1,2]]]]+=derUgen[[3;;4,3;;4]]/.{x[1]->0,x[2]->0,x[3]->Body[[Pruziny[[i,1,2]],1,1]],x[4]->Body[[Pruziny[[i,1,2]],1,2]],kk->Pruziny[[i,2]],l0->(Sqrt[((Body[[Pruziny[[i,1,2]],1,1]])^2)+((Body[[Pruziny[[i,1,2]],1,2]])^2)])},
Pruziny[[i,1,1]]==Pocet+2,
U[[2Pruziny[[i,1,2]]-1;;2Pruziny[[i,1,2]],2Pruziny[[i,1,2]]-1;;2Pruziny[[i,1,2]]]]+=derUgen[[3;;4,3;;4]]/.{x[1]->1,x[2]->0,x[3]->Body[[Pruziny[[i,1,2]],1,1]],x[4]->Body[[Pruziny[[i,1,2]],1,2]],kk->Pruziny[[i,2]],l0->(Sqrt[(((Body[[Pruziny[[i,1,2]],1,1]]-1))^2)+((Body[[Pruziny[[i,1,2]],1,2]])^2)])},
True,
For[j=1,j≤2,j++,
For[k=1,k≤2,k++,
U[[2Pruziny[[i,1,j]]-1;;2Pruziny[[i,1,j]],2Pruziny[[i,1,k]]-1;;2Pruziny[[i,1,k]]]]+=
derUgen[[2j-1;;2j,2k-1;;2k]]/.{x[1]->Body[[Pruziny[[i,1,1]],1,1]],x[2]->Body[[Pruziny[[i,1,1]],1,2]],x[3]->Body[[Pruziny[[i,1,2]],1,1]],x[4]->Body[[Pruziny[[i,1,2]],1,2]],kk->Pruziny[[i,2]],l0->√((((Body[[Pruziny[[i,1,1]],1,1]]-Body[[Pruziny[[i,1,2]],1,1]]))^2)+(((Body[[Pruziny[[i,1,1]],1,2]]-Body[[Pruziny[[i,1,2]],1,2]]))^2))}
(*zavorka For j*)]
(*zavorka For k*)]
(*zavorka Which*)]
(*zavorka For i*)]

(* vygenerování matice kinetické energie *)
T=Table[0,{2Pocet},{2Pocet}];
For[i=1,i≤Length[Pruziny],i++,
Which[
Pruziny[[i,1,1]]==Pocet+1,
T[[2Pruziny[[i,1,2]]-1;;2Pruziny[[i,1,2]],2Pruziny[[i,1,2]]-1;;2Pruziny[[i,1,2]]]]+=derTgen[[3;;4,3;;4]]/.{x[1]->0,x[2]->0,x[3]->Body[[Pruziny[[i,1,2]],1,1]],x[4]->Body[[Pruziny[[i,1,2]],1,2]],m->Pruziny[[i,2]]},
Pruziny[[i,1,1]]==Pocet+2,
T[[2Pruziny[[i,1,2]]-1;;2Pruziny[[i,1,2]],2Pruziny[[i,1,2]]-1;;2Pruziny[[i,1,2]]]]+=derTgen[[3;;4,3;;4]]/.{x[1]->1,x[2]->0,x[3]->Body[[Pruziny[[i,1,2]],1,1]],x[4]->Body[[Pruziny[[i,1,2]],1,2]],m->Pruziny[[i,2]]},
True,
For[j=1,j≤2,j++,
For[k=1,k≤2,k++,
T[[2Pruziny[[i,1,j]]-1;;2Pruziny[[i,1,j]],2Pruziny[[i,1,k]]-1;;2Pruziny[[i,1,k]]]]+=
derTgen[[2j-1;;2j,2k-1;;2k]]/.{x[1]->Body[[Pruziny[[i,1,1]],1,1]],x[2]->Body[[Pruziny[[i,1,1]],1,2]],x[3]->Body[[Pruziny[[i,1,2]],1,1]],x[4]->Body[[Pruziny[[i,1,2]],1,2]],m->Pruziny[[i,2]]}
(*zavorka For j*)]
(*zavorka For k*)]
(*zavorka Which*)]
(*zavorka For i*)]

(* Nalezení vlastních vektorů a vlastních hodnot *)
EigenSystem=Eigensystem[U.Inverse[T]];

(* sestavení módů a grafických prvků *)
amp=0.06; (* amplituda módů *)
(* tabulka souřadnic bodů pro jednotlivé módy s časovou proměnnou t, tzn. KmitajiciBody[mod] obsahuje polohy všech bodů v závislosti na čase t *)
KmitajiciBody=Table[Join[Table[{Body[[i,1,1]]+amp EigenSystem[[2,mod,2i-1]]Cos[(Sqrt[EigenSystem[[1,mod]]])t],Body[[i,1,2]]+amp EigenSystem[[2,mod,2i]]Cos[(Sqrt[EigenSystem[[1,mod]]])t]},{i,1,Pocet}],{{0,0},{1,0}}],{mod,1,2Pocet}];
(* tabulka pružin směřující na aktuální polohy bodů pro jednotlivé módy, tzn. KmitajiciPruziny[mod] obsahují polohy koncových bodů pružin v závislosti na čase t*)
KmitajiciPruziny=Table[Table[{KmitajiciBody[[mod,Pruziny[[i,1,1]]]],KmitajiciBody[[mod,Pruziny[[i,1,2]]]]},{i,1,Length[Pruziny]}],{mod,1,2Pocet}];


arrowMultiplier=0.25;
(* Grafické znázornění poměrů amplitud jednotlivých módů *)
groupBy=3; (* počet grafů na řádek *)
GraphicsGrid[Partition[Table[Graphics[{
Text[Style["ω = "<>ToString[Round[(Sqrt[EigenSystem[[1,mod]]]),0.01]]<>"(Subscript[ω, 0])",Medium],{0.09,0.41}],
{PointSize[Large],Point[Transpose[Body][[1]]]},
{Dotted,Line[Map[Body[[#,1]]&,Transpose[Pruziny][[1]]]]},
{Arrowheads[0.03],Red,Arrow[Transpose[{Transpose[Body][[1,1;;-3]],Transpose[Body][[1,1;;-3]]+arrowMultiplier Partition[EigenSystem[[2,mod]],2]}]]},
{Arrowheads[0.03],Blue,
Arrow[Transpose[{Transpose[Body][[1,1;;-3]],Transpose[Body][[1,1;;-3]]-arrowMultiplier Partition[EigenSystem[[2,mod]],2]}]]}
},PlotRange->{{-0.05,1.05},{-0.25,0.55}}],{mod,1,2Pocet}],groupBy,groupBy,1,{}],ImageSize->800]


(* Animace pohybu mostu pro jednotlivé módy (plus zakreslené poměry amplitud v daném módu *)
Manipulate[GraphicsGrid[{{
Graphics[{
{PointSize[Large],Point[Transpose[Body][[1]]]},
{Dotted,Line[Map[Body[[#,1]]&,Transpose[Pruziny][[1]]]]},
{Arrowheads[0.02],Red,Arrow[Transpose[{Transpose[Body][[1,1;;-3]],Transpose[Body][[1,1;;-3]]+arrowMultiplier Partition[EigenSystem[[2,mod]],2]}]]},
{Arrowheads[0.02],Blue,
Arrow[Transpose[{Transpose[Body][[1,1;;-3]],Transpose[Body][[1,1;;-3]]-arrowMultiplier Partition[EigenSystem[[2,mod]],2]}]]}
},PlotRange->{{-0.1,1.1},{-0.18,0.48}}],
Show[
Graphics[{PointSize[Large],Point[KmitajiciBody[[mod]]]}],
Graphics[Line[KmitajiciPruziny[[mod]]]],
PlotRange->{{-0.1,1.1},{-0.18,0.48}}
]}},ImageSize->800,PlotLabel->Style[Dynamic["Mód č. "<>ToString[mod]<>", ω = "<>ToString[Round[(Sqrt[EigenSystem[[1,mod]]]),0.01]]<>"(Subscript[ω, 0])"],16]]/.{t->tt},{{mod,1,Dynamic["Mód č. "<>ToString[mod]]},1,2Pocet,1},{{tt,0,"čas t"},0,∞},SaveDefinitions->True]

Export do .gif a .png

frames=20;
plotVerticalMargin=0.05;
arrowMultiplier=0.25;
imageWidth=750;
plotRange={{-0.1,1.1},{-plotVerticalMargin,((Sqrt[3])/(2.(pocetN+1)))+plotVerticalMargin}};

mods=Table[m,{m,1,2Pocet}];
(*zajimave mody pro pocetN=10*) (*mods = {42,41,39,38,40,37,34,32,16,15,14,13,11};*)
mods={22};
(*mods={1,10};*)

grafT[mod_]:=Show[
Graphics[{PointSize[Large],Point[KmitajiciBody[[mod]]]}],
Graphics[Line[KmitajiciPruziny[[mod]]]],
PlotRange->plotRange,ImageSize->imageWidth
];

grafA[mod_]:=Graphics[{
{PointSize[Large],Point[Transpose[Body][[1]]]},
{Dotted,Line[Map[Body[[#,1]]&,Transpose[Pruziny][[1]]]]},
{Arrowheads[0.02],Red,Arrow[Transpose[{Transpose[Body][[1,1;;-3]],Transpose[Body][[1,1;;-3]]+arrowMultiplier Partition[EigenSystem[[2,mod]],2]}]]},
{Arrowheads[0.02],Blue,
Arrow[Transpose[{Transpose[Body][[1,1;;-3]],Transpose[Body][[1,1;;-3]]-arrowMultiplier Partition[EigenSystem[[2,mod]],2]}]]}
},PlotRange->plotRange,ImageSize->imageWidth];


For[i=1,i≤Length[mods],i++,
tabulka=Table[grafT[mods[[i]]]/.{t->tt},{tt,0,((2π)/(Sqrt[EigenSystem[[1,mods[[i]]]]])),(1/frames)((2π)/(Sqrt[EigenSystem[[1,mods[[i]]]]]))}];
Export["mostN"<>ToString[pocetN]<>"m"<>ToString[mods[[i]]]<>".gif",tabulka,"GIF","DisplayDurations"->{0.05},"AnimationRepetitions"->∞];
Export["mostN"<>ToString[pocetN]<>"m"<>ToString[mods[[i]]]<>".png",grafA[mods[[i]]],"PNG"];
]