CoCalc Public FilesLS 30B W18 NOW / CO2-Model-Interactive-Demo.sagewsOpen with one click!
Authors: Samuel Aragones, Kuan-Ting Chen, Abraham Chorbajian, Timothy Hamilton, Kristin McCully, Mrinal Rath, Feifan Tang, Han Wang, Jeong Yoon
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[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]
# 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 "]))
Interact: please open in CoCalc