2018-04-19-000849.sagews

#Brenden Gregorio

%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.")

    # 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
    if delays:
        timestepcorrectionfactor = ceil(timestep / min(delays) * 2)
    else:
        timestepcorrectionfactor = 1
    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]

#FE 4.2
#7

var("N", "Ntau")
r=1.2
k=50
d=0.1
tau=2.8
Nprime=r*Ntau*(1-(Ntau/50))-d*N
sol=dde_solve(Nprime, N, {Ntau:2.8}, history=20, tmax=100, timestep=0.1)
t=srange(0,100,0.1)
list_plot(zip(t,sol), axes_labels=["Time", "N"], plotjoined=True)

(N, Ntau)

#Changing the tau value or r changes the oscillatory behavior of the system.

#FE 4.3
#6

var("S", "P")
Sprime=0.5-(0.23*S*P*P)
Pprime=(0.23*S*P*P)-0.4*P
sol=desolve_odeint([Sprime, Pprime], ics=[4,5], dvars=[S,P], times=t)
p=list_plot(zip(t, sol[:,0]))+list_plot(zip(t, sol[:,1]), plotjoined=True, color="black")
show(p)

(S, P)