Shared2018-04-19-000849.sagewsOpen in CoCalc
#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)