SIR model with recovered individuals may lose their immunity and become reinfected with the disease. But came with a failure about integration tolerances

15 次查看（过去 30 天）
Julian Pereira
Julian Pereira 2022-9-20
评论： Bjorn Gustavsson 2022-9-28
Hi everyone,
I am facing a very intrigating problem atm. I have this SIR model with individuals may lose their immunity and become reinfected with the disease. I tried to use a SIR simulater from matlab community but it seems it didn't work too. So I am trying to do myself. After I run the code, come with an warming:
Warning: Failure at t=4.135655e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.110223e-16) at time t.
In ode45 (line 352)
In Q10_SIR_model (line 15)
See bellow my code:
%% Q10_SIR_model
clc; close all ; clear;
% time
tspan = 0:0.01:100;
% Initial Conditions (Suceptible, Infected, Recovered)
q0 = [99; 1; 0];
% Parameters
beta = 0.5; gamma = 20; delta = 0.4;
%% Model SIR
% Initial conditions as vector
p = [beta; gamma; delta];
[t, y] = ode45(@(t,q)fc_ode_SIR(t,q,p),tspan,q0);
%% Plot
figure;
subplot(2,1,1)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]');
legend('Suceptible','Infected', 'Recovered', 'Location','Best');
y = y(:,3);%*N;
% figure(2);
subplot(2,1,2)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]');
legend('Infected', 'Location','Best');
And my script:
%% fc_ode_SIR
function dq = fc_ode_SIR(~,q,p) % q = initial conditions (Suceptible, Infected, Recovered)
% p = parameter (beta, gamma, delta, initial conditions)
beta = p(1); gamma = p(2); delta = p(3);
S = q(1); I = q(2); R = q(3);
% SIR Model
S = beta*S*I + delta*R;
I = beta*S*I - gamma*I;
R = gamma*I - delta*I;
dq = [-S;
I;
R];
end
I don't know where exactaly is my mistake as I run other SIR models and it was fine. I am quite sure it's something silly but if someone could find that I'd appreciate.
Cheers
  2 个评论
Julian Pereira
Julian Pereira 2022-9-28
Final version:
%% Q10_SIR_model
clc; close all ; clear;
% time
tspan = 0:0.01:100;
% Initial Conditions (Suceptible, Infected, Recovered)
q0 = [99; 1; 0];
% Parameters
beta = 0.5; gamma = 20; delta = 0.4;
%% Model SIR
% Initial conditions as vector
p = [beta; gamma; delta];
[t, y] = ode45(@(t,q)fc_ode_SIR(t,q,p),tspan,q0);
%% Plot
figure;
subplot(2,1,1)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]');
legend('Suceptible','Infected', 'Recovered', 'Location','Best');
y = y(:,3);%*N;
% figure(2);
subplot(2,1,2)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]');
legend('Infected', 'Location','Best');
And the function:
%% fc_ode_SIR
function dq = fc_ode_SIR(~,q,p) % q = initial conditions (Suceptible, Infected, Recovered)
% p = parameter (beta, mi, ni, initial conditions)
beta = p(1); gamma = p(2); delta = p(3);
S = q(1); I = q(2); R = q(3);
% SIR Model
dS = -beta*S*I + delta*R;
dI = beta*S*I - gamma*I;
dR = gamma*I - delta*R;
% dS = beta*S*I + delta*I; % Note changes here
% dI = beta*S*I - gamma*I - delta*I; % and here
% dR = gamma*I;
dq = [dS;
dI;
dR];
end
And the output:
Thanks all for your help!
Bjorn Gustavsson
Bjorn Gustavsson 2022-9-28
Doesn't seem unreasonable. Perhaps there's even some signs of what is necessary for recurring influenca outbreaks in the oscillations - however, the reality is obviously far more complicated.

请先登录，再进行评论。

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

回答（3 个）

