CoCalc Shared FilesLab 9 / Lab9-turnin.sagews
Author: Kelly Truong
Views : 71
# Lab 9: :(

# Name: Kelly
# 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 Exercise 1. By hand, find the equilibria of the system R0 = aR + bJ, J0 = cR + dJ.

#(0,0)

#Exercise 2. Declare a, b, c, d, R and J as symbolic variables.
a = 0.25
b=0.5
c=0.65
d=0.75
var("R,J")
Rprime(R)= a*R +b*J
Jprime(J)= c*R + d*J


(R, J)
#Exercise 3. Make a list of eight points, two in each quadrant of the R−J plane, to serve as initial conditions for simulations. Recall that points are defined using  parentheses: (R1, J1) is a point.
Rpoints=[-3,-1,1,3]
Jpoints=[-2,4,2,-4]
coordinates = zip(Rpoints,Jpoints)
coordinates

[(-3, -2), (-1, 4), (1, 2), (3, -4)]
#Exercise 4. Pick a positive value for a and d. (For simplicity, make them equal to each other.) Plot the system’s vector field, letting both R and J range from -10 to 10. Is the equilibrium stable or unstable?
a=0.25
d=0.25
Rprime(R)= a*R
Jprime(J)= d*J
plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10))

#IT IS UNSTABLE

#Exercise 5. Now, make a and d negative and plot the vector field. Is the equilibrium stable or unstable?
a=-0.50
d=-0.50
Rprime(R)= a*R
Jprime(J)= d*J

plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10))

#it is stable

#Exercise 6. Interpret the results of the previous two exercises in terms of the relationship between Romeo and Juliet. HINT: If you have trouble doing this from just the vector field, try running some simulations for various initial conditions using desolve_odeint.
# when a & d are both positive, as romeo's love for juliet increases, juliet's love for romeo increases as well . When a & d are negative, as romeo's love for juliet decreases, juliets love for him decreases as well

#Exercise 7. Keeping the absolute value of a and d the same, make one ofthem positive and the other negative. Plot the vector field as before. Does the equilibrium appear stable or unstable
a=-0.50
d=0.50
Rprime(R)= a*R
Jprime(J)= d*J

plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10))

#It is semi stable, a saddle point

#Exercise 8. Overlay a large red point at the equilibrium on top of the vector field (use the point command to plot points) and assign this plot to a variable. Then, use a for loop to simulate the dynamics of the system for each initial condition in the list you defined earlier, make a plot of the resulting trajectory and add it to the existing vector field. View the result. HINT: If your trajectories go outside the vector field, display the plot using show and specify xmin, xmax, ymin and ymax, as necessary.

p=plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10)) +  point([0,0], size =80, color  = "red")

for i in coordinates :
traj = desolve_odeint([Rprime,Jprime],ics=i,times=srange(0,50,0.1),dvars=[R,J])
p += list_plot(zip(traj[:,0],traj[:,1]), color = "purple", thickness = 2, plotjoined= True)

show(p,xmin=-10,xmax=10, ymin = -10, ymax=10)

#Exercise 9. Describe what you see in purely visual and dynamical terms, without referring to the meaning of the model. Is the equilibrium point stable or unstable?
#the trajectory moves away from the origin going up and down and towards it going left and right. that makes the equillibrium point unstable

#Exercise 10. Run simulations and plot time series for a couple of different initial conditions. Draw on the vector field, trajectories and time series to describe what is happening to Romeo and Juliet’s relationship in this model.
p=plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10)) +  point([0,0], size =80, color  = "red")

for i in coordinates :
traj = desolve_odeint([Rprime,Jprime],ics=i,times=srange(0,50,0.1),dvars=[R,J])
p += list_plot(zip(traj[:,0],traj[:,1]), color = "purple", thickness = 2, plotjoined= True)

show(p,xmin=-10,xmax=10, ymin = -10, ymax=10)



#Exercise 11. Pick a positive value for b and c, keeping them equal to each other. Plot the vector field. What kind of equilibrium point appears?
var("R","J")
b = 0.5
c = 0.5
Rprime = b*J
Jprime = c*R

(R, J)


plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10))

#a saddle node bifurcation

#Exercise 12. What happens if b and c are both negative?
b = -0.5
c = -0.5

Rprime = b*J
Jprime = c*R
plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10))

#a saddle node bifurcation in the other direction

