Wahadło matematyczne o długości ℓ z masą (Subscript[m, 2]), którego zaczep o masie (Subscript[m, 1]) może poruszać się bez tarcia po osi x (bez przybliżenia małych kątów). Ten układ ma dwa stopnie swobody. Aby opisać położenie elementów układu wystarczą dwie współrzędne uogólnione, na przykład (Subscript[x, 1]) oraz kąt ϕ.

y1[t]=0

x2[t]=x1[t] + l*Sin[ϕ[t]]

y2[t]=-l*Cos[ϕ[t]]

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

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

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

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

(* Energia kinetyczna dwóch punktowych mas *)

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

(* Energia potencjalna zależy wyłącznie od y2 *)

Epot=  FullSimplify[m2*g*y2[t]]

(* Funkcja Lagrange'a układu *)

L=FullSimplify[Ekin-Epot]

Etot=FullSimplify[Ekin+Epot]

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

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

d3=D[d1,t];

s1=Solve[d3==d2,ϕ''[t]];

FullSimplify[s1]

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

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

d6=D[d4,t];

s2=Solve[d6==d5,x1''[t]];

FullSimplify[s2]

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

l=1 ; (* m *)

m2=0.25; (* kg *)

m1=0.25; (* kg *)

g=9.81; (* m/s^2 *)

(* Warunki początkowe  *)

x10=0;

v10=0;

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

(Subscript[ω, 0])=0;

tmax=20; (* s *)

numerycznie=NDSolve[{ϕ''[t]==-((g Sin[ϕ[t]]+Cos[ϕ[t]] x1''[t])/l),x1''[t]==((l m2 (Sin[ϕ[t]] ((ϕ'[t])^2)-Cos[ϕ[t]] ϕ''[t]))/(m1+m2)),ϕ[0]==(Subscript[ϕ, 0]),ϕ'[0]==(Subscript[ω, 0]),x1[0]==x10,x1'[0]==v10},{x1[t],ϕ[t]},{t,0,tmax}]

numerycznie

Flatten[numerycznie]

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

Flatten[numerycznie][[1]]

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

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

(* Teraz mogę używać funkcji x1(t) i ϕ[t] tak, jak innych funkcji. Z ich pomocą mogę zdefiniować także inne potrzebne mi wielkości. *)

x1[2]

x1'[8.5]

ϕ[3]

ϕ'[11]

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

x2[1]

x2'[4]

y2[t_]=-l*Cos[ϕ[t]];

y2[4]

(* Całkowity pęd układu w kierunku osi x powinien być stałą wielkością, bo całkowita siła zewnętrzna ma zerową składową w tym kierunku  *)

px=m1*x1'[t]+m2*x2'[t];

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

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

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

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

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

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

Plot[{m1*x1'[t],m2*x2'[t],Evaluate[px]},{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","(Subscript[m, 1])(Subscript[x, 1]^('))(t), (Subscript[m, 2])(Subscript[x, 2]^('))(t), (Subscript[p, x])(t) [m/s]"}]

Modelowanie ruchu wahadła ze swobodnym zaczepem przy pomocy wbudowanych obiektów graficznych

pow=Line[{{-3/5,0},{6/5,0}}];

uchwyt[t_]:=Rectangle[{x1[t]-1/20,0},{x1[t]+1/20,1/10}];

odc[t_]:=Line[{{x1[t],0},{x2[t],y2[t]}}];

masa2[t_]:=Disk[{x2[t],y2[t]},1/20];

gpow=Graphics[{Thick,pow}];

guchwyt[t_]:=Graphics[uchwyt[t]];

godc[t_]:=Graphics[{Thick,odc[t]}]

gmasa2[t_]:=Graphics[masa2[t]]

Animate[Show[gpow,guchwyt[t],godc[t],gmasa2[t],PlotRange->{{-3/5,6/5},{-6/5,1/5}},AspectRatio->Automatic],{t,0,tmax}]