Séance 2 : Étudier les interactions écologiques à l'aide d'équations différentielles¶
# Commençons par importer les outils dont nous aurons besoin.
from functools import partial
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
%matplotlib inline
plt.rc('text', usetex=True)
Introduction¶
Il existe toute une classe de modèles permettant de rendre compte de l'évolution d'un certain nombre de grandeurs dans le temps et qui sont assez simple à écrire et à simuler : les systèmes d'équations différentielles du premier ordre (ODE).
- Ces modèles sont déterministes. Ils font partie de la famille des systèmes dynamiques
- Ces modèles sont en temps continu, les trajectoires sont définies sur un-sous ensemble de $\mathbb R$.
- Ces modèles sont non spatialisés, on parle de modèle de champ moyen. L'espace d'état est de la forme $\mathbb R^N$, les grandeurs ne sont pas entières parce que l'on regarde la densité moyenne (en individus par unité de volume).
- Ils peuvent être autonomes ou non. Nous allons commencer par un système autonome.
Dans le cas d'un modèle ODE du premier ordre autonome, on dispose de plusieurs outils pour représenter le comportement du modèle. Ils sont de plus en plus synthétiques (mais aussi de plus en plus couteux à établir). Dans l'ordre :
- Le tracé d'une trajectoire donne l'état du système au cours du temps pour certaines valeurs de paramètres et certaines conditions initiales. Pour le tracer il suffit d'intégrer l'ODE.
- Le tracé du diagramme de phase donne le comportement qualitatif du système pour certaines valeurs de paramètres, mais toutes conditions initiales. Pour le tracer, il faut chercher les isoclines zéro et les points d'équilibre.
- Le tracé du diagramme de bifurcation donne la position et la nature des équilibre du système pour (potentiellement) toutes valeurs de paramètres, et toutes conditions initiales.
Partie I : Croissance exponentielle - systèmes dynamiques linéaires en une dimension¶
Exercice 1 :¶
1) Résoudre analytiquement le système.
2) Simuler les trajectoire du système (S1) en utilisant sa solution analytique pour les paramètres et conditions initiales suivants :
a) $N_0 = 100$, $r = 0.5$
b) $N_0 = 100$, $r = -0.5$
(Indice : utiliser np.linspace et plt.plot).
Correction Exercice 1 - question 1
$$(S1) : \frac{dN}{dt} = rN $$Le système (S1) est linéaire, et on peut calculer sa solution analytique qvec $N_0$ la condition initiale :
$$N(t)=N_0 e^{r t}$$# Correction Exercice 1 - question 2
dt=0.1
T=5
tspace = np.linspace(0,T,int(T/dt))
N0=100
r=0.5
plt.plot(tspace, N0*np.exp(r*tspace))
N0=100
r=-0.5
plt.plot(tspace, N0*np.exp(r*tspace))
plt.legend(["r=0.5","r=-0.5"])
#plt.savefig('syst_lineaire_1D.pdf')
Intégration numérique¶
Exercice 2 :¶
1) Écrire une fonction euler_explicite(y0, r, dt, T) qui simule la trajectoire du système (S1) en utilisant l'aire du rectangle rouge.
2) Tracer la trajectoire pour les paramètres suivants :
a) $N_0 = 100$, $r = 0.5$, $T = 10$, $dt = 1$
b) $N_0 = 100$, $r = -0.5$, $T = 10$, $dt = 1$
c) $N_0 = 100$, $r = -0.5$, $T = 50$, $dt = 5$
Comparer avec la solution analytique.
3) Comparer l'écart entre les trajectoires estimées et la solution analytique pour $dt\in [0.001, 0.01,0.1,1,2]$.
# Correction Exercice 2 - Question 1
def euler_explicite(N0, r, dt, T):
'''Résolution numérique de l'ODE N' = r y sur [0,T].
En utilisant la méthode d'Euler explicite.
Args:
N0 (float): Condition initiale.
r (float): Taux de croissance
dt (float): Incrément de temps
T (float): Temps final.
Retour:
(iterable): Trajectoire du système
'''
N = np.zeros(int(T/dt))
N[0] = N0
for n in range(int(T/dt)-1):
N[n+1] = N[n] * (1+dt*r) # N[n] + le rectangle rouge.
return N
# Correction Exercice 2 - Question 2
# N0 r dt T
conditions = [(100, 0.5, 1, 10),
(100,-0.5, 1, 10),
(100,-0.5, 5, 50)]
fig, axes = plt.subplots(1,len(conditions), figsize=(len(conditions)*6,4))
for ax,(N0,r,dt,T) in zip(axes,conditions):
ax.set_title(r"$N_0={}, r={}, dt={}, T={}$".format(N0, r, dt, T))
tspace = np.linspace(0,T,int(T/dt))
ax.plot(tspace, euler_explicite(N0,r,dt,T))
ax.plot(tspace, N0*np.exp(r*tspace), color='k', label='Solution analytique')
ax.legend()
#plt.savefig('euler_explicte.pdf')
# Correction Exercice 2 - Question 3
dtspace = [0.001, 0.01,0.1,1,2]
T = 10
r = 1
N0 = 10
fig,ax = plt.subplots(1,2,figsize=(10,5))
for dt in dtspace:
tspace = np.linspace(0,T,int(T/dt))
N= euler_explicite(N0,r,dt,T)
ax[0].plot(tspace, N, label='dt={}'.format(dt))
ax[0].scatter(tspace, N, marker='.')
tspace = np.linspace(0,T,1000)
ax[0].plot(tspace,N0*np.exp(r*tspace), color='k', label='Solution analytique')
ax[0].legend()
ax[0].set(xlabel='Temps', ylabel='Densité de la population (N)')
# Plot l'erreur
for dt in dtspace:
tspace = np.linspace(0,T,int(T/dt))
y = euler_explicite(N0,r,dt,T)
ax[1].scatter(dt, abs(N0*np.exp(r*T)-y[-1]), marker='.', label='dt={}'.format(dt))
tspace = np.linspace(0,T,1000)
#plt.savefig('euler_explicte_dt.pdf')
Exercice 3 :¶
1) Ecrire une fonction euler_implicite(y0, r, dt, T) qui simule la trajectoire du système (S1)
en utilisant l'aire du rectangle vert.
2) Tracez la trajectoire pour les paramètres suivants :
a) y0 = 100, r = 0.5, T = 10, dt = 1
b) y0 = 100, r = -0.5, T = 10, dt = 1
c) y0 = 100, r = -0.5, T = 50, dt = 5
3) Comparer les résultats avec la méthode précédente.
# Correction Exercice 3 - Question 1
def euler_implicite(N0, r, dt, T):
'''Résolution numérique de l'ODE N' = r N sur [0,T].
En utilisant la méthode d'Euler implicite.
Args:
N0 (float): Condition initiale.
r (float): Taux de croissance
dt (float): Incrément de temps
T (float): Temps final.
Retour:
(iterable): Trajectoire du système
'''
N = np.zeros(T//dt)
N[0] = N0
for n in range(T//dt-1):
N[n+1] = N[n] * 1/(1-dt*r) #N[n] + le rectangle vert.
return N
# Correction Exercice 3 - Question 2
# N0 r dt T
conditions = [(100, 0.5, 1, 10),
(100,-0.5, 1, 10),
(100,-0.5, 5, 50)]
fig, axes = plt.subplots(1,len(conditions), figsize=(len(conditions)*6,4))
for ax,(N0,r,dt,T) in zip(axes,conditions):
ax.set_title(r"$N_0={}, r={}, dt={}, T={}$".format(N0, r, dt, T))
tspace = np.linspace(0,T,int(T/dt))
ax.plot(tspace, euler_implicite(N0,r,dt,T))
ax.plot(tspace, N0*np.exp(r*tspace), color='k', label='Solution analytique')
ax.legend()
#plt.savefig('euler_implicite.pdf')
# Correction Exercice 3 - Question 3
N0 = 100
T = 50
r = -.5
fig, axes = plt.subplots(1,2,figsize=(15,6))
for ax, label, integrator in ((axes[0],'implicite',euler_explicite), (axes[1], 'explicite',euler_implicite)):
for dt in [1,5]:
ax.plot(np.arange(0,T,dt),
integrator(N0,r,dt,T),
label='$dt = {}$'.format(dt))
ax.set(title = 'Euler {}, $r = {}$'.format(label,r))
ax.legend()
Partie II : Croissance logistique - systèmes dynamiques non-linéaires en une dimension¶
Exercice 4 :¶
1) Écrire la nouvelle forme du modèle, où $r$ représente le taux de naissance per capita et $K$ la capacité biologique du système (i.e. le nombre maximum d'individu que le système supporte à l'équilibre). Calculer analytiquement ses équilibres et déduire leur stabilité.
2) Simuler une trajectoire de ce modèle et représenter là sous forme de graphe pour $r=1$ et $K=10$. (Indice : utiliser partial et scipy.integrate.odeint).
3) Tester pour plusieurs conditions initiales.
Correction Exercice 4 - Question 1
Le taux de naissance est linéaire $b(N) = rN$ (c'est à dire constant per capita) et le taux de mort est $m(N) = r\frac{N^2}{K}$ (c'est à dire linéaire per capita). Ce qui donne:
$$(S2) : \frac{dN}{dt} = rN (1 - \frac{N}{K}) $$Les équilibres de (S2) correspondent à $\frac{dN}{dt}=rN (1 - \frac{N}{K})=0$, c'est à dire $N=0$ ou $N=K$.
Pour déterminer la stabilité des équilibres, on doit regarder le signe de la pente de la tangente à la courbe $N'=f(N)$, i.e le signe de la dérivée de $N'$ par rapport à $N$ calculée aux points d'équilibres :
$$\frac{dN'}{dN}= r - 2r\frac{N}{K}= r(1 - 2\frac{N}{K})$$Pour $N^*=0$, $\frac{dN'}{dN}(0)=r>0$, l'équilibre $N^*=0$ est instable.
Pour $N^*=K$, $\frac{dN'}{dN}(K)=-r<0$, l'équilibre $N^*=K$ est stable.
# Correction Exercice 4 - Question 2
def logistique(N, t, r, K):
"""Dérivée temporelle du modèle de croissance logistique.
Args:
N (float): Population de proies
t (float): Temps (non utilisé)
r (float): Taux d'accroissement.
K (float): Capacité biotique
Return: dx/dt
"""
return r * N * (1 - N/K)
# On crée une liste des temps pour lequel on va intégrer le système.
temps = np.linspace(0,20,1500)
# On utilise functools.partial pour définir l'application partielle de la fonction
# logistique où les paramètres sont fixés.
K = 10
r = 1
derivative = partial(logistique, r=r, K=K)
# On définit les conditions initiales.
conditions_initiales = 1
# Avec ces trois éléments on peut utiliser le solveur de scipy.
# Celui ci renvoie un np.array de taille len(time) par len(condition_initiale)
trajectoire = scipy.integrate.odeint(derivative,
y0=conditions_initiales,
t=temps)
# On affiche cette trajectoire à l'aide de matplotlib.
fig,ax = plt.subplots(1,1,figsize=(10,5))
ax.plot(temps,trajectoire)
ax.set(xlabel='Temps',ylabel='Densité de population', label='Trajectoire')
plt.show()
# Correction Exercice 4 - Question 3
# Différentes condtions initiales.
# Remarquez que les trajectoires ne se coupent jamais.
fig,ax = plt.subplots(1,1,figsize=(10,5))
for ci in (0,1,2,10,20):
trajectoire = scipy.integrate.odeint(derivative,
y0=ci,
t=temps)
ax.plot(temps,trajectoire)
ax.set(xlabel='Temps',ylabel='Densité de population', title='Trajectoires')
#plt.show()
#plt.savefig('syst_non_lineaire_1D.pdf')
# Correction Exercice 5
fig,ax = plt.subplots(1,1,figsize=(6,6))
Nspace = np.linspace(0,12)
plt.hlines(0,0,12,linestyles='--', alpha=.5, label="Espace d'état")
plt.plot(Nspace, derivative(Nspace,0), label='dN/dt = F(N) Valeur du flot')
# Placer les points d'équilibre
plt.scatter([K],[0], label='Équilibre stable', color='k')
plt.scatter([0],[0], label='Équilibre instable', color='w', edgecolors='k')
# Placer le flot
Ngrid = np.arange(12)
plt.quiver(Ngrid, [0]*len(Ngrid), derivative(Ngrid,0), [0]*len(Ngrid), label='Flot', color='r')
# Configuration du diagramme
ax.set(xlabel='Densité de population N', ylabel='dN/dt', title='Diagramme de phase (en 1D)')
ax.legend();
#plt.savefig('syst_non_lineaire_1D_phase.pdf')
### Correction Exercice 6
fig,ax = plt.subplots(1,1,figsize=(8,4))
plt.plot([0,10],[0,10], ls='-', color='k', label='Stable')
plt.hlines(0, 0,10, linestyles='--', label='Instable')
plt.plot([0,-10],[0,-10], ls='-', color='k')
plt.hlines(0, 0,-10, linestyles='--')
ax.set(ylabel="Densité de population à l'équilibre N*",
xlabel='Paramètre de contrôle K',
title='Diagramme de bifurcation sur K, r=1')
ax.legend();
#plt.savefig('syst_lineaire_1D_bifurcation.pdf')
Partie III : Modèles de Lotka-Volterra - systèmes dynamiques non-linéaires en deux dimensions¶
On ajoute une seconde espèce (sous forme d'une seconde dimension).
Stabilité des systèmes linéaires¶
Exercice 7 :¶
Étudier et représenter la stabilité des systèmes dynamiques linéaires suivants :
\begin{equation} A = \left [ \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right ] ,\ \ B = \left [ \begin{array}{cc} -1 & 1 \\ 1 & -1 \\ \end{array} \right ] ,\ \ C = \left [ \begin{array}{cc} -1 & 0 \\ 0 & 1 \\ \end{array} \right ], \ \ D = \left [ \begin{array}{cc} 0 & -1 \\ 1 & 0 \\ \end{array} \right ] ,\ \ E = \left [ \begin{array}{cc} -1 & -1 \\ 1 & 0 \\ \end{array} \right ] ,\ \ F = \left [ \begin{array}{cc} 1 & -1 \\ 1 & 0 \\ \end{array} \right ] \end{equation}# Correction Exercice 7
equilibria = {'Noeud instable': np.array([[1,0],[0,1]]),
'Noeud stable': np.array([[-1,0],[0,-1]]),
'Point selle': np.array([[-1,0],[0,1]]),
'Centre': np.array([[0,-1],[1,0]]),
'Spirale stable': np.array([[-1,-1],[1,0]]),
'Spirale instable': np.array([[1,-1],[1,0]]),
}
plt.rc('text', usetex=True)
fig, axes = plt.subplots(2,3,figsize=(15,10))
for ax,(k,A) in zip(axes.flat, equilibria.items()):
#### Trace le flot en 2D ####
# np.meshgrid permet de
X,Y = np.meshgrid(np.linspace(-1,1,10), np.linspace(-1,1,10))
# np.vectorize permet de donner à une fonction quelquonque (ici le lambda)
# des capacités de broadcasting (application à tous les éléments).
# Par défaut vectorize attend une fonction qui prend des scalaires et retourne des scalaires.
# Si ce n'est pas le cas il faut préciser la signature. Ici la fonction prend deux scalaires "(),()" et
# retourne un array de taille 2 "->(2)"
f = np.vectorize(lambda x,y: A@[x,y], signature='(),()->(2)')
U = f(X,Y)
ax.quiver(X, Y, U[:,:,0], U[:,:,1], color='grey') # champs de vecteur
ax.streamplot(X, Y, U[:,:,0], U[:,:,1], density=0.5) # trace le flot
# Calcule les valeurs propres et les vecteurs propres de A.
valp, vecp = np.linalg.eig(A)
# Teste la stabilité de l'équilibre.
stable = all(np.real(valp)<0)
# Représente l'équilibre.
ax.scatter(0, 0,
color='k' if stable else 'w',
edgecolor = 'k' if not stable else 'w',
label='Equilibrium')
# Légende en latex.
# Nécessite d'avoir plt.rc('text', usetex=True)
tex = (r'$A = \left [ \begin{array}{cc}'
+ r"{0[0][0]:2} & {0[0][1]:2} \\ {0[1][0]:2} & {0[1][1]:2}".format(A)
+ r' \end{array} \right ]\\'
+ r'\lambda_0 = {0[0]:.1f} \\ \lambda_1 = {0[1]:.1f} $'.format(valp))
ax.legend(title=tex, loc='upper right')
ax.set(title=k)
plt.tight_layout()
A) Modèle de Lotka-Volterra de compétition :¶
Exercice 8 :¶
1) Quelle est la nouvelle forme du système notée (S3) ?
2) Tracer les trajectoires comme des séries temporelles du modèle de LV pour les trois jeux de paramètes et les deux conditions initiales:
Jeux de paramètes :
$r_1=1$, $r_2=3$, $a_{1\to1}=1$, $a_{2\to1}=1.2$, $a_{1\to2}=1.2$ et $a_{2\to2}=1.5$
$r_1=1$, $r_2=3$, $a_{1\to1}=1$, $a_{2\to1}=1.2$, $a_{1\to2}=1.5$ et $a_{2\to2}=1$
$r_1=1$, $r_2=3$, $a_{1\to1}=1.2$, $a_{2\to1}=1$, $a_{1\to2}=1$ et $a_{2\to2}=1.5$
Conditions initiales : $(N_1(0)=0.3,N_2(0)=0.5)$ et $(N_1(0)=0.7,N_2(0)=0.2)$.
3) Tracer ces même trajectoires comme des courbes paramétrées dans l'espace des phases du modèle de LV. Décrire l'évolution dans chacun des cas.
Correction Exercice 8 - Question 1
Première étape, on change l'espace d'état à $\mathbb N^2$
$$(S3) = \begin{cases} \frac{dN_1}{dt} = r_1N_1 (1 - \frac{N_1}{K_1})\\ \frac{dN_2}{dt} = r_2N_2(1 - \frac{N_2}{K_2})\\ \end{cases}$$Deuxième étape, on ajoute le terme de compétition:
$$(S3) = \begin{cases} \frac{dN_1}{dt} = r_1N_1 (1 - \frac{N_1}{K_1} - a_{2\to1}N_2)\\ \frac{dN_2}{dt} = r_2N_2(1 - \frac{N_2}{K_2} - a_{1\to2}N_1)\\ \end{cases}$$En général on note $a_{1\to1} = K_1^{-1}$ et $a_{2\to2} = K_2^{-1}$, ainsi:
$$(S3) = \begin{cases} \frac{dN_1}{dt} = r_1N_1 (1 - a_{1\to1} N_1 - a_{2\to1}N_2)\\ \frac{dN_2}{dt} = r_2N_2(1 - a_{2\to2} N_2 - a_{1\to2}N_1)\\ \end{cases}$$# Correction Exercice 8 - Question 2
# Utilisez les jeux de paramètres suivants:
p1 = {'r':np.array([1,3], dtype=float), 'a':np.array([[1,1.2],[1.2,1.5]], dtype=float)} # fortes compétitions intra et interspécifique pour l'espèce 2
p2 = {'r':np.array([1,3], dtype=float), 'a':np.array([[1,2.1],[1.5,1]], dtype=float)} # fortes compétitions interspécifiques
p3 = {'r':np.array([1,3], dtype=float), 'a':np.array([[1.2,1],[1,1.5]], dtype=float)} # fortes compétitions intraspécifiques
plist = (p1,p2,p3)
# Conditions initiales:
ci_1 = [.3,.5]
ci_2 = [.7,.2]
ci_list = [ci_1, ci_2]
# Temps
tspace = np.linspace(0,25,100)
def lotka_volterra(y,t,r,a):
'''Flot du système de Lotka-Volterra
Args:
y (np.array): Densités de populations (vecteur),
t (float): Temps (non utilisé),
r (np.array): Taux de croissance,
a (np.array): Matrice d'interaction,
Retourne le flot du système de LV au point y.
'''
return r * y * (1-a@y)
# On peut afficher cette trajectoire en utilisant matplotlib.
def trace_trajectoire(trajectoire, temps, labels=("Espèce 1", "Espèce 2"),
title='Trajectoire du Modèle de Lotka-Volterra', ax=None):
"""Trace la trajectoire du système sur ax.
Args:
trajectoire (np.array)
temps (np,array)
labels (Iterable)
title (str)
ax (matplotlib.axis)
"""
# Si aucun axe n'est précisé on tracera sur l'axe "en cours".
if ax is None:
ax = plt.gca()
# Trace les lignes
ax.plot(temps, trajectoire[:,0], label=labels[0])
ax.plot(temps, trajectoire[:,1], label=labels[1])
# Ajoute une légende correspondant aux labels des lignes.
ax.legend()
# Ajoute un titre aux axes.
ax.set_xlabel("Temps")
ax.set_ylabel("Variables d'état")
ax.set_title(title)
# Calculer les trajectoires...
trajectoires = []
for p in plist:
flot = partial(lotka_volterra, **p)
for ci in ci_list:
tr = scipy.integrate.odeint(flot, ci, tspace)
trajectoires.append((tr,p, ci))
# Tracer les trajectoires au cours du temps...
fig, axes = plt.subplots(3,2,figsize=(15,20))
for ax,(tr,p, ci) in zip(axes.flat, trajectoires):
trace_trajectoire(tr, tspace,
ax=ax,
title="ci={}, r={}, a={}{}".format(ci, p['r'],p['a'][0],p['a'][1]))
# Correction Exercice 8 - Question 3
# Tracer les trajectoires dans l'espace des phases
fig, ax = plt.subplots(1,1,figsize=(4,4))
plt.grid()
for (tr,p, ci) in trajectoires:
l = ax.plot(tr[:,0],tr[:,1])
ax.scatter(tr[0,0], tr[0,1], color=l[0].get_color(), marker='o')
ax.scatter(tr[-1,0], tr[-1,1], color=l[0].get_color(), marker='x')
ax.set(xlabel='Espèce 1',
title="Espace d'état".format(ci, p['r'],p['a'][0],p['a'][1]),
ylabel='Espèce 2')
plt.tight_layout()
Exercice 9 :¶
1) Trouver les isoclines zéro et les équilibres du modèle de Lotka-Volterra competitif (S3).
2) Tracer le diagramme de phase pour les trois jeux de paramètres utilisés exercice 7. (Indice: utiliser les fonctions plt.plot, plt.scatter et plt.streamplot).
3) Etudier numériquement la stabilité des équilibres. En déduire les conditions de coexistance de deux espèces.
Correction Exercice 9 - Question 1 :
Les isoclines zéro de $N_1$ sont les droites : $N_1=0$ et $N_2 = \displaystyle{\frac{1-a_{11}N_1}{a_{21}}}$.
Les isoclines zéro de $N_2$ sont les droites : $N_2=0$ et $N_2=\displaystyle{\frac{1-a_{12}N_1}{a_{22}}}$.
Les équilibres sont donc :
- $N_1^*=0$ et $N_2^*=0$;
- $N_1^*=0$ et $N_2^*=\displaystyle{\frac{1}{a_{22}}}$;
- $N_2^*=0$ et $N_1^*=\displaystyle{\frac{1}{a_{11}}}$;
- $N_1^*=\displaystyle{\frac{a_{22}-a_{21}}{a_{22}a_{11}-a_{21}a_{12}}}$ et $N_2^*=\displaystyle{\frac{a_{11}-a_{12}}{a_{22}a_{11}-a_{21}a_{12}}}$.
Le dernier est "biologiquement" possible uniquement si $N_1^*>0$ et $N_2^*>0$.
La stabilité des équilibres peut être déterminée graphiquement en analysant les signes des dérivées temporelles de part et d'autre des isoclines.
# Correction Exercice 9 - Question 2 :
def plot_template_lv(param, ax=None, quiver=True, color_eq=False):
'''Plot the phase space of a LV system, with isoclines, equilibria and stability and
derivative vector field'''
if ax is None:
plt.figure(figsize=(15,15))
ax = plt.gca()
# Champ de vecteurs.
if quiver:
xmax = 1 if not (param['a']==0).any() else 10
x = np.linspace(0,xmax,25)
y = np.linspace(0,xmax,25)
X,Y = np.meshgrid(x,y, indexing='ij')
flot = lambda x,y:lotka_volterra(np.array([x,y]), 0 ,**param)
f = np.vectorize(flot, signature='(),()->(2)')
U = f(X,Y)
norm = np.sqrt(U[:,:,0]**2+U[:,:,1]**2)
ax.quiver(X, Y, U[:,:,0]/norm, U[:,:,1]/norm ,color='k', units='xy', angles='xy', alpha=.5)
if not (param['a']==0).any():
# Isoclines triviales
mx = np.max(1/param['a'])
ax.plot([0,mx],[0,0], color='r', ls='--',alpha=.5)
ax.plot([0,0],[mx,0], color='b', ls='--',alpha=.5)
# Isoclines non triviales
ax.plot([0,1/param['a'][0,0]], [1/param['a'][0,1],0], color='b', ls='--',alpha=.5)
ax.plot([0,1/param['a'][1,0]],[1/param['a'][1,1],0], color='r', ls='--',alpha=.5)
ax.text(1/param['a'][0,0], 0, r'$\frac{1}{a_{11}}$', color='b',
horizontalalignment='center', verticalalignment='top')
ax.text(0, 1/param['a'][0,1], r'$\frac{1}{a_{21}}$', color='b',
horizontalalignment='right', verticalalignment='center')
ax.text(1/param['a'][1,0], 0, r'$\frac{1}{a_{12}}$', color='r',
horizontalalignment='center', verticalalignment='top')
ax.text(0, 1/param['a'][1,1], r'$\frac{1}{a_{22}}$', color='r',
horizontalalignment='right', verticalalignment='center')
# Première bissectrice
ax.plot([0,1],[0,1], color='k', ls='-.', alpha=.1)
# Plot the equilibria.
det = np.linalg.det(param['a'])
internal = np.array([(param['a'][1,1]-param['a'][0,1])/det, (param['a'][0,0]-param['a'][1,0])/det])
equilibria = [[0,0], [0,1/param['a'][1,1]], [1/param['a'][0,0],0]]
stability = [False, param['a'][0,1]>param['a'][1,1], param['a'][1,0]>param['a'][0,0] ]
if all(internal>0):
equilibria.append(internal)
stability.append(not stability[-1])
ax.set_xlabel('N1')
ax.set_ylabel('N2')
colors = ('blue','orange','green')
stable_eq = [e for e,s in zip(equilibria, stability) if s]
for e,s in zip(equilibria, stability):
if not s:
ax.scatter(*e, color='w', edgecolor='k')
for i,e in enumerate(stable_eq):
if color_eq:
ax.scatter(*e, color=colors[i])
else:
ax.scatter(*e, color='k')
fig, axes = plt.subplots(1,3,figsize=(21,7))
for ax,p in zip(axes,plist):
plot_template_lv(p,ax=ax)
#plt.savefig('syst_lineaire_2D_competition.pdf')
# Les isoclines de N1 sont représentées en bleus, tandis que celles de N2 sont rouges.
Correction Exercice 9 - Question 3 :
Les 2 espèces coexistent et forment un équilibre stable uniquement si : $N_1^*=\displaystyle{\frac{a_{22}-a_{21}}{a_{22}a_{11}-a_{21}a_{12}}}>0$, $N_2^*=\displaystyle{\frac{a_{11}-a_{12}}{a_{22}a_{11}-a_{21}a_{12}}}>0$, $a_{22}-a_{21}>0$ et $a_{11}-a_{12}>0$, c'est à dire, lorsque la compétition intraspécifique est plus importante que la compétition interspécifique.
B) Modèle de Lotka-Volterra proie-prédateur :¶
Exercice 10 :¶
Soit le modèle proie (N) - prédateur (P) suivant : $$(S4) = \begin{cases} \displaystyle{\frac{dN}{dt} }= bN - aNP \\ \\ \displaystyle{\frac{dP}{dt} }= cNP - dP \\ \end{cases}$$
1) Déterminer analytiquement les équilibres et leur stabilité.
2) Simuler ce modèle, tracer ses trajectoires comme des séries temporelles, puis comme des courbes paramétrées dans l'espace d'état, pour les paramètres et conditions initiales suivantes : $(a=0.01, b=0.02, c=0.01, d=0.02)$ et $(a=0.01, b=0.01, c=0.01, d=0.06)$ et $(N_0=3.6,P_0=2.1)$, $(N_0=1.8,P_0=4)$.
3) Ajouter les équilibres, les isoclines et le flot à l'espace d'état pour former le diagramme de phase.
Correction Exercice 10 - Question 1
Les isoclines de $N$ sont les droites $N=0$ et $P=\displaystyle{\frac{b}{a}}$.
Les isoclines de $N$ sont les droites $P=0$ et $N=\displaystyle{\frac{d}{c}}$.
Les équilibres sont donc $(N^*,P^*)=(0,0)$ et $(N^*,P^*)=(\displaystyle{\frac{d}{c}},\displaystyle{\frac{b}{a}})$.
Pour déterminer la stabilité de ces équilibres on étudie la jacobienne du flot aux points d'équilibre :
$$J_f(N,P)=\left [ \begin{array}{cc} b-aP & -aN \\ cP & cN-d \\ \end{array} \right ]$$En $(0,0)$ on a : $J_f(0,0)=\left [ \begin{array}{cc} b & 0 \\ 0 & -d \\ \end{array} \right ]$ dont les valeurs propres $b$ et $-d$ sont réelles et de signes opposés : l'équilibre $(0,0)$ est donc un point selle.
En $(\displaystyle{\frac{d}{c}},\displaystyle{\frac{b}{a}})$ on a : $J_f(\displaystyle{\frac{d}{c}},\displaystyle{\frac{b}{a}})=\left [ \begin{array}{cc} 0 & \displaystyle{\frac{-ad}{c}} \\ \displaystyle{\frac{cb}{a}} & 0 \\ \end{array} \right ]$ dont les valeurs propres $i\sqrt{bd}$ et $-i\sqrt{bd}$ sont imaginaires pures : l'équilibre $(\displaystyle{\frac{d}{c}},\displaystyle{\frac{b}{a}})$ est donc un centre.
# Correction Exercice 10 - Question 2
params = [{"a":.01, "b":.02, "c":.01, "d":.02},
{"a":.01, "b":.01, "c":.01, "d":.06}]
initial_conditions = [[3.6,2.1], [1.8,4]]
temps = np.linspace(0,1000,1500)
def proie_predateur(b,d,a,c):
"""Convertit les paramètres du modèle proie-prédateur en paramètre du modèle général.
Args:
b: Taux de reproduction per capita des proies.
d: Taux de décès per capita des prédateurs
a: Taux de prédation des prédateurs sur les proies
c: Taux de conversion de biomasse de proies en biomasse de prédateurs"""
return {'r': np.array([b,-d]),
'a':np.array([[0,a/b],[-c/-d,0]])}
def trace_diagramme_de_phase_lv(a,b,c,d, xmax=None, ymax=None, quiver=True, ax=None):
"""Trace le squelette du diagramme de phase de Lotka-Volterra.
Args:
a,b,c,d (floats): Paramètres du modèle.
xmax, ymax (floats): Limites du graphe.
quiver (bool): Si vrai, affiche le champ de vecteur
ax (matplotlib axis): Axe sur lequel ajouter le diagramme.
"""
# Configuration du plot.
if ax is None:
fig = plt.figure(figsize=(5,5))
ax = plt.gca()
xmax = 2*d/c if xmax is None else xmax
ymax = 2*b/a if ymax is None else ymax
# Isoclines zéros.
ax.hlines([0, b/a], 0, xmax, color=['b','r'], linestyles='--', alpha=.4)
ax.vlines([0, d/c], 0, ymax, color=['r','b'], linestyles='--', alpha=.4)
# Points d'équilibre.
ax.scatter(d/c,b/a, color='w', edgecolor='k',)
ax.scatter(0,0, color='w', edgecolor='k',)
# Text
ax.text(d/c, 0, r'$\frac{d}{c}$', color='b',
horizontalalignment='center', verticalalignment='top')
ax.text(0, b/a, r'$\frac{b}{a}$', color='r',
horizontalalignment='right', verticalalignment='center')
# Axis labels
ax.set_xlabel('Proies')
ax.set_ylabel('Prédateurs')
ax.set_title('Diagramme de Phase')
# Champ de vecteur
if quiver:
x = np.linspace(0,xmax,25)
y = np.linspace(0,ymax,25)
X,Y = np.meshgrid(x,y, indexing='ij')
flot = lambda x,y:lotka_volterra(np.array([x,y]), 0 ,**proie_predateur(a=a, b=b, c=c, d=d))
f = np.vectorize(flot, signature='(),()->(2)')
U = f(X,Y)
norm = np.sqrt(U[:,:,0]**2+U[:,:,1]**2)
ax.quiver(X, Y, U[:,:,0]/norm, U[:,:,1]/norm ,color='k', units='xy', angles='xy', alpha=.5)
return ax
# Correction Exercice 10 - Question 3
for p in params:
plt.figure(figsize=(15,5))
ax_phase = plt.subplot(1,3,3)
for i, ic in enumerate(initial_conditions):
ax = plt.subplot(1,3,i+1)
trajectoire = scipy.integrate.odeint(partial(lotka_volterra, **proie_predateur(**p)),
y0=ic,
t=temps)
label = 'Trajectoire {}'.format(i+1)
trace_trajectoire(trajectoire,temps, ax=ax)
ax.set_title(label)
ax_phase.plot(*trajectoire.transpose(), color='C{}'.format(i+2), label=label)
trace_diagramme_de_phase_lv(quiver=True,ax=ax_phase,
xmax=trajectoire[:,0].max(),
ymax=trajectoire[:,1].max(), **p)
ax_phase.legend()
plt.show()
#plt.savefig('syst_lineaire_2D_competition_PP.pdf')
C) Variante du modèle de Lotka-Volterra proie-prédateur :¶
Exercice 11 :¶
Soit le modèle suivant :
\begin{equation} (S5) = \begin{cases} \frac{dN}{dt} = bN(1-\frac{N}{K}) - aNP \\ \frac{dP}{dt} = cNP - dP \\ \end{cases} \end{equation}1) Simuler quelques trajectoires. Quelles sont les différences avec le modèle (S4)?
2) Dessiner le diagramme de phase. Où sont les isoclines-zéro ? Où sont les équilibres ? Quelle est leur nature ?
3) Faire varier $c$. Faire un diagramme de bifurcation.
# Correction Exercice 11 - Question 1
params = [{"a":.01, "b":.02, "c":.01, "d":.02, "K":10},
{"a":.01, "b":.02, "c":.02, "d":.02, "K":10}]
initial_conditions = [[1,1],
[1.8,4]]
temps = np.linspace(0,3000,1500)
def proie_predateur_mod(b,d,a,c,K):
"""Convertit les paramètres du modèle proie-prédateur modifié en paramètre du modèle général.
Args:
b: Taux de reproduction per capita des proies.
d: Taux de décès per capita des prédateurs
a: Taux de prédation des prédateurs sur les proies
c: Taux de conversion de biomasse de proies en biomasse de prédateurs
K: capacité biotique des proies"""
return {'r': np.array([b,-d]),
'a':np.array([[1/K,a/b],
[-c/-d,0]])}
# Correction Exercice 11 - Question 2
def trace_diagramme_de_phase_lv_mod(a,b,c,d,K, xmax=None, ymax=None, quiver=True, ax=None):
"""Trace le squelette du diagramme de phase de Lotka-Volterra.
Args:
a,b,c,d (floats): Paramètres du modèle.
xmax, ymax (floats): Limites du graphe.
quiver (bool): Si vrai, affiche le champ de vecteur
ax (matplotlib axis): Axe sur lequel ajouter le diagramme.
"""
# Configuration du plot.
if ax is None:
fig = plt.figure(figsize=(5,5))
ax = plt.gca()
xmax = 4 if xmax is None else xmax
ymax = 4 if ymax is None else ymax
# Isoclines zéros.
ax.hlines(0, 0, xmax, color=['b','r'], linestyles='--', alpha=.4)
yspace = np.linspace(0,ymax)[1:]
ax.plot(K*(1-(a/b)*yspace), yspace,color='r',ls='--', alpha=.4)
ax.vlines([0, d/c], 0, ymax, color=['r','b'], linestyles='--', alpha=.4)
# Points d'équilibre.
ax.scatter(d/c,(b/a)*(1-(d/(c*K))),
color='k' if c>d/K else 'w',
edgecolor='w' if c>d/K else 'k',
s=100)
ax.scatter(0,0,
color='w' if c>d/K else 'k',
edgecolor='k' if c>d/K else 'w',
s=100)
# Text
ax.text(d/c, 0, r'$\frac{d}{c}$', color='b',
horizontalalignment='center', verticalalignment='top')
ax.text(0, b/a, r'$\frac{b}{a}$', color='r',
horizontalalignment='right', verticalalignment='center')
# Axis labels
ax.set_xlabel('Proies')
ax.set_ylabel('Prédateurs')
ax.set_title('Diagramme de Phase')
# Champ de vecteur
if quiver:
x = np.linspace(0,xmax,25)
y = np.linspace(0,ymax,25)
X,Y = np.meshgrid(x,y, indexing='ij')
flot = lambda x,y:lotka_volterra(np.array([x,y]), 0 ,**proie_predateur_mod(a=a, b=b, c=c, d=d, K=K))
f = np.vectorize(flot, signature='(),()->(2)')
U = f(X,Y)
norm = np.sqrt(U[:,:,0]**2+U[:,:,1]**2)
ax.quiver(X, Y, U[:,:,0]/norm, U[:,:,1]/norm ,color='k', units='xy', angles='xy', alpha=.5)
ax.set(xlim=(-0.1,xmax),ylim=(-0.1,ymax))
return ax
for p in params:
fig, axes = plt.subplots(1,3,figsize=(15,5))
for i, (ax,ic) in enumerate(zip(axes[:2],initial_conditions)):
trajectoire = scipy.integrate.odeint(partial(lotka_volterra, **proie_predateur_mod(**p)),
y0=ic,
t=temps)
label = 'Trajectoire {}'.format(i+1)
trace_trajectoire(trajectoire,temps, ax=ax)
ax.set_title(label)
axes[2].plot(*trajectoire.transpose(),
color='C{}'.format(i+2),
label=label, zorder=-99)
trace_diagramme_de_phase_lv_mod(quiver=True, ax=axes[2], **p)
axes[2].legend()
plt.show()
#plt.savefig('syst_lineaire_2D_competition_modPP.pdf')
# Correction Exercice 11 - Question 3
def trace_bifurcation_lvm(a,b,d,K):
low = np.linspace(0,d/K,100)+1e-5
high = np.linspace(d/K,2*d/K,100)+1e-5
position_eq = lambda c: (b/a)*(1-(d/(c*K)))
plt.plot(low,[0]*len(low),ls='-',color='k')
plt.plot(high,[0]*len(high),ls='--',color='k')
plt.plot(low,position_eq(low),ls='--',color='k')
plt.plot(high,position_eq(high),ls='-',color='k')
plt.ylim(-10,10)
plt.xlabel('c')
plt.ylabel('Prédateurs')
plt.title('Diagrame de bifurcation')
trace_bifurcation_lvm(a=.01, b=0.2, d=0.02, K=10)
