

In [1]:

    
from scipy.integrate import odeint

Solving 1st order ODEs

The logistic equation is a first-order non-linear differential equation that describes the evolution of a population as a function of the population size at a give generation. It can be written as:

${\displaystyle \tau {dp(t) \over dt} = p(t)(k-p(t))},$

where $p(t)$ y the population size at generation $t$, $k$ is the maximal size of the population, and $\tau$ a proportionality factor

Define the equation



In [2]:

    
def diff(p, generation):
    """
    Returns the as size of the population as a function of the generation
    defined in the following differential equation:
    
    dp/dg = p*(k-p)/tau,
    
    where p is the population size, g is the generation index, k is 
    the maximal population size (fixed to 1000) and tau a constant that describes the 
    number of individuals per generation (fixed to 1e4)
    
    p          -- (int) population size, dependent variable 
    generation -- (int) generation index, independent variable
    """
    
    tau = 1e4    # rate of individuals per generation
    pMax = 1000     # max population size
    return (p * (pMax-p) ) / tau



In [3]:

    
# define the independent variable (i.e., generations)
g = np.arange(200) # 200 generations

Numerical solution to the differential equation

It requieres the initial condition and the independent variable



In [4]:

    
# solve the differential equation 
p = odeint(diff, 2, g) # initial conditions is 2 individuals



In [5]:

    
# plot the solution
plt.plot(g,p);
plt.ylabel('Population size');
plt.xlabel('Generation');

Introducing optional arguments to the differential equation



In [6]:

    
def diff2(p, generation, tau, pMax):
    """
    Returns the as size of the population as a function of the generation
    defined in the following differential equation:
    
    dp/dg = p*(pMax-p/tau,
    
    where p is the population size, g is the generation index, pMax is 
    the maximal population size and tau a constant that describes the 
    number of individuals per generation
    
    p          -- (int) population size 
    generation -- (int) generation index
    pMax       -- (int) maximal number of individuals in a population
    tau        -- (float) rate of individuals per generation
    """
    return (p*(pMax-p))/tau



In [7]:

    
y = odeint(diff2, 2, g, args=(1e4, 1000)) # initial conditions is 2 individuals, tau = 1e4, pMax = 1000



In [8]:

    
# plot the solution
plt.plot(g,y);
plt.ylabel('Population size');
plt.xlabel('Generation');