In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline
plt.rcParams.update({'font.size': 14})

In [2]:
def predpreyeq(X, t, g=1):
x = X[0]; y = X[1];
dxdt = x - x*y
dydt = -g*y + x*y
return dxdt, dydt

def midpoint(tmax=30,h=0.01,x0=1,y0=1,g=1):
n = int(tmax/h)
ts = np.arange(n+1)*h
xs = np.zeros(n+1)
ys = np.zeros(n+1)
xs[0] = x0; ys[0] = y0; # set initial condition
for k in range(n):
dxdt, dydt = predpreyeq([xs[k], ys[k]], ts[k], g)
xhalf = xs[k] + (h/2)*dxdt
yhalf = ys[k] + (h/2)*dydt
thalf = ts[k] + h/2
dxdt, dydt = predpreyeq([xhalf, yhalf], thalf, g)
xs[k+1] = xs[k] + h*dxdt
ys[k+1] = ys[k] + h*dydt
return ts, xs, ys

g0 = 2.5
ts, xs, ys = midpoint(g=g0)
plt.plot(ts, xs, "-", label="Rabbits", linewidth=1.5)
plt.plot(ts, ys, "-", label="Foxes", linewidth=1.5)
plt.xlabel("Time")
plt.ylabel("Population")
plt.yticks(np.arange(0,6))
plt.legend();
plt.savefig("plot1a.pdf")

In [3]:
ts=np.linspace(0,20,400)
xis = np.linspace(0.5, g0, 9)
print(xis)
for xi in xis:
ts, xs, ys = midpoint(x0=xi, g=g0)
plt.plot(xs, ys, "-")
plt.plot([g0], [1], ".")  # the fixed point
plt.ylim(0,5)
plt.xlabel("Rabbits")
plt.ylabel("Foxes")
plt.savefig("plot1b.pdf")

[ 0.5   0.75  1.    1.25  1.5   1.75  2.    2.25  2.5 ]

In [4]:
def findperiod(ts,xs,ys, g=1):
"""This function estimates the period of the oscillations in the rabbit and fox populations.
I look for zeros of dx/dt, corresponding to maxima and minima of the functions
Rather than look for dx/dt=0 directly, I test for a change in the sign of dx/dt,
and then use linear interpolation to improve the time estimate.
The function returns the average of the periods found.
"""
Ts = []
for i in range(1,len(ts)):
pt0 = (xs[i-1],ys[i-1])
pt1 = (xs[i],ys[i])
dx0 = predpreyeq(pt0, ts[i-1], g)[0]
dx1 = predpreyeq(pt1, ts[i], g)[0]
if dx0*dx1 < 0:
# dx/dt has changed sign; use linear interpolation to estimate the 't' value.
t0, t1 = ts[i-1], ts[i]
dt = t1-t0
ddx = dx1-dx0
T = t0 - dx0*dt/ddx
Ts.append(T)

if len(Ts) >= 3:
if len(Ts) % 2 == 1:  # if an odd number of maxima and minima are found, do this ...
T = 2*(Ts[-1] - Ts[0]) / (len(Ts)-1)
else:
T = 2*(Ts[-2] - Ts[0]) / (len(Ts)-2)
else:
print("Can't estimate the period: more data needed. ", xs[0])
T = 0
return T

In [5]:
g0 = 1
npts=79
x0s = np.linspace(0.03,2.0,npts)
periods0 = np.zeros(len(x0s))
for i in range(len(x0s)):
ts, xs, ys = midpoint(x0=x0s[i], g=g0, tmax=60)
periods0[i] = findperiod(ts, xs, ys)

plt.xlabel("Initial condition $x_0$")
plt.ylabel("Period $T$")
plt.ylim(0,10)
plt.plot(x0s, periods0, '-', label="g = 1");
plt.plot([1], [2*np.pi], 'rs');
plt.legend()
plt.tight_layout()
plt.savefig("plot2b_1.pdf")

In [6]:
g0 = 2.5
npts=79
x0s = np.linspace(0.05,4.0,npts)
periods0 = np.zeros(len(x0s))
for i in range(len(x0s)):
ts, xs, ys = midpoint(x0=x0s[i], g=g0, tmax=60)
periods0[i] = findperiod(ts, xs, ys)

plt.xlabel("Initial condition $x_0$")
plt.ylabel("Period $T$")
plt.ylim(0,8)
plt.plot(x0s, periods0, '-', label='g = 2.5');
plt.plot([g0], [2*np.pi/np.sqrt(g0)], 'rs');
plt.legend()
plt.tight_layout()
plt.savefig("plot2b_2.pdf")

In [7]:
# this code can be used for both part 3(a) Harvesting and 3(b) Disease.
def predpreyeq(X, t, g=1, alpha=0, beta=0, tau=5):
x = X[0]; y = X[1];
dxdt =  x - x*y - alpha - beta*x*np.sin(t/tau)**2
dydt = -g*y + x*y
return dxdt, dydt

def midpoint(tmax=30,h=0.01,x0=1,y0=1,g=1,alpha=0,beta=0,tau=5):
n = int(tmax/h)
ts = np.arange(n+1)*h
xs = np.zeros(n+1)
ys = np.zeros(n+1)
xs[0] = x0; ys[0] = y0; # set initial condition
epsilon = 0.001
for k in range(n):
if xs[k] < epsilon or ys[k] < epsilon:
break
dxdt, dydt = predpreyeq([xs[k], ys[k]], ts[k], g, alpha, beta, tau)
xhalf = xs[k] + (h/2)*dxdt
yhalf = ys[k] + (h/2)*dydt
thalf = ts[k] + h/2
dxdt, dydt = predpreyeq([xhalf, yhalf], thalf, g, alpha, beta, tau)
xs[k+1] = xs[k] + h*dxdt
ys[k+1] = ys[k] + h*dydt
return ts, xs, ys, k

In [8]:
# In part 3, you could use either odeint or the midpoint method. Here I compare the output from both methods.
from scipy.integrate import odeint
n = 1000
ts = np.linspace(0,30,n)
# start near the fixed point, which is now a spiral source.
g0 = 1
alp = 0.1
x0, y0 = g0 +0.01, 1 - alp/g0 +0.01
Xs = odeint(predpreyeq, [x0,y0], ts, args=(1,0.1,0,5))
plt.xlabel("Rabbits"); plt.ylabel("Foxes")
plt.plot(Xs[:,0], Xs[:,1], 'r.')
ts,xs,ys,k = midpoint(alpha=0.1, x0=x0, y0=y0)
plt.plot(xs, ys, 'b-');

In [9]:
# An example plot, showing extinction.
alpha0=0.1
g0=1
ts, xs, ys, kmax = midpoint(alpha=alpha0, g=g0, tmax=60)
plt.figure(figsize=(6,6))
plt.xlabel('time')
plt.ylabel('population')
plt.plot(ts[:kmax], xs[:kmax], 'r-', label='rabbits')
plt.plot(ts[:kmax], ys[:kmax], 'b-', label='foxes')
plt.legend(loc=0);
plt.plot([0,55],[0,0],'-')
plt.savefig('plot3a_1.pdf')

In [10]:
# Same data, showing rabbits vs foxes.
plt.figure(figsize=(6,6))
plt.xlabel('rabbits')
plt.ylabel('foxes')
plt.plot(xs[:kmax], ys[:kmax], 'r-')
plt.plot([g0], [1-alpha0/g0], ".")  # the fixed point
plt.savefig('plot3a_3.pdf')

In [11]:
# Investigate survival time as a function of alpha
tmax = 200
npts = 150
alphas = np.linspace(0.05, 0.5, npts)
tsurvival = np.zeros(npts)
for i in range(npts):
ts, xs, ys, kmax = midpoint(alpha=alphas[i], x0=1.5, y0 = 1, tmax=tmax)
tsurvival[i] = ts[kmax]

plt.figure(figsize=(6,6))
plt.xlabel('$\\alpha$')
plt.ylabel('survival time')
plt.plot(alphas, tsurvival, 'r-');
plt.savefig('plot3a_2.pdf')

In [12]:
# This plot indicates that the survival time appears to diverge like 1/alpha for small alpha.
plt.plot(np.log(alphas), np.log(tsurvival), 'r-')
val=0.04
xs=np.linspace(val,1, 10)
ys=2*xs**(-1)
plt.plot(np.log(xs), np.log(ys), '-');

In [13]:
# An example plot with disease.
tau0 = 5
g0 = 1
print(g0)
ts, xs, ys, kmax = midpoint(beta=0.2,  g=g0, tau=tau0, tmax=400)
plt.xlabel('time')
plt.ylabel('population')
plt.plot(ts[:kmax], xs[:kmax], 'r-', label='rabbits')
plt.plot(ts[:kmax], ys[:kmax], 'b-', label='foxes')
plt.legend(loc=0);
plt.title("(a) $\\beta = 0.2, g = 1, \\tau = 5$")
plt.savefig('plot3b_1.pdf')

1

In [14]:
plt.figure(figsize=(5,5))
plt.plot(xs[:kmax], ys[:kmax], 'r-');
plt.xlabel('rabbits');
plt.ylabel('foxes');
plt.title("(b)")
plt.tight_layout()
plt.savefig('plot3b_2.pdf')

In [15]:
# Try the case where the period of the system without disease is the same as tau
tau0 = 5
g0 = (2*np.pi/5)**2
print(g0)
ts, xs, ys, kmax = midpoint(beta=0.2, x0=g0, g=g0, tau=tau0, tmax=400)
plt.xlabel('time')
plt.ylabel('population')
plt.plot(ts[:kmax], xs[:kmax], 'r-', label='rabbits')
plt.plot(ts[:kmax], ys[:kmax], 'b-', label='foxes')
plt.legend(loc=0);
plt.title("(c) $\\beta = 0.2, g = (2\\pi/\\tau)^2, \\tau = 5$")
plt.savefig('plot3b_3.pdf')

1.5791367041742972

In [16]:
from mpl_toolkits.mplot3d import Axes3D
tau = 5
zs = np.sin(ts/tau)**2
fig = plt.figure(figsize=(6,6))
ax.set_zlabel("$\\sin^2(t/\\tau)$");