How to calculate the difference error between Runge-Kutta Order 4 Method and euler method?

I want to calculate the difference error between Runge-Kutta order 4 method and Euler method. Because I know the Runge-Kutta order 4 method more than precision depend on euler. So, Can I calculate the difference error? This is my Runge-Kutta code. Thanks

tstart = 0;
tend = 180;
dt = 0.01;
T = (tstart:dt:tend).';
Y0 = [10 0 0 0 0 0 0 0];
f = @myode;
Y = fRK4(f,T,Y0);
M1 = Y(:,1);
M2 = Y(:,2);
M3 = Y(:,3);
M4 = Y(:,4);
M5 = Y(:,5);
M6 = Y(:,6);
O  = Y(:,7);
P  = Y(:,8);

figure
%subplot(3,1,1)
%hold on
plot(T,M1,'r','Linewidth',2) 
xlabel('Times (days)')
ylabel('M1 (gr/ml)')

figure
%subplot(3,1,2)
%hold on
plot(T,M2,'b','Linewidth',2)
xlabel('Times (days)')
ylabel('M2 (gr/ml)')

figure
%subplot(3,1,3)
%hold on
plot(T,M3,'g','Linewidth',2)
xlabel('Times (days)')
ylabel('M3 (gr/ml)')

figure
%subplot(3,1,3)
%hold on
plot(T,M4,'b','Linewidth',5)
xlabel('Times (days)')
ylabel('M4 (gr/ml)')

figure
%subplot(3,1,3)
%hold on
plot(T,M5,'r','Linewidth',5)
xlabel('Times (days)')
ylabel('M5 (gr/ml)')

figure
%subplot(3,1,3)
%hold on
plot(T,M6,'g','Linewidth',5)
xlabel('times (days)')
ylabel('M6 (gr/ml)')

figure
%subplot(2,1,1)
%hold on
plot(T,O,'k','Linewidth',2)
xlabel('Times (days)')
ylabel('O (gr/ml)')

figure
%subplot(2,1,2) 
%hold on
plot(T,P,'m','Linewidth',2)
xlabel('Times (days)')
ylabel('P (gr/ml)')

function Y = fRK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for j = 2:N
    t = T(j-1);
    y = Y(j-1,:);
    h  = T(j) - T(j-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(j,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end
function CM1 = myode (~,MM)
  M1 = MM(1);
  M2 = MM(2);
  M3 = MM(3);
  M4 = MM(4);
  M5 = MM(5);
  M6 = MM(6);
  O  = MM(7);
  P  = MM(8);
  
  delta=50;
  gamma=75;
  K1=10^-4;
  K2=5*10^-4;
  K3=10^-3;
  K4=5*10^-3;
  K5=10^-2;
  K6=5*10^-2;
  Ko=0.1;
  n=6;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  mu_2=10^-3;
  mu_3=10^-3;
  mu_4=10^-3;
  mu_5=10^-3;
  mu_6=10^-3;
  mu_o=10^-4;
  mu_p= 10^-5;
  
  sumter=K2*M2+K3*M3+K4*M4+K5*M5;
  
  CM1= zeros(1,8);
  CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*sumter-((Oa-n)*K6*M1*M6)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));
  CM1(2) = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
  CM1(3) = (K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);
  CM1(4) = (K3*M1*M3)-(K4*M1*M4)-(mu_4*M4);
  CM1(5) = (K4*M1*M4)-(K5*M1*M5)-(mu_5*M5);
  CM1(6) = (K5*M1*M5)-(K6*M1*M6)-(mu_6*M6);
  CM1(7) = (K6*M1*M6)-(Ko*M1*O)-(mu_o*O);
  CM1(8) = (Ko*M1*O)-(mu_p*P);
end

Answer 1

[Aquatris]

You can add a 'eul' function for euler integration and take the norm of the difference between the resulting outputs.

clear all ,clc

tstart = 0;
tend = 180;
dt = 0.001;
T = (tstart:dt:tend).';
Y0 = [10 0 0 0 0 0 0 0];
f = @myode;
Y = fRK4(f,T,Y0); % RK4 integration
Y2 = eul(f,T,Y0); % euler integration
 
norm((Y-Y2)./max(Y))

ans = 0.0280

%% here is the euler integration
function Y = eul(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for j = 2:N
    t = T(j-1);
    y = Y(j-1,:);
    h  = T(j) - T(j-1);

    Y(j,:) = y + f(t,y)*h; % y(n+1) = y(n) + dy*dt
  end
end
% here is your RK4
function Y = fRK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for j = 2:N
    t = T(j-1);
    y = Y(j-1,:);
    h  = T(j) - T(j-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(j,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end

function CM1 = myode (~,MM)
  M1 = MM(1);
  M2 = MM(2);
  M3 = MM(3);
  M4 = MM(4);
  M5 = MM(5);
  M6 = MM(6);
  O  = MM(7);
  P  = MM(8);
  
  delta=50;
  gamma=75;
  K1=10^-4;
  K2=5*10^-4;
  K3=10^-3;
  K4=5*10^-3;
  K5=10^-2;
  K6=5*10^-2;
  Ko=0.1;
  n=6;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  mu_2=10^-3;
  mu_3=10^-3;
  mu_4=10^-3;
  mu_5=10^-3;
  mu_6=10^-3;
  mu_o=10^-4;
  mu_p= 10^-5;
  
  sumter=K2*M2+K3*M3+K4*M4+K5*M5;
  
  CM1= zeros(1,8);
  CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*sumter-((Oa-n)*K6*M1*M6)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));
  CM1(2) = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
  CM1(3) = (K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);
  CM1(4) = (K3*M1*M3)-(K4*M1*M4)-(mu_4*M4);
  CM1(5) = (K4*M1*M4)-(K5*M1*M5)-(mu_5*M5);
  CM1(6) = (K5*M1*M5)-(K6*M1*M6)-(mu_6*M6);
  CM1(7) = (K6*M1*M6)-(Ko*M1*O)-(mu_o*O);
  CM1(8) = (Ko*M1*O)-(mu_p*P);
end

[Torsten]

Y(j,:) = y + f(t,y)*h; % y(n+1) = y(n) + dy*dt

instead of

Y(j,:) = y + f(h,y)*h; % y(n+1) = y(n) + dy*dt

in your code.

[Aquatris]

@Torsten Yes you are right, I changed it in the answer. Thank you for pointing it out.

I assumed the myode function was giving the dydt already, since it ignores the first function input. That's why I put h instead as a place holder.

[Torsten]

I assumed the myode function was giving the dydt already, since it ignores the first function input.

Yes, it doesn't matter in this case, but maybe someone else wants to use your function "eul" for a different system of equations that explicitly depends on t.

[cindyawati cindyawati]

thank you for your response @Torsten @Aquatris

[Sam Chak]

Hi @cindyawati cindyawati

If the integration step size is chosen to be extremely small, the error within the finite tspan could become negligible. In this case, the computed error is 10 times smaller. Based on this finding, what conclusions can you draw?

tstart  = 0;
tend    = 180;
f       = @myode;
Y0      = [10 0 0 0 0 0 0 0];

dt      = 0.00001;
T       = (tstart:dt:tend).';
Y       = fRK4(f,T,Y0); % RK4 integration
Y2      = eul(f,T,Y0); % euler integration

error   = norm((Y-Y2)./max(Y))

error = 0.0028