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()