SIR Model Blank.sagews

The following code finds the numerical solution to the following SIR model:

$S_{t+1}=S_{t}-aS_{t}I_{t}$

$I_{t+1}=I_{t}+aS_{t}I_{t}-bI_{t}$

$R_{t+1}=R_{t}+bI_{t}$

If we think of $S_{t+1}-S_{t}=\Delta S$ as the change in the susceptible population over one unit of time (a day), then $-aS_{t}I_{t}$ is the rate of change of susceptibles per day.  Let's drop the subscript notation and write $\Delta S=-aS_{t}I_{t}$  Let's consider looking at time increments of half a day.  Since $aS_{t}I_{t}$ is the rate of change of susceptibles per day, $0.5aS_{t}I_{t}$ is the rate of change of susceptibles per half day.  Let's think about this more generally... consider increments of $\Delta t$, then $aS_{t}I_{t}\Delta t$; t is the rate of change of susceptibles per whatever time increment you choose.

So, the code plots S,I,R over time using initial conditions and how the populations are changing over time.  

tdata=[] #Collecting data points for time
Sdata=[] #Collecting data pointsd for amount of susecpitible
Idata=[] #Data for amount of infected
Rdata=[] #Rdata amount of recovered
Ddata=[] #number of people dying
Datacollector=[['t', 'S', 'I', 'R','D']] #connects values that happen at the same time


#Inital conditions section AKA how much every pool/variable starts with
t=0 #start on day one
S=20000000 #Starting number of susceptible
I=90 #Starting number of infected
R=0 #starting number of recovered
D=0 #NO ONE DIED YET...

#Parameter section: constants that affect variables
a=.00000002 #likeyhood of getting infected
b=(1/31) #how long it takes to recover AKA two weeks
c=.01 #PEOPLE GET DEAD

#Euler's method
deltat=.01 

tdata.append(t) #Adds a time value to the list
Sdata.append(S) #Adds a susceptible value to the list
Idata.append(I) #adds a infect value to the list
Rdata.append(R) #adds a recovered value to the list
Ddata.append(D) #adds dead people LOL
Datacollector.append([t,S,I,R,D])

tfinal=120 #how long you are interested in AKA 100 days
finalstep=tfinal*(1/deltat) #finds total number of steps

for i in range(0,finalstep): #Calculuate tomorrow
    t=t+deltat #amount of time that's past in one step
    deltaS=-a*S*I     #calulates the change in populations
    deltaI=a*S*I-b*I-c*I 
    deltaR=b*I
    deltaD=c*I
    S=S+deltaS*deltat #Adds the change in population to the current population
    I=I+deltaI*deltat
    R=R+deltaR*deltat
    D=D+deltaD*deltat
    tdata.append(t) #Adds new values to the list
    Sdata.append(S)
    Idata.append(I)
    Rdata.append(R)
    Ddata.append(D)
    Datacollector.append([t,S,I,R,D])

#graph
Splot=list_plot(zip(tdata,Sdata), color='red', plotjoined=False)
Iplot=list_plot(zip(tdata,Idata), color='green',plotjoined=False)
Rplot=list_plot(zip(tdata,Rdata), color='blue', plotjoined=False)
Dplot=list_plot(zip(tdata,Ddata), color='pink', plotjoined=False)

Splot+Iplot+Rplot+Dplot

table(Datacollector)