MATLAB Answers

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

1 view (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
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
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!

Translated by