W wierzchołkach A,B,C i D kwadrata o boku 2a zaczepiono sprężyny o współczynnikach sprężystości kA, kB, kC i kD. Drugie końce sprężyn przymocowano do punktowej masy m, która może się poruszać w płaszczyźnie kwadrata. W położeniu równowagi  (środek kwadratu) żadna ze sprężyn nie jest ani rozciągnięta, ani ściśnięta, co oznacza,że długość swobodnej sprężyny wynosi l0=a*Sqrt[2]. Zakładając, że w chwili t=0 położenie masy i jej prędkość wynoszą odpowiednio r(t=0)=(x0,y0) oraz v(t=0)=(vx0,vy0), znajdujemy położenie i prędkość masy w dowolnej chwili t>0, r(t)= (x(t),y(t)) i v(t)=(x'(t),y'(t)).

(* współrzędne wierzchołków kwadratu *)

rA={-a,-a};
rB={a,-a};
rC={a,a};
rD={-a,a};

(* współrzędne masy punktowej m, czyli punktu P *)

rP={x[t],y[t]};

(* wektor predkosci masy punktowej m *)

vP={x'[t],y'[t]};

(* Odleglosci punktowej masy od wierzcholkow kwadratu, czyli dlugosci rozciagnietych linek *)

AP=Sqrt[ (rP-rA).(rP-rA)] ;
BP=Sqrt[ (rP-rB).(rP-rB)] ;
CP=Sqrt[ (rP-rC).(rP-rC)] ;
DP=Sqrt[ (rP-rD).(rP-rD)] ;

(* dlugosc swobodnej sprężyny, taka sama dla kazdej sprężyny *)

l0=a*Sqrt[2];

(* Energia potencjalna zalezy wylacznie od polozenia punktu P *)

Epot= (1/2)*kA*(AP-l0)^2 +  (1/2)*kB*(BP-l0)^2 + (1/2)*kC*(CP-l0)^2 + (1/2)*kD*(DP-l0)^2 ;

(* Energia kinetyczna to energia kinetyczna punktowej masy *)

Ekin=(1/2)*m*(x'[t]^2+y'[t]^2);

L=FullSimplify[Ekin-Epot];

Etot=FullSimplify[Ekin+Epot];

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

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

d3=D[L,x[t]];

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

d5=D[d1,t];

d6=D[d2,t];

s=FullSimplify[Solve[{d5==d3,d6==d4},{x''[t],y''[t]}]]

Eqs=s/. Rule->Equal

FullSimplify[Etot]

(* CForm[FullSimplify[s]]  *)

(* CForm[FullSimplify[Etot]] *)

(* Parametry układu *)

m=1/10;
a=10;
kA=1;kB=1;kC=1;kD=1;

(* Warunki początkowe: składowe wektora położenia i prędkości masy w chwili t=0 *)

x0=2;y0=4;
vx0=0;vy0=0;

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

(* Czas śledzenia ruchu masy *)

tmax=50;

numerycznie=NDSolve[Join[Flatten[Eqs],WP],{x[t],y[t]},{t,0,tmax}]

numerycznie

Flatten[numerycznie]

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

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

(* Teraz tych funkcji mogę używać tak, jak "zwykłych" funkcji typu Sin[x], Exp[x], itp. *)

x[1]

y[3]

x'[2]

y'[5]

x'''[3]

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

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

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

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

ParametricPlot[{x[t],
y[t]},{t,0,tmax},PlotRange->All,AxesLabel->{"x [m]","y [m]"}]

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

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

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

(* Pierwsza animacja bazuje wprost na ParametricPlot *)

Animate[ParametricPlot[{x[t],
y[t]},{t,0,s},PlotRange->{{-4,4},{-4,4}},AxesLabel->{"x [m]","y [m]"}],{s,0.0001,tmax}]

(* Druga animacja jest innego typu. Używam w niej prostych (wręcz najprostszych) wbudowanych obiektów graficznych takich jak jak linia i koło. *)

(* Przygotowanie sprężyny *)

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

skala=1/4;

kwadrat=Line[{{-a,-a},{a,-a},{a,a},{-a,a},{-a,-a}}];

gkw=Graphics[{Thick,Evaluate[kwadrat]}];

sA[t_]:=Sprezyna6[{-a,-a},{x[t],y[t]},skala]

sB[t_]:=Sprezyna6[{a,-a},{x[t],y[t]},skala]

sC[t_]:=Sprezyna6[{a,a},{x[t],y[t]},skala]

sD[t_]:=Sprezyna6[{-a,a},{x[t],y[t]},skala]

masa[t_]:=Disk[{x[t],y[t]},1/5]

gsA[t_]:=Graphics[sA[t]]

gsB[t_]:=Graphics[sB[t]]

gsC[t_]:=Graphics[sC[t]]

gsD[t_]:=Graphics[sD[t]]

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

Animate[Show[gkw,Evaluate[gsA[t]],Evaluate[gsB[t]],Evaluate[gsC[t]],Evaluate[gsD[t]],Evaluate[gmasa[t]],PlotRange->{{-11,11},{-11,11}},AspectRatio->Automatic],{t,0,tmax}]