Macm 316 Computing Assignment #6
Submit on Crowdmark by Wednesday, August 5, 2020, 11pm
Upload one...

Question

Macm 316 Computing Assignment #6
Submit on Crowdmark by Wednesday, August 5, 2020, 11pm
Upload one .pdf file with up to 4 pages: Pages 1-2 is your typed report (your discussions, data and
figures on two pages); Pages 3-4 is a listing of your code(s). The assignment is due at 11:00pm.
You will receive a Crowdmark link for the upload.
The deadline for this assignment is a Wednesday. I moved it to take into account the long weekend
Aug 1–3. You are free to upload sooner!
A (simplified) Covid-19 model. The so-called SEIR model is a standard way to model the
spread of infectious diseases. It is an example of a compartmental model. We will look at this
model in its most basic form, and solve it numerically for various parameters.
The population of N individuals is assigned to different compartments, and individuals can move
between compartments. The SEIR model uses four compartments, S,E, I,R, representing segments
of the population, so S+E+I+R = N (a simplification, meaning the total population is constant,
with some leeway in interpreting the term "total population").
S. Susceptible individuals. Those are healthy, and not immune individuals who may become
infected upon contact with an infected individual.
E. Exposed individuals. These are individuals who have been infected, but because of the
incubation period of the virus are not yet infectious themselves. They will transition to the
infected group;
I. Infected individuals. These are individuals who have been infected, and can pass on the
infection to susceptible individuals.
R. Removed individuals. (In an optimistic scenario referred to as recovered individuals.) These
include people who have recovered and are now immune, and people who have died.
Note that for Covid-19 we do not know yet whether recovered people are immune, and if so, for
how long. This model assumes they become immune. Interactions between people in compartments
are proportional to the number of people in each compartment – this gives a “quadratic” term.
Transfers from one compartment to another that do not involve interactions are assumed to occur
proportional to the number of people in one compartment, corresponding to linear terms.
The equations are as follows (the variable t is time, say measured in days):
dS
dt
= −βS I
N
susceptibles become exposed due to interaction with I, contact rate of β
dE
dt
= +βS I
N
− αE incoming S minus exposed moving to I, α = 1/incubation period
dI
dt
= αE − γI incoming exposed - “removed”, γ = 1/infectious period
dR
dt
= +γI
mrt 1
Macm 316 Computing Assignment #6
That was the “easy” part. The difficulty, even within the constraint of this model, is to get somewhat
realistic parameters – a major data collection problem in the real world. The other interesting part
is, that some of these parameters will also be time dependent – the effect of physical distancing
and lockdown or reopening. You probably have heard about the now infamous parameter R0,
the reproduction number. For our model we have R0 = β/γ. (Given that we have “removed”
individuals R, the R0 naming convention is somewhat unfortunate.)
A reasonable set of parameters is
1. α = 1/5.2, incubation period of roughly 5 days.
2. γ = 1/10, infectious period of 10 days. The actual period is probably longer, but this value
takes into account that sick people go to hospital or stay home, rather than continuing to
circulate among the general population).
3. R0 = 3.5, for now assumed constant (which would be bad news).
Note that the parameter β = R0γ.
When things become interesting is when we work with time-varying parameters, in particular
R0 = R0(t) (and, by association, β = β(t)). When you look at the model, β is the contact rate,
which can be controlled. If we all stayed in a remote log cabin by ourselves, β = R0 = 0; if
we practice physical distancing and wear masks when physical distancing is not possible, then R0
might stay below 1; if we hang out in bars or at large gatherings, then β and R0 will go through
the roof – leading to large numbers of Is, infected and infectious people.
We could model the effect of intervention with data for R0 as follows (three profiles are given):
Days (since outbreak XXXXXXXXXX XXXXXXXXXX XXXXXXXXXXafter 1000
R XXXXXXXXXX XXXXXXXXXX
R XXXXXXXXXX XXXXXXXXXX
R XXXXXXXXXX XXXXXXXXXX
What is missing? We need N , the number of people. For British Columbia we use N = 5 Million.
We also need initial conditions with I(0) and/or E(0) non-zero. So, “to get started” let us take
I(0) = 40; E(0) = 20 I(0); R(0) = 0, S(0) = N − I(0) − E(0) −R(0).
Once again, R(0) is the initial value of R(t), and is nor related to the reproductive number R0.
Last point: we run the simulation for 6 months, so take Tfinal = 180. Sample output:
mrt 2
Macm 316 Computing Assignment #6
Now you have all the information to become a member of the modeling team for the infectious
disease epidemic.
Your tasks
1. Use one of Matlab’s codes, for example, ode45 or ode15s, to solve the SEIR equations
numerically. This will require you to write a function that evaluates the right hand side of
the differential equations.
2. Run your codes with N = 5 Million people (BC scenario), and α = 1/5.2 and γ = 1/10.
3. At first, use constant R0 and run your code with different values of R0: 3.5, 2.5, 1.25, 0.9. For
the constant R0 case summarize your observations, but only report the results corresponding
to R0 = 3.5. R0 = 3.5 is your Scenario 1.
4. Your Scenarios 2, 3, and 4: Run your code with the time dependent R0 values given in the
table above. You may simply use the values as a piecewise constant function, or interpolate
the values using Matlab’s interp1 function. The first row corresponds to a scenario with
some lock-down measures; the second row shows an effective response; the third row shows
the consequences of allowing a lapse in prevention measures. You may want to plot your R0
values agains time.
5. Assuming a 4% death rate among the R(t) population, compare the number of deaths per 1
Million people for each of your four scenarios. A table will do!
6. Show your plots of number of active cases and total number of cases for your four scenarios.
7. BC has a total of about 5600 acute care beds, and 200 ICU (intensive care unit) beds.
Assume that 3500 of acute care beds, and 160 ICU beds are available for Covid cases, and
8% of infected people need to be hospitalized, and 1% of infected people must be admitted
to the ICU. Compare the four scenarios above in terms of hospital capacity: Plot the number
of Covid patients in hospital and in ICUs, and compare to the capacity,
Some programming notes – just some ideas, you are allowed to do things differently.
1. Rewrite in vector form. To get the system into a format for Matlab to process, we write
y(t) =

y1(t)
y2(t)
y3(t)
y4(t)
 =

S(t)
E(t)
I(t)
R(t)
 , dydt =

−βy1(t)
y3(t)
N
+βy1(t)
y3(t)
N
− αy2(t)
αy2(t) − γy3(t)
+γy3(t)

2. Assume your function describing the equation is
function yprime = seir(t,y)
mrt 3
Macm 316 Computing Assignment #6
The easiest way to get the parameters (α, γ and R0) into the function is to “hard
code” them. You then have to change your function every time you run your code with
different parameters.
Alternatively, you can pass parameters in a variable called par.
%Main program
par.alpha = 1/5.2;
par.gamma = 1/10;
par.rzero = 3.5;
par.N = 5.0e6;
function yprime = seir(t,y,par)
You can then refer to those values inside your function seir.m, as they are being passed
along. To call, for example, ode45 you would write
rtol = 1.e-6; atol=1.e-5;
options = odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] = ode45(@(t, y) seir(t,y, par) , [0,180], yinit, options)
atol and rtol are error tolerances for solving the differential equations.
yinit is the vector of initial values [S(0), E(0), I(0), R(0)]T . If you want Matlab to
give you more information, you could use the following options:
options = odeset('AbsTol', atol,'RelTol', rtol,'MaxOrder',5,'Stats','on');
3. A remaining issue is how to implement a time-varying R0 factor. One possibility is
to pass an array of length at least 180, that has an entry for each day. For the first
time-varying data set:
maxd=200; rtzero=zeros(1,maxd);
rtable = [ 1 3.5; 21 2.6; 71 1.9; 85 1.0; XXXXXXXXXX; XXXXXXXXXX];
for j=1:5,
constdays = min(rtable(j+1,1),maxd+1) - rtable(j,1);
rtzero( rtable(j,1):rtable(j,1)+constdays-1 ) = rtable(j,2)*ones(1,constdays);
end;
par.rtzero = rtzero;
4. You can test your code by setting γ = 0. Then you are essentially just solving an
exponential decay equation.
mrt 4

