Rozważamy klocek, który porusza się w jednym wymiarze pod wpływem dwóch sił:
(1) siły pochodzącej od sprężyny o stałej sprężystości k,
(2) siły tarcia o własnościach znanych z codziennego doświadczenia. Ta siła tarcia "dopasowuje się" do siły sprężyny i nie pozwala ruszyć klocka z miejsca, jeśli |k*x| ≤ Fs.  Gdy klocek się porusza, siła tarcia jest oczywiście skierowana przeciwnie do kierunku prędkości i ma wartość Fk < Fs.

(* I wersja, gdzie xp ma być dokładnie równe zer0 *)

Ft1[x_,xp_]:=Piecewise[{{k*x,Abs[k*x]≤Fs&&xp==0},{Fs*Sign[x],Abs[k*x]>Fs&&xp==0},{-Fk*Sign[xp],xp≠ 0}}];

(* II wersja, gdzie wystarczy, że xp jest odpowiednio (parametr eps) bliskie zeru *)

eps=10^(-10);

Ft[x_,xp_]:=Piecewise[{{k*x,Abs[k*x]≤Fs&&Abs[xp]≤ eps},{Fs*Sign[x],Abs[k*x]>Fs&&Abs[xp]≤ eps},{-Fk*Sign[xp],Abs[xp]>0}}];

(* Dane liczbowe *)

m=1/2; (* kg *)

Fs=1/5; (* N *)

Fk=1/7; (* N *)

k=1; (* N/cm *)

(* Dalsza część notatnika pochodzi z Przyklad_6_oscylator_harmoniczny_prosty_OK.nb *)

x0=10; (* cm *)

v0=1/1000;  (* cm/s *)

(* Wartość tmax została ustawiona metodą prób i błędów w taki sposób, by x'[tmax] była wciąż jeszcze różna od zera. Oczywiście tmax zależy od wartości Fs, Fk i k. Po czasie odpowiadającym zatrzymaniu oscylatora, Mathematica ma problem z numerycznym rozwiązaniem ! *)

tmax=77.75;

numerycznie=NDSolve[{m*x''[t]==-k*x[t]+Ft[x[t],x'[t]],x[0]==x0,x'[0]==v0},x[t],{t,0,tmax}]

Flatten[numerycznie]

(* Wyciągam funkcję x[t] z numerycznego rozwiązania *)

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

(* Teraz mogę używać funkcji x[t] tak,jak innych funkcji *)

x[0.1]

x'[0.2]

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

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

(* Modelowanie ruchu masy przy pomocy wbudowanych obiektów graficznych. Ruch sprężyny przygotowany został od podstaw w notebooku sprezyna.nb. Teraz korzystam jedynie z wyników. *)

(* pomocnicza funkcja do zdefiniowania łamanej udającej sprężynę *)

g[xx_,a_]:=If[xx≥a/18&&xx≤17*a/18,(1/3)*Sin[18*Pi/a*(xx-a/18)],0]

Plot[g[xx,9],{xx,0,9},Axes->True]

tab1[a_]:=Simplify[Table[{xx,g[xx,a]},{xx,0,a,a/36}],a>0]

spr1[a_]:=Line[tab1[a]];

tab2[xxp_,yyp_,xxk_]:=Transpose[{Transpose[tab1[xxk-xxp]][[1]]+xxp,Transpose[tab1[xxk-xxp]][[2]]+yyp}];

(* pozioma sprężyna o początku w punkcie (xxp,yyp) i końcu w punkcie (xxk,yyp)  *)

SprezynaPozioma[xxp_,yyp_,xxk_]:=Line[tab2[xxp,yyp,xxk]]

(* Długość swobodnej sprężyny *)

x1=15;

f[t_]:=x1+x[t]; (* x2: położenie końca sprężyny w czasie ruchu masy m; ta wielkość musi być dodatnia, by symulacja wyglądała rozsądnie  *)

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

spr2[xxp_,yyp_,t_]:=SprezynaPozioma[xxp,yyp,f[t]]

(* pozioma sprężyna o początku w punkcie (xxp,yyp) i końcu opisywanym przez funkcję f[t] *)

gspr2[xxp_,yyp_,t_]:=Graphics[spr2[xxp,yyp,t]]

masa[t_]:=Line[{{f[t],-1},{f[t]+1,-1},{f[t]+1,1},{f[t],1},{f[t],-1}}];

gmasa[t_]:=Graphics[{Thick,masa[t]}];

(* kreska pionowa w miejscu umocowania sprężyny *)

pion[xxp_]:=Line[{{xxp,-1},{xxp,1}}]

gpion[xxp_]:=Graphics[{Thick,pion[xxp]}]

(* pozioma powierzchna, po której przesuwa się masa m; zakres tej linii powinien być określony przez minimalną i maksymalną wartość f[t] dla t ∈ [0,tmax] *)

poziom=Line[{{-2,-1.1},{26,-1.1}}];

gpoziom=Graphics[{Thick,poziom}];

Animate[Show[gspr2[0,0,t],gpoziom,gpion[0],gmasa[t],PlotRange->{{-2,26},{-1.1,1}}],{t,0,tmax},AnimationRepetitions->1,AnimationRunning->False]