# Confused by Runge-Kutta method

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));

• I'll take a detailed look at this code after lunch. But one passing note, zeros(dim) generates a square matrix of dimensions dim by dim. If you want to create a column vector, you should use zeros(dim,1). – Emily Nov 5 '13 at 16:51
• $v^2=z^2+w^2\,$ of course. – Tony Piccolo Nov 5 '13 at 23:03
• You need $\cos\theta=z/v \;$ and $\; \sin\theta=w/v\,$.$\;$ So ... – Tony Piccolo Nov 6 '13 at 6:28
• My first code was wrong: try my final code. – Tony Piccolo Nov 16 '13 at 7:56
• At the end, if all is well, it should be fine to plot the resistive trajectory together with the unresistive one (C=0) in a uniform scale. – Tony Piccolo Nov 16 '13 at 11:49

The system is $$\begin {cases} x'=z \\ \\y'=w \\ \\ z'=-k\,z\sqrt{z^2+w^2} \\ \\w'=-k\,w\sqrt{z^2+w^2}-g \end {cases}$$with$$k=\frac 1{2m}C\rho S$$The initial values are $x(0)=0\;,\; y(0)=0\;,\;z(0)=v(0)\cos\theta(0)\;,\;w(0)=v(0)\sin\theta(0)\,$.

You must be clever in recoding your procedure.
The system now is $$\begin {cases} y'_1=f_1(t,y_1,y_2,y_3,y_4) \\ \\y'_2=f_2(t,y_1,y_2,y_3,y_4) \\ \\ y'_3=f_3(t,y_1,y_2,y_3,y_4) \\ \\y'_4=f_4(t,y_1,y_2,y_3,y_4) \end {cases}$$The previous one was$$\begin {cases} y'_1=f_1(t,y_1,y_2) \\ \\y'_2=f_2(t,y_1,y_2) \end {cases}$$You must recode completely without any reference to your particular system (this only conditions the choice of the functions).

Now $y_1=x\;$,$\;y_2=y\;$,$\;y_3=z\;$,$\;y_4=w$.

Edit

The choice of step value is not so important if one is not interested in saving resources (of time and space of memory). Do the changes in my comment and try.

If you don't get a range of 176 m, try this

final code

k2(1)= fz( y3+h*k1(3)/2);
k2(2)= fw( y4+h*k1(4)/2);
k2(3)=fx(k2(1),k2(2));
k2(4)=fy(k2(1),k2(2));

k3(1)= fz(y3-h*k1(3)+2*h*k2(3));
k3(2)= fw(y4-h*k1(4)+2*h*k2(4));
k3(3)=fx(k3(1),k3(2));
k3(4)=fy(k3(1),k3(2));

• ah, i guess I started to get it! but how I suppose to recode in that way? I mean in my test system $f_1$ and $f_2$ equations were use $y_1$ and $y_2$ parameters, thats why I made it like $y'_1 = f_1(t, y_1, y_2)$ – Danil Gholtsman Nov 6 '13 at 11:03
• and by the way is there 5 parameters necessery in the main system? I mean in main system there are no any correlation with $t$, isn't it depends on $w$ and $z$ only(or $y_3, y_4$)? – Danil Gholtsman Nov 6 '13 at 11:04
• I mean you add the code to your question. – Tony Piccolo Nov 6 '13 at 20:06
• Read this and you are done. – Tony Piccolo Nov 7 '13 at 9:42
• It is matlab not c. – Tony Piccolo Nov 7 '13 at 19:52