(* Aby przeliczyć automatycznie wybierz Menu/Evaluation/Evaluate notebook *)

Graphics[
{
Line[{{-0.45422694912576017`,0.9075932007059151`},{0.45422694912576006`,0.9075932007059151`},{0.`,1.6835338811872702`},{-0.45422694912576017`,0.9075932007059151`},{-1,0},{0,0},{-0.45422694912576017`,0.9075932007059151`},{0,0},{0.45422694912576006`,0.9075932007059151`},{1,0},{0,0}}],
{PointSize[Large],Point[{-0.45422694912576017`,0.9075932007059151`}],Point[{0.45422694912576006`,0.9075932007059151`}],Point[{0.`,1.6835338811872702`}]},
{PointSize[Large],Text["((Subscript[x, 1]),(Subscript[y, 1]))",{-0.62,0.9075932007059151`}],Text["((Subscript[x, 2]),(Subscript[y, 2]))",{0.65,0.9075932007059151`}],Text["((Subscript[x, 0]),(Subscript[y, 0]))",{0.`,1.7835338811872702`}]},
{EdgeForm[Thick],FaceForm[LightGray],
Rectangle[{-3/2,0},{3/2,-1/3}]}
}
]

(* Celem zadania jest zbadanie małych drgań układu z powyzszego rysunku. Trzy identyczne masy m są połaczone 7 sprężynami i długości spoczynkowej d i stałej sprężystości k. Cztery dolne sprężyny są przymocowane do ,,ściany''. *)

q={x0[t],y0[t],x1[t],y1[t],x2[t],y2[t]} (* Współrzędne uogólnione*)

