I need to solve flight of the shell problem (get XY trajectory using Runge-Kutta method) and at first check if method works on some simple test system.
So here's Runge-Kutta stuff
$k_1 = f(t_n, y_n)$
$k_2 = f(t_n + h/2, y_n + hk_1/2) $
$k_3 = f(t_n+h, y_n - hk_1 + 2hk_2)$
$y_{n+1} = y_n + h(k_1 + 4k_2 + k_3)/6$
where $h$ is step
Here's my test system
$y'_1 = -5y_1 - 10y_2 + 14e^{-x}$
$y'_2 = -10y_1 - 5y_2 + 14e^{-x}$
approximately it's $y1 = y2 = e^{-x}$
So I solve it on $[0;4]$ and checked it (by comparison with $y1 = y2 = e^{-x}$) (btw approximation is very good even with h=0.1). Here's the matlab code, if you interested in
Rnge-Kutta method
function [ res_y ] = RungeKutta(dim, size, grid, step, f1,f2,y1, y2)
k1=zeros(dim);
k2=zeros(dim);
k3=zeros(dim);
h = step;
res_y(1,1) = y1;
res_y(2,1) = y2;
for i=1: size
k1(1)= f1(grid(i),y1,y2);
k1(2)= f2(grid(i),y1,y2);
k2(1)= f1(grid(i)+h/2, y1+h*k1(1)/2, y2+h*k1(2)/2);
k2(2)= f2(grid(i)+h/2, y1+h*k1(1)/2, y2+h*k1(2)/2);
k3(1)= f1(grid(i)+h, y1-h*k1(1)+2*h*k2(1), y2-h*k1(2)+2*h*k2(2));
k3(2)= f2(grid(i)+h, y1-h*k1(1)+2*h*k2(1), y2-h*k1(2)+2*h*k2(2));
res_y(1,i+1) = y1 + h*(k1(1) + 4*k2(1) + k3(1))/6;
res_y(2,i+1) = y1 + h*(k1(2) + 4*k2(2) + k3(2))/6;
y1 = res_y(1,i+1);
y2 = res_y(2,i+1);
end
end
Main method
a = 0; b = 4;
h = 0.1; % step
t = a:h:b; %grid
n = 2;
m = size(t,2);
hold on;
plot(t, exp(-t),'b-')
plot(t, exp(-t),'r--')
hold off;
y1=1; y2 = 1;
f1_ptr = @f1;% out = -5 * y1 - 10 * y2 + (14)*exp(-x);
f2_ptr = @f2;% out = -10 * y1 - 5 * y2 + (14)*exp(-x);
res_y = RungeKutta(n,m-1,t,h,f1_ptr, f2_ptr,1,1);
hold on;
plot(t,res_y);
hold off;
f1 function
function [ out ] = f1( x, y1, y2, alpha, beta )
if nargin == 3
alpha = 5;
beta = 10;
end
out = -alpha * y1 - beta * y2 + (alpha + beta - 1)*exp(-x);
end
f2 function
function [ out ] = f2( x, y1, y2, alpha, beta )
if nargin == 3
alpha = 5;
beta = 10;
end
out = -beta * y1 - alpha * y2 + (alpha + beta - 1)*exp(-x);
end
But my main system, for the shell trajectory look like that
$x'' = - \frac{1}{2m} C \rho S \cos(\theta) v^2$
$y'' = - \frac{1}{2m} C \rho S \sin(\theta) v^2 - g$
Where's $v = \sqrt{x'^2 + y'^2}$ and angle $\theta = arctg(y'/x')$ and $C, \rho, S, m, g$ some constants.
Also I've got some initial conditions
$x(0) = 0, y(0) = 0, \theta(0) = 0.6, v(0) = 50 m/s$
So I've stuck here. I don't understand clearly how to solve this system using Runge-Kutta method. I mean, I don't understand how to code it in such situation.
I'm not good at math, I am a coder. Please, colud you help me?
UPD Thanks to Tony I think I began to understand problem. But method I wrote doesn't work :c
Here it is
function [ shell_result ] = RungeKuttaShell(step, y1,y2,y3,y4)
h = step; norma = 1;
i = 1;
shell_result(1,1) = y1;
shell_result(2,1) = y2;
shell_result(3,1) = y3;
shell_result(4,1) = y4;
while (norma > 0) %how much iterations?
k1(1)= fz(y3);
k1(2)= fw(y4);
k1(3)=fx(k1(1), k1(2));
k1(4)=fy(k1(1), k1(2));
k2(1)= fz( y3+h*k1(1)/2);
k2(2)= fw( y4+h*k1(2)/2);
k2(3)=fx(k2(1),k2(2));
k2(4)=fy(k2(1),k2(2));
k3(1)= fz(y3-h*k1(1)+2*h*k2(1));
k3(2)= fw(y4-h*k1(2)+2*h*k2(2));
k3(3)=fx(k3(1),k3(2));
k3(4)=fy(k3(1),k3(2));
shell_result(1,i+1) = y1 + h*(k1(1) + 4*k2(1) + k3(1))/6;
shell_result(2,i+1) = y2 + h*(k1(2) + 4*k2(2) + k3(2))/6;
shell_result(3,i+1) = y3 + h*(k1(3) + 4*k2(3) + k3(3))/6;
shell_result(4,i+1) = y4 + h*(k1(4) + 4*k2(4) + k3(4))/6;
y1 = shell_result(1,i+1);
y2 = shell_result(2,i+1);
y3 = shell_result(3,i+1);
y4 = shell_result(4,i+1);
i = i+1;
norma = y4; %height?
end
end
functions look like this
function [ output ] = fx( z,w)
m=15; C=0.2; rho=1.29; S=0.25; g = 9.81;
output = -((1/2*m)*C*rho*S)*z*sqrt(z*z + w*w);
end
function [ output ] = fy( z,w )
m=15; C=0.2; rho=1.29; S=0.25; g = 9.81;
output = -((1/2*m)*C*rho*S)*w*sqrt(z*z + w*w) - g;
end
function [ output_args ] = fw( input_args )
output_args = input_args;
end
function [ output_args ] = fz( input_args )
output_args = input_args;
end
And I use it just like that
shell = RungeKuttaShell(0.1,0,0,v0*cos(th0), v0*sin(th0));
zeros(dim)
generates a square matrix of dimensionsdim
bydim
. If you want to create a column vector, you should usezeros(dim,1)
. $\endgroup$