In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.spatial
import deplete2d
%matplotlib inline
For the well mixed model, we can just assume that the prey density is uniform throughout the entire environment. So the rate of prey encounters is just the rate of area swept by the predator per unit time. So if the radius of the predator is $r$ and his speed is $v$, then he covers $2rv\,dt$ area per unit time, in which there are $2rvu\, dt$ critters, where $u$ is the prey density. So the well-mixed equation is \begin{equation} \frac{du}{dt} = -2rvu, \end{equation} which has solution: \begin{equation} u(t) = u_0 \exp(-2rv t) \end{equation} Also useful later: time when well-mixed model as dropped to proportion $\epsilon$ of original value: \begin{equation} \frac{\log(\epsilon)}{-2rv} \end{equation}
In [2]:
def wellmixed(t, density=1000, radius=0.01, speed=1.):
return density * np.exp(-2*radius*speed*t)
deplete2d.py contains a class deplete2d which allows simulation of 2D static prey with a single moving predator. deplete2d specifies a 1x1 domain with periodic boundary conditions and a Poisson-distributed number of prey of specified, with a single predator situated at the centre.
The predator has a perceptive radius radius and a speed speed. The object has a method step that allows steps with no compoment larger than 0.5 to be taken. Each step is assumed to be at constant speed, and each prey that is within the predator's perceptive radius at some point during that step gets eaten (removed). The time at which it is eaten is defined to be the time that the predator passes closest to the prey point (it doesn't divert from its path to eat the prey). Each consumption event is logged, and get_output returns all output events.
In [45]:
# Test whether a step works
sim = deplete2d.deplete2d(density=10000, radius=0.02)
steps = np.array([[0., 0.3], [0., -0.15], [0.15, 0.], [0., 0.15], [0., -0.3],
[0.1, 0.], [0., 0.15], [0.03, 0.03], [-0.03, 0.03], [-0.03, -0.03], [0.03, -0.03]])
for step in steps:
sim.step(step)
# take a snapshot of the prey at time 1000 (enough time to have taken all steps)
prey = sim.get_prey(1000)
x = prey['position']
plt.scatter(x[:,0], x[:,1])
plt.axis('equal')
_ = plt.axis([0, 1, 0, 1])
In [46]:
t, pop = sim.get_output()
plt.plot(t, pop)
plt.xlabel('t')
plt.ylabel('population')
Out[46]:
There are two parameters that need to be fiddled if we hold the density fixed. One is the radius of the predator and one is the step length of the random walk. The way that the Brownian motion has been implemented is that x and y compoments for a step are drawn from Gaussian distributions with standard deviation $\sigma / \sqrt{2}$
Note that the timescales on the three graphs are clearly much different.
In [53]:
from sympy import *
init_printing()
Work out the mean interpoint spacing for fixed $\mu$. I'm sure this is in a book somewhere, but let's just figure it out. Assuming homogeneous Poisson process in 2D with rate $\mu$, probability of there being exactly one point in $[r, r+dr]$ and none in the disc of radius $r$ around another point is
\begin{align} P(\text{closest pt in r, r + dr}) &= P(\text{no pt closer than r}) P(\text{one pt in r, r + dr}) \\ &= \exp(-\mu \pi r^2) 2\pi r\mu\,dr \exp(-2\pi r\mu\,dr) \\ &= \exp(-\mu \pi r^2) 2\pi r \mu \, dr + O(dr^2) \end{align}So our nearest neighour distribution is \begin{equation} f(r) = 2 \pi r \mu \exp(-\mu \pi r^2) \end{equation} Let's first check this thing integrates to one
In [97]:
r = Symbol('r')
mu = Symbol('mu', positive=True)
f = 2*pi*r*mu *exp(-mu*pi*r**2)
integrate(f, (r, 0, oo))
Out[97]:
Now compute the mean interpoint spacing
In [87]:
integrate(r*f, (r, 0, oo))
Out[87]:
So if our density is 10000, then our mean interpoint spacing is 1/200 = 0.005. So we should look at values of $r$ and $\sigma$ relative to this
In [98]:
density = 1000.
radii = [0.1, 1, 10.] * 1 / (2*np.sqrt(density))
sigmas = radii.copy()
for sigma in sigmas:
for r in radii:
params = {'density': density, 'radius': r, 'speed': 1}
sim = deplete2d.deplete2d(**params)
deplete2d.brownianwalk(sim, T=3*log(0.001)/(-2*params['radius']*params['speed']), sigma=sigma)
[t, pop] = sim.get_output()
plt.figure()
plt.plot(t, pop, t, wellmixed(t, **params))
plt.legend(('random walk', 'well-mixed'))
plt.title('$\mu$ = %g, r = %.3g, $\sigma$ = %.3g'% (params['density'], r, sigma))
plt.xlabel('t')
plt.ylabel('population')
Next up: same as above, but with averages over lots of realisations so we can see what's going on in the large radius case. Need to interpolate output onto fixed time grid to do this