(OverDot[q])=D[q,t] (* Prędkości uogólnione; kropka to tylko symbol graficzny OverDot[], operator różniczkowania to ', czyli Derivative[1]*)

num=Length[q] (* Ilość stopni swobody *)

T=1/2 m (OverDot[q]).(OverDot[q]) (* Energia kinetyczna *)

U=
1/2 k (Sqrt[ (x1[t]+a)^2+y1[t]^2]-d)^2+
1/2 k (Sqrt[(x2[t]-a)^2+y2[t]^2]-d)^2+
1/2 k (Sqrt[x1[t]^2+y1[t]^2]-d)^2+
1/2 k (Sqrt[x2[t]^2+y2[t]^2]-d)^2+
1/2 k (Sqrt[ (x1[t]-x2[t])^2+(y1[t]-y2[t])^2]-d)^2+
1/2 k (Sqrt[ (x1[t]-x0[t])^2+(y1[t]-y0[t])^2]-d)^2+
1/2 k (Sqrt[ (x0[t]-x2[t])^2+(y0[t]-y2[t])^2]-d)^2

(* Energia potencjalna układu 7 sprężyn *)

sysMIN=Table[D[U,q[[i]]]==0,{i,1,num}] (* Warunek konieczny istnienia minimum potencjału, z punktu widzenia czysto matematycznego pewnej funkcji wielu zmiennych *)

min={x0->(0&),y0->(H&),x1->(-b&),x2->(b&),y1->(h&),y2->(h&)} (* szukamy minimum w postaci pewnego stanu symetrycznego (funkcji stałych ) *)

a=d (* Dla uproszczenia zakładam, że podstawa ,,choinki'' to 2 długości sprężyny *)

sysMIN2=Simplify[sysMIN/.min,b>0] (* Konwersja układu r. funkcyjnych na algebraiczne *)

trojkat={h->d Sqrt[3]/2,H->d Sqrt[3],b->d/2} (* ,,Odgaduję, że minimum (prawdopodobnie globalne) to ,,trójkąt równoboczny''.
*)

sysMIN2/.{h->d Sqrt[3]/2,H->d Sqrt[3],b->d/2} (* Sprawdzam powyższy ,,ansatz''*)

Simplify[%,d>0]

trojkatFun=min/.{h->d Sqrt[3]/2,H->d Sqrt[3],b->d/2} (* minimum lokalne (tr. równoboczny) zapisane jako funkcje *)

jac=Simplify[Table[D[U,q[[i]],q[[j]]],{i,1,num},{j,1,num}]/.trojkatFun,d>0];
(* Obliczam macierz drugich pochodnych, tzw. jakobian *)

Det[jac] (* Wyznacznik jakobianu powinien być dodatni, jeżeli w danym punkcie faktycznie jest minimum *)

k MatrixForm[jac/k] (* Jakobian, pierwsza z macierzy niezbędna do wyznaczenia małych drgań *)

Collect[1/2q.jac.q,k] (* Energia potencjalna w przybliżeniu małych drgań *)

kin=Table[D[T,D[q[[i]],t],D[q[[j]],t]],{i,1,num},{j,1,num}]/.trojkatFun;(* Analogicznie obliczam macierz ,,energii kinetycznej'' *)

Collect[1/2(OverDot[q]).kin.(OverDot[q]),m](* W tym wypadku energia kinetyczna nie zmienia się *)

m MatrixForm[kin/m] (* Macierz energii kinetycznej; w badanym przypadku proporcjonalna do macierzy jednostkowej *)

(* Ze względu na fakt, że macierz energii kinetycznej jest ,,jednostkowa'', uogólniony problem własny dla dwóch macierzy ,,kin'' i ,,jac'', sprowadza się do zwykłego problemu własnego macierzy ,,jac''. Co więcej, widać że, częstości własne są proporcjonalne do Sqrt[k/m]. Dlatego bez utraty ogólności można przyjąć k=1 i m=1.
*)

Det[-ω^2 kin+jac]==0 (* Formalne równanie charakterystyczne dla problemu małych drgań *)

CharacteristicPolynomial[jac/k,λ] (* Wielomian charakterystyczny uproszczonego problemu własnego ( λ=ω^2 ) *)

Solve[%==0,λ] (* rozwiązanie daje nam kwadraty częstości własnych *)

N[%]

Chop[%] (* liczby te są oczywiście rzeczywiste, pojawiające się błędne wartości urojone  są na poziomie precyzji procesora *)

{wartWlasne, wektoryWlasne}=Eigensystem[jac/k]; (* Kompletne rozwiązanie problemu własnego dla macierzy jac; aby je wyswietlić nalezy usunać średnik *)

wartWlasne//N (* Wartości własne są juz posortowane *)

Table[wektoryWlasne[[i]].wektoryWlasne[[j]]//FullSimplify,{i,1,num},{j,1,num}]//MatrixForm
(*Wektory własne są ortogonalne, ale nie są znormalizowane*)

N[%]

d=1 (* Ustalam długość sprężyny jako jednostkę, aby pokazać drgania własne graficznie *)

A=0.1 (* Dowolna (ale ,,mała'')  amplituda małych drgań *)

point=
Function[{np,ii,t},Partition[(q/.trojkatFun)+A wektoryWlasne[[ii]]/Sqrt[ wektoryWlasne[[ii]]. wektoryWlasne[[ii]]]Cos[Sqrt[wartWlasne[[ii]]]t],2][[np]]
]

point[1,1,1]

Manipulate[
ListPlot[Partition[(q/.trojkatFun)+A wektoryWlasne[[i]]/Sqrt[ wektoryWlasne[[i]]. wektoryWlasne[[i]]]Cos[Sqrt[wartWlasne[[i]]]t],2],PlotStyle->{PointSize[Large]},PlotRange->{{-d,d},{-d,2d}},AspectRatio->1,
Epilog->Line[
{
point[2,i,t],
point[3,i,t],
point[1,i,t],
point[2,i,t],
{-d,0},
{0,0},
point[2,i,t],
{0,0},
point[3,i,t],
{d,0},
{0,0}
}
](*Line neds here*)
],
{i,{1,2,3,4,5,6}},
{t,0,20}]

Manipulate[
ListPlot[Partition[(q/.trojkatFun)+A wektoryWlasne[[i]]/Sqrt[ wektoryWlasne[[i]]. wektoryWlasne[[i]]]Cos[Sqrt[wartWlasne[[i]]]t],2],PlotStyle->{PointSize[Large]},PlotRange->{{-d,d},{-d,2d}},AspectRatio->1],
{i,{1,2,3,4,5,6}},
{t,0,20}]

A=0.2

Manipulate[
TableForm@{Table[
ListPlot[Partition[(q/.trojkatFun)+A wektoryWlasne[[i]]/Sqrt[ wektoryWlasne[[i]]. wektoryWlasne[[i]]]Cos[Sqrt[wartWlasne[[i]]]t],2],PlotStyle->{PointSize[Large]},PlotRange->{{-d,d},{-d,2d}},AspectRatio->1,
Axes->False,
Epilog->Line[
{
point[2,i,t],
point[3,i,t],
point[1,i,t],
point[2,i,t],
{-d,0},
{0,0},
point[2,i,t],
{0,0},
point[3,i,t],
{d,0},
{0,0}
}
](*Line neds here*)
],
{i,3}
],
Table[
ListPlot[Partition[(q/.trojkatFun)+A wektoryWlasne[[i]]/Sqrt[ wektoryWlasne[[i]]. wektoryWlasne[[i]]]Cos[Sqrt[wartWlasne[[i]]]t],2],PlotStyle->{PointSize[Large]},PlotRange->{{-d,d},{-d,2d}},AspectRatio->1,
Axes->False,
Epilog->Line[
{
point[2,i,t],
point[3,i,t],
point[1,i,t],
point[2,i,t],
{-d,0},
{0,0},
point[2,i,t],
{0,0},
point[3,i,t],
{d,0},
{0,0}
}
](*Line neds here*)
],
{i,4,6}
]},
{t,0,120}]