Runge Kutta method computational cost

I have written a Matlab script that performs 2nd and 4th order Runge-Kutta method for different values of the step h.

h = [0.5, 0.2, 0.05, 0.01]

It could be expected that, the lower the step, the greater the computational time needed to perform the method. However, the code I have written does not follow that law and I cannot seem to find the reason behind it. It could be my code or my own computer. Thank you in advance.

The script for the 2nd order method is the following:

i     = 1;
x_end = zeros(1, 4);
ct    = zeros(1, 4);

for h = [0.5, 0.2, 0.05, 0.01] 
    
    tic
    t    = t0:h:tf;
    x    = zeros(1, length(t));
    x(1) = x0;
    k    = 1;
    while k < length(t)
        
        xkp1P  = x(k) + h  * f(x(k),t(k));
        x(k+1) = x(k) + h/2*(f(x(k),t(k)) + f(xkp1P,t(k+1)));
        k      = k + 1;

    end
    ct(i)    = toc;                                                        % Computational time [s]
    x_end(i) = x(end);                                                            
    figure(1)
    plot(t, x)
    hold on
    
    i = i + 1;    
end

The script for the 4th order method is the following:

i     = 1;
x_end = zeros(1, 4);
ct    = zeros(1, 4);

for h = [0.5, 0.2, 0.05, 0.01] 
   
    t    = t0:h:tf;
    x    = zeros(1, length(t));
    x(1) = x0;   
    k    = 1;
    tic
    while k < length(t)
        
        K1 = f(x(k)         , t(k)      );
        K2 = f(x(k)+0.5*h*K1, t(k)+0.5*h);
        K3 = f(x(k)+0.5*h*K2, t(k)+0.5*h);
        K4 = f(x(k)+    h*K3, t(k)+    h);
        
        x(k+1) = x(k) + h/6*(K1 + 2*K2 + 2*K3 + K4);
        k  = k + 1;

    end
    ct(i)    = toc;         % Computational time [s]  
    x_end(i) = x(end);                                                         
    figure(2)
    plot(t, x)
    hold on
    
    i = i + 1;    
end

Answer 1

[John D'Errico]

tic and toc are really bad ways to estimate the time required for (almost) anything, unless the time is pretty long, and you only want a poor estimate.

So use timeit instead. That is, you would need to wrap your code inside a function, calling timeit on that function.

timeit worries about warmup problems with the code, since the first time you call any function, it gets cached. That takes extra time. It also calls your code many times, taking an average of the calls.

[Rafael Félix Soriano]

It does work better, but it still does not work for h = 0.2. It does however work perfectly when taking a longer time interval (from 0 to 8 for example).

The modified code for RK2 with timeit is as follows:

i     = 1;
x_end = zeros(1, 4);
ct    = zeros(1, 4);

for h = [0.5, 0.2, 0.05, 0.01] 
    
    t    = t0:h:tf;
%     x    = zeros(1, length(t));
%     x(1) = x0;
%     k    = 1;
%     while k < length(t)
%         
%         xkp1P  = x(k) + h  * f(x(k),t(k));
%         x(k+1) = x(k) + h/2*(f(x(k),t(k)) + f(xkp1P,t(k+1)));
%         k      = k + 1;
% 
%     end
h
    [x]   = RK2(x0, t0, tf, h);
    f = @() RK2(x0, t0, tf, h);
    ct(i) = timeit(f)
    

%     ct(i)    = toc;                   % Computational time [s]
    x_end(i) = x(end);                                                            
    figure(4)
    plot(t, x)
    hold on
    
    i = i + 1;    
end

function [x] = RK2(x0, t0, tf, h)
    
    f    = @(x,t) x - t^2 + 1;
    t    = t0:h:tf;
    x    = zeros(1, length(t));
    x(1) = x0;
    k    = 1;
    while k < length(t)
        
        xkp1P  = x(k) + h  * f(x(k),t(k));
        x(k+1) = x(k) + h/2*(f(x(k),t(k)) + f(xkp1P,t(k+1)));
        k      = k + 1;

    end

end

Answer 2

[James Tursa]

Side Note, You should not be calculating f(x(k),t(k)) twice in your 2nd order method. You should be doing it like your 4th order method. E.g., this

        xkp1P  = x(k) + h  * f(x(k),t(k));
        x(k+1) = x(k) + h/2*(f(x(k),t(k)) + f(xkp1P,t(k+1)));

should be something like this instead:

        K1 = f(x(k),t(k));
        xkp1P  = x(k) + h  * K1;
        x(k+1) = x(k) + h/2*(K1 + f(xkp1P,t(k+1)));

Answer 3

[Meysam Mahooti]

https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle

https://www.mathworks.com/matlabcentral/fileexchange/61130-runge-kutta-fehlberg-rkf78?s_tid=srchtitle