#Exercise 13. Keeping the absolute value of b and c the same, make one of them positive and the other negative. Plot the vector field as before.
b = -0.5
c = 0.5

Rprime = b*J
Jprime = c*R
plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10))

#it's a neutral center

#Exercise 14. Overlay a large red point at the equilibrium on top of the vector field (use the point command to plot points) and assign this plot to a variable. Then, use a for loop to simulate the dynamics of the system for each initial condition in the list you defined earlier, make a plot of the resulting trajectory and add it to the existing vector field. View the result.

p=plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10)) +  point([0,0], size =80, color  = "red")

for i in coordinates :
traj = desolve_odeint([Rprime,Jprime],ics=i,times=srange(0,50,0.1),dvars=[R,J])
p += list_plot(zip(traj[:,0],traj[:,1]), color = "purple", thickness = 2, plotjoined= True)

show(p,xmin=-10,xmax=10, ymin = -10, ymax=10)

#Exercise 15. Describe what you see in purely visual and dynamical terms, without referring to the meaning of the model. Is the equilibrium point stable or unstable?
#stable since all trajectories spiral towards the origin

#Exercise 16. Run simulations and plot time series for a couple of different initial conditions. Draw on the vector field, trajectories and time series to describe what is happening to Romeo and Juliet’s relationship in this model.
p=plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10)) +  point([0,0], size =80, color  = "red")

for i in coordinates :
traj = desolve_odeint([Rprime,Jprime],ics=i,times=srange(0,50,0.1),dvars=[R,J])
p += list_plot(zip(traj[:,0],traj[:,1]), color = "purple", thickness = 2, plotjoined= True)

show(p,xmin=-10,xmax=10, ymin = -10, ymax=10)

#Exercise 17. Plot the vector field for this system. Since it will be hard to tell exactly what’s going on from the vector field alone, pick an initial condition, run a simulation and overlay the trajectory on top of the vector field.
var("R","J")
Rprime(R)= J
Jprime(J)=-R-0.05*J
p=plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10)) +  point([0,0], size =80, color  = "red")

traj = desolve_odeint([Rprime,Jprime],ics=coordinates[0],times=srange(0,50,0.1),dvars=[R,J])
p += list_plot(zip(traj[:,0],traj[:,1]), color = "purple", thickness = 2, plotjoined= True)

show(p,xmin=-10,xmax=10, ymin = -10, ymax=10)

(R, J)
#Exercise 18. Plot time series for R and J. Describe what is happening in terms of Romeo and Juliet’s relationship
for i in coordinates :
traj = desolve_odeint([Rprime,Jprime],ics=i,times=srange(0,50,0.1),dvars=[R,J])
p += list_plot(zip(traj[:,0],traj[:,1]), color = "purple", thickness = 2, plotjoined= True)

show(p,xmin=-10,xmax=10, ymin = -10, ymax=10)


#Exercise 19. This kind of equilibrium is called a stable spiral. Why is the word “stable” appropriate?
#because the trajectory is spiraling towards the eq point

#Exercise 20. Plot the vector field for this system and overlay a trajectory ontop of it.
var("R","J")
Rprime(R)= J
Jprime(J)=-R+0.05*J
p=plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10)) +  point([0,0], size =80, color  = "red")

traj = desolve_odeint([Rprime,Jprime],ics=coordinates[2],times=srange(0,50,0.1),dvars=[R,J])
p += list_plot(zip(traj[:,0],traj[:,1]), color = "purple", thickness = 2, plotjoined= True)

show(p,xmin=-10,xmax=10, ymin = -10, ymax=10)

(R, J)
#Exercise 21. Plot time series for R and J. Describe what is happening in terms of Romeo and Juliet’s relationship.
var("R","J")
Rprime(R)= J
Jprime(J)=-R+0.05*J
p=plot_vector_field((Rprime,Jprime),(R,-10,10),(J,-10,10)) +  point([0,0], size =80, color  = "red")

for i in coordinates :
traj = desolve_odeint([Rprime,Jprime],ics=i,times=srange(0,50,0.1),dvars=[R,J])
p += list_plot(zip(traj[:,0],traj[:,1]), color = "purple", thickness = 2, plotjoined= True)

show(p,xmin=-10,xmax=10, ymin = -10, ymax=10)

(R, J)
#As romeo's love for Juliet increases, Juliet's love for him decreases with some increase in proportion to her own love

#Exercise 22. What might this kind of equilibrium be called?
#unstalbe spiral, a source