Assignment 3
Due Wed, Mar 23, 23:59pm, commited to your repo on the Abacus GitLab server.
1.
The SIR model (W. O. Kermack and A. G. McKendrick, 1927, see Wikipedia page on this topic for references and additional detail) describes the outbreak and evolution of an infectious desease.
The model describes the evolution of a population with 𝑁 individuals, each of which belongs at any given time to one of three _compartments_:
𝑆(𝑡) is used to represent the number of individuals not yet infected with the disease at time t, and which are susceptible to the disease.
𝐼(𝑡) is the number of individuals that are infected and contagious.
𝑅(𝑡) represents individuals who have been infected and then removed from the disease, either due to immunization or due to death. Individuals in this category are not able to be infected again or to transmit the infection to others.
The total population 𝑁=𝑆(𝑡)+𝐼(𝑡)+𝑅(𝑡) is constant. In this model there is a general progression of individuals:
→→
The differential equations are interestingly quite similar to reaction networks in chemistry or nuclear physics, where a member of compartment has to interact with a member of compartment to create a new member of . This is equivalent to a reaction. The transition from to is equivalent to a radioactive decay.
The system of differential equations is:
𝑑𝑆𝑑𝑡=−𝛽𝑆𝐼𝑁
𝑑𝐼𝑑𝑡=𝛽𝑆𝐼𝑁−𝛾𝐼
𝑑𝑅𝑑𝑡=𝛾𝐼
𝛾 represents the mean recovery/death rate in units of d−1, and 1/𝛾 is the mean infective period 𝑃inf=1/𝛾. The infective period is the asymptomatic incubation period plus the symptomatic time until the individual transitions to compartment .
𝛽 represents the infection rate. If 𝑛inf is the number of members of each member of infects, then 𝛽=𝑛inf/𝑃inf=𝑛inf𝛾.
The final input parameters are the population size 𝑁 and the initial values 𝑆0, 𝐼0 and 𝑅0.
1.1
Solve the model equations for 𝛾=1/10d−1 which implies an infective period of 𝑃inf=1/𝛾=10d. Use as initial values (𝑆0,𝐼0,𝑅0)=(𝑁−𝐼0,574,0) and 𝑁=1.1×107. These are the parameters that describe approximately the situation of the Corona virus infection in China, focussing on the city of Wuhan, on January 22 2020:
Cvirus
This data could at the time be viewed and exlored at the WHO web page.
By Feb 4 (13 days later) the number of infected individuals had increased to 24.4×103 individuals.
The goal is to estimate 𝑛inf. Start with 𝑛inf=2 implying that each infected person will infect on average 2 uninfected individuals. Make a simulation for the duration of 13d and repeat while adjusting the 𝑛inf parameter (to within one decimal) until it approximately matches the Wuhan data, i.e. after 13 days the compartment has risen from its initial value to approximately 25k. What is the best value of 𝑛inf?
Save the plot for this part to the file Fig1.1.png.
Note: The parameter 𝑛inf is now widely known as the reproductive number 𝑅0, not to be confused with the variable 𝑅 we are using here for the removed population.
%pylab ipympl
1.2
Using these parameters make a simulation for 100𝑑, and plot the three compartments 𝑆, 𝐼 and 𝑅. We are also interested in the rate of hospital admissions per day. These are 5% of the rate of those who transition from 𝐼 to 𝑅, i.e. 5% of 𝑑𝑅𝑑𝑡. Add as well a line for the 1% fraction of those in 𝑅 who die. Finally, add a horizontal line for the capacity of hospital beds available in terms of admission capacity per day. China has on average 4.34 beds per 1000 people. Assume that each patient stays 4 days in hospital, and that 20% of the available hospital beds can be used for Corona virus patients.
1.2.1
After how many days from the beginning of the simulations on Jan 22nd does the hospital admission per day according to the simulation with the 𝑛inf parameter from 1.1 exceed the available capacity? How many people would have died at this point? (Find these numbers approximately by interacting with your graph zooming in.) Further, in this scenario, by how much would the daily admission exceed capacity when it peaks, and what is the asymptotic number of people who would have died?
1.2.2
Obviously this is not the scenario that unfolded, because 𝑛inf was dramatically reduced through physical distancing. Esimate how large 𝑛inf may at most be so that the hospital capacity is not overwhelmed.
You can see that this scenario is also not what unfolded, at least not in most cases. We now understand why this is also not what happened, and why in most places 𝑛inf was further reduced and kept lower.
Problem 2
𝑓(𝑥)=sin(𝑥0.9)
and
𝑔(𝑥)=cos(2𝑥+𝑥⎯⎯√)
Plot 𝑔(𝑥) vs. 𝑓(𝑥) for 𝑥∈[0,…,50] using 1000 equal intervals.
Make a snake scatter plot of the first 40 points (40 is then the snake length_) of the 𝑥 array from part 2.1. Make the size and color proportional to 𝑥 within the plotted range (snake length) so that the _head of the snake is on color and large and the tail is small and the other color.
Create a multi-processing script and add it to the assignment commit, that creates 960 frames, each of which shows a snake scatter as in 2.2, but each with a different, subsequent starting point in the 𝑥 array from its 1st to 960th value. The multi-processing script should use 6 threads. Each frame is written to a png image file with the number of the x value in the name. Finally, use the ffmpeg program to combine the frames into a movie called trig_snake.mp4. Use -framerate 30 to generate a swift movie experience.
Problem 3
A network of species x1, x2 and x3 is defined by the following transmutations:
x1 + x2 -> 2 x3
x3 + x1 -> 2 x2
x3 + x2 -> 2 x1
For example, one particle of type x1 reacts with one particle of type x2 to make 2 particles of type x3. The transmutations are preserving the total number of species. The rate at which each of these transmutations take place is r12 for transmutation x1 + x2 -> 2 x3 and so on. Then, the following coupled system of ODEs describe the time evolution of each species. Note that negative trems on the RHS correspond to destruction terms and positive terms are oroduction terms.
𝑑𝑥1𝑑𝑡=−𝑟12𝑥1𝑥2−𝑟31𝑥3𝑥1+2𝑟32𝑥3𝑥2𝑑𝑥2𝑑𝑡=−𝑟12𝑥1𝑥2+2𝑟31𝑥3𝑥1−𝑟32𝑥3𝑥2𝑑𝑥3𝑑𝑡=+2𝑟12𝑥1𝑥2−𝑟31𝑥3𝑥1−𝑟32𝑥3𝑥2
3.1
Integrate this system for the following parameters:
# parameters, initial conditions, integration range
r12 = 0.15; r31 = 0.05; r32 = 0.25
state0 = [0.4, 1.5, 0.01]
tmin,tmax = (0,25)
Make a line plot of the time evolution of each of the three species as well as their sum (to check total number conservation).
3.2
Find the equilibrium solutions, which are obtained by setting the LHS of the above given ODE to zero. You will have to find a SciPy library that can find the roots of multivariable functions. In most cases for any root finding algorithm to work it is key to have a good initial guess for the solution. Take your initial guess for x1, x2 and x3 from the ODE solution from 3.1 for t=10.
Add markers for the equilibrium solution to the plot made in 3.1, at t=25.
Note: It is instructive to use this problem to explore different methods offered by the root finding package, the dependence of the found solution on the initial guess and on tolerance.
Problem 4
Create a python function func1 that calculates
𝑓(𝑥)=𝑠𝑖𝑛(𝑥3)𝑥3𝑒𝑥𝑝(−𝑥)
Make a plot for 𝑥∈[−2.5,1] using an x array representing 100 intervals.
Calculate the derivative 𝑑𝑓𝑑𝑥 numerically and make a plot of the derivative.
Calculate the integral of 𝑓(𝑥) for 𝑥∈[−2.5,1] using
a library that uses a conventional discrete integration
a Monte-Carlo method
Solve the ODE
𝑑𝑥𝑑𝑡=𝑓(𝑥)
for the initial value 𝑥01=−1.465, and a second time for 𝑥02=−1.463 for 𝑡∈[0,20]. Make a plot of 𝑥(𝑡) for both initial value and explain briefly why the two trajectories is so different althoug the initial values are almost the same.