 CoCalc Public FilesLab 2 / Lab2-turnin.ipynb
Author: Ryan Kien
Views : 64
Compute Environment: Ubuntu 20.04 (Default)
In [ ]:
# Name: Ryan Kien
# I worked on this code with: Fujia Guo, Nushrat Esha, Xavier Herrera

# 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


In :
#1
def dde_solve(dde, statevar, delayedvars, history, tmax, timestep):
# Check validity of delays.
if min(delayedvars.values()) < 0:
raise ValueError("This function will not work with negative delays. "
"Consider consulting a fortune teller instead.")
if any(val<timestep for val in delayedvars.values()):
raise ValueError("Time step should not be larger than delay.")

# Set up variables and delays.

delayedvars = delayedvars.items()
dde = dde.subs({v: statevar for v, delay in delayedvars if delay == 0})
delayedvars = [(v, delay) for v, delay in delayedvars if delay != 0]
allvars = [str(statevar)] + [str(v) for v, delay in delayedvars]
delays = [delay for v, delay in delayedvars]

# Set up fast functions.
dde_func = fast_float(dde, *allvars)
history_func = fast_float(history, "t")

# Adjust the timestep if necessary
mindelay = min(delays) if delays else timestep
timestepcorrectionfactor = ceil(timestep / mindelay)
timestep /= timestepcorrectionfactor

# A function to perform history lookups.
def lookup(t):
"""Does a history lookup at each delay from t, stores result in allvars[1:]"""
for i, delay in enumerate(delays):
if t - delay <= 0:
allvars[i+1] = history_func(t - delay)
else:
r = (t - delay) / timestep
n = floor(r)
r -= n
allvars[i+1] = result[n]*(1 - r) + result[n + 1]*r

# Set up for the first iteration.
result = [history_func(0)]
lookup(0)
for t in sxrange(0, tmax - timestep, timestep):
# Compute k1. Note history lookup has already been done.
allvars = result[-1]
k1 = dde_func(*allvars)
# Compute k2.
lookup(t + timestep/2)
allvars += timestep/2 * k1
k2 = dde_func(*allvars)
# Compute k3. Note history lookup has already been done.
allvars = result[-1] + timestep/2 * k2
k3 = dde_func(*allvars)
# Compute k4.
lookup(t + timestep)
allvars = result[-1] + timestep * k3
k4 = dde_func(*allvars)
# Finally, compute the RK4 weighted average.
result.append(result[-1] + (k1 + 2*k2 + 2*k3 + k4)/6 * timestep)
return result[::timestepcorrectionfactor]

In :
#2

2
In :
L = 6
Vmax = 80
a = 0.2
n = 5
h = 400
tau = 0.2

In :
#3
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: tau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
list_plot(list(zip(time, sol)), plotjoined=True) In :
#4
output=[]
for x in sol:
output.append((Vmax*x^n)/(h+x^n)*x)

In :
#5
list_plot(list(zip(time, output)), plotjoined=True)
#oscillations occur at a very high rate, and stabilize after a few minutes In :
#6a
tau = .11
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: tau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
list_plot(list(zip(time, sol)), plotjoined=True) In :
#6b
tau = 1
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: tau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
list_plot(list(zip(time, sol)), plotjoined=True) In :
#6c
tau = .5
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: tau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
list_plot(list(zip(time, sol)), plotjoined=True) In :
#7
@interact
def equation(newtau=(0.11,1,.01)):
L = 6
Vmax = 80
a = 0.2
n = 5
h = 400
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: newtau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))
#the minimum delay for persisten osciallations is when tau = .20

In :
#8
@interact
def equation(steepness=(1,25)):
L = 6
Vmax = 80
a = 0.2
n = steepness
h = 400
tau =.2
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: tau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))
#n controls the steepness

In :
#9
@interact
def equation(newtau=(0.11,1,.01), steepness=(1,25)):
L = 6
Vmax = 80
a = 0.2
n = steepness
h = 400
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: newtau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))

In :
#10
@interact
def equation(newtau=(0.11,1,.01), steepness=(1,25)):
L = 6
Vmax = 80
a = 0.2
n = steepness
h = 400
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: newtau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))
#5 is the minimum value for where persistent oscillations occur

In :
#11a
@interact
def equation(newtau=(0.11,1,.01), steepness=(1,25)):
L = 6
Vmax = 80
a = 0.2
n = steepness
h = 400
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: newtau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))
#when tau = .3, the minimum steepness value for oscillations is 4

In :
#11b
@interact
def equation(newtau=(0.11,1,.01), steepness=(1,25)):
L = 6
Vmax = 80
a = 0.2
n = steepness
h = 400
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: newtau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))
#when tau = .5, the minimum steepness value for persistent oscillations is 3

In :
#11c
@interact
def equation(newtau=(0.11,1,.01), steepness=(1,25)):
L = 6
Vmax = 80
a = 0.2
n = steepness
h = 400
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: newtau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))
#when tau = .7, the minimum steepness value for persistent oscillations is 3

In :
#11d
@interact
def equation(newtau=(0.11,1,.01), steepness=(1,25)):
L = 6
Vmax = 80
a = 0.2
n = steepness
h = 400
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: newtau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))
#when tau = 1, the minimum steepness value for persistent oscillations is 3

In :
#11e
@interact
def equation(newtau=(0.11,1,.01), steepness=(1,25)):
L = 6
Vmax = 80
a = 0.2
n = steepness
h = 400
var("x","xtau")
time = srange(0,100,0.1)
delay={xtau: newtau}
xprime = L-((Vmax * xtau^n)/(h + xtau^n))*(a*x)
sol=dde_solve(xprime, x, delay, history=0, tmax=100, timestep=0.1)
return(list_plot(list(zip(time, sol)), plotjoined=True))
#when tau = .11, the minimum steepness value for persistent oscillations is 9

In [ ]: