CoCalc Public FilesProject.sagews
Authors: Jessica Booth, Laine Nowak
Views : 64
Compute Environment: Ubuntu 20.04 (Default)
#Plotting All Populations separately
var ("P, H, G, C")
K=20
r=0.27

a1=1
b1=1
a2=1
b2=1
a3=1
b3=1
a4=1
b4=1
d1=0.3
d2=0.8
d3=.3
k2=.2
#manipulating k1 so we do not define it here
#P_Prime=r*P*(1-(P/K)) - ((a1*P)/(1+b1*P))*H -((a4*P)/(1+b4*P))*G
#H_Prime=c1*((a1*P)/(1+b1*P))*H - (d1*H) - ((a2*H)/(1+(b2*H)))*C -k1*H*G
#G_prime =c3*((a3*P)/(1+b3*P))*G - (d3*G) - ((a4*G)/(1+(b4*G)))*C -k2*H*G
#C_Prime=c2*((a2*H)/(1+(b2*H)))*C - (d2*C) + c4*((a4*G)/(1+(b4*G)))*C
t=srange(0,4000)
P_prime=0.1*P*(1-P/500) - (P/(100+P))*H - (P/(100+P))*G
H_prime=0.3*(P/(100+P))*H -.09*H - (H/(100+H))*C -(50*H*G)
G_prime=0.6*(P/(100+P))*G -.3*G -(G/(100+G))*C -(30*H*G)
C_prime=0.7*(H/(100+H))*C + 1.8*(G/(100+G))*C -.1*C

#solving the four differential equations
#Model=desolve_odeint([(r*P*(1-(P/K)) - ((a1*P)/(1+b1*P))*H -((a4*P)/(1+b4*P))*G,c1*((a1*P)/(1+b1*P))*H - (d1*H) - ((a2*H)/(1+(b2*H)))*C -(k1*H*G),(c3*((a3*P)/(1+b3*P))*G - (d3*G) - ((a4*G)/(1+(b4*G)))*C -(k2*H*)G), c2*((a2*H)/(1+(b2*H)))*C - (d2*C) + c4*((a4*G)/(1+(b4*G)))*C)], dvars=[P,H, G, C], times=t, ics=[10,20, 15, 5])
Model=desolve_odeint([P_prime,H_prime,G_prime,C_prime], dvars=[P,H,G,C], times=t,ics=[3,2,2,1])
#plotting the four equations for P,H,G,C
list_plot(list(zip(t,Model[:,0])), color="green", legend_label="Plants", plotjoined=True, axes_labels=["$Time$", "$Population Sizes$"], thickness=2)
list_plot(list(zip(t,Model[:,1])), color="brown", legend_label="Herbivore 1", plotjoined=True, thickness=2)
list_plot(list(zip(t,Model[:,2])), color="darkblue", legend_label="Herbivore 2", plotjoined=True, thickness=2)
list_plot(list(zip(t,Model[:,3])), color="violet", legend_label="Carnivore", plotjoined=True, thickness=2)

#try:
#increase interspecies competition
#make half densities higher
#make growth rates lower

(P, H, G, C)
#Possible constants of interaction between herbivores
#decreasing sigmoid
f(x)=(1/(1+x^2))
plot(f, xmin=0,xmax=30, ymin=0)

#increasing sigmoid
h(x)=(x/(1+x^3))
plot(h, xmin=0, ymin=0)

#inverted parabola
g(x)=(-(x^2)-15)+35
plot(g)