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

조회 수: 5(최근 30일)
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
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일
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
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

Community Treasure Hunt

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

Start Hunting!

Translated by