Kshitij · Accepted Answer

renewode/Question2Solution.m
% At first, use constant R0 and run your code with different values of R0: 3.5, 2.5, 1.25, 0.9. For
% the constant R0 case summarize your observations, but only report the results corresponding
% to R0 = 3.5. R0 = 3.5 is your Scenario 1.
% Test the function 
% Run your codes with N=5 Million people 
% (BC scenario), and ? = 1/5.2 and
% gamma=1/10. 
N=5e6; 
a=1/5.2;
gamma=1/10;
R0=[3.5, 2.5, 1.25, 0.9];
for ii=1:length(R0) 
b=gamma*R0(ii); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(ii);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(ii)]'  % Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[0,180],yinit,options)
figure(),
subplot(4,1,1),plot(t,y(:,1),'b','linewidth',2); grid on
xlabel('Time'); ylabel('S(t)');
str=['The initial condtion of R0=', num2str(R0(ii))];
title(str)
subplot(4,1,2),plot(t,y(:,2),'b','linewidth',2); grid on
xlabel('Time'); ylabel('E(t)');
subplot(4,1,3),plot(t,y(:,3),'b','linewidth',2); grid on
xlabel('Time'); ylabel('I(t)');
subplot(4,1,4),plot(t,y(:,4),'b','linewidth',2); grid on
xlabel('Time'); ylabel('R(t)');
end
renewode/Question4Solution.m
% Your Scenarios 2, 3, and 4: Run your code with the time dependent R0 values given in the
% table above. You may simply use the values as a piecewise constant function, or interpolate
% the values using Matlab’s interp1 function. The first row corresponds to a scenario with
% some lock-down measures; the second row shows an effective response; the third row shows
% the consequences of allowing a lapse in prevention measures. You may want to plot your R0
% values agains time.
%Scenario 4
clear all
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3.5,2.6,1.9,1,0.55,0.55];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
figure(i)
subplot(4,1,1),plot(t,y(:,1),'b','linewidth',2); grid on
xlabel('Time'); ylabel('S(t)');
str=['The initial condtion of R0=', num2str(R0(i))];
title(str)
hold off
hold on
subplot(4,1,2),plot(t,y(:,2),'b','linewidth',2); grid on
xlabel('Time'); ylabel('E(t)');
hold off
hold on
subplot(4,1,3),plot(t,y(:,3),'b','linewidth',2); grid on
xlabel('Time'); ylabel('I(t)');
hold off
subplot(4,1,4),plot(t,y(:,4),'b','linewidth',2); grid on
xlabel('Time'); ylabel('R(t)');
hold off
end
%% Scenario 3
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3,2.2,0.7,0.8,1,0.9];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
figure(i)
subplot(4,1,1),plot(t,y(:,1),'b','linewidth',2); grid on
xlabel('Time'); ylabel('S(t)');
str=['The initial condtion of R0=', num2str(R0(i))];
title(str)
hold off
hold on
subplot(4,1,2),plot(t,y(:,2),'b','linewidth',2); grid on
xlabel('Time'); ylabel('E(t)');
hold off
hold on
subplot(4,1,3),plot(t,y(:,3),'b','linewidth',2); grid on
xlabel('Time'); ylabel('I(t)');
hold off
subplot(4,1,4),plot(t,y(:,4),'b','linewidth',2); grid on
xlabel('Time'); ylabel('R(t)');
hold off
end
%% Scenario 4
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3,2.2,0.9,2.5,3.2,0.85];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
figure(i)
subplot(4,1,1),plot(t,y(:,1),'b','linewidth',2); grid on
xlabel('Time'); ylabel('S(t)');
str=['The initial condtion of R0=', num2str(R0(i))];
title(str)
hold off
hold on
subplot(4,1,2),plot(t,y(:,2),'b','linewidth',2); grid on
xlabel('Time'); ylabel('E(t)');
hold off
hold on
subplot(4,1,3),plot(t,y(:,3),'b','linewidth',2); grid on
xlabel('Time'); ylabel('I(t)');
hold off
subplot(4,1,4),plot(t,y(:,4),'b','linewidth',2); grid on
xlabel('Time'); ylabel('R(t)');
hold off
end
renewode/question_5.asv
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3.5,2.6,1.9,1,0.55,0.55];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
end
% Death amaong R(t) is 4%
Total_rt1=max(y(:,4));
death1=.04*Total_rt1;
death_per_mil2=death/5;
%%
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3,2.2,0.7,0.8,1,0.9];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
end
% Death amaong R(t) is 4%
Total_rt2=max(y(:,4));
death2=.04*Total_rt2;
death_per_mil3=death/5;
%%
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3,2.2,0.9,2.5,3.2,0.85];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
end
% Death amaong R(t) is 4%
Total_rt3=max(y(:,4));
death=.04*Total_rt3;
death_per_mil4=death/5;
%%
renewode/question_5.m
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3.5,2.6,1.9,1,0.55,0.55];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
end
% Death amaong R(t) is 4%
Total_rt1=max(y(:,4));
death1=.04*Total_rt1;
death_per_mil2=death1/5;
Death_scenerio2=death_per_mil2;
%%
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3,2.2,0.7,0.8,1,0.9];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
end
% Death amaong R(t) is 4%
Total_rt2=max(y(:,4));
death2=.04*Total_rt2;
death_per_mil3=death2/5;
Death_scenerio3=death_per_mil3;
%%
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3,2.2,0.9,2.5,3.2,0.85];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
end
% Death amaong R(t) is 4%
Total_rt3=max(y(:,4));
death3=.04*Total_rt3;
death_per_mil4=death3/5;
Death_scenerio4=death_per_mil4;
%%
T = table(Death_scenerio2,Death_scenerio3,Death_scenerio4)
renewode/seir.m
function yprime=seir(t,y,a,b,gamma,N)
yprime=zeros(4,1);
% The following are the differential equations 
yprime(1)=-b*y(1)*y(3)/N;
yprime(2)=b*y(1)*y(3)/N-a*y(2);
yprime(3)=a*y(2)-gamma*y(3); 
yprime(4)=gamma*y(3);
end
renewode/sixth.m
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3.5,2.6,1.9,1,0.55,0.55];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
total_infection2= y(:,3)+y(:,4);
total2{i}=total_infection2;
total_time2{i}=t;
end
%%
Total_time2 = cat(1, total_time2{:})
Total_active_case2 = cat(1, total2{:})
figure;
plot(Total_time2,Total_active_case2)
grid on
xlabel('Time in days'); ylabel('Total active cases');
title('scenerio 2')
%%
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3,2.2,0.7,0.8,1,0.9];
total=[];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
total_infection= y(:,3)+y(:,4);
total{i}=total_infection;
total_time{i}=t;
end
%%
Total_time = cat(1, total_time{:})
Total_active_case = cat(1, total{:})
figure;
plot(Total_time,Total_active_case)
grid on
xlabel('Time in days'); ylabel('Total active cases');
title('scenerio 3')
%%
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3,2.2,0.9,2.5,3.2,0.85];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]' ;% Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options);
total_infection1= y(:,3)+y(:,4);
total1{i}=total_infection1;
total_time1{i}=t;
end
%%
Total_time1 = cat(1, total_time1{:})
Total_active_case1 = cat(1, total1{:})
figure;
plot(Total_time1,Total_active_case1)
grid on
xlabel('Time in days'); ylabel('Total active cases');
title('scenerio 4')
renewode/TestQ1.m
% Test the function 
% Run your codes with N=5 Million people 
% (BC scenario), and ? = 1/5.2 and
% gamma=1/10. 
N=5e6; 
a=1/5.2;
gamma=1/10;
R0=0;   % We chose R0=0 at first 
b=gamma*R0; 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0;
% The simulation of question 1: 
yinit=[S0,E0,I0,R0]'  % Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[0,180],yinit,options)
figure(),
subplot(4,1,1),plot(t,y(:,1),'b','linewidth',2); grid on
xlabel('Time'); ylabel('S(t)');
subplot(4,1,2),plot(t,y(:,2),'b','linewidth',2); grid on
xlabel('Time'); ylabel('E(t)');
subplot(4,1,3),plot(t,y(:,3),'b','linewidth',2); grid on
xlabel('Time'); ylabel('I(t)');
subplot(4,1,4),plot(t,y(:,4),'b','linewidth',2); grid on
xlabel('Time'); ylabel('R(t)');
renewode/Untitled.mlx
[Content_Types].xml
        
_rels/.rels
       
matlab/document.xml
       % Test the function 
% Run your codes with N=5 Million people 
% (BC scenario), and ? = 1/5.2 and
% gamma=1/10. 
N=5e6; 
a=1/5.2;
gamma=1/10;
R0=0;   % We chose R0=0 at first 
b=gamma*R0; 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0;
% The simulation of question 1: 
yinit=[S0,E0,I0,R0]'  % Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[0,180],yinit,options)
figure(),
subplot(4,1,1),plot(t,y(:,1),'b','linewidth',2); grid on
xlabel('Time'); ylabel('S(t)');
subplot(4,1,2),plot(t,y(:,2),'b','linewidth',2); grid on
xlabel('Time'); ylabel('E(t)');
subplot(4,1,3),plot(t,y(:,3),'b','linewidth',2); grid on
xlabel('Time'); ylabel('I(t)');
subplot(4,1,4),plot(t,y(:,4),'b','linewidth',2); grid on
xlabel('Time'); ylabel('R(t)');        % At first, use constant R0 and run your code with different values of R0: 3.5, 2.5, 1.25, 0.9. For
% the constant R0 case summarize your observations, but only report the results corresponding
% to R0 = 3.5. R0 = 3.5 is your Scenario 1.
% Test the function 
% Run your codes with N=5 Million people 
% (BC scenario), and ? = 1/5.2 and
% gamma=1/10. 
N=5e6; 
a=1/5.2;
gamma=1/10;
R0=[3.5, 2.5, 1.25, 0.9];
for ii=1:length(R0) 
b=gamma*R0(ii); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(ii);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(ii)]'  % Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[0,180],yinit,options)
figure(),
subplot(4,1,1),plot(t,y(:,1),'b','linewidth',2); grid on
xlabel('Time'); ylabel('S(t)');
str=['The initial condtion of R0=', num2str(R0(ii))];
title(str)
subplot(4,1,2),plot(t,y(:,2),'b','linewidth',2); grid on
xlabel('Time'); ylabel('E(t)');
subplot(4,1,3),plot(t,y(:,3),'b','linewidth',2); grid on
xlabel('Time'); ylabel('I(t)');
subplot(4,1,4),plot(t,y(:,4),'b','linewidth',2); grid on
xlabel('Time'); ylabel('R(t)');
end
        % Your Scenarios 2, 3, and 4: Run your code with the time dependent R0 values given in the
% table above. You may simply use the values as a piecewise constant function, or interpolate
% the values using Matlab’s interp1 function. The first row corresponds to a scenario with
% some lock-down measures; the second row shows an effective response; the third row shows
% the consequences of allowing a lapse in prevention measures. You may want to plot your R0
% values agains time.
clear all
N=5e6; 
a=1/5.2;
gamma=1/10;
T=[0 21 71 85 91 200];
R0=[3.5,2.6,1.9,1,0.55];
for i=1:5
b=gamma*R0(i); 
I0=40; 
E0=20*I0;
S0=N-I0-E0-R0(i);
% The simulation of question 1: 
yinit=[S0,E0,I0,R0(i)]'  % Initial Conditions
rtol =1.e-6; atol=1.e-5;
options=odeset('AbsTol', atol,'RelTol',rtol,'MaxOrder',5);
[t,y] =ode45(@(t,y)seir(t,y,a,b,gamma,N),[T(i),T(i+1)],yinit,options)
hold on
subplot(4,1,1),plot(t,y(:,1),'b','linewidth',2); grid on
xlabel('Time'); ylabel('S(t)');
str=['The initial condtion of R0=', num2str(R0(i))];
title(str)
hold off
hold on
subplot(4,1,2),plot(t,y(:,2),'b','linewidth',2); grid on
xlabel('Time'); ylabel('E(t)');
hold off
hold on
subplot(4,1,3),plot(t,y(:,3),'b','linewidth',2); grid on
xlabel('Time'); ylabel('I(t)');
hold off
subplot(4,1,4),plot(t,y(:,4),'b','linewidth',2); grid on
xlabel('Time'); ylabel('R(t)');
hold off
end        function yprime=seir(t,y,a,b,gamma,N)
yprime=zeros(4,1);
% The following are the differential equations 
yprime(1)=-b*y(1)*y(3)/N;
yprime(2)=b*y(1)*y(3)/N-a*y(2);
yprime(3)=a*y(2)-gamma*y(3); 
yprime(4)=gamma*y(3);
end
matlab/output.xml
   manual code ready  0.4   matrix   yinit 4×1 4 1 double      4999160
         800
          40
           0
 double double [[{"value":"{\"value\":\"4999160\"}"},{"value":"{\"value\":\"800\"}"},{"value":"{\"value\":\"40\"}"},{"value":"{\"value\":\"0\"}"}]]  16  matrix   t 285×1 285 1 double          0
    0.0034
    0.0067
    0.0101
    0.0135
    0.0303
    0.0472
    0.0640
    0.0808
    0.1651
 double double [[{"value":"{\"value\":\"0\"}"},{"value":"{\"value\":\"0.0034\"}"},{"value":"{\"value\":\"0.0067\"}"},{"value":"{\"value\":\"0.0101\"}"},{"value":"{\"value\":\"0.0135\"}"},{"value":"{\"value\":\"0.0303\"}"},{"value":"{\"value\":\"0.0472\"}"},{"value":"{\"value\":\"0.0640\"}"},{"value":"{\"value\":\"0.0808\"}"},{"value":"{\"value\":\"0.1651\"}"}]]  19  matrix   y 285×4 285 4 double 	1.0e+06 *
    4.9992    0.0008    0.0000         0
    4.9992    0.0008    0.0000    0.0000
    4.9992    0.0008    0.0000    0.0000
    4.9992    0.0008    0.0000    0.0000
    4.9992    0.0008    0.0000    0.0000
    4.9992    0.0008    0.0000    0.0000
    4.9992    0.0008    0.0000    0.0000
    4.9992    0.0008    0.0000    0.0000
    4.9992    0.0008    0.0001    0.0000
    4.9992    0.0008    0.0001    0.0000
 double 6 double [[{"value":"{\"value\":\"4.9992\"}"},{"value":"{\"value\":\"0.0008\"}"},{"value":"{\"value\":\"0.0000\"}"},{"value":"{\"value\":\"0\"}"}]]  19  figure  ada9e41e-efcc-4eaa-a062-feeb60a7eb6f data:image/png;base64,

Macm 316 Computing Assignment #6 Submit on Crowdmark by Wednesday, August 5, 2020, 11pm Upload one .pdf file with up to 4 pages: Pages 1-2 is your typed report (your discussions, data and figures on...

Solution

Answer To This Question Is Available To Download

Related Questions & Answers

Submit New Assignment