Przykład 2: Ruch punktu materialnego P o masie m w jednorodnym polu grawitacyjnym po płaskiej krzywej stopnia trzeciego o równaniu 
y=a*(x-c)^2+b*(x-c)^3.
Jest to po prostu ruch koralika po gładkim drucie o zadanym kształcie, przy czym wszystkie punkty “drutu” należą do płaszczyzny xy.

W tym przypadku posłużymy się formalizmem Lagrange’a, by uzyskać równanie różniczkowe opisujące ruch koralika. Wykorzystanie II zasady dynamiki nie jest praktyczne, bo siły reakcji drutu działające na koralik są bardzo skomplikowane.

Mamy układ o jednym stopniu swobody: aby znaleźć punkt materialny (koralik na drucie), wystarczy podać współrzędną x, bo y policzymy z równania krzywej.

(* Równanie krzywej, po której porusza się punkt materialny P.  *)

y[t]=a*(x[t]-c)^2+b*(x[t]-c)^3;

y'[t]=Simplify[D[y[t],t]]

(* Energia potencjalna zależy wyłącznie od wysokości, na której znajduje się punkt P *)

Epot=  m*g*y[t];

(* Energia kinetyczna to energia kinetyczna punktowej masy *)

Ekin=Simplify[(1/2)*m*(x'[t]^2+y'[t]^2)]

(* Funkcja Lagrange'a L jest różnicą energii kinetycznej i energii potencjalnej *)

L=FullSimplify[Ekin-Epot]

(* Energia całkowita Etot jest sumą energii kinetycznej i energii potencjalnej *)

Etot=FullSimplify[Ekin+Epot];

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

d2=FullSimplify[D[d1,t]]

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

s=FullSimplify[Solve[d2==d3,x''[t]]]

(* Wyciągamy równanie różniczkowe *)

eq1=s[[1]]/. Rule->Equal

(* Tworzymy listę warunków początkowych, na razie w zapisie ogólnym *)

WP={x[0]==x0,x'[0]==vx0};

(* Numeryczne znajdowanie x[t] wymaga kompletu informacji: równania (eq1) oraz warunków początkowych (WP)  *)

eqs=Join[eq1,WP]

(* Parametry liczbowe krzywej trzeciego stopnia oraz wartości masy m i przyspieszenia ziemskiego *)

a=3/4;
b=-1/4;
c=1;
m=1;
g=981/100;

(* Warunki początkowe powinny być tak dobrane, by punkt materialny nie wyskoczył poza maksima krzywej *)

x0=1/100;
vx0=0;

(* czas śledzenia ruchu po krzywej *)
tmax=10;

numerycznie=NDSolve[eqs,x[t],{t,0,tmax}]

numerycznie

Flatten[numerycznie]

(* Wyciągam funkcje x[t] i y[t] z rozwiązania numerycznego *)

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

y[t_]:=Evaluate[(a*(x[t]-c)^2+b*(x[t]-c)^3)]

(* Mogę teraz używać x[t] i y[t] tak, jak innych funkcji. *)

x[1]

y'[4]

(* Mogę też rysować wykresy różnych wielkości *)

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

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

(* Przy okazji mogę znaleźć przybliżoną wartość okresu ruchu koralika po drucie *)

FindRoot[x'[t],{t,3}]

T=t/.FindRoot[x'[t],{t,3}]

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]"}]

(* Zasada zachowania całkowitej energii mechanicznej przy założonym braku tarcia daje możliwość niezależnego sprawdzenia rozwiązania numerycznego *)

Plot[Evaluate[Etot],{t,0,tmax},PlotRange->{9.590,9.5915},AxesLabel->{"t [s]","Etot(t) [J]"}]

(* Teraz możemy bardzo łatwo stworzyć animację ruchu koralika po drucie *)

Animate[Show[Plot[  a*(x-c)^2+b*(x-c)^3, {x,0,3}  ,PlotRange->{{0,3},{0,1}},AxesLabel->{"x [m]","y [m]"}],ListPlot[{{x[t],
y[t]}}, PlotStyle->{PointSize[0.03],Orange}]],{t,0,T}]

aaa=Table[Show[Plot[  a*(x-c)^2+b*(x-c)^3, {x,0,3}  ,PlotRange->{{0,3},{0,1}},AxesLabel->{"x [m]","y [m]"}],ListPlot[{{x[t],
y[t]}}, PlotStyle->{PointSize[0.03],Orange}]],{t,0,T,0.03}];

Export[
 NotebookDirectory[] <> "krzywa_trzeciego_stopnia_2020.gif", aaa]

Uzupełnimy teraz ten notatnik, aby pokazać, jak wygląda siła reakcji więzów działająca na koralik. Aby znaleźć ruch koralika sił nie musieliśmy rozważać, ale II zasada dynamiki oczywiście obowiązuje, więc z niej możemy otrzymać informację o sile reakcji drutu !

(* siła wypadkowa *)

Fwyp[t_]=m*{x''[t],y''[t]};

Fwyp[3]

(* siła grawitacji nie zależy od czasu i jest skierowana pionowo w dół *)

Fgraw= m*{0,-g};

(* siła reakcji więzów *)

Freakcji[t_]= Fwyp[t]-Fgraw;

Freakcji[1]

Plot[Norm[Freakcji[t]],{t,0,T},AxesLabel->{"t [s]"," |(Subscript[(OverVector[F]), reakcji])| [N]"},PlotRange->All]

skala=1/50;

VFreakcji[t_]=Arrow[{{x[t],y[t]},{x[t],y[t]}+skala*Freakcji[t]}];

GVFreakcji[t_]:=Graphics[{Red,Arrowheads[0.02],VFreakcji[t]}]

Animate[Show[Plot[  a*(x-c)^2+b*(x-c)^3, {x,0,3}  ,PlotRange->{{0,3},{-0.04,1.2}},AxesLabel->{"x [m]","y [m]"}],ListPlot[{{x[t],
y[t]}}, PlotStyle->{PointSize[0.02],Orange}],GVFreakcji[t]],{t,0,T}]

Bezpośrednie wykorzystanie II zasady dynamiki, gdy działa skomplikowana siła reakcji jest bardzo trudne !