Podwójne wahadło matematyczne płaskie

x1[t]=l1*Sin[phi1[t]];
y1[t]=l1*Cos[phi1[t]];

x2[t]=x1[t]+l2*Sin[phi2[t]];
y2[t]=y1[t]+l2*Cos[phi2[t]];

T1=(1/2)*m1*(D[x1[t],t]^2+D[y1[t],t]^2);
T2=(1/2)*m2*(D[x2[t],t]^2+D[y2[t],t]^2);
T=T1+T2;

U1=-m1*g*y1[t];
U2=-m2*g*y2[t];
U=U1+U2;

L=Simplify[T - U];

d1=D[L,phi1'[t]];
d2=D[L,phi2'[t]];
d3=D[L,phi1[t]];
d4=D[L,phi2[t]];

d5=D[d1,t];
d6=D[d2,t];

s=FullSimplify[Solve[{d5==d3,d6==d4},{phi1''[t],phi2''[t]}]];

eq1=s/. Rule->Equal;

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

l1=1/2 ; (* m *)

l2=3/4;  (* m *)

m1=1/4; (* kg *)

m2=1/5; (* kg *)

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

(* Warunki początkowe  *)

wp={phi1[0]==Pi/3,phi1'[0]==0,phi2[0]==Pi/4,phi2'[0]==0}

tmax=30; (* s *)

numerycznie=NDSolve[{eq1[[1]],wp},{phi1[t],phi2[t]},{t,0,tmax},WorkingPrecision->24];

(* Wyciągam funkcje phi1[t] i phi2[t] z numerycznego rozwiązania *)

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

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

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

phi1[3]

phi1'[11]

phi2[2]

phi2'[8]

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

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

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

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

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

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

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

(* Modelowanie ruchu wahadła podwójnego przy pomocy wbudowanych obiektów graficznych;
Uwaga na znak "-" przy zmiennych y1 i y2 ze względu na układ współrzędnych, którego używa Mathematica.
 *)

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

odc1[t_]=Line[{{0,0},{x1[t],-y1[t]}}];

masa1[t_]=Disk[{x1[t],-y1[t]},1/35];

odc2[t_]=Line[{{x1[t],-y1[t]},{x2[t],-y2[t]}}];

masa2[t_]=Disk[{x2[t],-y2[t]},1/35];

gsufit=Graphics[sufit];

godc1[t_]:=Graphics[{Thick,odc1[t]}]

gmasa1[t_]:=Graphics[masa1[t]]

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

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

Animate[Show[gsufit,godc1[t],gmasa1[t],godc2[t],gmasa2[t],PlotRange->{{-6/5,6/5},{-7/5,1/10}},AspectRatio->Automatic],{t,0,2.1}]

aaa=Table[Show[gsufit,godc1[t],gmasa1[t],godc2[t],gmasa2[t],PlotRange->{{-6/5,6/5},{-7/5,1/10}},AspectRatio->Automatic],{t,0,2.1,0.04}];

tmax

Export[
 NotebookDirectory[] <> "podwojne_wahadlo_plaskie.gif", aaa]