Time-dependent immunity and drug resistant SIR model
2 views (last 30 days)
Show older comments
I have to solve this SIR model with Euler's forward. But my teacher said that it is wrong with my graph. Can anyone help me ?
Here is my code
%initial values
S(1)= 1;
I_1(1)= 0;
I_2 (1) = 0;
R(1)= 0;
tMax= 50;
dt = 0.01;
time_vector = 0:dt:tMax;
nIterations = length(time_vector);
t= 0:dt:tMax;
v= 0.0024;
tau = 0.6;
u= 0.8;
h= 0.5;
rho = 0.8;
R = 0.2;
r= re(R,v,t);
b = BetaL(h,tau);
n=0 ;
for i = time_vector
n = n+1;
S(n+1)= S(n)-dt*h*S(n)+dt*rho*I_1(n)+dt*b.*R(n);
I_1(n+1) =I_1(n)+ dt*(1-u)*S(n) - dt*rho*I_1(n) - dt*u*h*I_1(n);
I_2(n+1) =I_2(n)+ dt*u*h*S(n) + dt*u*h*I_1(n) - dt*r(n)*I_2(n);
R(n+1)= R(n) + dt*r(n)*I_2(n) - dt*b*R(n);
end
hold on
plot (t,S(1:end-1),'-r');
plot (t,I_1(1:end-1),'-b');
plot (t,I_2(1:end-1),'.b');
plot (t,R(1:end-1),'-g');
xlabel ('t');
ylabel ('population');
title ('Drug resistant model');
legend ('S','I_1','I_2','R');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/245222/image.png)
1 Comment
Aquatris
on 28 Oct 2019
What is BetaL function?
In the line with "S(n+1) = ..." you do not use the b(h)
In the line "I2(n+1) = ..." you are missing (1-exp(-vt^2).
In the line "R(n+1) = ..." you are missing (1-exp(-vt^2)
Accepted Answer
ME
on 28 Oct 2019
Edited: ME
on 28 Oct 2019
Assuming you have your BetaL and re functions setup correctly then there are the errors listed in the comment from Aquatris but you are also missing a multiple "h" in the I_1(n+1) equation.
Your equations at the bottom of your question also state that
=0.2 but you appear to have set R=0.2 instead and then defined a function r(n) instead of using
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/245236/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/245236/image.png)
If fixing those items doesn't help then I think you might need to give some more information as to what the BetaL and re functions are supposed to be.
More Answers (0)
See Also
Categories
Find more on 2-D and 3-D Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!