Układ zawiera dwa cienkie pręty połączone przegubem, umieszczone w płaszczyźnie xy. Lewy koniec pierwszego pręta o długości R umieszczony jest w początku układu, a pręt może obracać się wokół osi z. Prawy koniec pierwszego pręta jest połączony przegubem z lewym końcem drugiego preta o długości L. Prawy koniec drugiego pręta ślizga się bez tarcia po osi x. Układ ma tylko jeden stopien swobody - do jego opisu można na przykład wybrać kąt, jaki tworzy pierwszy pręt z osią x.

(* Współrzędne punktów A i B *)

xA=R*Cos[ϕ[t]];
yA=R*Sin[ϕ[t]];

xB=R*Cos[ϕ[t]]+Sqrt[(R*Cos[ϕ[t]])^2+L^2-R^2];
yB=0;

(* Aby sprawdzić poprawność tych wzorów parametryzuje i rysuje odcinki OA i AB. Chwilowo przyjmuję ϕ[t]=alfa *)

(* Korzystam z ogólnego wzoru:
x=xA + t*(xB-xA);
y=yA + t*(yB-yA);
t należy do przedziału [0,1] *)

x1=l*Cos[alfa] ;
y1=l*Sin[alfa]; (* l w przdziale [0,R] *)

x2=R*Cos[alfa]  + (l/L)*(R*Cos[alfa]+Sqrt[(R*Cos[alfa])^2+L^2-R^2]-R*Cos[alfa]);
y2=R*Sin[alfa] + (l/L)*(-R*Sin[alfa]);

R=1;L=2;alfa=Pi/4;
P1=ParametricPlot[{x1,y1},{l,0,R},PlotRange->All];
P2=ParametricPlot[{x2,y2},{l,0,L},PlotRange->All];
Show[P1,P2]

R=1;L=2;alfa=Pi/2;
P1=ParametricPlot[{x1,y1},{l,0,R},PlotRange->All];
P2=ParametricPlot[{x2,y2},{l,0,L},PlotRange->All];
Show[P1,P2]

R=1;L=2;alfa=3*Pi/4;
P1=ParametricPlot[{x1,y1},{l,0,R},PlotRange->All];
P2=ParametricPlot[{x2,y2},{l,0,L},PlotRange->All];
Show[P1,P2]

R=1;L=2;alfa=Pi;
P1=ParametricPlot[{x1,y1},{l,0,R},PlotRange->All];
P2=ParametricPlot[{x2,y2},{l,0,L},PlotRange->All];
Show[P1,P2]

R=1;L=2;alfa=5*Pi/4;
P1=ParametricPlot[{x1,y1},{l,0,R},PlotRange->All];
P2=ParametricPlot[{x2,y2},{l,0,L},PlotRange->All];
Show[P1,P2]

R=1;L=2;alfa=3*Pi/2;
P1=ParametricPlot[{x1,y1},{l,0,R},PlotRange->All];
P2=ParametricPlot[{x2,y2},{l,0,L},PlotRange->All];
Show[P1,P2]

R=1;L=2;alfa=7*Pi/4;
P1=ParametricPlot[{x1,y1},{l,0,R},PlotRange->All];
P2=ParametricPlot[{x2,y2},{l,0,L},PlotRange->All];
Show[P1,P2]

R=1;L=2;alfa=2*Pi;
P1=ParametricPlot[{x1,y1},{l,0,R},PlotRange->All];
P2=ParametricPlot[{x2,y2},{l,0,L},PlotRange->All];
Show[P1,P2]

R=.;L=.;alfa=.;

x1=x1/.{alfa->ϕ[t]}

y1=y1/.{alfa->ϕ[t]}

x2=x2/.{alfa->ϕ[t]}

y2=y2/.{alfa->ϕ[t]}

(* Teraz licze energie kinetyczna T1 pierwszego preta WPROST Z DEFINICJI tzn. calkujac (sumujac) energie poszczegolnych elementow *)

rho1=m1/R;

dl1podl=Simplify[Sqrt[D[x1,l]^2+D[y1,l]^2],R>0&&L>0];

dm1=rho1*dl1podl;

x1p=D[x1,t]; y1p=D[y1,t];
v1p2=x1p*x1p+y1p*y1p;

Simplify[v1p2]

T1=Simplify[Integrate[(1/2)*dm1*v1p2,{l,0,R}]]

