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

6 vues (au cours des 30 derniers jours)
James Metz
James Metz le 26 Nov 2020
Commenté : Alan Stevens le 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 commentaires
KSSV
KSSV le 26 Nov 2020
Check the values......are you getting NaN's after some point of the loop.
James Metz
James Metz le 26 Nov 2020
Used isnan() function, but nothing yeilded a 'true'

Connectez-vous pour commenter.

Réponse acceptée

Alan Stevens
Alan Stevens le 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 commentaires
James Metz
James Metz le 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 le 26 Nov 2020
It wouldn't be necessary if your equations were linear, but terms like S*I make it so.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Problem-Based Optimization Setup dans Help Center et File Exchange

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by