Proies et prédateurs

Plus il y a de lapins, plus les renards sont bien nourris, donc plus ils se reproduisent et mangent des lapins et finalement moins il y a de lapins. Moins il y a de lapin, plus les renards ont faim et moins ils se reproduisent et finalement plus il y a de lapins.

Cette jolie animation montre ce cycle.

Les équations de Lotka-Volterra permette de simuler l'évolution de populations de proies et de prédateurs. Elles sont :

$$\frac{dx}{dt} = x(\alpha - \beta \; y)$$

$$\frac{dy}{dt} = - y(\gamma - \delta \; x)$$

Où $x$ représente le nombre de proies et $y$ celui des prédateurs. On note que :

• $\alpha$ représente la capacité de reprodution des proies
• $\beta$ représente l'appétit des prédateurs
• $\gamma$ représente le taux de morts de faim chez les prédateurs
• $\delta$ représente le taux de reproduction des prédateurs bien nourris

a = 3   # on choisit les paramètres et on les stocke en variables globales
b = 2
c = 3
d = 0.3

Exercice

• Tracer la courbe qui montre l'evolution du nombre de proies et de prédateurs dans le temps
• Tracer la courbe paramétrique proies/prédateur
• Faire varier le ratio proies / prédateurs et regarder l'impact sur la courbe paramétrique
• Trouver la valeur d'équilibre qui fait que les deux populations restent constates

Solution

import scipy.integrate as sci
import numpy as np
import matplotlib.pylab as plt

%matplotlib inline

def f(Y,t,a,b,c,d):
    return Y[0]*(a-b*Y[1]), -Y[1]*(c-d*Y[0]) # attention à la porté des variables

t = np.linspace(0,10,1000)
res = sci.odeint(f, (10,4), t, args=(a,b,c,d)) # 10000  lapins et 4000 renards
plt.figure(figsize=(13,3), dpi=100)
plt.plot(t,res)
plt.grid()

On voit qu'il s'agit de courbe périodiques ce qui implique que dans un repères nombre de proies / nombre de prédateur, on repasse toujours par les mêmes points et donc qu'on tourne en rond.

Pour tracer ce cycle, regardons res afin de comprendre sa structure :

res

array([[10.        ,  4.        ],
       [ 9.51200644,  3.99704519],
       [ 9.04888295,  3.98839291],
       ...,
       [ 2.37116055,  1.30247454],
       [ 2.38126123,  1.27299365],
       [ 2.39279978,  1.24422046]])

On dessine donc la seconde colonne en fonction de la première soit res[:,1] en fonction de res[:,0] :

t = plt.linspace(0,3,100)           # pas trop de cycles et pas pour voir les accélérations
for init_pred in range(1,7):   # on fait varier le ratio proie/predateur initial
    res = sci.odeint(f, (10,init_pred), t, args=(a,b,c,d)) 
    plt.plot(res[:,0],res[:,1],'.',label="%2.1f"%(10/init_pred)) # '.' demande de mettre un point par valeur. 
plt.xlabel(u'proies')
plt.ylabel(u'prédateurs')
plt.legend()

<matplotlib.legend.Legend at 0x7ff0bc81e610>

Surprise La courbe 5.0 est incluse dans la courbe 10.0 alors qu'ensuite chaque courbe est incluse dans celle de ratio plus petit !

Zoomons en choisissant des ratios proie/prédateur dans la zone d'inversion :

t = np.linspace(0,2.5,100)   # temps pas assez long pour avoir des cycles complets
for init_pred in np.arange(0.5,3,0.5):   # arange pour avoir des réels
    res = sci.odeint(f, (10,init_pred), t, args=(a,b,c,d)) 
    plt.plot(res[:,0],res[:,1],'.',label="%2.1f"% (10/init_pred)) 
plt.xlabel(u'proies')
plt.ylabel(u'prédateurs')
plt.legend()
plt.grid()  # on commence les cycles en x = 10

Le cycle de ration 6.7 semble un points ! Un point d'équilibre.

Optimisons

Calculons la ratio qui mène au cercle le plus petit c.a.d à l'équilibre naturel. On note qu'il s'agit du cycle le plus court, on va donc cherche à minimiser le cycle et pour cela on va utiliser minimize de scipy.optimize.

Pour commencer définissons la fonction à optimiser :

def ecart(ratio):

    t = np.linspace(0,9,1000)
    res = sci.odeint(f, (10., 10./ratio), t, args=(a,b,c,d))
    return max(res[:,0]) - min(res[:,0])

import scipy.optimize as sco

mini = sco.minimize(ecart,[3],method='TNC',bounds=[(1,10)])
print("Minimum pour le ratio = %f" % mini.x[0])
t = np.arange(0,5,0.01)
res = sci.odeint(f, (10., 10./mini.x[0]), t, args=(a,b,c,d))
plt.plot(res[:,0],res[:,1]);

Minimum pour le ratio = 6.666669

La ration de 20/3 donne la position d'équilibre. On voit sur la figure qui suit que dans ce cas les populations restent constantes.

plt.plot(t,res);