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.