Davide Masiello
Davide Masiello 2022-9-20
Ran in:
I guess the problem is in this three lines
S = beta*S*I + delta*R;
I = beta*S*I - gamma*I;
R = gamma*I - delta*I;
where you have redefined the variables, which seem to depend on the variables values themselves.
If I had to guess, those expressions are rather the values of the time derivatives of the variables.
Below, I tried to code it that way and it works, but please check the math carefully.
clear,clc
tspan = [0,100];
q0 = [99; 1; 0];
beta = 0.5; gamma = 20; delta = 0.4;
p = [beta; gamma; delta];
[t, y] = ode45(@(t,q)fc_ode_SIR(t,q,p),tspan,q0);
%% Plot
figure;
subplot(2,1,1)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]');
legend('Suceptible','Infected', 'Recovered', 'Location','Best');
y = y(:,3);%*N;
% figure(2);
subplot(2,1,2)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]');
legend('Infected', 'Location','Best');
function dq = fc_ode_SIR(~,q,p)
beta = p(1); gamma = p(2); delta = p(3);
S = q(1); I = q(2); R = q(3);
dS = beta*S*I + delta*R;
dI = beta*S*I - gamma*I;
dR = gamma*I - delta*I;
dq = [ -dS;
dI;
dR
];
end
  2 个评论
Julian Pereira
Julian Pereira 2022-9-20
Thank you so much! Oh it was a silly mistake indeed! Basically I was using the same variable for two diferente meanings.... Thank you!!!
Bjorn Gustavsson
Bjorn Gustavsson 2022-9-20
@Julian Pereira, no it is not. The loss of R due to transition back to S couldn't be proportional to the amount of I. This makes it possible, at least in principle that the loss of R (delta*I) becomes larger than R which would lead to a negatve number of R...

请先登录，再进行评论。

David Goodmanson
David Goodmanson 2022-9-20
编辑：David Goodmanson 2022-9-20
Hi Julian,
I have doubts about the equations, because with dq(1) = -S as you have, the equations are equivalent to
S = -beta*S*I - delta*R;
I = beta*S*I - gamma*I;
R = gamma*I - delta*I;
and it doesn't seem like recovered but reinfectible people should decrease the number of susceptible people. Should the upper right sign be positive? Also, should the lower right term be -delta*R instead of -delta*I?
Let's say whatever changes need to be done, or not, are done. The real problem is that the second line above redefines I, and then the new value, rather than the original value of I, is plugged into the third line. And the same thing is true for S redefined on line 1 and used on line 2, actually.
Bjorn Gustavsson
Bjorn Gustavsson 2022-9-20
You have an error in the SIR-equations. They should look something like this:
%% fc_ode_SIR
function dq = fc_ode_SIR(~,q,p) % q = initial conditions (Suceptible, Infected, Recovered)
% p = parameter (beta, gamma, delta, initial conditions)
beta = p(1); gamma = p(2); delta = p(3);
S = q(1); I = q(2); R = q(3);
% SIR Model
S = beta*S*I + delta*R;
I = beta*S*I - gamma*I;
R = gamma*I - delta*R; % NOTE THE CHANGE HERE!!!!!!
dq = [-S;
I;
R];
end
The source of susceptibles should either be some fraction of the recovered as I've changed it above, that is a loss of the recovered should be proportional to the number of recovered. Or some fraction of the infected should transfer directly to the susceptible:
%% fc_ode_SIR
function dq = fc_ode_SIR(~,q,p) % q = initial conditions (Suceptible, Infected, Recovered)
% p = parameter (beta, gamma, delta, initial conditions)
beta = p(1); gamma = p(2); delta = p(3);
S = q(1); I = q(2); R = q(3);
% SIR Model
S = beta*S*I + delta*I; % Note changes here
I = beta*S*I - gamma*I - delta*I; % and here
R = gamma*I;
dq = [-S;
I;
R];
end
That was your goof.
Also worth to not is that these simple SIR-models are decidedly un-biological. The constant transition-rates means that there will be a time-independent recovery-rate - if we at time zero turn off the infection-rate the reduction in infected will be a simple exponential decay - similar to radioactive decay. This is clearly not sensible.
HTH

类别

在 Help Center 和 File Exchange 中查找有关 Variables 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by