Z Skrypty dla studentów Ekonofizyki UPGOW
Modelowanie: Armata
Mamy: \[\vec F = m \vec a\] więc: \[\vec F = m \vec \ddot x\] w kartezjańskim układzie współrzędnych mamy: \[\begin{cases} F_x = m \ddot x \\ F_y = m \ddot y \end{cases}\] Siły w ogólności zależą od prędkości, czasu i położenia. W przypadku lotu pocisku możemy założyć, że działa na niego siła tarcia oraz siła ciężkości.
\[\vec F = \vec T + \vec Q \] Gdzie Q to ciężar \[\vec Q = \begin{cases} 0 \\ -mg \end{cases}\]
Siła tarcia zależy od prędkości ruchu posicku względem powietrza. Jeżeli będzie wiał wiatr \[\vec u(t) = \begin{cases} u(t) \\0 \end{cases}\] to trzeba jego prękość odjąć od prędkości ruchu pocisku względem ziemi (wiatr pozorny). \[\vec T = - \gamma (\vec v-\vec{u(t)}) = \begin{cases} -\gamma (v_x-u(t) ) \\ -\gamma v_y \end{cases} \]
Czyli mamy układ równań:
\[\begin{cases} \dot x=v_x \\ \dot y=v_y\\ \dot v_x = -\gamma/m \; (v_x-u(t)) \\ \dot v_y=-\gamma/m\; v_y-g \end{cases}\]
Kod źródłowy
clear all close all disp('This is a gun!') function dx = gun(X,T) g=10; m=1; gama=1.5; u=(T>0.1)*10; dx = zeros(4,1); dx(1) = X(3); % x dx(2) = X(4); % y dx(3) = -gama/m*( X(3)-u ); % vx dx(4) = -gama/m*X(4) - g; % vy return end n=125; t = linspace(0,1.22,n); i=1; subplot(2,1,1) hold on alfas=linspace(0.01,pi/2,10); v0=1.; for alfa=alfas X0 = [0,0,v0*cos(alfa),v0*sin(alfa)]; [X,T,MSG]=lsode(@gun,X0,t); plot(X(:,1),X(:,2),'r-') i0=find(X(:,2)<0)(1)-1; i1=find(X(:,2)<0)(1); x_interp(i)=X(i0,1)+( 0 -X(i0,2))/(X(i1,2)-X(i0,2))*(X(i1,1)-X(i0,1)); plot(x_interp(i), 0 ,'g*') i+=1; endfor hold off grid on ylim([0,.05]) xlim([0,.1]) subplot(2,1,2) plot(alfas*180/pi,x_interp,'bo') alfa_fine=linspace(0,90,1000); x_fine=interp1(alfas*180/pi,x_interp,alfa_fine,'cubic'); [zasieg,alfa_i]=max(x_fine); hold on plot(alfa_fine,x_fine) plot(alfa_fine(alfa_i),zasieg,'ro') printf("Najwiekszy zasieg dla strzalu pod katem %f stopni\n",alfa_fine(alfa_i))