Prosta vrtavka

Eulerjeve dinamične enečbe

eq={(Subscript[j, 1])(Subscript[ω, 1])'[t] - (Subscript[ω, 2])[t]  (Subscript[ω, 3])[t]((Subscript[j, 2])-(Subscript[j, 3])) == 0,(Subscript[j, 2])(Subscript[ω, 2])'[t] - (Subscript[ω, 3])[t]  (Subscript[ω, 1])[t]((Subscript[j, 3])-(Subscript[j, 1])) == 0,(Subscript[j, 3])(Subscript[ω, 3])'[t] - (Subscript[ω, 1])[t]  (Subscript[ω, 2])[t]((Subscript[j, 1])-(Subscript[j, 2])) == 0}

Primer: rešitev za začetni pogoj (Subscript[ω, 1])(t=0)=0.

Nastavek za rešitev zapišemo z Jacobijevimi eliptičnimi funkcijami JacobiSN[α t,m], JacobiCN[α t,m] in JacobiDN[α t,m], saj je odvod ene enak produktu ostalih dveh, funkcija JacobiSN[α t,m] pa tudi zadošča začetnemu pogoju.

Odvodi Jacobijevih eliptičnih funkcij

D[JacobiSN[α t,m],t]
D[JacobiCN[α t,m],t]
D[JacobiDN[α t,m],t]

Res, odvod ene je produkt drugih dveh

Reševanje sistema

Rešitev iščemo v obliki

rl ={
(Subscript[ω, 1])-> Function[t, (Subscript[a, 1])JacobiSN[α t,m]],
(Subscript[ω, 2])-> Function[t, (Subscript[a, 2])JacobiCN[α t,m]],
(Subscript[ω, 3])-> Function[t, (Subscript[a, 3])JacobiDN[α t,m]]};

Simplify[eq/.rl]

Dobimo algebraični sistem za neznane parametre (Subscript[a, i])

eqeq=Simplify[eq/.rl]/.{JacobiSN-> Function[{t,m},1],JacobiCN-> Function[{t,m},1],JacobiDN-> Function[{t,m},1]}

Privzeli bomo, da so vztrajnostni momenti med seboj različni. V primeru, če so enaki(simetrična prosta vrtavka), znamo rešitev poiskati s trigonometričnimi  funkcijami.

as = (Subscript[j, 1])<(Subscript[j, 2])<(Subscript[j, 3])

slv=Simplify[Solve[eqeq,{(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3])}],as]

Length[slv]

Prva rešitev je trivialna, ostale so si podobne do predznaka in vrstnega reda. Poglejmo na primeru.

slv/.{(Subscript[j, 1])-> 1,(Subscript[j, 2])-> 2,(Subscript[j, 3])-> 3}

Odločimo se za drugo rešitev. S pomočjo nje izrazimo iz začetnega pogoja vrednosti parametrov α in m.

slvαm=Solve[Thread[{(Subscript[ω, 2])[0],(Subscript[ω, 3])[0]}== {(Subscript[a, 2]),(Subscript[a, 3])}/.slv[[2]]],{α,m}]

Splošna rešitev

Dobili smo rešitev

vω=Simplify[{(Subscript[ω, 1])[t],(Subscript[ω, 2])[t],(Subscript[ω, 3])[t]}/.rl/.slv[[2]]/.slvαm[[1]],as]

Primer

primer = {(Subscript[j, 1])-> 1,(Subscript[j, 2])-> 2,(Subscript[j, 3])-> 3,(Subscript[ω, 1])[0]-> 0,(Subscript[ω, 2])[0]-> 1,(Subscript[ω, 3])[0]-> -1}

Plot[Evaluate[vω/.primer],{t,0,10}]

Gibanje je periodično. Period je

perioda =4EllipticK[m]/.slvαm[[1]]
perioda=perioda/.primer

Telesna vrtilna količina rešitve

Prostorska vrtilna količina je seveda konstantna, telesna pa je funkcija časa.

vk=ParametricPlot3D[Evaluate[DiagonalMatrix[{(Subscript[j, 1]),(Subscript[j, 2]),(Subscript[j, 3])}].vω/.primer],{t,0,perioda},PlotStyle-> Directive[Red,Thick]]

