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

1 次查看（过去 30 天）
James Metz
James Metz 2020-11-26
评论： Alan Stevens 2020-11-26
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:

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

采纳的回答

Alan Stevens
Alan Stevens 2020-11-26
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 个评论
显示 2更早的评论
James Metz
James Metz 2020-11-26
Oh that makes sense. If you don't mind me asking, why is it necessary in a code?
Alan Stevens
Alan Stevens 2020-11-26
It wouldn't be necessary if your equations were linear, but terms like S*I make it so.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Biological and Health Sciences 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by