20200402.sir-epidemic.sage

repnum = 3

var('s i r t')
p = desolve_system_rk4([-repnum*i*s, repnum*i*s-i, i], [s,i,r], ics=[0, N(1/repnum), N((repnum-log(repnum)-1)/repnum), N(log(repnum)/repnum)], ivar=t, end_points=[-5,8], step=0.01)
list_plot([(tmp[0], tmp[1]) for tmp in p], color="green", plotjoined=True) + list_plot([(tmp[0], tmp[2]) for tmp in p], color="red", plotjoined=True) + list_plot([(tmp[0], tmp[3]) for tmp in p], color="blue", plotjoined=True)
list_plot([(tmp[0], tmp[1]) for tmp in p], color="green", plotjoined=True, ymin=1e-5, ymax=1, scale="semilogy") + list_plot([(tmp[0], tmp[2]) for tmp in p], color="red", plotjoined=True, ymin=1e-5, ymax=1, scale="semilogy") + list_plot([(tmp[0], tmp[3]) for tmp in p], color="blue", plotjoined=True, ymin=1e-5, ymax=1, scale="semilogy")

G = N(-lambert_w(-repnum*exp(-repnum))/repnum)
qx = [(N(j/100), N(((1-G)^2/sqrt(G))*exp(repnum*(1-G)*(j/100)))) for j in range(-500,801)]
q = [(t, ((1-G)^2 + G*x)/((1-G)^2 + x), ((1-G)^4*x)/(((1-G)^2 + G*x)*((1-G)^2 + x)), ((1-G)*G*x)/((1-G)^2 + G*x)) for (t,x) in qx]
list_plot([(tmp[0], tmp[1]) for tmp in q], color="green", plotjoined=True) + list_plot([(tmp[0], tmp[2]) for tmp in q], color="red", plotjoined=True) + list_plot([(tmp[0], tmp[3]) for tmp in q], color="blue", plotjoined=True)
list_plot([(tmp[0], tmp[1]) for tmp in q], color="green", plotjoined=True, ymin=1e-5, ymax=1, scale="semilogy") + list_plot([(tmp[0], tmp[2]) for tmp in q], color="red", plotjoined=True, ymin=1e-5, ymax=1, scale="semilogy") + list_plot([(tmp[0], tmp[3]) for tmp in q], color="blue", plotjoined=True, ymin=1e-5, ymax=1, scale="semilogy")