$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: $S'=-a SI/P_0\quad I'=a SI/P_0 - b I\qquad R'=b I$ We will assume that our total population is $P_0=330$ million people. We start with 100 infected people, $I(0)=\frac{100}{1,000,000}=0.0001$, $S(0)=P_0-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.
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 plot show(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$.
var("S,I") a=1/6 b=1/20 P0=330 Sprime(S,I)=(-a*S*I)/P0 Iprime(S,I)=(a*S*I)/P0-b*I vec1=plot_vector_field([Sprime,Iprime], (S,0,350), (I,0,200), color="green", axes_labels=["S","I"], frame=False, aspect_ratio=1) show(vec1)
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.
var("S,I") a=1/6 b=1/20 P0=330 Sprime(S,I)=(-a*S*I)/P0 Iprime(S,I)=(a*S*I)/P0-b*I vec1=plot_vector_field([Sprime,Iprime], (S,0,350), (I,0,200), color="green", axes_labels=["S","I"], frame=False, aspect_ratio=1)+point([300,25], color="blue", size=40)+text("Outbreak Start", [300,25-10])+point([110,130], color="blue", size=40)+text("Max I", [110, 130-10])+point([28,0], color="blue", size=40)+text("Outbreak End", [32, 10]) show(vec1)
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)=P_0-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.
var("S,I") a=1/6 b=1/20 P0=330 Sprime(S,I)=(-a*S*I)/P0 Iprime(S,I)=(a*S*I)/P0-b*I vec1=plot_vector_field([Sprime,Iprime], (S,0,350), (I,0,200), color="green", axes_labels=["S","I"], frame=False, aspect_ratio=1)+point([300,25], color="blue", size=40)+text("Outbreak Start", [300,25+10])+point([110,130], color="blue", size=40)+text("Max I", [110, 130-10])+point([28,0], color="blue", size=40)+text("Outbreak End", [32, 10])+point([329.999,.0001], color="blue", size=40)+text("Start", [329.999, .0001+10]) show(vec1)
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:
def Euler2D(Xprime,Yprime,X0,Y0,t_end,dt): Xlist=[X0] # starts list for X values Ylist=[Y0] # starts list for Tvalues X=X0 # initial value for X Y=Y0 # initial value for Y for t in srange(dt,t_end,dt,include_endpoint=true): # loops thru time values after tstart Xp=Xprime(X,Y) # gives derivative X' at (X,Y) Yp=Yprime(X,Y) # gives derivative Y' at (X,Y) dX=dt*Xp # computes delta X dY=dt*Yp # computes delta Y X=X+dX # finds new X value Y=Y+dY # finds new Y value Xlist.append(X) # apends new X value to Xlist Ylist.append(Y) # apends new Y value to Ylist XYlist=[Xlist,Ylist]# Combines Xlist and Ylist to make a list of lists return XYlist #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 variables Jprime(J,R)=R # formula for J' Rprime(J,R)=-J # formula for R' tend=7. dt=0.001 JRlist=Euler2D(Jprime,Rprime,1.,0.,tend,dt) Jvals=JRlist[0] # extracts the J values from JRlist Rvals=JRlist[1] # extracts the R values from JRlist # make a vector field plot pvf=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 plot ptr=list_plot(list(zip(Jvals,Rvals)),plotjoined=True,axes_labels=["J","R"],aspect_ratio=1) ## combine vector field and trajectory and show show(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.
var("S,I") a=1/6 b=1/20 P0=330 I0=.0001 S0=330-I0 Sprime(S,I)=(-a*S*I)/P0 Iprime(S,I)=(a*S*I)/P0-b*I tend=365. dt=1. SIlist=Euler2D(Sprime,Iprime,330,.0001,tend,dt) Svals=SIlist[0] Ivals=SIlist[1] pvf1=plot_vector_field([Sprime,Iprime],(S,0,350),(I,0,200), color="green", axes_labels=["S","I"], frame=False, aspect_ratio=1) ptr1=list_plot(list(zip(Svals,Ivals)), plotjoined=True, axes_labels=["S","I"], aspect_ratio=1) show(pvf1+ptr1)
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 step tend=2.*pi JRlist=Euler2D(Jprime,Rprime,1.,0.,tend,dt) tvals=srange(0.,tend,dt,include_endpoint=True) # list of time values Jvals=JRlist[0] # extracts the J values from JRlist Rvals=JRlist[1] # extracts the R values from JRlist pvf=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 ## trajectory ptr=list_plot(list(zip(Jvals,Rvals)),plotjoined=True,axes_labels=["J","R"],aspect_ratio=1) list_frames=[] # creates a list for frames to go into for j in srange(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 value p=pvf+pt+ptr # combine into one plot list_frames.append(p) # appends this plot to list_frames ani=animate(list_frames) # creates an animation from the list of plots show(ani,iterations=3) # displays the animation only repeat 3 times