S is the number (in millions) of suceptible people, I is the number (in millions) of infectious people, and R is the number (in millions) of recovered and immune people. If we ignore births and deaths, the equations that describe the changes in these populations are:
We will assume that our total population is P0=330 million people. We start with 100 infected people, I(0)=1,000,000100=0.0001, S(0)=P0−I(0)=329.9999 million people and R(0)=0. (No one has recovered yet.)
Note that if we do not choose to track R the number of recovered people, that the change equations for S and I are defined only in terms of each other.
Part 1: Trajectories
Throughout this lab, examples will be given from the Romeo and Juliet model. We start by plotting a vector field:
var("J,R")Jprime(J,R)=R# formula for J'Rprime(J,R)=-J# formula for R'pvf=plot_vector_field([Jprime,Rprime],(J,-2,2),(R,-2,2),color="green",axes_labels=["J","R"],frame=False,aspect_ratio=1)# makes a vector field plotshow(pvf)# displays vector field plot
Exercise 1. For the SIR model, with a=1/6 and b=1/20, plot a vector field in the SI-plane with S values going from 0 to 350 and I values going from 0 to 200.
Exercise 2. Using the CoCalc functions point() and text() which were introduced in Lab 5, add the following labeled points to your vector field:
A point near the beginning of the disease where the number of susceptible people is high and the number of infectious people is low. Label this point "Outbreak Start".
Follow the arrows of the vector field and estimate the highest point of the trajectory starting from your first point. Add this high point to the vector field and label it "Max I". This represents the maximum number of people infected by the disease at a single time.
Following the trajectory you are imagining down to the S-axis and add a point near where you think it would intersect the S-axis. Label this point "Outbreak Ends". This point represents the end of the disease, when there are no more infected people.
From the above exercise you should have gotten an idea of how the number of susceptible and number of infected people change over the course of the disease (i.e. you have an idea of what the trajectory looks like). In the next few exercises we will plot a trajectory to see if it matches your intuition.
Exercise 3. Add the point I(0)=0.0001, S(0)=P0−I(0)=329.9999 to your vector field and label it "Start". These correspond to 100 initial infected people and just under 330 million susceptible people.
We are now going to create a trajectory starting from the point in the previous exercise. To do this we will use our Euler2D function from Lab 3:
defEuler2D(Xprime,Yprime,X0,Y0,t_end,dt):Xlist=[X0]# starts list for X valuesYlist=[Y0]# starts list for TvaluesX=X0# initial value for XY=Y0# initial value for Yfortinsrange(dt,t_end,dt,include_endpoint=true):# loops thru time values after tstartXp=Xprime(X,Y)# gives derivative X' at (X,Y)Yp=Yprime(X,Y)# gives derivative Y' at (X,Y)dX=dt*Xp# computes delta XdY=dt*Yp# computes delta YX=X+dX# finds new X valueY=Y+dY# finds new Y valueXlist.append(X)# apends new X value to XlistYlist.append(Y)# apends new Y value to YlistXYlist=[Xlist,Ylist]# Combines Xlist and Ylist to make a list of listsreturnXYlist#returns a list which is made up of the list of X values and the list of Y values
Euler2D allows us to solve systems of differential equations (such as our Romeo and Juliet model or our SIR model). Here is an example of how we can plot a trajectory using the output of Euler2D:
var("J,R")# make J and R symbolic variablesJprime(J,R)=R# formula for J'Rprime(J,R)=-J# formula for R'tend=7.dt=0.001JRlist=Euler2D(Jprime,Rprime,1.,0.,tend,dt)Jvals=JRlist# extracts the J values from JRlistRvals=JRlist# extracts the R values from JRlist# make a vector field plotpvf=plot_vector_field([Jprime,Rprime],(J,-2.,2.),(R,-2.,2.),color="green",axes_labels=["J","R"],frame=False,aspect_ratio=1)## make a trajectory plotptr=list_plot(list(zip(Jvals,Rvals)),plotjoined=True,axes_labels=["J","R"],aspect_ratio=1)## combine vector field and trajectory and showshow(pvf+ptr)
Exercise 4. Run the code for Euler2D above, then use it to plot a trajectory (for our SIR model) starting at the point from the previous exercise. Your time values should start at 0 days and go to 365 days in 1 day increments.
We can bring trajectories to life using animations. Here is an example with Romeo and Juliet:
var("J,R")Jprime(J,R)=R# formula for J'Rprime(J,R)=-J# formula for R'dt=0.001# time steptend=2.*piJRlist=Euler2D(Jprime,Rprime,1.,0.,tend,dt)tvals=srange(0.,tend,dt,include_endpoint=True)# list of time valuesJvals=JRlist# extracts the J values from JRlistRvals=JRlist# extracts the R values from JRlistpvf=plot_vector_field([Jprime(J,R),Rprime(J,R)],(J,-2.,2.),(R,-2.,2.),color="green",axes_labels=["J","R"],frame=False,aspect_ratio=1)# makes a vector field plot## trajectoryptr=list_plot(list(zip(Jvals,Rvals)),plotjoined=True,axes_labels=["J","R"],aspect_ratio=1)list_frames=# creates a list for frames to go intoforjinsrange(0,len(tvals),200):# loops going from 0 to the number of time values in tvals skipping by 200 each time.pt=point([Jvals[j],Rvals[j]],color="red",size=40)# plots the point at the j-th time valuep=pvf+pt+ptr# combine into one plotlist_frames.append(p)# appends this plot to list_framesani=animate(list_frames)# creates an animation from the list of plotsshow(ani,iterations=3)# displays the animation only repeat 3 times