CoCalc Public FilesProject.sagewsOpen with one click!
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)