2020-01-16-160601.sagews

#4.2.1
var("H,P,G")
n=5
k1=.2
k2=.2
k3=.2
Hprime= ((1)/(1+G^n))-k1*H
Pprime= H- k2*P
Gprime= P-k3*G
t=srange(0,100,.1)
results=desolve_odeint([Hprime,Pprime,Gprime], ics=[1,5,8], times=t, dvars= [H,P,G])
list_plot(zip(t,results[:,0]), plotjoined=True, color="red") + list_plot(zip(t,results[:,1]), plotjoined=True) + list_plot(zip(t,results[:,2]), plotjoined=True, color= "green")

(H, P, G)

n=8
k1=.2
k2=.2
k3=.2
Hprime= ((1)/(1+G^n))-k1*H
Pprime= H- k2*P
Gprime= P-k3*G
t=srange(0,100,.1)
results=desolve_odeint([Hprime,Pprime,Gprime], ics=[1,5,8], times=t, dvars= [H,P,G])
list_plot(zip(t,results[:,0]), plotjoined=True, color="red") + list_plot(zip(t,results[:,1]), plotjoined=True) + list_plot(zip(t,results[:,2]), plotjoined=True, color="green")

n=10
k1=.2
k2=.2
k3=.2
Hprime= ((1)/(1+G^n))-k1*H
Pprime= H- k2*P
Gprime= P-k3*G
t=srange(0,100,.1)
results=desolve_odeint([Hprime,Pprime,Gprime], ics=[1,5,8], times=t, dvars= [H,P,G])
list_plot(zip(t,results[:,0]), plotjoined=True, color="red") + list_plot(zip(t,results[:,1]), plotjoined=True) + list_plot(zip(t,results[:,2]), plotjoined=True, color="green")

#4.2.2
k1=.2
k2=.2
k3=.2
Hprime= ((1)/(1+G^n))-k1*H
Pprime= H- k2*P
Gprime= P-k3*G
t=srange(0,100,.1)
results=desolve_odeint([Hprime,Pprime,Gprime], ics=[1,5,8], times=t, dvars= [H,P,G])
list_plot(zip(results[:,0], results[:,1], results[:,2]), plotjoined=True)

3D rendering not yet implemented

#4.2.3
n=10
Hprime= ((1)/(1+G^n))-k1*H
Gprime=H- k3*G
t=srange(0,10,.1)
results=desolve_odeint([Hprime, Gprime], ics=[2,3], times=t, dvars=[H,G])
list_plot(zip(t,results[:,0]),plotjoined=True, axes_labels = ["time","H values"]) + list_plot(zip(t,results[:,1]),plotjoined=True, axes_labels = ["time","G values"])

n=60
Hprime=((1)/(1+G^n))-k1*H
Gprime=H- k3*G
t=srange(0,10,.1)
results=desolve_odeint([Hprime, Gprime], ics=[2,3], times=t, dvars=[H,G])
list_plot(zip(t,results[:,0]),plotjoined=True, axes_labels = ["time","H values"]) + list_plot(zip(t,results[:,1]),plotjoined=True, axes_labels = ["time","G values"])

n=189
Hprime=((1)/(1+G^n))-k1*H
Gprime=H- k3*G
t=srange(0,10,.1)
results=desolve_odeint([Hprime, Gprime], ics=[2,3], times=t, dvars=[H,G])
list_plot(zip(t,results[:,0]),plotjoined=True, axes_labels = ["time","H values"]) + list_plot(zip(t,results[:,1]),plotjoined=True, axes_labels = ["time","G values"])