 CoCalc Public FilesLS 30B W18 NOW / CO2-Model-Interactive-Demo.sagews
Views : 92
Compute Environment: Ubuntu 18.04 (Deprecated)
%auto
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]

# Interactive for Mackey-Glass Model
@interact
def respiration (tau = (0.1,0.2,0.001), n = (0.1, 10, 0.01)):
var("x", "x_tau")
L = 6
Vmax = 80
a = 0.2
h = 400
max_time = 200
tstep = 0.05
xprime = L - (Vmax * x_tau^n)/(h + x_tau^n) * a * x
delay = {x_tau: tau}
sol = dde_solve(xprime,x,delay, history = 0, tmax = max_time, timestep = tstep)
time = srange(0,max_time,tstep)
show(list_plot(zip(time, sol), plotjoined = true, axes_labels = ["Time (minutes)","Concentration of Carbon Dioxide "]))