Wahadło matematyczne o długości ℓ i masie 𝓂, którego punkt zawieszenia porusza się ze stałą prędkością kątową ω po okręgu leżącym w pionowej płaszczyźnie 

Dla tego układu energia całkowita nie jest stała !

x[t]=r*Cos[w*t]+l*Sin[ϕ[t]];

y[t]=r*Sin[w*t]-l*Cos[ϕ[t]];

x'[t]=D[x[t],t];

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

(* Energia kinetyczna punktowej masy *)

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

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

Epot=  FullSimplify[m*g*y[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];

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

FullSimplify[s]

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

r=1/3; (* m *)

w=2*Pi;(* rad/s *)

l=1 ; (* m *)

m=1/4; (* kg *)

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

(* Warunki początkowe  *)

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

(Subscript[w, 0])=0;

tmax=80; (* s *)

numerycznie=NDSolve[{ϕ''[t]==((r (w^2) Cos[t w-ϕ[t]]-g Sin[ϕ[t]])/l),ϕ[0]==(Subscript[ϕ, 0]),ϕ'[0]==(Subscript[w, 0])},ϕ[t],{t,0,tmax},WorkingPrecision->24,MaxSteps->50000]

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]

x[t_]:=r*Cos[w*t]+l*Sin[ϕ[t]];

y[t_]:=r*Sin[w*t]-l*Cos[ϕ[t]];

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 ruchu wahadła przy pomocy wbudowanych obiektów graficznych *)

kolo=Circle[{0,0},r];

s=Text[StyleForm["S",FontSize->12,FontWeight->"Bold"],{-1/20,-1/20}];

promien[t_]=Line[{{0,0},{r*Cos[w*t],r*Sin[w*t]}}];

punkt[t_]=Disk[ {r*Cos[w*t],r*Sin[w*t]},1/50];

odc[t_]=Line[{{r*Cos[w*t],r*Sin[w*t]},{x[t],y[t]}}];

masa[t_]=Disk[{x[t],y[t]},1/20];

gs=Graphics[s];

gkolo=Graphics[{Thick,kolo}];

gpromien[t_]:=Graphics[promien[t]];

gpunkt[t_]=Graphics[punkt[t]];

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

gmasa[t_]:=Graphics[masa[t]];

Animate[Show[gkolo,gs,gpromien[t],gpunkt[t],godc[t],gmasa[t],PlotRange->{{-r-l-1/20*l,r+l+1/20*l},{-r-l-1/20*l,r+l+1/20*l}},AspectRatio->Automatic],{t,0,tmax}]