Na gładkim stole znajdują się dwie masy, m1 i m2,
połączone sprężyną o współczynniku sprężystości k
i długości l0. Opisać ruch tego układu, rozwiązując numerycznie równanie Newtona dla ruchu względnego.

r={x[t],y[t]}

V=(1/2)*k*(Sqrt[r.r]-l0)^2

Fx=-D[V,x[t]]

Fy=-D[V,y[t]]

eq1=μ*D[r,{t,2}]=={Fx,Fy};

(* Parametry ruchu jednostajnego srodka masy *)

R0={X0,Y0};

V0={VX0,VY0};

R[t_]=V0*t+R0;

X0=2/30;

Y0=-1/30;

VX0=1/30;

VY0=-1/30;

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

m1=1/2; (* kg *)

m2=1/4; (* kg *)

μ=(m1*m2)/(m1+m2);

M=m1+m2;

l0=5/10 ; (* m *)

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

tmax=50; (* s *)

x0=1/5;

y0=-1/10;

vx0=1/10;

vy0=-1/10;

WP={x[0]==x0,x'[0]==vx0,y[0]==y0,y'[0]==vy0}

eqs={μ* x''[t]==Fx,μ* y''[t]==Fy,x[0]==x0,x'[0]==vx0,y[0]==y0,y'[0]==vy0};

numerycznie=NDSolve[eqs,{x[t],y[t]},{t,0,tmax}]

Flatten[numerycznie]

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

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

r1[t_]=R[t]-(m2/M)*{x[t],y[t]};

r2[t_]=R[t]+(m1/M)*{x[t],y[t]};

Ekin=(1/2)*m1*D[r1[t],t].D[r1[t],t]+(1/2)*m2*D[r2[t],t].D[r2[t],t];

Epot=(1/2)*k*(Sqrt[{x[t],y[t]}.{x[t],y[t]}]-l0)^2;

Etot=Ekin+Epot;

x1[t_]:=r1[t][[1]]

y1[t_]:=r1[t][[2]]

x2[t_]:=r2[t][[1]]

y2[t_]:=r2[t][[2]]

CM[t_]:=R[t]

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]"}]

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 środek masy porusza się ruchem jednostajnym po linii prostej, gdyż nie działają żadne siły zewnętrzne !