Na powierzchni równi, która ma masę M i może poruszać się
bez tarcia po poziomym stole, położono na wysokości h klocek o masie m. Opisać ruch równi i klocka 
(a) przy pomocy zasad dynamiki Newtona,
(b) korzystając z formalizmu Lagrange’a

Najpierw wypisujemy równania wynikające z II zasady dynamiki Newtona

Na równię i klocek działają sily przyciągania ziemskiego, siły nacisku i reakcji na nacisk. Możemy także załozyć, ze równia porusza się tylko w kierunku x. Uwaga: w tym przypadku 
R=|R|=|N| ≠m*g*Cos[α] ! Wiemy tylko, ze nacisk jest prostopadły do powierzchni równi, a siły N i R to akcja i reakcja (III zasada dynamiki Newtona). Oczywiście na równię też działa siła ciężkości i reakcja na nacisk ze strony stołu, ale tych sił nie musimy rozpatrywać.

eq1=m*x1''[t]==R*Sin[α];

eq2=m*y1''[t]==-m*g+R*Cos[α];

eq3=M*x2''[t]==-R*Sin[α];

Niewiadomymi w tych równaniach są x1’’[t], y1’’[t], x2’’[t] oraz wartość siły R. Potrzebujemy jeszcze jednego równania, które ma charakter czysto geometryczny i pokazuje, że klocek znajduje się na równi:

eq4=(h-y1[t])==(x1[t]-x2[t])*Tan[α];

Liczę drugą pochodną po czasie obu stron tego równania

eq4n=D[eq4,{t,2}]

Teraz rozwiązujemy układ równań na x1’’[t], y1’’[t], x2’’[t] i R

sol=Simplify[Solve[{eq1,eq2,eq3,eq4n},{x1''[t],y1''[t],x2''[t] , R}]]

Teraz zapiszemy rozwiązania w postaci układu równań różniczkowych (bardzo prostego, bo wszystkie przyspieszenia są stałe) i rozwiążemy go dla zadanych warunków początkowych.

diffeqs=sol/. Rule->Equal

WP={x1[0]==0,y1[0]==h,x2[0]==0,x1'[0]==0,y1'[0]==0,x2'[0]==0};

soldiffeqs=Simplify[DSolve[{x1''[t]==((g M Sin[α])/(M Cos[α]+(m+M) Sin[α] Tan[α])),y1''[t]==-((g (m+M) Sin[α] Tan[α])/(M Cos[α]+(m+M) Sin[α] Tan[α])),x2''[t]==-((g m Sin[α])/(M Cos[α]+(m+M) Sin[α] Tan[α])),x1[0]==0,y1[0]==h,x2[0]==0,x1'[0]==0,y1'[0]==0,x2'[0]==0},{x1[t],y1[t],x2[t]},t]]

Z rozwiązania w postaci listy dostaję szukane funkcje opisujące ruch klocka i równi

xs1[t_]=x1[t]/.soldiffeqs[[1]]

ys1[t_]=y1[t]/.soldiffeqs[[1]]

xs2[t_]=x2[t]/.soldiffeqs[[1]]

Warto sprawdzić, czy spełnione jest prawo zachowania energii mechanicznej.

Etot=FullSimplify[(1/2)*m*(xs1'[t]^2+ys1'[t]^2)+(1/2)*M*xs2'[t]^2+m*g*ys1[t],m>0&&M>0&&g>0&&0<α<Pi/2]

Siły ciężkości działające na równię oraz klocek mają kierunek pionowy. Pęd w kierunku poziomym też powinien być stały w czasie i wynosić zero.

Px=FullSimplify[m*xs1'[t]+M*xs2'[t],m>0&&M>0&&g>0&&0<α<Pi/2]

Py=FullSimplify[m*ys1'[t],m>0&&M>0&&g>0&&0<α<Pi/2]

(* Skladowa pozioma calkowitej sily dzialajacej na uklad *)

Fx=FullSimplify[m*xs1''[t]+M*xs2''[t],m>0&&M>0&&g>0&&0<α<Pi/2]

(* Składowa pionowa calkowitej sily dzialajacej na uklad:
Fy = -(m+M)*g + R, gdzie R to sila reakcji podloza -> R = Fy + (m+M)*g *)

Fy=FullSimplify[m*ys1''[t],m>0&&M>0&&g>0&&0<α<Pi/2]

R=FullSimplify[Fy + (m+M)*g]

ClearAll["Global`*"];

Teraz ten sam problem rozwiążemy, stosując formalizm Lagrange’a

Musimy tylko zapisać energię kinetyczną i potencjalną układu ! W ogóle nie rozpatrujemy działających sił.
Dlatego na rysunku wystarczy zaznaczyć położenie elementów układu.
Musimy zacząć od stwierdzenia, że nasz układ ma dwa stopnie swobody, zakładając, że klocek znajduje się na równi. Do opisu układu wystarczają x1[t] i x2[t], bo y1[t] dostajemy z warunku  eq4=(h-y1[t])=(x1[t]-x2[t])*Tan[α];

eq4=(h-y1[t])==(x1[t]-x2[t])*Tan[α]; (* warunek z podobienstwa trojkatow *)

s1=Simplify[Solve[eq4,y1[t]]]

s1[[1]]

y1[t_]=y1[t]/.s1[[1]]

y1[t]

y1'[t]

Liczę energię kinetyczną, energię potencjalną, funkcję Lagrange’a i całkowitą energię mechaniczną układu

Ekin=FullSimplify[(1/2)*m*(x1'[t]^2+y1'[t]^2)+(1/2)*M*x2'[t]^2,m>0&&M>0&&g>0&&0<α<Pi/2]

Epot=FullSimplify[m*g*y1[t],m>0&&M>0&&g>0&&0<α<Pi/2]

L=FullSimplify[Ekin-Epot ,m>0&&M>0&&g>0&&0<α<Pi/2]

Etot=FullSimplify[Ekin+Epot ,m>0&&M>0&&g>0&&0<α<Pi/2]

Równania Lagrange’a

d1=FullSimplify[D[L,x1'[t]]]

d2=FullSimplify[D[d1,t]]

d3=FullSimplify[D[L,x1[t]]]

d4=FullSimplify[D[L,x2'[t]]]

d5=FullSimplify[D[d4,t]]

d6=FullSimplify[D[L,x2[t]]]

Rozwiązujemy układ równań na x1’’[t] i  y1’’[t]

sol=Simplify[Solve[{d2==d3,d5==d6},{x1''[t],x2''[t]}]]

Teraz zapiszemy rozwiązania w postaci układu równań różniczkowych (bardzo prostego, bo wszystkie przyspieszenia są stałe) i rozwiążemy go dla zadanych warunków początkowych

diffeqs=sol/. Rule->Equal

WP={x1[0]==0,x2[0]==0,x1'[0]==0,x2'[0]==0};

soldiffeqs=Simplify[DSolve[Join[Flatten[diffeqs],WP],{x1[t],x2[t]},t]]

Z rozwiązania w postaci listy dostaję szukane funkcje opisujące ruch klocka i równi

xs1[t_]=x1[t]/.soldiffeqs[[1]]

xs2[t_]=x2[t]/.soldiffeqs[[1]]

ys1[t_]=Simplify[y1[t]/.{x1[t]->xs1[t],x2[t]->xs2[t]}]

Warto sprawdzić, czy spełnione jest prawo zachowania energii mechanicznej.

Etot=FullSimplify[(1/2)*m*(xs1'[t]^2+ys1'[t]^2)+(1/2)*M*xs2'[t]^2+m*g*ys1[t],m>0&&M>0&&g>0&&0<α<Pi/2]

Siły ciężkości działające na równię oraz klocek mają kierunek pionowy. Pęd w kierunku poziomym też powinien być stały w czasie i wynosić zero.

Px=FullSimplify[m*xs1'[t]+M*xs2'[t],m>0&&M>0&&g>0&&0<α<Pi/2]

Na koniec zobaczmy jeszcze, pod jakim kątem względem powierzchni stołu (β) porusza się klocek. Nie jest to kąt α, ale kąt wynikający ze stosunku składowych wektora prędkości klocka  !

tangensbeta=FullSimplify[-ys1'[t]/xs1'[t],m>0&&M>0&&g>0&&0<α<Pi/2]

Dla bardzo ciężkiej równi (M ⟶∞) powinniśmy dostać wzory słuszne dla “zwykłego” zsuwania się klocka po równi bez udziału tarcia

Limit[xs1[t],M->∞]

Limit[ys1[t],M->∞]

Limit[xs2[t],M->∞]

Limit[tangensbeta,M->∞]

stmax=Solve[ys1[t]==0,t]

stmax[[2]]

tmax=t/.stmax[[2]]

xs2min=Simplify[xs2[tmax]]

Teraz pobawimy się elementami przyszłej animacji. Zaczynamy od samej równi ...

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

M=1; (* kg *)

m= 1; (* kg *)

h=2;

a=6;

α=ArcTan[h/a];

β=ArcTan[tangensbeta];

N[β]

rownia=Line[{{0,h},{0,0},{a,0},{0,h}}];

gr=Graphics[{Thick,rownia}]

x1=a/3;

y1=h-x1*Tan[α];

P1={x1,y1};

dx=dy=1/5; (* rozmiary klocka *)

P2=P1+{dx*Cos[α],-dy*Sin[α]};

P3=P2+{dy*Sin[α],dy*Cos[α]};

P4=P1+{dy*Sin[α],dy*Cos[α]};

klocek=Line[{P1,P2,P3,P4,P1}];

gk=Graphics[{Thick,klocek}];

Show[gr,gk]

(* lekkie przesunięcie klocka w górę o eps *)

eps=1/50;

P1[[2]]=P1[[2]]+eps;

P2[[2]]=P2[[2]]+eps;

P3[[2]]=P3[[2]]+eps;

P4[[2]]=P4[[2]]+eps;

klocek=Line[{P1,P2,P3,P4,P1}];

gk=Graphics[{Thick,klocek}];

Show[gr,gk]

(* teraz chcemy ożywić ten obrazek -> współrzędne równi i klocka stają się funkcjami czasu *)

p1[t_]:={xs1[t],eps+ys1[t]}

p2[t_]:=p1[t]+{dx*Cos[α],eps-dy*Sin[α]}

p3[t_]:=p2[t]+{dy*Sin[α],dy*Cos[α]}

p4[t_]:=p1[t]+{dy*Sin[α],dy*Cos[α]}

kloc[t_]:=Line[{p1[t],p2[t],p3[t],p4[t],p1[t]}]

gkloc[t_]:=Graphics[{Thick,kloc[t]}]

rown[t_]:=Line[{{xs2[t],h},{xs2[t],0},{xs2[t]+a,0},{xs2[t],h}}];

grown[t_]:=Graphics[{Thick,rown[t]}];

stol=Line[{{xs2min-a/10,-eps},{a+a/10,-eps}}];

gstol=Graphics[{Thick,stol}];

N[tmax]

tor[t_]:=Line[{{xs1[0],ys1[0]},{xs1[t],ys1[t]}}];

gtor[t_]:=Graphics[{Red,Thin,tor[t]}];

Animate[Show[grown[t],gkloc[t],gstol,gtor[t],PlotRange->{{xs2min-a/10,a+a/10},{-h/10,h+h/10}}],{t,0,tmax}]