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 = fig.add_subplot(111, projection='3d')
ts, xs, ys, kmax = midpoint(beta=0.2, x0=g0, g=g0, tau=tau0, tmax=400)
ax.plot(xs,ys, np.sin(ts/tau)**2, 'r-')
ts, xs, ys, kmax = midpoint(beta=0.2, x0=g0, g=g0, tau=tau0, tmax=20)
ax.plot(xs,ys, np.sin(ts/tau)**2, 'b-')
plt.xlabel("rabbits"); plt.ylabel("foxes"); 
plt.title("(d)")
ax.set_zlabel("$\\sin^2(t/\\tau)$");
#plt.plot(xs, ys, zs)
plt.savefig("plot3b_4.pdf")