Homework Coding.sagews

#Exercise 4.2.1
#Verify for H/P/G model that values less than 8 result in stable equilibrium, but as n passes 8, equilibrium becomes unstable and a stable oscillation is created. 
#Hprime= 1/(1+G^n)-0.2*H
#Pprime= H-0.2*P
#Gprime= P-0.2*G
var("H,P,G")
t=srange(0,300,1)
#n value:4
sol=desolve_odeint([1/(1+G^4)-0.2*H,H-0.2*P,P-0.2*G],ics=[1,2,3],dvars=[H,P,G],times=t)
list_plot(zip(t,sol[:,0]),plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)+list_plot(zip(t,sol[:,2]),plotjoined=True)
#n value:6
sol=desolve_odeint([1/(1+G^6)-0.2*H,H-0.2*P,P-0.2*G],ics=[1,2,3],dvars=[H,P,G],times=t)
list_plot(zip(t,sol[:,0]),plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)+list_plot(zip(t,sol[:,2]),plotjoined=True)
#n value:7
sol=desolve_odeint([1/(1+G^7)-0.2*H,H-0.2*P,P-0.2*G],ics=[1,2,3],dvars=[H,P,G],times=t)
list_plot(zip(t,sol[:,0]),plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)+list_plot(zip(t,sol[:,2]),plotjoined=True)
#n value:8
sol=desolve_odeint([1/(1+G^8)-0.2*H,H-0.2*P,P-0.2*G],ics=[1,2,3],dvars=[H,P,G],times=t)
list_plot(zip(t,sol[:,0]),plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)+list_plot(zip(t,sol[:,2]),plotjoined=True)
#n value:10
sol=desolve_odeint([1/(1+G^10)-0.2*H,H-0.2*P,P-0.2*G],ics=[1,2,3],dvars=[H,P,G],times=t)
list_plot(zip(t,sol[:,0]),plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)+list_plot(zip(t,sol[:,2]),plotjoined=True)

(H, P, G)

#Exercise 4.2.2
t=srange(0,100,1)
#initial conditions = [1,2,3]
sol=desolve_odeint([1/(1+G^8)-0.2*H,H-0.2*P,P-0.2*G],ics=[1,2,3],dvars=[H,P,G],times=t)
list_plot(sol)
#initial conditions = [0.5,1,2]
sol=desolve_odeint([1/(1+G^8)-0.2*H,H-0.2*P,P-0.2*G],ics=[0.5,1,2],dvars=[H,P,G],times=t)
list_plot(sol)
#initial conditions = [0.1,0.5,1]
sol=desolve_odeint([1/(1+G^8)-0.2*H,H-0.2*P,P-0.2*G],ics=[0.1,0.5,1],dvars=[H,P,G],times=t)
list_plot(sol)

3D rendering not yet implemented

3D rendering not yet implemented

3D rendering not yet implemented

#Exercise 4.2.3
#Without the middleman, the negative feedback model will not oscillate
#n value:8
t=srange(0,100,1)
sol=desolve_odeint([1/(1+G^8)-0.2*H, H-0.2*G], ics=[1,2], dvars=[H,G], times=t)
list_plot(zip(t,sol[:,0]),plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)
#n value:10
sol=desolve_odeint([1/(1+G^10)-0.2*H, H-0.2*G], ics=[1,2], dvars=[H,G], times=t)
list_plot(zip(t,sol[:,0]),plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)
#n value:100
sol=desolve_odeint([1/(1+G^100)-0.2*H, H-0.2*G], ics=[1,2], dvars=[H,G], times=t)
list_plot(zip(t,sol[:,0]),plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)

#Exercise 4.1.7
#a=2
var("X,V")
t=srange(0,10)
sol=desolve_odeint([V,-X-(2*V^3-V)],ics=[2,4],dvars=[X,V],times=t)
list_plot(zip(t,sol[:,0]),axes_labels=["Time","H/P/G"], plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)#a=2
#a=0.5
var("X,V")
t=srange(0,100)
sol=desolve_odeint([V,-X-(0.5*V^3-V)],ics=[2,4],dvars=[X,V],times=t)
list_plot(zip(t,sol[:,0]),axes_labels=["Time","H/P/G"], plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)#a=2
#a=50
var("X,V")
t=srange(0,100)
sol=desolve_odeint([V,-X-(50*V^3-V)],ics=[2,4],dvars=[X,V],times=t)
list_plot(zip(t,sol[:,0]),axes_labels=["Time","H/P/G"], plotjoined=True)+list_plot(zip(t,sol[:,1]),plotjoined=True)

(X, V)

(X, V)

(X, V)

#Exercise 6.5.1
list1=[(1,1)]
M=matrix(RDF, [[1.2,0], [0,1.25]])
for i in srange (0,100):
    current_state=list1[-1]
    next_state=M*vector(current_state)
    list1.append(next_state)
list_plot(list1,plotjoined=True)

#Exercise 6.5.1
list2=[(6,3)]
M=matrix(RDF, [[1.2,0], [0,0.9]])
for i in srange (0,50):
    current_state=list2[-1]
    next_state=M*vector(current_state)
    list1.append(next_state)
list_plot(list2,plotjoined=True)

#Exercise 6.5.3
var("X,Y")
Xand1=1.2*X
Yand1=0.9*Y
t=srange(0,100)
#sol=desolve_odeint([Xand1,Yand1],ics=[6,3],dvars[X,Y],times=t)
#plot(sol)

die=zip(Xand1,t)
plot(die)

(X, Y)

Error in lines 5-5 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1013, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> TypeError: zip argument #1 must support iteration

#Lab 9 
# Lab 9: Linear approximations to surfaces and vector fields

# Name: Tasha Tsao
# 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 put each exercise in a new cell with an exercise number at the top.
# Use enough comments that you and the grader can understand your code.
# Label axes on all graphs.

# 1
var("y")
plot3d(x^2+y^2,(x,-4,4),(y,-4,4))

y

3D rendering not yet implemented

#2
plot(x^2)+plot(0, color="red", axes_labels=("x","f(x,y)"))

#3
plot(y^2, axes_labels=("y","f(x,y)"))+plot(0, color="red")

#4
plot3d(x^2+y^2,(x,-4,4),(y,-4,4))+plot3d(0,(x,-4,4),(y,-4,4), color="violet")

3D rendering not yet implemented

#5
plot3d(2*(x-2)-(y-3)-7,(x,-3,3),(y,-3,3))+point3d([2,3,-7], color="red")

3D rendering not yet implemented

#6
#plot plane tangent to the point (1,-1,2)
#Tangent to X at 1 is 2; tangent to Y at -1 is -2 
plot3d(x^2+y^2,(x,-3,3),(y,-3,3),opacity=0.5)+point3d([1,-1,2], color="red")+plot3d(2*(x-1)-2*(y+1)+2,(x,-3,3),(y,-3,3))

3D rendering not yet implemented

#7
plot3d(2*x+y^2,(x,-5,5),(y,-5,5),opacity=0.8)+point3d([2,1,5],color="red")+plot3d(2*(x-2)+2*(y-1)+5,(x,-5,5),(y,-5,5),opacity=0.5)

3D rendering not yet implemented

#8
#This is a stable spiral equilibrium 
sol=desolve_odeint([0.2*X*(1-X/7)-(X/(1+X))*Y,1/10*Y*(1-Y/X)],ics=[1,1],dvars=[X,Y],times=srange(0,100))
var("X,Y")
#Regular magnification
plot_vector_field((0.2*X*(1-X/7)-(X/(1+X))*Y,1/10*Y*(1-Y/X)),(X,-1,1),(Y,-1,1))+point((0.239,0.24),size=30,color="red")+list_plot(sol,plotjoined=True)

(X, Y)

#Zoomed in
plot_vector_field((0.2*X*(1-X/7)-(X/(1+X))*Y,1/10*Y*(1-Y/X)),(X,0,1),(Y,0,1))+point((0.239,0.24),size=30,color="red")+list_plot(sol,plotjoined=True)

#9
var("x,y")
g=0.2*x*(1-x/7)-(x/(1+x))*y
h=1/10*y*(1-y/x)

(x, y)

jacM=jacobian((g,h),(x,y))
show(jacM)

$\displaystyle \left(\begin{array}{rr} -0.0571428571428571 \, x - \frac{y}{x + 1} + \frac{x y}{{\left(x + 1\right)}^{2}} + 0.200000000000000 & -\frac{x}{x + 1} \\ \frac{y^{2}}{10 \, x^{2}} & -\frac{y}{5 \, x} + \frac{1}{10} \end{array}\right)$

#10
jacPointM=jacM.subs({x:1, y:2})
jacPointM

[-0.357142857142857 -1/2] [ 2/5 -3/10]

jacPointM.eigenvalues()

[-2/35*I*sqrt(61) - 23/70, 2/35*I*sqrt(61) - 23/70]

#Because the eigenvalues contain an imaginary component, the trajectories will contain rotation. Because the real component of the eigenvalues are negative, the point is stable. Together, the equilibirum is a stable spiral. 

#11
sol2=desolve_odeint([-0.357142857142857*x-(1/2)*y, (2/5)*x-(3/10)*y],dvars=[x,y],times=srange(0,100),ics=[0.5,0.5])
#This point is a stable spiral, just like the equilibrium point of the nonlinear system. 
plot_vector_field((-0.357142857142857*x-(1/2)*y, (2/5)*x-(3/10)*y),(x,0,1),(y,0,1))+point((0.239,0.24),size=30,color="red")+list_plot(sol, plotjoined=True)

︠ef343490-1cb2-4d3c-b542-263630f57a55︠