Open with one click!
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 # use enough comments that you and the grader can understand your code.
In [2]:
#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[0] = result[-1] k1 = dde_func(*allvars) # Compute k2. lookup(t + timestep/2) allvars[0] += timestep/2 * k1 k2 = dde_func(*allvars) # Compute k3. Note history lookup has already been done. allvars[0] = result[-1] + timestep/2 * k2 k3 = dde_func(*allvars) # Compute k4. lookup(t + timestep) allvars[0] = 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 [12]:
#2 menu = {"Ice Cream":4.95, "Cookies":2, "Candy":1} menu["Cookies"]
2
In [95]:
L = 6 Vmax = 80 a = 0.2 n = 5 h = 400 tau = 0.2
In [96]:
#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 [97]:
#4 output=[] for x in sol: output.append((Vmax*x^n)/(h+x^n)*x)
In [98]:
#5 list_plot(list(zip(time, output)), plotjoined=True) #oscillations occur at a very high rate, and stabilize after a few minutes
In [90]:
#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 [86]:
#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 [88]:
#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 [112]:
#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 [113]:
#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 [115]:
#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 [116]:
#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 [117]:
#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 [118]:
#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 [119]:
#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 [120]:
#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 [121]:
#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 [ ]: