# SIR program, for studying an epidemic using Euler’s method
# First, specify the starting and ending points, stepsize, and total number of observation points
tstart=0
tfin=110
stepsize=0.5
length=((tfin-tstart)/stepsize)+1
# Next, specify values of parameters, and initial values of variables
a=0.000005
b=1/14
S=49600
I=400
R=0
t=tstart
# Set up empty lists for the values we're about to compute
Svalues=[]
Ivalues=[]
Rvalues=[]
tvalues=[]
# The following loop does three things:
# (1) stores the current values of S, I, R, and t into the lists created above;
# (2) computes the next values of S, I, R using Euler's method;
# (3) increases t by the stepsize
for i in range(length):
# Store current values
Svalues.append(S)
Ivalues.append(I)
Rvalues.append(R)
tvalues.append(t)
# Compute rates of change using SIR equations
Sprime=-a*S*I
Iprime=a*S*I-b*I
Rprime=b*I
# Net change equals rate of change times stepsize
DeltaS=Sprime*stepsize
DeltaI=Iprime*stepsize
DeltaR=Rprime*stepsize
# New values equal current values plus net change
S=S+DeltaS
I=I+DeltaI
R=R+DeltaR
t=t+stepsize
# Next time through the loop, the above new values play the role of current values
# Zip the t values with the S/I/R values into lists of ordered pairs, and create plots of these
Splot=list_plot(list(zip(tvalues,Svalues)),marker='o',color='blue')
Iplot=list_plot(list(zip(tvalues,Ivalues)),marker='o',color='red')
Rplot=list_plot(list(zip(tvalues,Rvalues)),marker='o',color='green')
# Now plot the computed S,I,R values together on a single graph, with axes labelled appropriately
SIRgraph=Splot+Iplot+Rplot
show(SIRgraph,axes_labels=['$t$ (days)','$S,I,R$ (individuals)'])