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.
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.
Pour plus de lisibilité du code on désencapsule les paramètres et la valeur des paramètres dans le corps de la fonction en créant des valeurs locales :
...
...
@@ -241,7 +241,7 @@ Il reste à mettre un peu mieux en forme 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.
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 :
...
...
@@ -278,7 +278,7 @@ 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)
Un script complet rassemblant les différentes étapes que nous avons vu plus haut est disponible sur la page de [corrections](https://forgemia.inra.fr/ludovic.mailleret/biomaths-mam4/-/blob/master/biomaths_mam4_malthus.ipynb).
***
...
...
@@ -300,13 +300,13 @@ Pour rappel, le modèle s'écrit :
#### Dynamiques par rapport au temps, $`E`$ constant
Nous allons tout d'abord intégrer simplement ce modèle, comme nous l'avons fait plus haut, avec $`r=1,\ K_a=2,\ K=10,\ E=0.2`$ ou $`E=0.85`$ et $`x(0)=10`$.
Nous allons tout d'abord intégrer simplement ce modèle, comme nous l'avons fait plus haut, avec $r=1,\ K_a=2,\ K=10,\ E=0.2$ ou $E=0.85$ et $x(0)=10$.
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
#### 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 :
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.
...
...
@@ -320,40 +320,40 @@ Le principe de la simulation est maintenant de faire varier $`E`$ au cours du te
On choisira pour valeurs de paramètres par exemple :
- $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.
L'attendu théorique est un stock qui s'éteint si $T_{overF}$ est trop important.
***
#### Représentation dans l'espace $`(E,\ x)`$
#### 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`$.
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})`$
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`$*
*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.
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.
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](https://forgemia.inra.fr/ludovic.mailleret/biomaths-mam4/-/blob/master/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 :
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)`$.
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)$.
Dans un second temps, nous supposerons que la taille de la population d'oiseaux varie lentement, en fonction des captures de chenilles et d'un taux de mortalité. Nous observerons une situation d’enchaînement de bifurcations dynamiques qui crée un cycle d'hysteresis que nous représenterons en fonction du temps et dans l'espace $`(y,x)`$.
Dans un second temps, nous supposerons que la taille de la population d'oiseaux varie lentement, en fonction des captures de chenilles et d'un taux de mortalité. Nous observerons une situation d’enchaînement de bifurcations dynamiques qui crée un cycle d'hysteresis que nous représenterons en fonction du temps et dans l'espace $(y,x)$.
***
...
...
@@ -361,11 +361,11 @@ Dans un second temps, nous supposerons que la taille de la population d'oiseaux
On choisira pour valeurs de paramètres par exemple :
- $`r=5,\ K=10,\ h=0.5,\ \alpha = 1`$ et la densité de population d'oiseaux $`y=7`$
- $r=5,\ K=10,\ h=0.5,\ \alpha = 1$ et la densité de population d'oiseaux $y=7$
La simulation ne pose pas de difficulté particulière par rapport à ce qu'on a fait précédemment (si ce n'est de se rappeler que le carré en Python s'écrit `**2`). Il peut être intéressant de représenter différentes dynamiques afin d'illustrer la multi-stabilité dans ce modèle.
Le tracé du diagramme de bifurcation se fait par une méthode très similaire au modèle sur la pêche via un array `x_plot` et le calcul des $`y`$ qui correspondent aux valeurs des équilibres.
Le tracé du diagramme de bifurcation se fait par une méthode très similaire au modèle sur la pêche via un array `x_plot` et le calcul des $y$ qui correspondent aux valeurs des équilibres.
***
...
...
@@ -384,18 +384,18 @@ Le modèle prend maintenant la forme :
\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`$.
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`.
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.
Vous pouvez vous rapporter à la page de [corrections](https://forgemia.inra.fr/ludovic.mailleret/biomaths-mam4/-/blob/master/biomaths_mam4_tordeuse.ipynb) pour confronter vos résultats aux attendus.
***
...
...
@@ -414,11 +414,11 @@ Nous considérons le modèle de dynamique de populations (classique) de Lotka et
\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.
Nous utiliserons par ailleurs les méthodes `streamplot()` et `quiver()` des systèmes d'axes pour représenter l'orientation du champs de vecteur et quelques échantillons de trajectoires approximées. Ces deux méthodes calculent le champs de vecteur ou une approximation des trajectoires sur une grille dans le plan $`(x,y)`$. Il convient d'abord de définir cette grille :
Nous utiliserons par ailleurs les méthodes `streamplot()` et `quiver()` des systèmes d'axes pour représenter l'orientation du champs de vecteur et quelques échantillons de trajectoires approximées. Ces deux méthodes calculent le champs de vecteur ou une approximation des trajectoires sur une grille dans le plan $(x,y)$. Il convient d'abord de définir cette grille :
```python
# définition de l'échantillonnage selon $x$ et $y$
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.*
*Pour cette fin de partie vous pouvez vous aider du fichier de [correction](https://forgemia.inra.fr/ludovic.mailleret/biomaths-mam4/-/blob/master/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:
...
...
@@ -483,7 +483,7 @@ On peut facilement représenter l'évolution de cette fonction au cours du temps
- Le calcul d'une fonction de deux variables à valeurs dans $`\mathbb{R}`$ sur une grille (qui est cependant similaire à ce que l'on a fait plus haut pour la représentation du champs de vecteur):
- Le calcul d'une fonction de deux variables à valeurs dans $\mathbb{R}$ sur une grille (qui est cependant similaire à ce que l'on a fait plus haut pour la représentation du champs de vecteur):
```python
# définition de l'échantillonnage selon $x$ et $y$
...
...
@@ -497,7 +497,7 @@ X, Y = np.meshgrid(x_grid, y_grid)
int_premiere_grid=int_premiere([X,Y],params_LV)
```
- La représentation de cette fonction en fonction de $`x`$ et $`y`$ :
- La représentation de cette fonction en fonction de $x$ et $y$ :
fig4.suptitle("Intégrale première\n modèle de Lotka Volterra",va='top',fontsize='18')
```
Finalement nous pouvons surimposer sur cette figure le plan correspondant à la valeur initiale de l'intégrale première en fonction de $`x`$ et $`y`$, ainsi que l'évolution de cette quantité au cours du temps le long de la trajectoire simulée:
Finalement nous pouvons surimposer sur cette figure le plan correspondant à la valeur initiale de l'intégrale première en fonction de $x$ et $y$, ainsi que l'évolution de cette quantité au cours du temps le long de la trajectoire simulée:
```python
# tracé du plan correspondant a la valeur initiale de l'intégrale première
...
...
@@ -541,7 +541,7 @@ Nous considérons le modèle de dynamique de populations attribué à Rosenzweig
\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
...
...
@@ -553,9 +553,9 @@ On prendra comme valeurs de paramètres :
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`$
#### Diagramme de bifurcations $y_\infty$ en fonction de $K$
Il y a de nombreuses manières de faire des diagrammes de bifurcations en fonction de paramètres, notamment :
...
...
@@ -578,4 +578,4 @@ Une fois l'espace de paramètres exploré au travers de cette boucle, on pourra
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).
Il est fort possible que nous manquions de temps pour faire cette partie, je vous invite à consulter [la correction](https://forgemia.inra.fr/ludovic.mailleret/biomaths-mam4/-/blob/master/biomaths_mam4_RMA.ipynb).
