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...

1 answer below »







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) 1..20 21..70 71..84 85..90 91..110 111..1000 after 1000 R0 3.5 2.6 1.9 1.0 0.55 0.55 0.5 R0 3 2.2 0.7 0.8 1.00 0.90 0.5 R0 3 2.2 0.9 2.5 3.20 0.85 0.5 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; 91 0.55; 1001 0.5]; 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
d6e65cee71e43b426a5fee7af5ad6891f98b618e6872c06f7efb81e4f12392a6-us0qvygf.mlx a316comp06-vnwfsg35.pdf
Answered Same DayAug 05, 2021

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

Kshitij answered on Aug 06 2021
148 Votes
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...
SOLUTION.PDF
SOLUTION.PDF

Answer To This Question Is Available To Download

Related Questions & Answers

More Questions »

Submit New Assignment

Copy and Paste Your Assignment Here
Have files?