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

Question

cindyawati cindyawati 2024-3-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2089871-how-to-calculate-the-difference-error-between-runge-kutta-order-4-method-and-euler-method

评论： Sam Chak 2024-3-4

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

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Aquatris 2024-3-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2089871-how-to-calculate-the-difference-error-between-runge-kutta-order-4-method-and-euler-method#answer_1420451

编辑：Aquatris 2024-3-4

在 MATLAB Online 中打开

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

cindyawati cindyawati 2024-3-4

thank you for your response @Torsten @Aquatris

Sam Chak 2024-3-4

在 MATLAB Online 中打开

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