Dwa wahadła o tej samej długości l z masami m1 i m2 zawieszono na tej samej wysokości w odległości a. Masy m1 i m2 połączono nieważką sprężyną o stałej sprężystości k. Zakładamy, że swobodna sprężyna ma długość a.

(* Położenie masy m1 *)

x1[t]=l*Sin[ϕ1[t]];

y1[t]=-l*Cos[ϕ1[t]];

(* Położenie masy m2 *)

x2[t]=a+l*Sin[ϕ2[t]];

y2[t]=-l*Cos[ϕ2[t]];

(* Prędkości obu mas *)

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

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

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

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

Ekin=FullSimplify[(1/2)*m1*(x1'[t]^2+y1'[t]^2)+(1/2)*m2*(x2'[t]^2+y2'[t]^2)]

Epot=  FullSimplify[m1*g*y1[t]+m2*g*y2[t]+(1/2)*k*(Sqrt[(x1[t]-x2[t])^2+(y1[t]-y2[t])^2]-a)^2]

L=Ekin-Epot

Etot=Ekin+Epot

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

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

d3=D[d2,t];

d4=D[L,ϕ2[t]];

d5=D[L,ϕ2'[t]];

d6=D[d5,t];

s=Simplify[Solve[{d1==d3,d4==d6},{ϕ1''[t],ϕ2''[t]}]]

eq1=(s/. Rule->Equal)[[1]]

(* Parametry układu: długość i masy wahadeł matematycznych, wartość przyspieszenia ziemskiego, odległość między punktami zawieszenia wahadeł i stała sprężystości sprężyny *)

l=1/2 ; (* m *)

m1=1/4; (* kg *)

m2=1/4; (* kg *)

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

a=1/2; (* m *)

k=1/5; (* N/m *)

(* Warunki początkowe  *)

(Subscript[ϕ, 10])=-Pi/6;

(Subscript[ω, 10])=0;

(Subscript[ϕ, 20])=0;

(Subscript[ω, 20])=0;

WP={ϕ1[0]==(Subscript[ϕ, 10]),ϕ1'[0]==(Subscript[ω, 10]),ϕ2[0]==(Subscript[ϕ, 20]),ϕ2'[0]==(Subscript[ω, 20])};

eq2=Flatten[Append[eq1,WP]]

tmax=50; (* s *)

numerycznie=NDSolve[eq2,{ϕ1[t],ϕ2[t]},{t,0,tmax},WorkingPrecision->24]

Flatten[numerycznie]

(* Wyciągam funkcje z numerycznego rozwiązania *)

ϕ1[t_]:=Evaluate[ϕ1[t]/.Flatten[numerycznie]]

ϕ2[t_]:=Evaluate[ϕ2[t]/.Flatten[numerycznie]]

Plot[ϕ1[t],{t,0,tmax}]

Plot[ϕ2[t],{t,0,tmax}]

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

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

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

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

Plot[x1[t],{t,0,tmax}]

Plot[y1[t],{t,0,tmax}]

Plot[Sqrt[x1[t]^2+y1[t]^2],{t,0,tmax}]

Plot[x2[t],{t,0,tmax}]

Plot[y2[t],{t,0,tmax}]

Plot[Sqrt[x2[t]^2+y2[t]^2],{t,0,tmax}]

Plot[Evaluate[Ekin],{t,0,tmax}]

Plot[Evaluate[Epot],{t,0,tmax}]

Plot[Evaluate[Ekin+Epot],{t,0,tmax}]

(* Tu używam modułu dla sprężyny *)

gsprezyna[Pp0_,Pk0_,skala0_]:=Module[{Pp=Pp0,Pk=Pk0,skala=skala0,h,t,a,tab1,d,theta,i,tab5,tab6,Sprezyna6},h[t_,a_,skala_]:=If[t≥a/18&&t≤17*a/18,skala*Sin[18*Pi/a*(t-a/18)],0];
tab1[a_,skala_]:=Simplify[Table[{t,h[t,a,skala]},{t,0,a,a/36}],a>0];
d[Pp_,Pk_]:=Sqrt[(Pk-Pp).(Pk-Pp)];
theta[Pp_,Pk_]:=ArcTan[Pk[[1]]-Pp[[1]],Pk[[2]]-Pp[[2]]];
tab5[Pp_,Pk_,skala_]:=Table[RotationMatrix[theta[Pp,Pk]].tab1[d[Pp,Pk],skala][[i]],{i,1,37}];
tab6[Pp_,Pk_,skala_]:=Transpose[{Transpose[tab5[Pp,Pk,skala]][[1]]+Pp[[1]],Transpose[tab5[Pp,Pk,skala]][[2]]+Pp[[2]]}];
Sprezyna6[Pp_,Pk_,skala_]:=Line[tab6[Pp,Pk,skala]];
Graphics[Sprezyna6[Pp,Pk,skala]]]

(* Modelowanie *)

sufit=Rectangle[{-3,0},{3,1/30}];

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

odc2[t_]=Line[{{a,0},{x2[t],y2[t]}}];

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

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

gsufit=Graphics[sufit];

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

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

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

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

gmojasprezyna[t_,skala_]:=gsprezyna[{x1[t],y1[t]},{x2[t],y2[t]},skala];

Animate[Show[gsufit,godc1[t],godc2[t],gmasa1[t],gmasa2[t],gmojasprezyna[t,1/25],PlotRange->{{-1/2,1},{-6/10,1/10}}],{t,0,tmax}]

tmax

aaa = Table[Show[gsufit,godc1[t],godc2[t],gmasa1[t],gmasa2[t],gmojasprezyna[t,1/25],PlotRange->{{-1/2,1},{-6/10,1/10}}],{t,0,2 , 0.05}];

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