Commit a27311f8 authored by Ludovic Mailleret's avatar Ludovic Mailleret
Browse files

update

parent 6a7b4c1c
......@@ -417,7 +417,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
......@@ -431,7 +431,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
"version": "3.8.12"
}
},
"nbformat": 4,
......
%% Cell type:markdown id: tags:
%% Cell type:markdown id:e0a4436e tags:
 
# Notebook d'aide/correction pour la séance EPU MAM4 Biomaths (suite 3)
 
*séance du 21/03/2022*
 
*Ludovic Mailleret, Mars 2022*
 
## Populations en interaction : le modèle proies-prédateurs de Lotka et Volterra
 
### Dynamiques et plan de phase
 
Nous considérons le modèle de dynamique de populations (classique) de Lotka et Volterra
 
$$
\left\{\begin{array}{l}
\dot x = rx - c xy\\
\dot y = bxy - m y
\end{array}\right.
$$
 
Il n'y a pas de difficulté particulière à la simulation par rapport au modèle de la tordeuxe du bourgeon de l'épinette avec population d'oiseaux variables.
Il n'y a pas de difficulté particulière à la simulation par rapport au modèle de la tordeuse du bourgeon de l'épinette avec population d'oiseaux variables.
 
%% Cell type:code id: tags:
%% Cell type:code id:44329d14 tags:
 
``` python
# on nettoie l'espace de travail et on reload les modules
%reset -f
 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
 
 
# script de simulation
# densités initiales des populations
x0 = 1
y0 = 2.5
 
# encapsulation de la densité initiale
etat0_LV = np.array([x0, y0])
 
# définition d'un vecteur tspan pour tracer E_varie()
t_0 = 0 # temps initial
t_fin = 30.0 # temps final, ici on prend un temps final long pour observer les dynamiques lentes
pas_t = 0.01 # pas de temps de récupération des variables entre t_0 et t_fin
 
# définition du tspan
tspan = np.arange(t_0, t_fin, pas_t)
 
# paramètres du modèle
r = 1.0
c = 1.0
b = 1.0
m = 1.0
 
# encapsulation des paramètres dans un array
params_LV = np.array([r, c, b, m])
 
 
# définition du modèle de Lotka Volterra
def modele_LV(etat, t, params):
x, y = etat # recupere les variables d'etat
r, c, b, m = params # recupere les parametres
etatdot = [r*x - c*x*y, # dot x
b*x*y - m*y] # dot y
return etatdot # renvoie la derivee
 
 
# simulation proprement dite
int_LV = odeint(modele_LV, etat0_LV, tspan, args=(params_LV,), hmax=pas_t)
```
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:90bbe5eb tags:
 
A priori, tout est OK:
 
%% Cell type:code id: tags:
%% Cell type:code id:e19e8e6c tags:
 
``` python
int_LV
```
 
%% Output
 
array([[1. , 2.5 ],
[0.98511256, 2.4998135 ],
[0.97045034, 2.4992576 ],
...,
[0.30945175, 0.61472634],
[0.31065286, 0.61049964],
[0.31187176, 0.60630933]])
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:c93c1956 tags:
 
Nous pouvons passer à la représentation graphique, d'abord en fonction du temps:
 
%% Cell type:code id: tags:
%% Cell type:code id:ab1ab902 tags:
 
``` python
# création d'une figure, et d'un système d'axe
fig1, ax1 = plt.subplots(1, 1, figsize=(14, 8))
 
# titre de la figure
fig1.suptitle("Dynamiques Proies - Prédateurs\n modèle de Lotka Volterra", va='top', fontsize='18')
 
# tracé de x et y contre le temps
ax1.plot(tspan, int_LV[:, 0], color = 'C0', label = "proies $x$")
ax1.plot(tspan, int_LV[:, 1], color = 'C1', label = "prédateurs $y$")
 
# labellisation des axes
ax1.set_xlabel('temps', fontsize='16')
ax1.set_ylabel('densités de populations', fontsize='16')
 
# légende
ax1.legend(fontsize='14')
 
# modification éventuelle des bornes des axes
ax1.set_ylim(bottom=None, top=None)
 
# ajout d'une grille
ax1.grid()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:6c3349b4 tags:
 
Il est intéressant aussi de représenter les trajectoires dans l'espace d'état $(x,y)$, en combinaison avec une analyse qualitative du plan de phase (isoclines nulles, équilibres), une représentation du champs de vecteurs et de quelques morceaux de trajectoires sur le plan.
 
Commençons par calculer les isoclines nulles et les équilibres.
 
%% Cell type:code id: tags:
%% Cell type:code id:12a59ede tags:
 
``` python
# array annexes pour le calcul et la représentation des isoclines nulles
xplot = np.arange(0, 3, .1)
yplot = np.arange(0, 3, .1)
 
# isoclines nulles de xdot
null_x_x = np.zeros_like(yplot) # x = 0 isocline nulle de xdot
null_x_y = np.ones_like(xplot)*(r/c) # y = r/c isocline nulle de xdot
 
# isoclines nulles de ydot
null_y_y = np.zeros_like(xplot) # y = 0 isocline nulle de ydot
null_y_x = np.ones_like(yplot)*(m/b) # x = m/b isocline nulle de ydot
 
# équilibres
eq_extinct = [0, 0]
eq_coex = [r/c, m/b]
```
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:5a23929c tags:
 
Puis on commence à tracer le plan de phase:
 
%% Cell type:code id: tags:
%% Cell type:code id:3ac182f1 tags:
 
``` python
# création d'une figure, et d'un système d'axe
fig2, ax2 = plt.subplots(1, 1, figsize=(14, 8))
 
# titre de la figure
fig2.suptitle("Dynamiques Proies - Prédateurs\n modèle de Lotka Volterra", va='top', fontsize='18')
 
# tracé des isoclines nulles, des équilibres, et de la trajectoire simulée
# isoclines nulles de x
ax2.plot(null_x_x, yplot, color = 'C2')
ax2.plot(xplot, null_x_y, color = 'C2', label = "isoclines nulles de $\dot x$")
# isoclines nulles de y
ax2.plot(xplot, null_y_y, color = 'C1')
ax2.plot(null_y_x, yplot, color = 'C1', label = "isoclines nulles de $\dot y$")
 
# équilibres
ax2.plot(eq_extinct[0], eq_extinct[1], marker ='.', color = 'C3', markersize = 16)
ax2.plot(eq_coex[0], eq_coex[1], marker ='.', color = 'C0', markersize = 16)
# trajectoires
ax2.plot(int_LV[:, 0], int_LV[:, 1], color = 'C0', linewidth = 3, label = "trajectoire")
 
# labellisation des axes
ax2.set_xlabel('Proies $x$', fontsize='16')
ax2.set_ylabel('Prédateurs $y$', fontsize='16')
 
# légende
ax2.legend(fontsize='14', loc = "upper right")
 
# modification éventuelle des bornes des axes
ax2.set_ylim(bottom=-.1, top=None)
ax2.set_xlim(left=-.1, right=None);
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:3ca475d4 tags:
 
Il est intéressant (et joli) de surimposer sur cette figure une représentation du champs de vecteur et quelques idées de trajectoire.
 
Pour cela il nous faut définir une grille de valeurs dans le plan $(x,y)$ et utiliser la méthode `quiver()` des systèmes d'axes (champs de vecteurs) et la méthode `streamplot()` (échantillons de trajectoires).
 
La définition de la grille se fait par la fonction `meshgrid()` de `numpy` sur la base d'`array` définissant l'amplitude et le pas des variations de $x$ et $y$.
 
%% Cell type:code id: tags:
%% Cell type:code id:dacfee71 tags:
 
``` python
# définition de l'échantillonnage selon $x$ et $y$
x_grid = np.linspace(0.1, 3.0, 10) # au passage on change un peu de np.arange()
y_grid = np.linspace(0.1, 3.0, 10)
 
# grille X,Y selon x_grid et y_grid
X, Y = np.meshgrid(x_grid, y_grid)
```
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:caf0d0be tags:
 
On calcule ensuite les valeurs des dérivées $\dot x$ et $\dot y$ sur cette grille:
 
%% Cell type:code id: tags:
%% Cell type:code id:8e31bb68 tags:
 
``` python
# dérivées dot_x et dot_y sur la grille
dx, dy = modele_LV([X, Y], 0, params_LV)
```
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:56cac5d5 tags:
 
A partir de là on peut tracer le champs de vecteur et l'échantillon de trajectoires dans le système d'axes `ax2`:
 
%% Cell type:code id: tags:
%% Cell type:code id:0f6c4549 tags:
 
``` python
# tracé du champs de vecteurs
# Attention à l'option angles ='xy' les fleches sont tracees avec orientation en unité naturelle de l'ecran
# et pas en unité naturelle de la figure
ax2.quiver(X, Y, dx, dy, angles = 'xy', color = 'grey', width = .002)
 
# tracé des échantillons de trajectoires
ax2.streamplot(X, Y, dx, dy, density = 0.4, maxlength = 0.25, color = "purple")
 
# on réaffiche la figure
display(fig2)
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:2e0e29eb tags:
 
 
### Intégrale première
 
Finalement, nous représentons l'intégrale première de ce système.
 
%% Cell type:code id: tags:
%% Cell type:code id:c687a568 tags:
 
``` python
# définition de la fonction conservée au cours du temps dans le modèle de Lotka Volterra
def int_premiere(etat, params):
x, y = etat # recupere les variables d'etat
r_l, c_l, b_l, m_l = params # recupere les parametres
H_xy = -r_l*np.log(y) + c_l*y - m_l*np.log(x) + b_l*x # calcule la valeur de la fonction
return H_xy
```
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:30756115 tags:
 
On peut vérifier la constance de cette quantité au cours du temps:
 
%% Cell type:code id: tags:
%% Cell type:code id:b56a898f tags:
 
``` python
# création d'une figure, et d'un système d'axe
fig3, ax3 = plt.subplots(1, 1, figsize=(14, 8))
 
# titre de la figure
fig3.suptitle("Intégrale première au cours du temps\n modèle de Lotka Volterra", va='top', fontsize='18')
 
# tracé des isoclines nulles, des équilibres, et de la trajectoire simulée
# isoclines nulles de x
ax3.plot(tspan, int_premiere([int_LV[:,0], int_LV[:,1]], params_LV))
 
# limites des axes
ax3.set_ylim(bottom=0, top=3)
 
# labellisation des axes
ax3.set_xlabel('Intégrale première', fontsize='16')
ax3.set_ylabel('Temps', fontsize='16')
 
# ajout d'une grille
ax3.grid()
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:108b373e tags:
 
Et finalement représenter ces quantités en 3D, en fonction des valeurs de $x$ et $y$. Pour cela nous utilisons une grille, comme plus haut pour le plan de phase, et evaluons la fonction sur cette grille.
 
%% Cell type:code id: tags:
%% Cell type:code id:a1701cfa tags:
 
``` python
# définition de l'échantillonnage selon $x$ et $y$
x_grid = np.linspace(0.15, 3.0, 30) # au passage on change un peu de np.arange()
y_grid = np.linspace(0.15, 3.0, 30)
 
# grille X,Y selon x_grid et y_grid
X, Y = np.meshgrid(x_grid, y_grid)
 
# calcul de la fonction conservee sur la grille
int_premiere_grid = int_premiere([X, Y], params_LV)
```
 
%% Cell type:markdown id: tags:
%% Cell type:markdown id:bd7c3be5 tags:
 
Nous créons une figure avec des paramètres adequats pour un graphique 3D et tracons la fonction en fonction des valeurs de $x$ et $y$ sur la grille.
 
%% Cell type:code id: tags:
%% Cell type:code id:38ffbda5 tags:
 
``` python
# création d'une figure, et d'un système d'axe un peu particulier puisqu'il est en 3D
fig4, ax4 = plt.subplots(1, 1, subplot_kw={"projection": "3d"}, figsize=(14, 8))
 
# nous importons des colormaps pour représenter la fonction integrale_premiere
from matplotlib import cm
 
# nous tracons l'integrale première sur la grille
ax4.plot_surface(X, Y, int_premiere_grid, cmap=cm.viridis, antialiased=True, alpha =.7)
 
# réglage de l'angle de vision en fonction de l'élévation et de l'azimut
ax4.view_init(elev=10, azim= 30)
 
# labellisation des axes
ax4.set_xlabel('Proies $x$', fontsize='16')
ax4.set_ylabel('Prédateurs $y$', fontsize='16')
ax4.set_zlabel("intégrale première", fontsize='16')
 
# titre de la figure
fig4.suptitle("Intégrale première\n modèle de Lotka Volterra", va='top', fontsize='18');
```
 
%% Output