CoCalc Public FilesSIR Model Blank.sagewsOpen with one click!
Author: Mitchell Krieger
Views : 103
Compute Environment: Ubuntu 18.04 (Deprecated)

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




If we think of St+1St=ΔSS_{t+1}-S_{t}=\Delta S as the change in the susceptible population over one unit of time (a day), then aStIt-aS_{t}I_{t} is the rate of change of susceptibles per day.  Let's drop the subscript notation and write ΔS=aStIt\Delta S=-aS_{t}I_{t}  Let's consider looking at time increments of half a day.  Since aStItaS_{t}I_{t} is the rate of change of susceptibles per day, 0.5aStIt0.5aS_{t}I_{t} is the rate of change of susceptibles per half day.  Let's think about this more generally... consider increments of Δt\Delta t, then aStItΔtaS_{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