CoCalc Public FilesLab 6 / Lab6-turnin.sagews
Author: Kelly Truong
Views : 374
# Lab 6: Simulating Differential Equations in SageMath

# Name: Kelly
# I worked on this code with: Haylee

# 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
def eulers3(diff,initial,steps,time) :
range = srange(0,time+steps,steps)
values = [initial]
for i in range :
a = diff(values[-1])*steps
values.append(values[-1]+a)
return list_plot(zip(range,values), axes_labels=["time","values"],plotjoined=true)


xPrime(x)=-2.1*x
eulers3(xPrime,1,1,10)
eulers3(xPrime,1,0.5,10)
eulers3(xPrime,1,0.1,10)

# the graph when the step size is 0.1 looks most similar to how the function of derivative should be

#2 Exercise 2. Run the simulation above and plot the results with list plot.
t = srange(0,10,0.1)
solution = desolve_odeint(-2.1*x,1,t,x)
list_plot(zip(t,solution), axes_labels=["time","solution"])
#nice

#3 Exercise 3. Look at the variable sol itself. What does the output of desolve odeint consist of?
#the derivative/slope of the function -2.1x at a certain point

#4Exercise 4. Suppose one rabbit population starts off with 1000 individuals and grows at a per-capita rate of 0.03 per year while another starts off with 2000 individuals grows at a per-capita rate of 0.02 per year. Use simulation and plotting to find the approximate time at which the first population becomes larger than the second population.

time =  srange(0,100)

rabs = desolve_odeint(0.03*x,ics=1000,times=time,dvars=x)
bits = desolve_odeint(0.02*x,ics=2000,times=time,dvars=x)

list_plot(rabs) +list_plot(bits, color = "purple")
#around time 70 the first rabbit population exceeds the 2nd population

#5Exercise 5. Simulate the logistic equation,
time = srange(0,100.05,0.05)
logistic = desolve_odeint(0.1*x*(1-(x/1000)),ics=500,times=time,dvars=x)
list_plot(logistic)

#6 Exercise 6. Run the predator-prey simulation above and plot time series of the results
var("N, P")
t = srange(0,100,0.1)
sol=desolve_odeint([0.5*N - 0.01*N*P, 0.5*0.01*N*P - 0.2*P], ics=[50,75],
dvars=[N,P], times=t)
list_plot(zip(t, sol[:,0])) + list_plot(zip(t,sol[:,1]), color="red")
#syntax [:,0] accesses all elements with the 2nd coordinate 0

(N, P)
#7ook at the variable sol itself. What does the output of desolve odeint consist of for a two-variable model? What does sol[:,0] output? What does sol[:,1] output?

#for a two variable model it contains the changes of both of the populations in arrays. the notation sol[:,0] accesses all elements in the column with 0 as the second coordinate. likely sol[:,1] accesses all elements in the colum with 1 as the second coordinate

#8 Exercise 8. Plot the trajectory of the predator-prey simulation above.
var("N, P")
t = srange(0,100,0.1)
sol=desolve_odeint([0.5*N - 0.01*N*P, 0.5*0.01*N*P - 0.2*P], ics=[50,75],
dvars=[N,P], times=t)
list_plot(zip(sol[:,1], sol[:,0]),plotjoined=true)

(N, P)
#9 Simulate the predator-prey model for two more sets of initial values and plot the results as time series and trajectories.
var("N, P")
t = srange(0,100,0.1)
sol=desolve_odeint([0.5*N - 0.01*N*P, 0.5*0.01*N*P - 0.2*P], ics=[100,200],
dvars=[N,P], times=t)
list_plot(zip(sol[:,1], sol[:,0]),plotjoined=true)
list_plot(zip(t,sol[:,1]), legend_label ="N") + list_plot(zip(t,sol[:,0]), color = "red", legend_label = "P")

