Bąk Lagrange’a

Wyprowadzenie modelu

Needs["VariationalMethods`"]
(* współrzędne konfiguracyjne *)
variables={α[t],β[t],γ[t]};
(* lagranżian i równania *)
lagrangian=1/2 mI1 (((β'[t])^2)+((α'[t])^2) ((Sin[β[t]])^2))+1/2 mI3(((γ'[t]+α'[t] Cos[β[t]]))^2)- m g r Cos[β[t]];
eqs=EulerEquations[lagrangian,variables,t](* położenia środka ciężkości i innych plus wizualizacja bąka *)
rotZ[α_]=RotationMatrix[α,UnitVector[3,3]];  (* macierze rotacji elementarnych *)rotY[α_]=RotationMatrix[α,UnitVector[3,2]];
rotStoB=rotZ[α[t]-Pi/2].rotY[β[t]].rotZ[γ[t]+Pi/2]; (* transformacja ukł. S w B *)
rPosB={0,0,r};  (* środek ciężkości *)
rPosS[t_]=rotStoB.rPosB;
base={0,0,0};  (* wierzchołek bąka *)
pBaseB={rBase,0,r};  (* punkt na krawędzi podstawy - rBase promień podstawy *)
pBaseS=rotStoB.pBaseB;
lagrCone={Cone[{rPosS[t],base},rBase],Line[{base,3/2rPosS[t]}],Lighter[Yellow],Line[{pBaseS,base}],Line[{pBaseS/.rBase->-rBase,pBaseS}]};

Symulacja

(* parametry, warunki początkowe, czas symulacji *)
params={r->1,rBase->1/5,mI1->1,mI3->0.8,m->2,g->9.81};
initial={α[0]==0,β[0]==0.2,γ[0]==0,α'[0]==0,β'[0]==0,γ'[0]==0};
tMax=10;
(* rozwiązujemy i wizualizujemy *)
whichVar={1,2,3};
sol=First[NDSolve[{eqs,initial}/.params,Head/@variables,{t,0,tMax}]];
Plot[Evaluate[variables[[whichVar]]/.sol],{t,0,tMax},PlotLegends->variables[[whichVar]]]

nFrames=500;trace[t_]:={Green,Line[Map[Function[Evaluate[3/2rPosS[#]/.params/.sol]],Range[0,t,0.025]]]} 
frames=Table[Graphics3D[{Gray,Thick,lagrCone,trace[t]}/.params/.sol,PlotRange->{{-1.6,1.6},{-1.6,1.6},{0,2}},ImageSize->600, Background->Lighter[White],Boxed->False,Axes->True,AxesLabel->{"x","y","z"}],{t,tMax/nFrames,tMax,tMax/nFrames}];
ListAnimate[frames]

Dodajmy tarcie tylko na wirowaniu

(* W trzecim równaniu proporcjonalne do prędkości wirowania *)
eqs[[3]]
eqs[[3,2]]=k γ'[t];
eqs[[3]]

Dodajmy tarcie na pozostałych dwóch ,,osiach”

(* Podobnie jak poprzednio w dwóch pozostałych *)
eqs[[1,2]]=k α'[t];
eqs[[2,2]]=k β'[t];

Symulacja

(* parametry, warunki początkowe, czas symulacji *)
params={r->1,rBase->1/5,mI1->1,mI3->0.8,m->2,g->9.81,k->2};
initial={α[0]==0,β[0]==0.2,γ[0]==0,α'[0]==0,β'[0]==0,γ'[0]==20};
tMax=10;
(* rozwiązujemy i wizualizujemy *)
whichVar={1,2,3};
sol=First[NDSolve[{eqs,initial}/.params,Head/@variables,{t,0,tMax}]];
Plot[Evaluate[variables[[whichVar]]/.sol],{t,0,tMax},PlotLegends->variables[[whichVar]]]

nFrames=500;
frames=Table[Graphics3D[{Gray,Thick,lagrCone,trace[t]}/.params/.sol,PlotRange->{{-1.6,1.6},{-1.6,1.6},{-1.6,2}},ImageSize->600, Background->Lighter[White],Boxed->False,Axes->True,AxesLabel->{"x","y","z"}],{t,tMax/nFrames,tMax,tMax/nFrames}];
ListAnimate[frames]