Détermination de paramètres pour un modèle d'évolution de population



In [1]:

    
%matplotlib inline
pylab.rcParams['figure.figsize'] = (10.0, 8.0)

from numpy import *
from scipy.optimize import root

La population et ses dérivées en fonction des paramètres

$$p_{p_0, m, n} : t \mapsto \frac{m p_0 e^{mt}}{n p_0 e^{m t} + m - n p_0}$$



In [2]:

    
def population(p0, m, n, t):
    return m*p0 * exp(m*t) / (n*p0*exp(m*t) + m - n*p0)

def ddm_population(p0, m, n, t):
    """Dérivée partielle par rapport à m
    """
    return p0*exp(m*t)*(n*p0*(exp(m*t)-1)+t*m*(m-n*p0))/(n*p0*exp(m*t)+m-n*p0)**2
    
def ddn_population(p0, m, n, t):
    """Dérivée partielle par rapport à n
    """
    return m*p0*p0*exp(m*t)*(exp(m*t)-1)/(n*p0*exp(m*t) + m - n*p0)**2

Le système d'équations à résoudre



In [3]:

    
def solveur(p0, t1, p1, t2, p2, facteur=1, method="hybr"):
    
    def residu(x):
        """Fonction donc on cherche le zéro
        """
        m, n = x[0], x[1]
        return array([
                population(p0, m, n, t1) - p1,
                population(p0, m, n, t2) - p2
            ])

    def jacobien(x):
        m, n = x[0], x[1]
        return array([
                [ddm_population(p0, m, n, t1), ddm_population(p0, m, n, t2)],
                [ddn_population(p0, m, n, t1), ddn_population(p0, m, n, t2)]
            ])
    
    def solveur_parametre(x0):
        resultat = root(residu, x0, jac=jacobien, method=method)
        return resultat
    
    return solveur_parametre

Résolution

Définition des paramètres



In [4]:

    
pop_0  = 29981336 # Population initiale
pop_7  = 31359340 # Population après 7 ans
pop_13 = 32402940 # Population après 13 ans

Numériquement, le problème n'est pas très bien posé. On utilise donc le changement de paramètres suivant $$ p_{p_0, m, n} = \frac{1}{k}p_{k\, p_0, m, \frac{n}{k}}$$

On choisit un point initial et on résout



In [7]:

    
facteur = 1e-7

x0 = array([1e-5, 1e-12/facteur]) # Point initial

s = solveur(pop_0*facteur, 7, pop_7*facteur, 13, pop_13*facteur, method="hybr")

resultat = s(x0)

m_sol = resultat.x[0]
n_sol = resultat.x[1]*facteur

print "m : {0}\nn : {1}\npopulation limite : {2}".format(m_sol, n_sol, m_sol/n_sol)









    



m : 0.0309178139908
n : 7.98381577349e-10
population limite : 38725610.4949

Graphe de l'évolution de population avec les paramètres solution



In [8]:

    
population_parametree = vectorize(lambda t:population(pop_0, m_sol, n_sol, t))
t = range(300)
plot(t, population_parametree(t))
plot([0, 7, 13], [pop_0, pop_7, pop_13], "o")
xlabel(u"Années")
ylabel("Population")
axhline(m_sol/n_sol)









    Out[8]:





<matplotlib.lines.Line2D at 0x7f6bda4313d0>



In [ ]: