ClearAll["Global`*"]

(* Pręt ślizgający sie po osiach układu współrzędnych *)

L=2;alfa=Pi/6;

(* rysunek sprawdza parametryzacje dowolnego punktu preta *)

x1=(L-l)*Cos[alfa];y1=l*Sin[alfa];
P1=ParametricPlot[{x1,y1},{l,0,L},PlotRange->{{0,2.1},{0,1.1}}]

L=.;alfa=.;

x=(L-l)*Cos[ϕ[t]];y=l*Sin[ϕ[t]];
xp=D[x,t];yp=D[y,t];
vp2=xp*xp+yp*yp;

vp2

Simplify[vp2]

(* T to calkowita energia kinetyczna preta *)

T=Integrate[(1/2)*(m/L)*vp2,{l,0,L}]

(* U to calkowita energia potencjalna preta *)

U=(1/2)*m*g*L*Sin[ϕ[t]]

(* funkcja Lagrange'a preta *)

LL=T-U

d1=D[LL,ϕ'[t]];
d2=D[LL,ϕ[t]];d3=D[d1,t];S1=Solve[d3==d2,ϕ''[t]]

(* teraz policzymy energie kinetyczna, szukajac predkosci katowej *)

xA=L*Cos[ϕ[t]];

yA=0;

rA={xA,yA,0}

vA=D[rA,t]

xB=0;

yB=L*Sin[ϕ[t]];

rB={xB,yB,0}

vB=D[rB,t]

rCM=(rA+rB)/2

vCM=D[rCM,t]

(* wyrazenie ponizej zawiera predkosc katowa i powinno byc rowne predkosci liniowej punktu A. Predkosc katowa ma tylko skladowa w kierunku z *)

rCMA=rA-rCM

vCM+{0,0,w3}\[Cross]rCMA

Solve[vA==vCM+{0,0,w3}\[Cross]rCMA,w3]

(* Mamy prosty zwiazek miedzy pochodna ϕ po czasie a predkoscia katowa. Do jego znalezienia uzylismy punktu A. Ten sam wynik powinnismy otrzymac dla punktu B ! *)

rCMB=rB-rCM

Solve[vB==vCM+{0,0,w3}\[Cross]rCMB,w3]

(* inny sposob na policzenie energii kinetycznej preta *)

Tnew=(1/2)*m*vCM.vCM + (1/2)*(1/12)*m*L^2*(ϕ'[t])^2;

Simplify[Tnew==T]

(* dostajemy to samo *)

(* Parametry układu: długość i masa wahadła matematycznego oraz wartość przyspieszenia ziemskiego *)

L=1 ; (* m *)

m=1; (* kg *)

g=981/100; (* m/s^2 *)

(* Warunki początkowe  *)

(Subscript[ϕ, 0])=Pi/3;

(Subscript[ω, 0])=0;

tmax=20; (* s *)

numerycznie=NDSolve[{ϕ''[t]==-((3*g *Cos[ϕ[t]])/(2*L)),ϕ[0]==(Subscript[ϕ, 0]),ϕ'[0]==(Subscript[ω, 0])},ϕ[t],{t,0,tmax},WorkingPrecision->24]

numerycznie

Flatten[numerycznie]

(* Wyciągam funkcję ϕ[t[ z numerycznego rozwiązania *)

ϕ[t_]:=Evaluate[ϕ[t]/.Flatten[numerycznie]]

(* Teraz mogę używać funkcji ϕ[t] tak, jak innych funkcji *)

ϕ[3]

ϕ'[11]

Plot[ϕ[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","ϕ(t) [rad]"}]

Plot[ϕ'[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","ϕ'(t) [rad/s]"}]

Plot[Evaluate[T],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Ekin(t) [J]"}]

Plot[Evaluate[U],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Epot(t) [J]"}]

Plot[Evaluate[T+U],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Etot(t) [J]"}]

(* Modelowanie ruchu pręta przy pomocy wbudowanych obiektów graficznych *)

xaxis=Line[{{-1.1*L,0},{1.1*L,0}}];

yaxis=Line[{{0,-1.1*L},{0,1.1*L}}];

gxaxis:=Graphics[{Thin,xaxis}]

gyaxis:=Graphics[{Thin,yaxis}]

odcAB[t_]=Line[{{L*Cos[ϕ[t]],0},{0,L*Sin[ϕ[t]]}}];

godcAB[t_]:=Graphics[{Thick,odcAB[t]}]

Animate[Show[gxaxis,gyaxis,godcAB[t],PlotRange->{{-1.1*L,1.1*L},{-1.1*L,1.1*L}},AspectRatio->Automatic],{t,0,tmax}]