Potrójne wahadło matematyczne

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

x3[t]=x2[t]+l3*Sin[phi3[t]];

y3[t]=y2[t]+l3*Cos[phi3[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);

T3=(1/2)*m3*(D[x3[t],t]^2+D[y3[t],t]^2);

T=T1+T2+T3;

U1=-m1*g*y1[t];

U2=-m2*g*y2[t];

U3=-m3*g*y3[t];

U=U1+U2+U3;

Calkowita=U+T;

L=Simplify[T-U];

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

d2=D[L,phi2'[t]];

d3=D[L,phi3'[t]];

d4=D[L, phi1[t]];

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

d6=D[L,phi3[t]];

d7=D[d1,t];

d8=D[d2,t];

d9=D[d3,t];

s=Simplify[Solve[{d7==d4, d8==d5,d9==d6},{phi1''[t],phi2''[t],phi3''[t]}]];

eq2=s/. Rule->Equal;

eqs=Flatten[eq2];

(* Parametry układu podane DOKŁADNIE *)

l1=1/2;

l2=1/3;

l3=1/4;

m1=1/2;

m2=1/4;

m3=1;

g=981/100;

phi10=Pi/3;

phi20=Pi/4;

phi30=Pi/5;

tmax=20

om0=0;

WP={phi1[0]==phi10,phi1'[0]==om0,phi2[0]==phi20,phi2'[0]==om0,phi3[0]==phi30,phi3'[0]==om0};

wynik=NDSolve[{eqs,WP},{phi1[t],phi2[t],phi3[t]},{t,0,tmax},WorkingPrecision->24,MaxSteps->50000];

Flatten[wynik]

phi1[t_]:=Evaluate[phi1[t]/.Flatten[wynik][[1]]]

phi1[2]

phi2[t_]:=Evaluate[phi2[t]/.Flatten[wynik][[2]]]

phi2[2]

phi3[t_]:=Evaluate[phi3[t]/.Flatten[wynik][[3]]]

Plot[phi1[t],{t,0,tmax}, PlotRange->All]

Plot[phi2[t],{t,0,tmax}, PlotRange->All]

Plot[phi3[t],{t,0,tmax}, PlotRange->All]

Plot[Evaluate[Calkowita],{t,0,tmax}, PlotRange->All]

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

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

(* Modelowanie zachowania układu przy pomocy wbudowanych obiektów graficznych *)

masa[t_]=Disk[{x1[t],-y1[t]},1/20];

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

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

odc3[t_]=Line[{{x2[t],-y2[t]},{x3[t],-y3[t]}}];

masa3[t_]=Disk[{x3[t],-y3[t]},1/20];

gsufit=Graphics[sufit];

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

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

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

godc3[t_]:=Graphics[{Thick,odc3[t]}]

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

gmasa3[t_]:=Graphics[masa3[t]]

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