Konstanti gibanja sta kinetična energija in velikost vrtilne količine.

T=(1/2)((Subscript[j, 1])(((Subscript[ω, 1])[0])^2)+(Subscript[j, 2])(((Subscript[ω, 2])[0])^2)+(Subscript[j, 3])(((Subscript[ω, 3])[0])^2))/.primer

L = (Sqrt[((((Subscript[j, 1])(Subscript[ω, 1])[0]))^2)+((((Subscript[j, 2])(Subscript[ω, 2])[0]))^2)+((((Subscript[j, 3])(Subscript[ω, 3])[0]))^2)])/.primer

Animate[Show[vk,Graphics3D[{PointSize[Large],Point[Evaluate[DiagonalMatrix[{(Subscript[j, 1]),(Subscript[j, 2]),(Subscript[j, 3])}].vω/.primer]]},
PlotRange->{{-(Sqrt[2 (Subscript[j, 1])T]),(Sqrt[2 (Subscript[j, 1])T])},{-(Sqrt[2 (Subscript[j, 2])T]),(Sqrt[2 (Subscript[j, 2])T])},{-(Sqrt[2 (Subscript[j, 3])T]),(Sqrt[2 (Subscript[j, 3])T])}}/.primer,Boxed->False]/.t-> $t],{$t,0,perioda}]

Tir leži na preseku sfere s polmerom L in Binetovega elipsoda s polosmi (Sqrt[2 T (Subscript[j, i])])

elipsoid=Graphics3D[{Ellipsoid[{0,0,0},{(Sqrt[2 T (Subscript[j, 1])]),(Sqrt[2 T (Subscript[j, 2])]),(Sqrt[2 T (Subscript[j, 3])])}/.primer],Sphere[{0,0,0},L]},Boxed->False];

Show[elipsoid,vk]

Integracija vektorja kotne hitrosti

(mw = Table[-Sum[LeviCivitaTensor[3][[i,j,k]](Subscript[ω, k])[t],{k,3}],{i,3},{j,3}])//MatrixForm

(mq = Table[(Subscript[q, i,j])[t],{i,3},{j,3}])//MatrixForm

eqeqeq=Thread[ Flatten[D[mq,t]- mq.mw]== 0]

Splošne rešitve Mathematica ne more poiskati. Izračunajmo numerično rešitev za primer.

Numerična rešitev

eqeqeq=eqeqeq/.rl/.slv[[2]]/.slvαm[[1]]/.primer

Začetni pogoj na rotacijo

zp=Thread[Flatten[mq- IdentityMatrix[3]/.t-> 0]== 0]

np = 10;
nslv=NDSolve[Join[eqeqeq,zp],Flatten[mq],{t,0,np perioda},AccuracyGoal->12,PrecisionGoal->12]

Test natančnosti rešitve

Plot[Evaluate[Flatten[Transpose[mq].mq-IdentityMatrix[3]]/.nslv[[1]]],{t,0, np perioda},PlotRange->All]

Ali je rešitev periodična ?

Plot[Evaluate[Flatten[mq]/.nslv[[1]]],{t,0, perioda}]

Plot[Evaluate[Norm[Flatten[mq-IdentityMatrix[3]]]/.nslv[[1]]],{t,0,10perioda}]

Približna perioda je

Flatten[mq]/.nslv[[1]]/.t->5 perioda

Predstavitev rešitve

Rešitev predstavimo kot rotacijo telesnega koordinatnega sistema

e1=mq.{1,0,0}/.nslv[[1]];
e2=mq.{0,1,0}/.nslv[[1]];
e3=mq.{0,0,1}/.nslv[[1]];

d=1.5;
ks = Graphics3D[{Red,Arrow[{{0,0,0},{d,0,0}}],Arrow[{{0,0,0},{0,d,0}}],Arrow[{{0,0,0},{0,0,d}}]},PlotRange->{{-1.5,1.5},{-1.5,1.5},{-1.5,1.5}},Boxed->False];

Animate[Show[ks,Graphics3D[{Thick,Line[{{0,0,0},e1}],Line[{{0,0,0},e2}],Line[{{0,0,0},e3}]},PlotRange->{{-d,d},{-d,d},{-d,d}},Boxed->False]/.t-> $t],{$t,0,5perioda}]