SharedLab 10 / Lab 10-turnin.sagewsOpen in CoCalc
Author: Patrick Ku
Views : 12
# Lab 9 # Name: Patrick Ku # I worked on this code with: # Please do all of your work for this week's lab in this worksheet. If # you wish to create other worksheets for scratch work, you can, but # this is the one that will be graded. You do not need to do anything # to turn in your lab. It will be collected by your TA at the beginning # of (or right before) next week’s lab. # Be sure to clearly label which question you are answering as you go and to # use enough comments that you and the grader can understand your code.
#1 # x prime represents plants because they have logistic growth, and they have a saturation death rate of being eaten by the herivore. # y prime is the herbivores becasue their birth rate is proportional to the amount of plants they consume, and they have a natural death rate and a saturation function descibing predation of the carnivores. # z prime is the carnivores because their birth rate is proportional to the amount of herbivores they consume, and they have a natural death rate and no predation.
#2 _=var('X,Y,Z') a1=5 b1=3 a2=0.1 b2=2 d1=0.4 d2=0.01 xprime=X*(1-X)-(a1*X*Y)/(1+b1*X) yprime=(a1*X*Y)/(1+b1*X)-d1*Y-(a2*Y*Z)/(1+b2*Y) zprime=(a2*Y*Z)/(1+b2*Y)-d2*Z time=srange(0,1000,1) sol=desolve_odeint([xprime,yprime,zprime],dvars = [X,Y,Z],times = time,ics = [1,0.5,8]) list_plot(zip(time,sol[:,0]),color="green",axes_labels=["time","population"],legend_label="plants",plotjoined=True)+list_plot(zip(time,sol[:,1]),color="blue",legend_label="herbivores",plotjoined=True)+list_plot(zip(time,sol[:,2]),color="red",legend_label="carnivores",plotjoined=True) list_plot(sol,axes_labels=["X'","Y'","Z'"],plotjoined=True)
3D rendering not yet implemented
#3 sol1=desolve_odeint([xprime,yprime,zprime],dvars = [X,Y,Z],times = time,ics = [1.5,1,4]) list_plot(zip(time,sol[:,0]),color="green",axes_labels=["time","population"],legend_label="plants A",plotjoined=True)+list_plot(zip(time,sol[:,1]),color="blue",legend_label="herbivores A",plotjoined=True)+list_plot(zip(time,sol[:,2]),color="red",legend_label="carnivores A",plotjoined=True)+list_plot(zip(time,sol1[:,0]),color="turquoise",axes_labels=["time","population"],legend_label="plants B",plotjoined=True)+list_plot(zip(time,sol1[:,1]),color="indigo",legend_label="herbivores B",plotjoined=True)+list_plot(zip(time,sol1[:,2]),color="magenta",legend_label="carnivores B",plotjoined=True) list_plot(sol,axes_labels=["X'","Y'","Z'"],plotjoined=True,legend_label="trajectory A")+list_plot(sol1,axes_labels=["X'","Y'","Z'"],plotjoined=True, legend_label="trajectory B",color="red")
3D rendering not yet implemented
#4 #the time series look more similar because it has the same bounded region and most of the plot is roughly the same #the trajectories look more different becasue the red trajectory started at a much lower state
#5 @interact def f(b1=(2,3,0.1)): _=var('X,Y,Z') a1=5 a2=0.1 b2=2 d1=0.4 d2=0.01 xprime=X*(1-X)-(a1*X*Y)/(1+b1*X) yprime=(a1*X*Y)/(1+b1*X)-d1*Y-(a2*Y*Z)/(1+b2*Y) zprime=(a2*Y*Z)/(1+b2*Y)-d2*Z time=srange(0,1000,1) solb1=desolve_odeint([xprime,yprime,zprime],dvars = [X,Y,Z],times = time,ics = [2,1.5,7]) show(list_plot(zip(time,solb1[:,0]),color="green",axes_labels=["time","population"],legend_label="plants",plotjoined=True,ymax=12)+list_plot(zip(time,solb1[:,1]),color="blue",legend_label="herbivores",plotjoined=True)+list_plot(zip(time,solb1[:,2]),color="red",legend_label="carnivores",plotjoined=True)) show(list_plot(solb1,axes_labels=["X'","Y'","Z'"],plotjoined=True))
Interact: please open in CoCalc
#6 #when you increase b1, the plot starts its chaotic and non-periodic oscillatory behavior earlier.
#7 t=srange(0,21,1) iterate=[0.2] r=2.5 for i in t[1:]: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) list_plot(zip(t,iterate),plotjoined=True)
t=srange(0,21,1) iterate=[0.2] r=3.1 for i in t[1:]: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) list_plot(zip(t,iterate),plotjoined=True)
t=srange(0,21,1) iterate=[0.2] r=3.5 for i in t[1:]: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) list_plot(zip(t,iterate),plotjoined=True)
t=srange(0,21,1) iterate=[0.2] r=3.9 for i in t[1:]: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) list_plot(zip(t,iterate),plotjoined=True)
# For initial conditions 2.5, 3.1, 3.5, there seems to be either periodic oscillations or no oscillations. For initial conditon 3.9, the oscillation is completely chaotic, showing chaotic behavior.
#8 t=srange(0,21,1) iteratea=[0.4] r=3.9 for i in t[1:]: f=r*iteratea[-1]*(1-iteratea[-1]) iteratea.append(f) list_plot(zip(t,iteratea),plotjoined=True,axes_labels=["time","Xn+1"])
t=srange(0,21,1) iterateb=[0.401] r=3.9 for i in t[1:]: f=r*iterateb[-1]*(1-iterateb[-1]) iterateb.append(f) list_plot(zip(t,iteratea),plotjoined=True,axes_labels=["time","Xn+1"])+list_plot(zip(t,iterateb),plotjoined=True, color="turquoise")
#the plot changes completely even though the initial condition is only changed slightly
#9 t=srange(0,21,1) iterate=[0.35] r=3.8 for i in t[1:]: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) iterate1=[0.351] for i in t[1:]: f=r*iterate1[-1]*(1-iterate1[-1]) iterate1.append(f) list_plot(zip(t,iterate),plotjoined=True,axes_labels=["time","Xn+1"])+list_plot(zip(t,iterate1),plotjoined=True, color="turquoise")
t=srange(0,21,1) iterate=[0.33] r=4 for i in t[1:]: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) iterate1=[0.329] for i in t[1:]: f=r*iterate1[-1]*(1-iterate1[-1]) iterate1.append(f) list_plot(zip(t,iterate),plotjoined=True,axes_labels=["time","Xn+1"])+list_plot(zip(t,iterate1),plotjoined=True, color="turquoise")
t=srange(0,21,1) iterate=[0.1] r=3.75 for i in t[1:]: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) iterate1=[0.101] for i in t[1:]: f=r*iterate1[-1]*(1-iterate1[-1]) iterate1.append(f) list_plot(zip(t,iterate),plotjoined=True,axes_labels=["time","Xn+1"])+list_plot(zip(t,iterate1),plotjoined=True, color="turquoise")
#the function is very sensitive to different initial conditions. The plots behave very differently (or eventually) when the initial conditions are changed
#10 time=srange(0,8,1) difference=[] for i in time: s=abs(iteratea[i]-iterateb[i]) difference.append(s) difference
[0.00100000000000000, 0.000776100000000057, 0.00264170997171922, 0.00551595115092895, 0.00841197403009952, 0.0213825233808967, 0.0101309498031747, 0.0363628324247411]
#extra scratch work # \/ \/ \/ \/ var('x,y') dx=2*x dy=3*(y^2)-6 plot_vector_field((dx,dy),(x,-5,5),(y,-5,5))+point((0,sqrt(2)),size=40)+point((0,-sqrt(2)),size=40)
(x, y)
dx=6*x^2+y dy=x-6*y plot_vector_field((dx,dy),(x,-0.05,0.05),(y,-0.05,0.05))+point((0,0),size=40)+point((-1/36,-1/216),size=40)
dx=2*x+3 dy=-4*y+4 plot_vector_field((dx,dy),(x,-1.5,1.5),(y,-1.5,1.5))+point((-3/2,1),size=40)
t=srange(0,101,1) iterate=[0.01] r=4 for i in t: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) list_plot(zip(t,iterate),plotjoined=True)
t=srange(0,21,1) iterate=[0.6] r=4 for i in t: f=r*iterate[-1]*(1-iterate[-1]) iterate.append(f) list_plot(zip(t,iterate),plotjoined=True) iterate
[0.600000000000000, 0.960000000000000, 0.153600000000000, 0.520028160000000, 0.998395491228058, 0.00640773729417265, 0.0254667127877661, 0.0992726373102061, 0.357670323166729, 0.918969052370147, 0.297859732624244, 0.836557249221031, 0.546916871987090, 0.991195228491788, 0.0349089900276012, 0.134761409771416, 0.466403088831347, 0.995484990239702, 0.0179784977886481, 0.0706210856236468, 0.262534991555938, 0.774441479058645]