Wahadło matematyczne o długości ℓ i masie 𝓂 bez przybliżenia małych kątów

(* 


 *)

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

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

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

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

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

(* Energia kinetyczna punktowej masy *)

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

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

l=1 ; (* m *)

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]])/l),ϕ[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]"}]

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

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/20];

gsufit=Graphics[sufit];

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

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

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

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