Na gladkim stole znajduja sie dwie masy, m1 i m2,
polaczone sprezyna o wspolczynniku sprezystosci k
i dlugosci l0. Opisac ruch tego ukladu o czterech stopniach swobody.

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

Epot=(1/2)*k*(Sqrt[(x1[t]-x2[t])^2+(y1[t]-y2[t])^2]-l0)^2;

L=FullSimplify[Ekin-Epot];

Etot=FullSimplify[Ekin+Epot];

d1=FullSimplify[D[L,x1'[t]]]

d2=FullSimplify[D[d1,t]]

d3=FullSimplify[D[L,x1[t]]]

d4=FullSimplify[D[L,y1'[t]]]

d5=FullSimplify[D[d4,t]]

d6=FullSimplify[D[L,y1[t]]]

d7=FullSimplify[D[L,x2'[t]]]

d8=FullSimplify[D[d7,t]]

d9=FullSimplify[D[L,x2[t]]]

d10=FullSimplify[D[L,y2'[t]]]

d11=FullSimplify[D[d10,t]]

d12=FullSimplify[D[L,y2[t]]]

eq1={d2==d3,d5==d6,d8==d9,d11==d12};

eq2=FullSimplify[Solve[eq1,{x1''[t],y1''[t],x2''[t],y2''[t]}]]

eq3=eq2/.Rule->Equal

(* Parametry układu: wartosc mas m1 i m2, dlugosc swobodnej sprezyny, stala sprężystości sprezyny *)

m1=0.5; (* kg *)

m2=0.25; (* kg *)

l0=5/10 ; (* m *)

k= 0.5; (* N/m *)

tmax=50; (* s *)

WP={x1[0]==0,y1[0]==0,x2[0]==2/10,y2[0]==-1/10,x1'[0]==0,y1'[0]==0,x2'[0]==1/10,y2'[0]==-1/10};

eqs=FullSimplify[Flatten[Append[eq3,WP]]]

numerycznie=NDSolve[eqs,{x1[t],y1[t],x2[t],y2[t]},{t,0,tmax}]

Flatten[numerycznie]

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

Plot[x2[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","(Subscript[x, 2])(t) [m]"}]

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

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

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

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

CM[t_]:={m1*x1[t]+m2*x2[t],m1*y1[t]+m2*y2[t]}/(m1+m2);

ParametricPlot[CM[t],{t,0,tmax}]

Modelowanie ruchu układu przy pomocy wbudowanych obiektów graficznych

hlp[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,hlp[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]]

gs[t_,skala_]:=Graphics[Sprezyna6[{x1[t],y1[t]},{x2[t],y2[t]},skala]]

masa1[t_]=Disk[{x1[t],y1[t]},50/1000];

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

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

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

cm[t_]:=Line[{CM[0],CM[t]}]

gcm[t_]:=Graphics[{Pink,cm[t]}]

Animate[Show[gcm[t],gs[t,1/20],gmasa1[t],gmasa2[t],PlotRange->{{-0.2,2},{-2,0.2}},AspectRatio->Automatic],{t,0,tmax}]

Linia pokazuje, ze srodek masy porusza sie ruchem jednostajnym po linii prostej, gdyz nie dzialaja zadne sily zewnetrzne !