# Runge-Kutta method accuracy

I got Runge-Kutta method here and I solve this system using it.

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}$

UPD: Initial conditions here is $x0_1 = x0_2 = 1$

So I need to solve it on $[0;4]$. Well, I was thought I solve it right, becasue of I checked how exact solution and these approximate solution plots looks like(on the left it's I zoomed plot until saw difference)

Also I checked how looks plot of difference between exact solution and approximate solution depending on step (let's call it e/h) looks like.

So e/h it looks like this

But when I checked e/h^4 dependence it looked like these

I show it to the teacher and she said that my solution is wrong, it's not suppose to be like these! I show my code to her asked for help but she said that she don't understands matlab :c

Am I really done something worng? And if yes what I've done wrong? And if not how to prove that I'm right?

Here's my code btw

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;

%e/h and e/h^4 plots

fig_a = figure;
set(fig_a,'name','e/h','numbertitle','off')
hold on;

counter = 0;

for h=0.001:0.01:0.1
y1=1; y2 = 1;
t = a:h:b;
m = size(t,2);
counter = counter + 1;

result_appr = RungeKutta(n,m-1,t,h,f1_ptr, f2_ptr,y1,y2);
result_exact = exp(-t);

result_difference = abs(result_appr(1, :) - result_exact);

e1(counter) = max(result_difference);
e2(counter) = max(result_difference);

hh = h*h*h*h;

ehh1(counter)=e1(counter)/hh;
ehh2(counter)=e2(counter)/hh;

end;

h=0.001:0.01:0.1;
plot(h,e1,'c');
plot(h,e2,'c');
hold off;

fig_b = figure;
set(fig_b ,'name','e/h^4','numbertitle','off')

hold on;
plot(h,ehh1,'r')
plot(h,ehh2,'b')
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

• Where are the initial conditions? – Amzoti Nov 23 '13 at 14:49
• @Amzoti oops sorry, I'll update question – Danil Gholtsman Nov 23 '13 at 14:57
• @Amzoti I updated it! – Danil Gholtsman Nov 23 '13 at 15:07
• @DanilGholtsman. When you have such a problem, after a few hours, you are no more able to see anything. In this case, ask a friend to look at it. Trust me, I starting in computing sciences 53 years ago. Does this change your results ? Please reply and post. Glad to have been able to help (hoping that this is the case). – Claude Leibovici Nov 23 '13 at 16:00
• I need to know where the test about $e/h^4$ comes from for a third order Runge-Kutta method: can your teacher give you some written reference ? – Tony Piccolo Nov 25 '13 at 17:27