from IPython.html.widgets import interact
from scipy.integrate import odeint
from sympy import symbols,Symbol,pprint
import os.path
if not os.path.exists("sde.py"):
import urllib
print urllib.urlretrieve(r"https://raw.githubusercontent.com/cchuang2009/PySDE/master/sde.py", "sde.py")
from sde import SDE_solver, Milstein, Euler
The deterministic logistic ordinary differential equation is: $$ \frac{dN}{dt} = r \cdot N \cdot \Big(1- \frac{N}{K}\Big) $$ where $N$ is the population size. $r$ is the density-independent growth rate, $K$ is the carrying capacity (maximum population size) and $t$ is time.
def ode(N,t,r=1.,K=100.):
return r * N * (1 - N/K)
def runode(r,K,T,N0=1):
t = linspace(0, T)
N = N0 + odeint(ode, 1, t, (r,K))
return t,N
@interact
def f(r=1.,K=100.,T=10.):
t,N = runode(r,K,T)
plot(t,N,'k-')
xlabel('t')
ylabel('N')
ylim(-0.1*max(N),max(N)*1.1)
return t,N
This ODE has a solution: $$ N(t) = \frac{K}{1 + \frac{K}{N(0)} \cdot e^{-r t}} $$
def logistic(t,r,K,N0):
Q = K/N0 - 1
return K / (1 + Q * exp(-r * t))
@interact
def f(r=1.,K=100.,N0=1,T=10.):
t = linspace(0, T)
N = logistic(t, r, K, N0)
plot(t, N, 'k-')
xlabel('t')
ylabel('N')
ylim(-0.1*max(N),max(N)*1.1)
return t,N
The stochastic logistic function is (Campillo et al. 2013): $$ dN = r \cdot N \cdot \Big(1 - \frac{N}{K}\Big)N \cdot dt + \rho \sqrt{\Big(\lambda + \mu + r \cdot \frac{N}{K}\Big)N} \cdot dB_t = \\ $$ where $\lambda$ is the birth rate, $\mu$ is the death rate ($r = \lambda - \mu$), $\rho$ is the level of stochasticity, and $B_t$ is a standard Brownian motion.
We use an SDE solver by cchuang2009.
def runsde(lam=1.5,mu=0.5,K=100.,rho=0.1,N0=1.,T=20.):
x,dx = symbols('x dx')
r = lam - mu
drift = r * x * (1 - x/K)
diffusion = rho * ((lam + mu + r * x/K) * x)**0.5
t = linspace(0, T)
_,N = Milstein(drift,diffusion,N0,0,T,len(t)-1)
return t,N
@interact
def f(lam=0.5,mu=0.1,K=100.,rho=0.5,N0=1,T=20.):
reps = 5
for _ in range(reps):
t,N_sde = runsde(lam,mu,K,rho,N0,T)
plot(t, N_sde, 'k-')
N = logistic(t,lam-mu,K,N0)
plot(t, N, "b--", lw=3)
ylim(-0.1*max(N),max(N)*1.2)
xlabel('t')
ylabel('N')
return T,N_sde,N
I'm Yoav Ram.
This notebook is released under a CC-BY-SA 3.0 license.
The notebook can be found online at http://ipython.yoavram.com.