(* Teraz liczę tę samą energię, wykorzystując wzór na moment bezwładności pręta wokół osi prostopadłej do pręta i przechodzącej przez koniec pręta. W przypadku pierwszego pręta łatwo jest stwierdzić, że jego prędkość kątowa wynosi w1=ϕ'[t]. *)

w1=ϕ'[t];
I1=(1/3)*m1*R^2;
T1N=(1/2)*I1*w1^2;
Simplify[T1N==T1]

(* Teraz liczę energię kinetyczną T2 drugiego pręta WPROST Z DEFINICJI tzn. całkując (sumując) energię poszczególnych elementów *)

rho2=m2/L;

dl2podl=Simplify[Sqrt[D[x2,l]^2+D[y2,l]^2],R>0&&L>0];

dm2=rho2*dl2podl;

x2p=D[x2,t]; y2p=D[y2,t];
v2p2=x2p*x2p+y2p*y2p;

Simplify[v2p2,R>0&&L>0&&L>R]

T2=Simplify[Integrate[(1/2)*dm2*v2p2,{l,0,L}],R>0&&L>0&&L>R]

xA=R*Cos[ϕ[t]];
yA=R*Sin[ϕ[t]];

xB=R*Cos[ϕ[t]]+Sqrt[(R*Cos[ϕ[t]])^2+L^2-R^2];
yB=0;

(* środek masy pierwszego pręta *)

xs1=(1/2)*(0+xA);

ys1=(1/2)*(0+yA);

(* środek masy drugiego pręta *)

xs2=(1/2)*(xA+xB);

ys2=(1/2)*(yA+yB);

(* prędkość punktu B *)

xBp=D[xB,t]; yBp=D[yB,t];
vB={xBp,yBp,0};

(* położenie punktu B względem środka masy drugiego pręta *)

rBs2={xB-xs2,yB-ys2,0};

(* prędkość środka masy drugiego pręta *)

vs2={D[xs2,t],D[ys2,t],0};

(* prędkość kątowa drugiego pręta *)

omega2m={0,0,omega2};

solw2=Solve[vB==vs2+omega2m\[Cross]rBs2,omega2]

w2=omega2/.solw2[[1]]

(* Teraz liczę energię kinetyczną T2 drugiego preta, jako sumę energii kinetycznej jego środka masy i energii kinetycznej względem środka masy *)

T2N=(1/2)*m2*vs2.vs2 + (1/2)*(1/12)*m2*L^2*w2^2;

Simplify[T2N==T2,R>0&&L>0&&L>R]

(* Do znalezienia prędkości kątowej drugiego pręta wykorzystaliśmy ruch punktu B. Oczywiście równie dobrze możemy wykorzystać punkt A, bo związek między prędkością liniową punktu i prędkością katową jest taki sam dla wszystkich punktów bryły sztywnej. *)

(* prędkość punktu A *)

xAp=D[xA,t]; yAp=D[yA,t];
vA={xAp,yAp,0};

(* położenie punktu A względem środka masy drugiego pręta *)

rAs2={xA-xs2,yA-ys2,0};

solw2n=Solve[vA==vs2+omega2m\[Cross]rAs2,omega2]

w2n=omega2/.solw2n[[1]]

Simplify[w2n==w2]

(* Teraz liczymy energię potencjalną obu prętów *)

U1=m1*g*ys1;

U2=m2*g*ys2;

EPOT=Simplify[(U1+U2),(L>0)&&(R>0)&&(L>R)]

EKIN=Simplify[(T1+T2),(L>0)&&(R>0)&&(L>R)];

ETOT=Simplify[(T1+T2)+(U1+U2),(L>0)&&(R>0)&&(L>R)];

LL=Simplify[T1+T2-(U1+U2),(L>0)&&(R>0)&&(L>R)];

d1=D[LL,ϕ'[t]];

d2=D[LL,ϕ[t]];

d3=D[d1,t];

S1=Solve[d3==d2,ϕ''[t]]

S1N=Simplify[S1,R>0&&L>0&&L>R]

(* CForm[Simplify[S1]] *)

S2N=S1N/. Rule->Equal

(* Parametry układu: długości i masy prętów oraz wartość przyspieszenia ziemskiego *)

R=1; (* m *)

m1=1; (* kg *)

L=2 ; (* m *)

m2=1; (* kg *)

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

(* Warunki początkowe  *)

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

(Subscript[ω, 0])=0;

tmax=20; (* s *)

numerycznie=NDSolve[{S2N[[1]][[1]],ϕ[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]

ClearAll[xA,yA,xB,yB];

xA[t_]=R*Cos[ϕ[t]];
yA[t_]=R*Sin[ϕ[t]];

xB[t_]=R*Cos[ϕ[t]]+Sqrt[(R*Cos[ϕ[t]])^2+L^2-R^2];
yB[t_]=0;

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[xA[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","xA(t) [m]"}]

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

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

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

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

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

Plot[Evaluate[ETOT],{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+R),0},{1.1*(L+R),0}}];

yaxis=Line[{{0,-1.1*(L+R)},{0,1.1*(L+R)}}];

gxaxis:=Graphics[{Thin,xaxis}]

gyaxis:=Graphics[{Thin,yaxis}]

odcOA[t_]=Line[{{0,0},{xA[t],yA[t]}}];

odcAB[t_]=Line[{{xA[t],yA[t]},{xB[t],yB[t]}}];

godcOA[t_]:=Graphics[{Thick,odcOA[t]}]

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

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