Wahadło matematyczne o zmiennej długości ℓ(t) i masie 𝓂 bez przybliżenia małych kątów jako przykład więzów zależnych od czasu

l[t_]=l0*(2+Sin[w*t]);

x[t]=l[t]*Sin[ϕ[t]]

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

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

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

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

Epot=  FullSimplify[m*g*y[t]]

(* Energia kinetyczna punktowej masy *)

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

(* Funkcja Lagrange'a układu *)

L=FullSimplify[Ekin-Epot]

Etot=FullSimplify[Ekin+Epot](* energia całkowita nie będzie zachowana ! *)

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

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

d3=D[d1,t];

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

s1=FullSimplify[s]

(* Parametry układu: długość (l0), częstość zmian długości (w),  masa (m) wahadła matematycznego oraz wartość przyspieszenia ziemskiego (g) *)

l0=1 ; (* m *)

w=3/10; (* 1/s *)

m=1/4; (* kg *)

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

(* Warunki początkowe  *)

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

(Subscript[ω, 0])=0;

tmax=20; (* s *)

numerycznie=NDSolve[{ϕ''[t]==-((g Sin[ϕ[t]]+2 l0 w Cos[t w] ϕ'[t])/(2 l0+l0 Sin[t w])),ϕ[0]==(Subscript[ϕ, 0]),ϕ'[0]==(Subscript[ω, 0])},ϕ[t],{t,0,tmax},WorkingPrecision->24]

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]

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

x[t_]=l[t]*Sin[ϕ[t]];

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

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[x[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","x(t) [m]"}]

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

Plot[Sqrt[x[t]^2+y[t]^2],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","l(t) [m]"}]

(* Modelowanie ruchu wahadła przy pomocy wbudowanych obiektów graficznych *)

sufit=Rectangle[{-1/2,0},{1/2,1/20}];

odc[t_]=Line[{{0,0},{x[t],y[t]}}];

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

gsufit=Graphics[sufit];

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

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

Animate[Show[gsufit,godc[t],gmasa[t],PlotRange->{{-5/2,5/2},{-7/2,1/2}},AspectRatio->Automatic],{t,0,tmax}]

T0=N[2*Pi*Sqrt[l/g]]