MATLAB Answers

How can I store values in empty matrix to plot later?

2 views (last 30 days)
James Metz
James Metz on 26 Nov 2020
Commented: Alan Stevens on 26 Nov 2020
I am trying to simulate an epidemic model based on the SIR model. I am having trouble with my code. My graphs are not showing up like i think they should. It is not an equation error. I think there may be a problem with the storage of my values or something. Please help.
%Epidemic Simulation
%Author: James Metz
%Date: Nov 15, 2020
%Define Natural Paramteres:
p = 0.0144; %Natural Birth Rate (Davidson County)
u = 0.008; %Natural Death Rate (Davidson County)
f = 0.03; %Infection Rate
r = 0.075; %Recovery Rate
m = 0.0001; %Death due to Infection Rate
v = 0.092; %Vaccination Rate
%Define Evaluation Time Paramteres:
dt = 1; %Time increments (days)
tEnd = 365; %Simulation length (days)
t = 1:dt:tEnd;
%Initialize Variables:
I = 10; %Initial infected population
R = 0; %Initial recovered population
S = 692587; %Initial susceptible population
I_1d = zeros(1, tEnd);
R_1d = zeros(1, tEnd);
S_1d = zeros(1, tEnd);
S_1d(1) = 692587;
I_1d(1) = 10;
R_1d(1) = 0;
%Loop through times:
for idx = 1:dt:tEnd-1
%Initialization:
S = S_1d(idx);
R = R_1d(idx);
I = I_1d(idx);
%Find changes in population numbers
dS_dt = -f*S*I - S*u - S*v + S*p; %Drop in unifected population
dI_dt = f*S*I - r*I - m*I - u*I; %Drop in infected population
dR_dt = r*I + S*v - u*R; %Gain in Recovered population
%Store new values:
S_1d(idx+1) = S_1d(idx) + dS_dt;
R_1d(idx+1) = R_1d(idx) + dR_dt;
I_1d(idx+1) = I_1d(idx) + dI_dt;
end
figure(1)
subplot(2,2,1)
plot(t, S_1d)
xlabel('Time (days)')
ylabel('Susceptible Population')
title('Drop in Susceptible Population due to Infection')
subplot(2,2,2)
plot(t, R_1d)
xlabel('Time (days)')
ylabel('Recovered Population')
title('Recovery Rate of Infected Persons')
subplot(2,2,3)
plot(t, I_1d)
xlabel('Time (days)')
ylabel('Infected Population')
title('Infection Rate of Susceptible Persons')
subplot(2, 2, 4)
plot(t, S_1d, 'green')
plot(t, R_1d, 'Blue')
plot(t, I_1d, 'red')
This is how my graphs are turining out:

  2 Comments

KSSV
KSSV on 26 Nov 2020
Check the values......are you getting NaN's after some point of the loop.
James Metz
James Metz on 26 Nov 2020
Used isnan() function, but nothing yeilded a 'true'

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 26 Nov 2020
You need to wotk with normalized population numbers. You can renormalize at the end.
%Epidemic Simulation
%Author: James Metz
%Date: Nov 15, 2020
%Define Natural Paramteres:
p = 0.0144; %Natural Birth Rate (Davidson County)
u = 0.008; %Natural Death Rate (Davidson County)
f = 0.03; %Infection Rate
r = 0.075; %Recovery Rate
m = 0.0001; %Death due to Infection Rate
v = 0.092; %Vaccination Rate
%Define Evaluation Time Paramteres:
dt = 1; %Time increments (days)
%tEnd = 365; %Simulation length (days)
tEnd =365;
t = 0:dt:tEnd;
elems = numel(t);
%Initialize Variables:
I0 = 10; %Initial infected population
R0 = 0; %Initial recovered population
S0 = 692587; %Initial susceptible population
I_1d = zeros(1, elems);
R_1d = zeros(1, elems);
S_1d = zeros(1, elems);
S_1d(1) = 1; % Normalize the numbers
I_1d(1) = I0/S0;
R_1d(1) = R0/S0;
%Loop through times:
for idx = 1:elems-1
%Initialization:
S = S_1d(idx);
R = R_1d(idx);
I = I_1d(idx);
%Find changes in population numbers
dS_dt = -f*S*I - S*u - S*v + S*p; %Drop in unifected population
dI_dt = f*S*I - r*I - m*I - u*I; %Drop in infected population
dR_dt = r*I + S*v - u*R; %Gain in Recovered population
%Store new values:
S_1d(idx+1) = S + dS_dt*dt;
R_1d(idx+1) = R + dR_dt*dt;
I_1d(idx+1) = I + dI_dt*dt;
end
S_1d = S_1d*S0; % Renormalize the numbers
R_1d = R_1d*S0;
I_1d = I_1d*S0;
figure(1)
subplot(2,2,1)
plot(t, S_1d)
xlabel('Time (days)')
ylabel('Susceptible Population')
title('Drop in Susceptible Population due to Infection')
subplot(2,2,2)
plot(t, R_1d)
xlabel('Time (days)')
ylabel('Recovered Population')
title('Recovery Rate of Infected Persons')
subplot(2,2,3)
plot(t, I_1d)
xlabel('Time (days)')
ylabel('Infected Population')
title('Infection Rate of Susceptible Persons')
subplot(2, 2, 4)
plot(t, S_1d, 'green')
plot(t, R_1d, 'Blue')
plot(t, I_1d, 'red')
This results in

  4 Comments

Show 1 older comment
Alan Stevens
Alan Stevens on 26 Nov 2020
Normalizing the population numbers means dividing the values of S, I and R by the total number of people in the population. Then your equations work on fractions of the population. So a fraction of 1 means the whole population, a fraction of 0.9 means 90% of the population and so on.
James Metz
James Metz on 26 Nov 2020
Oh that makes sense. If you don't mind me asking, why is it necessary in a code?
Alan Stevens
Alan Stevens on 26 Nov 2020
It wouldn't be necessary if your equations were linear, but terms like S*I make it so.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!