CO2-Model-Interactive-Demo.sagews

%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