sol=desolve_odeint([0.5*N - 0.01*N*P, 0.5*0.01*N*P - 0.2*P], ics=[1234,2345],
dvars=[N,P], times=t)
list_plot(zip(sol[:,1], sol[:,0]),plotjoined=true)
list_plot(zip(t,sol[:,1]), legend_label ="N") + list_plot(zip(t,sol[:,0]), color = "red", legend_label = "P")

(N, P)
#Exercise 10. How does the prey birth rate affect predator-prey dynamics? Simulate the model for three values of the birth rate and plot the results as time series. (Keep your initial conditions the same across simulations.) What seems to happen as the birth rate increases?

#I didn't want to copy paste it over and over again so I made a function where the only thing that changes is the birthrate of the N population
def laziness(birthrate):
var("N, P")
t = srange(0,100,0.1)
sol=desolve_odeint([birthrate*N - 0.01*N*P, 0.5*0.01*N*P - 0.2*P], ics=[100,200],
dvars=[N,P], times=t)
return list_plot(zip(t,sol[:,1]), legend_label ="N") + list_plot(zip(t,sol[:,0]), color = "red", legend_label = "P")

laziness(0.7)
laziness(0.9)
laziness(1.1)

#increased birthrate of the prey, the prey population begins to surpass the predator population, the predator population also decreases

#11 Run a parameter scan to test the effect of the prey birth rate on the amplitude of the prey oscillations. Use birth rate values ranging from 0.1 to 1.5 in steps of 0.1. Plot your results and describe any trend you see. HINT: You can use the functions max and min, which return the maximum and minimum values of a list, to measure oscillation amplitude.

#We can do many numerical experiments to see the effects of changing a parameter. This is called a parameter scan.

#amplitude is max-min/ 2
# sol[:,1] prey    &   sol[:,0]predator

birthrate = srange(0.1,1.5,0.1)
var("N, P")
amplitudes = []
for b in birthrate :
t = srange(0,100,0.1)
sol=desolve_odeint([b*N - 0.01*N*P, 0.5*0.01*N*P - 0.2*P], ics=[100,200], dvars=[N,P], times=t)
amplitudes.append((max(sol[:,1])-min(sol[:,1]))/2.0)

list_plot(zip(birthrate,amplitudes), plotjoined =true)


(N, P)
# as the prey birthrate increases the amplitude decreases

#Exercise 12. Use a similar parameter scan to find the effect of another parameter or initial condition on the dynamics of the Lotka-Volterra model. HINT: Run some simulations manually first to establish a reasonable range for your parameter or initial condition

predation = srange(0.1,1.5,0.1)
var("N, P")
amplitudes = []
for r in predation :
t = srange(0,100,0.1)
sol=desolve_odeint([0.5*N - r*N*P, 0.5*r*N*P - 0.2*P], ics=[100,200], dvars=[N,P], times=t)
amplitudes.append((max(sol[:,1])-min(sol[:,1]))/2.0)

list_plot(zip(predation,amplitudes), plotjoined =true)

(N, P)

#13 Exercise 13. Change one of the parameters in the predator-prey model above and plot the vector field. You can change the plotting range as necessary to get a good look. Describe your output.
v(R,F) = [0.5*R - 0.01*R*F, 0.005*R*F - 0.2*F]
plot_vector_field(v, (R,0,100), (F,0,100), axes_labels=["R", "F"])
#changed parameter for rabbit birth rate, the trajectory loops, higher rabbit population increases the fox population which decreases the rabbit population (negative feedback loop?)


#Exercise 14. In another version of the Romeo-Juliet system you saw in class, both Romeo and Juliet are completely self-absorbed. Romeo’s love changes at a rate proportional to itself, with a proportionality constant of 0.5. Juliet’s love also changes at a rate proportional to itself, with a proportionality constant of -0.25. Write the model for this system and plot its vector field, showing both positive and negative values of R and J. Describe your output.
v(R,J)=[0.5*R,-0.25*J]
plot_vector_field(v, (R,0,100),(J,0,100), axes_labels=["R","J"])

#as juliets love decreases romeo's love increases