CoCalc -- Collaborative Calculation in the Cloud
SharedLab 10 / Lab 10-turnin.sagewsOpen in CoCalc
# 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]