## przykład numerycznego  rozwiązania równania różniczkowego  II rzedu 

kill(all)$ reset()$

### Wzory

 P(t)   siła zewntrzna powodująca ruch masy  (zależna od czasu)
 F(t)    siła tarcia  ( zależna od czasu) 

eq0:m*'diff(x(t),t,2)+P+F(t)*signum(diff(x(t),t))+k*x(t)=0;

Powyższe równanie jest trudne do rozwiązania analitycznie (medofą ode3(...))
dlatego  rozwiązuje je numerycznie metodą Ruge Kutty
Równanie eq0 ( rzędu II - bo ma drugą pochodną  x(t)) musi być  zastąpione przez układ dwóch równań różniczkowych rzedu I ( z pierwszą pochodna x(t))

poszukiwane wielkości funkcyjne to:
    położenie x(t)  
      oraz 
    prędkość v(t)= diff(x(t),t)

całość rozwiązywana jest numerycznie metodą Rungego Kutty patrz nizej

Równanie różniczkowe pokazane wyżej  jest zastąpione układem dwu równań 
rózniczkowych rzedu I to jest eq1 i eq2
gdzie zamiast pierwszej pochodnej  x(t)' wprowadzona wprost prędkość v(t),
stąd  druga pochodna x(t) ( przyśpieszenie) wyrażono w postaci pierwszej pochodnej
z v(t)' 

Poniżej wykorzytano Maximę TYLKO do pokazania jakie podstawienia zastosowano ( te opisane w tekście powyżej)
prosze zauważyć, że równanie przekształcono tak aby druga pochodna po czasie z x(t) była po lewej stornie a reszta  po prawej stronie i  przyjęto ogóne oznaczenie  drugiej pochodnej x(t) po czasie jako pierwszej pochodnej  z  pierwszej pochodnej   v(t),    

eq1:diff(x(t),t)=v(t);
eq2:diff(v(t),t)=-P(t)/m-signum(v)*F(t)-k*x(t);

gdzie

diff(x(t),t,2)=diff(v(t),t);

stad  możemy zbudować układ dwóch równań rzedu I
uklad=[eq1,eq2]

### Numeryka 

kill(all);reset();

      k - sztywność 
      m - masa  
      t1 - czas obliczeń
      F=F(t)=  funkcja opisująca zmianę siły  tarcia
      P sila wymuszająca 
      x(t)  poszukiwana funkcja opisująca przemieszczenie
aby  dalsze obliczenia były możliwe muszą być określone dane liczbowe!!

k:1; 
P:0.1;  
m:1.0; 
t1:16; 
F1:1;
x0:2; 
f:0.1;
ω:2*%pi*f
;

Przebieg zmian funkcji opisujacej  F(t)   tu zadana w postaci wzoru wariantowego
ale może byc dowolna  np:
F(t):=F1*sin(t) 

P(t):=P*sin(ω*t);

F(t):=F1*cos(ω*t);

wxplot2d([F(t)],[t, 0,t1], [xlabel,"t"],[ylabel,"F(t)"]);

 poniżej jest układ równań rzedu I zastepujący równanie rzedu II

uklad:[v,-P(t)/m-signum(v)*F(t)-k*x/m];

Warunek początkowy - [ połozenie_pocz, predkość_pocz]

war_pocz:[0.1,0.0]$

tu startuje metoda Rungego Kutty  rk(.......)

[t,tp,tk,krok czasu] - okresla czas dla jakiego są wykonywane obliczenia,
 1/1024  to założony krok czasowy obliczeń 

sol:rk(uklad,[x,v],war_pocz,[t,0,t1,0.01])$

Ponieważ rozwiazanie jest numeryczne i ma postać macierzy
należy wydobyć interesujące nas wyniki 
w pierwszej  kolumnie jest czas
w drugiej położenie x(t)
w trzeciej prędkość v(t)_

wykreslenie zależności x(t) 

sol_pts:makelist([first(sol[i]),second(sol[i])],i,1,length(sol))$

wxplot2d([discrete,sol_pts], [xlabel,"t"],[ylabel,"x(t)"])$

predkośc  v(t)

sol_vel:makelist([first(sol[i]),third(sol[i])],i,1,length(sol))$

wxplot2d([discrete,sol_vel], [xlabel,"t"],[ylabel,"v(t)"])$

portret Fazowy

sol_Mphase:makelist([second(sol[i]),third(sol[i])],i,1,length(sol))$

wxplot2d([discrete,sol_Mphase], [xlabel,"x(t)"],[ylabel,"v(t)"]);