Do masy (Subscript[m, 1]) poruszającej się bez tarcia po równi pochyłej, przymocowano wahadło matematyczne o długości l i masie (Subscript[m, 2]).

(* Zanim zajmiemy się rozwiązaniem równania różniczkowego, pobawimy się elementami przyszłej animacji. Zaczynamy od samej równi. *)

h=2;

a=6;

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

gr=Graphics[{Thick,rownia}]

(* 
Patrz notatki:
*)

Clear[a,h]

prostaAH= y==-(h/a)*x+h;

Solve[y1==(a/h)*x1 + b1,b1];

y1=-(h/a)*x1+h;

prostal1=Simplify[y==(a/h)*x +  ((-a x1+h y1)/h)];

(* prosta l3 *)

x2=x1+dx;

y2=-(h/a)*x2+h;

prostal2=FullSimplify[y==(a/h)*(x-x2)+y2];

prostal3=y==-(h/a)*x+h+dy;

h=2;a=6;dx=1/5; dy=1/5;x1=2;eps=1/20;

plot=ContourPlot[Evaluate[{prostaAH,prostal1,prostal2,prostal3}],{x,0,6},{y,0,5/2},AspectRatio->Automatic]

Clear[h,a,dx, dy,x1,eps];

(* Punkty P1, P2, P3, P4 stanowią wierzchołki prostokąta, który przedstawia masę m. Są one punktami wspólnymi odpowiednich par prostych  *)

(* P1 *)

sp1=Solve[{prostaAH,prostal1},{x,y}];

P1={x,y}/.sp1[[1]];

(* Wprowadzam mały parametr eps, by dla potrzeb rysunku oddzielić klocek od równi *)

P1=P1+{0,eps};

sp2=Solve[{prostaAH,prostal2},{x,y}];

P2={x,y}/.sp2[[1]];

P2=P2+{0,eps};

sp3=FullSimplify[Solve[{prostal3,prostal2},{x,y}]];

P3={x,y}/.sp3[[1]];

P3=P3+{0,eps};

sp4=FullSimplify[Solve[{prostal3,prostal1},{x,y}]];

P4={x,y}/.sp4[[1]];

P4=P4+{0,eps};

h=2;a=6;dx=1/3; dy=1/3;x1=2;eps=1/100;

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

gk=Graphics[{Thick,klocek}]

Show[gr,gk]

srodekklocka=(1/2)*(P1+P2);

l=3;

φ=Pi/4;

xk=srodekklocka[[1]]+l*Sin[φ];

yk=srodekklocka[[2]]-l*Cos[φ];

odc=Line[{srodekklocka,{xk,yk}}];

masa2=Disk[{xk,yk},1/10];

godc:=Graphics[{Thick,odc}];

gmasa2:=Graphics[masa2];

Show[gr,gk,godc,gmasa2]

ClearAll["Global`*"];

(* Teraz zajmiemy się opisem ruchu układu, korzystając z formalizmu Lagrange'a. Nasz układ ma dwa stopnie swobody. Do opisu położenia obu mas wystarczy podanie współrzędnej x1 klocka oraz kąta wahadła φ. (rysunek) *)

(* (h-y1)/x1 = h/a *)

y1[t]=h-(h/a)*x1[t];

y1'[t]=D[y1[t],t]

x1'[t]=D[x1[t],t]

x2[t]=x1[t]+l*Sin[φ[t]];

x2'[t]=D[x2[t],t]

y2[t]=y1[t]-l*Cos[φ[t]];

y2'[t]=D[y2[t],t]

(* Energia potencjalna zależy wyłącznie od składowych y wektorów położenia mas m1 i m2. *)

Epot=  FullSimplify[m1*g*y1[t] + m2*g*y2[t]]

(* Energia kinetyczna jest sumą dwóch energii kinetycznych *)

Ekin=FullSimplify[(1/2)*m1*(x1'[t]^2+y1'[t]^2)+(1/2)*m2*(x2'[t]^2+y2'[t]^2)]

(* Funkcja Lagrange'a układu *)

L=FullSimplify[Ekin-Epot]

Etot=FullSimplify[Ekin+Epot]

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

d2=D[L,x1[t]];

d3=D[d1,t];

d4=D[L,φ'[t]];

d5=D[L,φ[t]];

d6=D[d4,t];

s=FullSimplify[Solve[{d3==d2,d6==d5},{x1''[t],φ''[t]}]]

eq=s/. Rule->Equal

(* Parametry układu: wymiary równi pochyłej (a i h), masa klocka (m1),długość wahadła (l), masa wahadła (m2) oraz wartość przyspieszenia ziemskiego (g) *)

a=12 ; (* m *)

h=2; (* m *)

m1=1/4; (* kg *)

l=3; (* m *)

m2=1/5; (* kg *)

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

(* Warunki początkowe  *)

x10=1/10;

v10=0;

φ0=Pi/6;

ω0=0;

tmax=4; (* s *)

WP={x1[0]==x10,x1'[0]==v10,φ[0]==φ0,φ'[0]==ω0};

eq2=Flatten[Append[eq,WP]]

numerycznie=NDSolve[eq2,{x1[t],φ[t]},{t,0,tmax},WorkingPrecision->20]

Flatten[numerycznie]

(* Wyciągam funkcje x1[t[ i φ[t] z numerycznego rozwiązania *)

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

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

x1[2]

x1'[1]

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

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

φ[1]

φ'[1]

y1[t_]=h-(h/a)*x1[t];

x2[t_]=x1[t]+l*Sin[φ[t]];

y2[t_]=y1[t]-l*Cos[φ[t]];

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

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

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

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

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[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 *)

eps=1/20;

h=2;

a=15;

dx=3/10;

dy=3/10;

p1[t_]:={x1[t],eps+((h (a-x1[t]))/a)}

p2[t_]:={dx+x1[t],eps+((h (a-dx-x1[t]))/a)}

p3[t_]:={dx+((a dy h)/((a^2)+(h^2)))+x1[t],eps+(1/(a ((a^2)+(h^2))))(a (h^3)+(a^3) (dy+h)-(a^2) h (dx+x1[t])-(h^3) (dx+x1[t]))}

p4[t_]:={((a dy h)/((a^2)+(h^2)))+x1[t],eps+h+(((a^2) dy)/((a^2)+(h^2)))-((h x1[t])/a)}

(* początek wahadła w srodku odcinka P1P2 *)

pk[t_]:=(1/2)*(p1[t]+p2[t]);

(* koniec wahadła *)

pp[t_]:=pk[t]+{l*Sin[φ[t]],-l*Cos[φ[t]]}

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

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

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

gr=Graphics[{Thick,rownia}];

odc2[t_]:=Line[{pk[t],pp[t]}];

masa22[t_]:=Disk[pp[t],2/10];

godc2[t_]:=Graphics[{Thick,odc2[t]}]

gmasa22[t_]:=Graphics[masa22[t]]

Animate[Show[gr,gkloc[t],godc2[t],gmasa22[t],PlotRange->{{0,15.5},{-3,2.5}}],{t,0,tmax}]