Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.