Comparing results of "ODE45" function versus 4th order Runge-Kutta for Influenza disease in MATLAB

Question

I have written a code for a system of nonlinear differential equations which models the Influenza disease. I was surprisingly shocked when I wanted to simulate this system with 4th order Runge-Kutta and ode45function. The problem is that the results of these 2 methods are not the same. In fact, when I write a code for 4th order Runge-Kutta algorithm, the states will become larger and larger in each step time till they reach infinity. But the results are reasonable and rational when I simulate the system with ode45function. I do not think the algorithm I wrote is incorrect. Can someone explain to me where I'm doing wrong In case if something is not valid in the coding?.

The codes I have written are as follows :

clc;clear;close all;

t = linspace(0,100,101)';
global alpha eta p f ep delta q beta kessi

kessi = 0.526;                % Transition rate for the exposed
alpha = 0.244;                % Recovery rate for the (symptomatic)
eta = alpha;                  % Infected
p = 0.667;                    % Recovery rate for the asymptomatic
f = 0.98;                     % 1-fatallity rate(one minus fatality rate)
ep = 0;                       % Infectivity reduction factor for the exposed
delta = 1;                    % Infectivity reduction factor for the asymptomatic
q = 0.5;                      % Contact reduction by isolation
beta = 1;                     % From CALCULATIONS

Influenza_Fun = @(t,x,u) [-beta*x(1)*(ep*x(2)+(1-q)*x(3)+delta*x(4))-u(1)*x(1)
                           beta*x(1)*(ep*x(2)+(1-q)*x(3)+delta*x(4))-kessi*x(2)
                            p*kessi*x(2)-alpha*x(3)-u(2)*x(3)
                            (1-p)*kessi*x(2)-eta*x(4)
                            f*alpha*x(3)+u(2)*x(3)+eta*x(4)+u(1)*x(1)];


% Initial Condtions
Ics = [15000 200 500 300 0]';                     % [S0 E0 I0 A0 R0]';
u = [0 0]';

[~,States] = ode45(@(t,x)Influenza_Fun(t,x,u),t,Ics);

for i=1:5
   plot(t,States(:,i),'LineWidth',2)
   hold on
   grid on

end
legend('S','E','I','A','R')
legend show

The above code will simulate the system based on ode45function. The following code is my RK4 algorithm :

clc;clear;close all;

% { Control of Influenza disease based on Genetic Algorithm using desired Objective FUnction}

global alpha eta p f ep delta q beta kessi


kessi = 0.526;                % Transition rate for the exposed
alpha = 0.244;                % Recovery rate for the (symptomatic)
eta = alpha;                  % Infected
p = 0.667;                    % Recovery rate for the asymptomatic
f = 0.98;                     % 1-fatallity rate(one minus fatality rate)
ep = 0;                       % Infectivity reduction factor for the exposed
delta = 1;                    % Infectivity reduction factor for the asymptomatic
q = 0.5;                      % Contact reduction by isolation
beta = 1;                     % From CALCULATIONS

% Initial Condtions
Ics = [15000 200 500 300 0]';                     % [S0 E0 I0 A0 R0]';

max_days = 100;
y(:,1) = Ics;
t = zeros(1,max_days);

h = 1;
t(1) = 0;
LB = [0 0]';
UB = [1 1]';
g = @yMatrix; 
u = [0 0]';

Influenza_Fun = @(t,x,u) [-beta*x(1)*(ep*x(2)+(1-q)*x(3)+delta*x(4))-u(1)*x(1)
                           beta*x(1)*(ep*x(2)+(1-q)*x(3)+delta*x(4))-kessi*x(2)
                            p*kessi*x(2)-alpha*x(3)-u(2)*x(3)
                            (1-p)*kessi*x(2)-eta*x(4)
                            f*alpha*x(3)+u(2)*x(3)+eta*x(4)+u(1)*x(1)];


%% 2. mAiN Loop

for time = 1:max_days


    k1 = Influenza_Fun(t(time),y(:,time),u);
    k2 = Influenza_Fun(t(time)+(h/2),y(:,time)+k1,u);
    k3 = Influenza_Fun(t(time)+(h/2),y(:,time)+k2,u);
    k4 = Influenza_Fun(t(time)+h,y(:,time)+k3,u);

    % y(i+1) = y(i) + (h/6)*(k1+2k2+2k3+k4)
    y(:,time+1) = y(:,time)+(h/6)*(k1+2*k2+2*k3+k4);
    t(time+1) = t(time)+h;

end

y = y';

for i=1:size(y,2)

    plot(t,y(:,i),'LineWidth',2)
    hold on
    grid on


end

score 0 · Answer 1 · answered Aug 21 '21 at 07:44

0

It should be

    k2 = Influenza_Fun(t(time)+(h/2),y(:,time)+(h/2)*k1,u);

etc.

Also use some way to exhibit the time steps that are actually used by ode45, for instance using the [ T,Y] = output format, discounting the extra points from the Refine option. Or use the function evaluation count, at least the last call variant in the documentation should exhibit this. Use a similar step size or function call count in RK4 to get comparable results.

answered Aug 21 '21 at 07:44

Lutz Lehmann

25,219
2
22
51

That is correct. We can change the time step for better results. What I'm asking is that why the results of "ode45" and RK4 should be different from each other when the step time for both is 1? – MMd.NrC Aug 21 '21 at 10:24
1

You need to study closer how ode45 works. The time step there is most surely not 1, this is only the time step for the points where the output is produced. You do not get any indication of how many internal steps were taken. To simulate this with RK4 you could compute the solution with time step 0.1 and only include every 10th step in the output array. – Lutz Lehmann Aug 21 '21 at 11:03
So You mean that for example when I use t = linspace(0,100,101) with step time 1, ode45 will solve the solution just for this range or it takes an arbitrary step time. Does this time vector with its step time even will be included in solving after all? – MMd.NrC Aug 21 '21 at 11:17
No, this output time array will not influence the internal solver steps. After each internal step the solver checks if it passed the next output time, and interpolates the corresponding value from the "dense output" polynomial. – Lutz Lehmann Aug 21 '21 at 11:20
Thank You so much for the Explanations. So what do You suggest me to do if I want to get same results for each method? – MMd.NrC Aug 21 '21 at 15:53

Comparing results of "ODE45" function versus 4th order Runge-Kutta for Influenza disease in MATLAB

1 Answers1