Cześć,
Mam program do numerycznego obliczania układów równań różniczkowych metodą Rungego Kutty 4 rzędu.
Równania te opisują lot pocisku artyleryjskiego, więc można mniej więcej przewidzieć prawidłowe przebiegi wykresów f(X,Y) czy f(kąt nachylenia, t).
Oto dwa programy, do których chciałabym się odnieść (różnice występują jedynie w pętli for, reszta bez zmian):
clc;
clear all;
%PARAMETRY OBLICZEŃ
tp=0;
tk=100;
h=0.1;
%WARUNKI POCZĄTKOWE
x=0;
y=0;
v=250;
k=40; %kąt nachylenia lufy podany w stoniach
kat=(k*pi)/180; %przeliczenie kąta na radiany
t=tp;
j=1;
for i=tp:h:tk
v1=row1(v,kat);
v2=row1(v+0.5*h , kat+0.5*v1);
v3=row1(v+0.5*h , kat+0.5*v2);
v4=row1(v+h , kat+v3);
kat1=row2(v,kat);
kat2=row2(v+0.5*h , kat+0.5*kat1);
kat3=row2(v+0.5*h , kat+0.5*kat2);
kat4=row2(v+h , kat+kat3);
x1=row3(v,kat);
x2=row3(v+0.5*h , kat+0.5*x1);
x3=row3(v+0.5*h , kat+0.5*x2);
x4=row3(v+h , kat+x3);
y1=row4(v,kat);
y2=row4(v+0.5*h , kat+0.5*y1);
y3=row4(v+0.5*h , kat+0.5*y2);
y4=row4(v+h , kat+y3);
v=v+(h/6)*(v1+2*v2+2*v3+v4);
kat=kat+(h/6)*(kat1+2*kat2+2*kat3+kat4);
x=x+(h/6)*(x1+2*x2+2*x3+x4);
y=y+(h/6)*(y1+2*y2+2*y3+y4);
if y <= 0
v=0;
else
wektorx(j)=x;
wektory(j)=y;
wektorv(j)=v;
wektorkat(j)=(kat*180)/3.14;
wektort(j)=t;
t=t+h;
j=j+1;
end
end
zasieg=max(wektorx)
wysokoscmaksymalna=max(wektory)
subplot(2,1,1)
plot(wektorx,wektory);
xlabel('Przemieszczenie pocisku względem osi X [m]')
ylabel('Przemieszczenie pocisku względem osi Y [m]')
title('Wykres obrazujący tor lotu pocisku')
grid on
subplot(2,1,2)
plot(wektort,wektorv);
xlabel('Czas lotu pocisku [s]')
ylabel('Prędkość pocisku [m/s]')
title('Wykres obrazujący zmiany prędkosci pocisku względem czasu')
grid on
for i=tp:h:tk
v1=row1(v,kat);
kat1=row2(v,kat);
x1=row3(v,kat);
y1=row4(v,kat);
v2=row1(v+(h/2),kat+(h/2)*kat1);
kat2=row2(v+(h/2),kat+(h/2)*kat1);
x2=row3(v+(h/2),kat+(h/2)*kat1);
y2=row4(v+(h/2),kat+(h/2)*kat1);
v3=row1(v+(h/2),kat+(h/2)*kat2);
kat3=row2(v+(h/2),kat+(h/2)*kat2);
x3=row3(v+(h/2),kat+(h/2)*kat2);
y3=row4(v+(h/2),kat+(h/2)*kat2);
v4=row1(v+h,kat+h*kat3);
kat4=row2(v+h,kat+h*kat3);
x4=row3(v+h,kat+h*kat3);
y4=row4(v+h,kat+h*kat3);
x=x+(h/6)*(x1+2*x2+2*x3+x4);
y=y+(h/6)*(y1+2*y2+2*y3+y4);
v=v+(h/6)*(v1+2*v2+2*v3+v4);
kat=kat+(h/6)*(kat1+2*kat2+2*kat3+kat4);
if y <= 0
v=0;
else
wektorx(j)=x;
wektory(j)=y;
wektorv(j)=v;
wektorkat(j)=(kat*180)/3.14;
wektort(j)=t;
t=t+h;
j=j+1;
end
end
Pierwszy program nie działa (generuje wykresy których w żaden sposób nie można porównać z rzeczywistą trajektorią lotu) - mimo że sposób wyznaczania współczynników k1,k2,k3,k4 został wiernie odtworzony ze skryptu :)
Mianowicie:
k1=f(x,y)
k2=f(x +0.5h, y+0,5k1)
k3=f(x +0.5h, y+0.5k2)
k4=f(x+h, y+k3)
f(i+1)=f(i)+(h/6)(k1 + 2k2 + 2*k3 + k4)
gdzie h- krok całkowania, f- rozpatrywana funkcja.
Druga wersja działa natomiast bez problemów, wykresy układają się wiarygodnie, tak jakbyśmy tego oczekiwali.
Nie rozumiem jednak dlaczego współczynniki dobrano w taki sposób. W żaden sposób nie odzwierciedlają one wzoru który znalazłam w literaturze.
Będę bardzo wdzięczna, jeżeli ktoś znajdzie chwilę czasu aby mi to wytłumaczyć.
Pozdrawiam.