Commit 63c0b476 authored by Ludovic Mailleret's avatar Ludovic Mailleret
Browse files

update simulations_python.md

parent 52d617fb
......@@ -12,32 +12,22 @@ Support de cours pour la séance de biomaths du 22/03 - EPU MAM4 - 2020-2021
***
*Rédaction : L. Mailleret, Mars 2021*
*Rédaction : L. Mailleret, Mars 2022*
## Généralités
## Généralités
### Aides et solutions
Ce support de cours est accompagné d'un [répertoire gitlab](https://forgemia.inra.fr/ludovic.mailleret/biomaths-mam4) qui rassemble plusieurs notebooks jupyter détaillant différents éléments de correction pour chaque partie, quelques aides pour la rédaction de certains éléments de code, ainsi que les figures de résultats que vous devez pouvoir générer.
Dans la mesure du possible, notamment au début de la séance, n'utilisez pas ces fichiers et essayez de faire les choses par vous-mêmes. Vous apprendrez bien plus efficacement ainsi.
### Début de script
De manière générale, quel que soit le langage, j'ai pour habitude de toujours commencer dans un environnement "propre", avec toutes les variables définies précédemment effacées. Cela permet de ne pas récupérer d'anciennes valeurs de paramètres ou variables issues de précédentes exécutions de script et rend en général le code plus facile à débugger.
Si vous utilisez Spyder, par le menu `Exécution > Configuration par fichier` (ou `Ctrl+F6` par défaut) cocher la case `Effacer toutes les variables avant exécution`.
Si vous utilisez un notebook Jupyter, entrez la commande magique suivante dans la première cellule :
```python
......@@ -63,8 +53,6 @@ import numpy as np
import matplotlib.pyplot as plt
```
Si vous utilisez Spyder, je préfère fermer toutes les fenêtres graphiques en début de script, afin d'être sûr que ce qui est visualisé est bien ce qui vient d'être calculé.
```python
......@@ -72,8 +60,6 @@ Si vous utilisez Spyder, je préfère fermer toutes les fenêtres graphiques en
plt.close('all')
```
Si vous utilisez Jupyter les graphiques sont inclus dans le notebook et cette étape est inutile.
***
......@@ -83,17 +69,20 @@ Si vous utilisez Jupyter les graphiques sont inclus dans le notebook et cette é
La fonction `odeint()` permet d'intégrer une équation différentielle ordinaire. Elle prend pour arguments essentiels :
- une fonction renvoyant la valeur de la dérivée de l'état en fonction de la valeur de l'état `etat`, du temps `t` et de paramètres `params`:
```python
votre_modele(etat, t, params)
```
```python
votre_modele(etat, t, params)
```
- un `array` de l'état au temps initial `etat0` :
```python
etat0 = np.array([x0, y0])
```
```python
etat0 = np.array([x0, y0])
```
- un `array`, `tspan`, spécifiant les valeurs du temps auxquelles vous voulez récupérer les valeurs des variables d'état issues de l'intégration numérique :
```python
tspan = np.arange(t0, tfin, pas_t)
```
```python
tspan = np.arange(t0, tfin, pas_t)
```
Nous allons voir plus spécifiquement sur le modèle de Malthus les différentes syntaxes employées pour définir, intégrer et représenter graphiquement les données simulées.
......@@ -106,9 +95,11 @@ Le reste de la séance vous laissera une plus grande part d'improvisation, l'ess
### Modèle de Malthus
Le modèle de Malthus, ou à croissance exponentielle, s'écrit pour une population de densité $`x`$ comme nous l'avons vu en cours :
```math
\dot x = (n-m)\ x
```
avec $`n`$ et $`m`$ les taux de natalité et de mortalité, respectivement.
#### Paramétrisation
......@@ -124,7 +115,8 @@ etat0_malthus = np.array([x0])
```
ainsi que des paramètres du modèle :
```python
```python
# paramètres du modèle
# taux de natalité
n = 3.0
......@@ -134,10 +126,12 @@ m = 2.0
```
que l'on encapsule dans un `array` :
```python
# encapsulation des paramètres dans un array
params_malthus = np.array([n, m])
```
On définit finalement les paramètres de l'intervalle de temps sur lequel on souhaite intégrer le modèle :
```python
......@@ -154,8 +148,6 @@ qui permettent de définir un `array` partant de `t_0` jusque `t_fin` tous les `
tspan = np.arange(t_0, t_fin, pas_t)
```
#### Définition du modèle
On définit une fonction qui renvoie la dérivée $`\dot x`$ en fonction de la valeur de l'état $`x`$, du temps et des paramètres.
......@@ -165,14 +157,12 @@ Pour plus de lisibilité du code on désencapsule les paramètres et la valeur d
```python
# définition du modèle de Malthus
def modele_malthus(etat, t, params):
x = etat # on recupere l'etat
x = etat # on recupere l'etat
n_l, m_l = params # on fera bien attention à l'ordre des paramètres défini plus haut
xdot = (n_l-m_l)*x # on calcule la derivee de l'etat
xdot = (n_l-m_l)*x # on calcule la derivee de l'etat
return xdot
```
#### Intégration numérique
On fait un simple appel à la fonction `odeint()` et on en stocke la sortie dans un `array` (qui est *en colonne*) qui représente les valeurs numériques prises par la densité de la population à chaque temps défini dans `tspan` :
......@@ -192,8 +182,6 @@ L'appel à `odeint()` nécessite :
*NB: d'autres arguments, par exemple de tolérance, peuvent être passés à `odeint()`, je vous renvoie sur l'aide de la fonction, mais les paramètres par défaut sont déjà satisfaisants*
A ce stade, l'exécution du script montre que l'intégration numérique est faite :
```python
......@@ -206,21 +194,14 @@ array([[ 1.00000000e-01],
[ 4.70826724e+07],
[ 4.75558612e+07],
[ 4.80338056e+07]])
```
Le travail restant est essentiellement lié à la représentation graphique des données simulées ; c'est une tâche qui n'en est pas moins importante.
#### Représentation des trajectoires en fonction du temps
Les représentations graphiques utilisent le module `matplotlib.pyplot` que l'on a importé en début de script sous l'alias `plt`.
On commence par créer un objet fenêtre graphique `fig1` associé à un objet système d'axes `ax1` qui permettra la représentation des données :
```python
......@@ -229,8 +210,6 @@ On commence par créer un objet fenêtre graphique `fig1` associé à un objet
fig1, ax1 = plt.subplots()
```
On trace ensuite le résultat de l'intégration numérique `int_malthus` en fonction du temps `tspan` en utilisant la méthode `plot()` de l'objet `ax1` :
```python
......@@ -243,12 +222,10 @@ De nombreuses options peuvent être passées à la méthode `plot()` . Ici j'ai
*NB : Je vous renvoie à l'aide de `plt.subplots` ou `plt.plot` ainsi qu'aux nombreux tutoriels sur internet pour différentes options graphiques.*
Il reste à mettre un peu mieux en forme la figure :
- labellisation des axes :
```python
# labellisation des axes
ax1.set_xlabel('temps')
......@@ -256,23 +233,23 @@ Il reste à mettre un peu mieux en forme la figure :
```
- mise en forme d'un titre :
```python
# titre de la figure
fig1.suptitle('Simulation du modèle de Malthus\n $n = {}, m = {}$'.format(n, m), fontsize = '12')
```
notez l'utilisation de la méthode `format()` des chaînes de caractères qui permet de renseigner les valeurs des termes accolades par les valeurs du taux de natalité $`n`$ et du taux de mortalité $`m`$ utilisés dans la simulation.
- modifier le minimum et le maximum des axes :
```python
# modification des bornes des axes
ax1.set_ylim(bottom=None, top=np.max(int_malthus))
```
- ajouter une grille :
```python
# ajout d'une grille
ax1.grid()
......@@ -288,14 +265,10 @@ ax1.set_ylim(bottom=np.min(int_malthus), top=np.max(int_malthus))
mais il convient de recalculer alors les bornes des axes.
#### Résultats de la simulation
![Simulations du modèle de Malthus](./img/Malthus.png)
Petites variations par rapport à ce que nous avons vu plus haut :
- création d'une figure avec deux subplots en lignes (modification de l'appel à `plt.subplots`), dans lesquels j'ai tracé les solutions en échelle linéaire et logarithmique
......@@ -303,12 +276,8 @@ Petites variations par rapport à ce que nous avons vu plus haut :
*NB : Vous pouvez modifier les paramètres, conditions initiales, et observer les modifications des simulations*
Un script complet rassemblant les différentes étapes que nous avons vu plus haut est disponible sur la page de [corrections](./biomaths_mam4_malthus.ipynb)
***
## Populations exploitées
......@@ -320,6 +289,7 @@ L'idée est ici de reprendre l'étude faite en cours en simulant une pêcherie d
Nous simulerons notamment l'effet de la durée de la période de surpêche sur la survie du stock, et illustrerons ces résultats avec des figures représentant la densité de population en fonction du temps, et la densité de population en fonction de l'effort de pêche.
Pour rappel, le modèle s'écrit :
```math
\dot x = rx\left(\frac{x}{K_a}-1\right)\left(1-\frac{x}{K}\right) - E x
```
......@@ -332,70 +302,52 @@ Nous allons tout d'abord intégrer simplement ce modèle, comme nous l'avons fai
La simulation n'est pas plus compliquée que dans les exemples précédents. Il en va de même pour la représentation graphique de la densité de population par rapport au temps.
#### Dynamiques par rapport au temps, $`E`$ varie au cours du temps
Le principe de la simulation est maintenant de faire varier $`E`$ au cours du temps pour étudier les effets de la surpêche. Pour cela :
- créer une fonction `E_varie()`, prenant pour argument le temps $`t`$ et les paramètres $`E_{sust},\ E_{overF},\ T_{sust},\ T_{overF},`$ qui selon la valeur du temps $`t`$ renvoie la valeur $`E_{sust}`$ (effort de pêche durable) ou $`E_{overF}`$ (surpêche) comme décrit dans la figure ci-dessous.
![effort de peche](./img/effort_peche.png)
*NB : un simple `if` renvoyant la valeur $`E_{sust}`$ ou $`E_{overF}`$ selon la valeur du temps $`t`$ fera l'affaire à ce stade*
*NB : un simple `if` renvoyant la valeur $`E_{sust}`$ ou $`E_{overF}`$ selon la valeur du temps $`t`$ fera l'affaire à ce stade*
- répliquer la figure ci-dessus
- utiliser la fonction `E_varie()` pour simuler le modèle
On choisira pour valeurs de paramètres par exemple :
- $`r=1,\ K_a=2,\ K=10,\ E_{sust}=0.2,\ E_{overF}=0.85`$
- $`T_{sust} = 10`$, et : $`T_{overF }=9.2`$ ou : $`T_{overF} =9`$
L'attendu théorique est un stock qui s'éteint si $`T_{overF}`$ est trop important.
***
#### Représentation dans l'espace $`(E,\ x)`$
On veut ici représenter le diagramme de bifurcation du modèle en fonction de $`E`$ et surajouter dans ce graphique les valeurs de $`x`$ simulées plus haut en fonction de $`E`$.
La partie bifurcation nécessite de tracer le lieu des points d'équilibre dans le plan $`(E, x)`$ en générant un `array` $`x_{plot}`$ entre $`K_a`$ et $`K`$ et, comme dans le cours, en calculant la valeur $`E_{plot}`$ définissant les équilibres et finalement en traçant la courbe définie par $`(E_{plot}, x_{plot})`$
*NB : on n'oubliera pas de tracer $`x=0`$*
Pour la représentation des simulations numériques, on a déjà la valeur de la densité de population $`x`$ au long du temps `tspan`. Nous avons aussi du calculer `E_varie()` au cours du temps pour tracer la figure ci-dessus. Il suffit alors de tracer la valeur de l'état $`x`$ au cours du temps en fonction de la valeur de `E_varie()` au cours du temps.
Un aperçu des résultats attendus et le cas échéant quelques aides à l'écriture du programme sont disponibles dans le fichier la page de [corrections](./biomaths_mam4_AE.ipynb), sur laquelle vous pourrez aussi confronter vos résultats aux attendus.
***
### Tordeuse du bourgeon de l'épinette
Nous reprenons l'étude faite en cours sur le modèle de la tordeuse du bourgeon de l'épinette [(Ludwig, Jones et Holling, 1978)](https://jmahaffy.sdsu.edu/courses/f09/math636/lectures/bifurcation_ode/ludwig_ecology_78.pdf). Le modèle que nous allons étudier ici est légèrement différent de celui vu en cours, car nous allons expliciter la population d'oiseaux prédateurs. Avec $`x`$ la densité de tordeuses et $`y`$ la densité d'oiseaux, le modèle s'écrit :
```math
\dot x =rx\left(1-\frac{x}{K}\right) - \frac{\alpha x^2}{h^2+x^2}\ y
```
$`\alpha`$ est le nombre maximal de chenilles que les oiseaux peuvent consommer *per capita* ; les autres paramètres ont été vus en cours.
$`\alpha`$ est le nombre maximal de chenilles que les oiseaux peuvent consommer *per capita* ; les autres paramètres ont été vus en cours.
Nous allons dans un premier temps étudier le modèle à population d'oiseaux constante, représenter plusieurs dynamiques en fonction du temps et tracer le diagramme de bifurcations correspondant dans l'espace $`(y,x)`$.
......@@ -420,6 +372,7 @@ Le tracé du diagramme de bifurcation se fait par une méthode très similaire a
Lorsque l'on prend en compte que la taille de la population d'oiseaux varie au cours du temps en fonction des captures, le modèle est directement en 2 dimensions. En supposant que la dynamique des oiseaux est lente par rapport à celle des insectes, on peut cependant déduire le comportement en projetant la dynamique des tordeuses à l'équilibre et en étudiant la dynamique résultante sur les oiseaux (approche dite "lent-rapide"). Comme ici, selon la taille de la population d'oiseaux, différents équilibres peuvent exister et être globalement stable, cela va donner un lieu a un comportement dynamique intéressant que l'on pourra entièrement appréhender sur le diagramme de bifurcation tracé précédemment.
Le modèle prend maintenant la forme :
```math
\left\{
\begin{array}{l}
......@@ -428,27 +381,22 @@ Le modèle prend maintenant la forme :
\end{array}
\right.
```
avec $`0<\varepsilon\ll 1`$, et $`n`$ et $`m`$ les taux associés à la natalité et à la mortalité des prédateurs. En pratique fixons $`\varepsilon=0.01`$, et les paramètres $`n=5`$ et $`m=3`$.
La seule difficulté ici est qu'il ne s'agit plus d'intégrer une équation différentielle de dimension 1, mais un système de deux équations différentielles.
Pour cela, il faut fournir à `odeint()` un modèle qui renvoie les valeurs des dérivées de l'état $`\dot x`$ et $`\dot y`$ sous forme d'un array ligne de type `[xdot, ydot]` . `odeint()` renvoie ensuite les valeurs de l'intégration numérique sous forme d'un tableau en dimension 2 dont la $`k`$-ième colonne représente la valeur de la $`k`$-ième variable le long de `tspan`.
Une fois l'intégration numérique effectuée vous tracerez :
- la densité de population de tordeuses en fonction du temps
- la densité de population de tordeuses en fonction de la densité de population d'oiseaux, au cours du temps, surimposée au diagramme de bifurcation du modèle avec la population d'oiseaux constante.
Vous pouvez vous rapporter à la page de [corrections](./biomaths_mam4_tordeuse.ipynb) pour confronter vos résultats aux attendus.
***
## Populations en interactions
### Le modèle proies-prédateurs de Lotka-Volterra
......@@ -463,9 +411,8 @@ Nous considérons le modèle de dynamique de populations (classique) de Lotka et
\dot y = bxy - m y
\end{array}\right.
```
Nous allons simuler ce modèle pour $`r=1, c=1, b=1, m=1`$, représenter les dynamiques de populations en fonction du temps et dans le plan de phase $`(x,y)`$. Cette partie est calquée sur ce que nous avons fait précédemment pour le modèle de tordeuse avec population d'oiseaux variable.
Nous allons simuler ce modèle pour $`r=1, c=1, b=1, m=1`$, représenter les dynamiques de populations en fonction du temps et dans le plan de phase $`(x,y)`$. Cette partie est calquée sur ce que nous avons fait précédemment pour le modèle de tordeuse avec population d'oiseaux variable.
Les subtilités concernent la représentation dans le plan de phase. Pour cela nous calculerons les isoclines nulles et représenterons les équilibres, et la trajectoire simulée.
......@@ -502,14 +449,14 @@ et les échantillons de trajectoires :
ax2.streamplot(X, Y, dx, dy, density = 0.4, maxlength = 0.25, color = "purple")
```
#### Intégrale première
Finalement, nous représentons l'intégrale première (quantité conservée au cours du temps dans ce système) en fonction du temps, puis dans un graphique en 3D en fonction de $x$ et de $y$, ainsi qu'une trajectoire du système dans cet espace. La quantité conservée au cours du temps dans le modèle de Lotka Volterra est la suivante:
```math
H(x,y) = b\ x-m\ln(x) + c\ y- r\ln(y)
```
Nous commençons par définir cette fonction, ce qui devrait être assez direct:
```python
......@@ -523,8 +470,6 @@ def int_premiere(etat, params):
On peut facilement représenter l'évolution de cette fonction au cours du temps en l'évaluant le long de la trajectoire simulée plus haut et en la traçant en fonction de `tspan`.
*Pour cette fin de partie vous pouvez vous aider du fichier de [correction](./biomaths_mam4_LV.ipynb), mais l'essentiel des commandes utiles est repris ci-dessous.*
La partie 3D est plus technique et utilise différentes fonctions spécifiques aux graphiques 3D, à commencer par:
......@@ -583,32 +528,30 @@ ax4.plot(int_LV[:,0], int_LV[:,1], int_premiere([int_LV[:,0], int_LV[:,1]], para
display(fig4)
```
### Le modèle proies-prédateurs de Rosenzweig-MacArthur
Nous considérons le modèle de dynamique de populations attribué à Rosenzweig et MacArthur (voir [Rosenzweig & MacArthur (1963)](http://www.jstor.com/stable/2458702), [Turchin (2003)](http://www.jstor.com/stable/2458702), [Smith (2008)](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.710.95&rep=rep1&type=pdf)).
```math
\left\{\begin{array}{l}
\dot x = \displaystyle rx\left(1-\frac{x}{K}\right) - c \frac{x}{h+x} y\\
\dot y = b\displaystyle \frac{x}{h+x} y - m y
\end{array}\right.
```
avec $`x`$ et $`y`$ les densités de proies et prédateurs, respectivement.
avec $`x`$ et $`y`$ les densités de proies et prédateurs, respectivement.
#### Simulation des dynamiques et plan de phase
La simulation et la représentation graphique est vraiment dans la droite ligne de ce qui a été fait plus haut pour le modèle de Lotka Volterra.
On prendra comme valeurs de paramètres :
```math
r=1,\ K=10,\ c= 1,\ h=2,\ b=2,\ m=1.
```
Pour ces valeurs, il existe un cycle limite attractif. Vous pouvez diminuer $`K`$ pour voir l'effondrement du cycle en un équilibre positif asymptotiquement stable, voire l'extinction des prédateurs.
Pour ces valeurs, il existe un cycle limite attractif. Vous pouvez diminuer $`K`$ pour voir l'effondrement du cycle en un équilibre positif asymptotiquement stable, voire l'extinction des prédateurs.
#### Diagramme de bifurcations $`y_\infty`$ en fonction de $`K`$
......@@ -618,12 +561,8 @@ Il y a de nombreuses manières de faire des diagrammes de bifurcations en foncti
- en suivant numériquement les équilibres (ou des invariants) et en caractérisant numériquement les changements de stabilité portant la signature des bifurcations (e.g [Blyth *et alii*, 2020](https://arxiv.org/pdf/2008.05226.pdf)),
- par force brute numérique, en simulant le système jusqu'à atteindre les états asymptotiques.
C'est cette dernière méthode que nous allons utiliser ici, en combinaison avec l'approche analytique. Ce n'est pas la méthode la plus élégante, mais elle permettra d'obtenir le diagramme recherché.
Le principe de l’algorithme est le suivant :
- fixer une valeur de paramètre
......@@ -635,23 +574,6 @@ Le principe de l’algorithme est le suivant :
Une fois l'espace de paramètre exploré au travers de cette boucle, on pourra tracer les valeurs asymptotiques calculées en fonction du paramètre, et compléter la figure avec les équilibres instables calculés à la main.
En pratique on réservera l'algorithme présenté ci-dessus au calcul des extrema du cycle limite, et on tracera directement les équilibres que nous savons caractériser sans simulation.
Il est fort possible que nous manquions de temps pour faire cette partie, je vous invite à consulter [la correction](./biomaths_mam4_RMA.ipynb).
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment