using DifferentialEquations
using Plots
Populations isolées
Le modèle de Malthus
Nous considérons le modèle proposé par Malthus (1798) : \dot x = (n-m)x, \tag{1} avec n le taux de natalité, et m le taux de mortalité.
Préliminaires
Pour simuler ce modèle, c’est à dire intégrer numériquement les solutions au problème de Cauchy correspondant à l’Equation 1 avec x(0)=x_0\geq0, nous allons utiliser les routines de DifferentialEquations.jl
avec la méthode par défaut. Nous utiliserons aussi le package Plots.jl
pour les représentations graphiques1.
1 si vous n’avez jamais installé un package, il faudra commencer par le faire avant de pouvoir l’utiliser…
Conditions initiales, paramètres et temps
Nous définissions la condition initiale de la simulation :
= 0.1 x0
Nous définissons les paramètres, et les encapsulons dans un vecteur de paramètres :
= 3.0 # natalité
n = 2.0 # mortalité
m = [n, m] # packing in an array par_malthus
Enfin, les propriétés du temps d’intégration :
= (0.0, 10.0)
tspan = 0.1 tstep
Système dynamique
On définit le modèle comme une fonction renvoyant la dérivée de x en fonction de l’état (u
est la notation pour DifferentialEquations.jl
), de paramètres (p
), et du temps (t
). L’ordre des arguments est important, c’est sous cette forme que les routines d’intégration attendent le modèle.
function malthus(u, p, t)
= p # unpacking
n, m = u # use x notation
x return dx = (n-m)x # return derivative
end
L’ordre des arguments u, p, t
est important pour le solver de DifferentialEquations.jl
.
Intégration
On commence par définir le problème de Cauchy à intégrer, comme un ODEProblem
avec arguments: le modèle, la condition initiale, les bornes d’intégration, les paramètres ainsi que l’option saveat = tstep
permettant de récupérer la solution tout les tstep
pas de temps (il y a plein d’options de sortie de l’intégration, dont même une fonction du temps !)
= ODEProblem(
prob_malthus # modèle
malthus, # condition initiale
x0, # tspan
tspan, # paramètres
par_malthus; = tstep, # keyword argument, option de sortie
saveat )
On intègre le modèle via solve
, défini par DifferentialEquations.jl
pour des struct
de type ODEProblem
:
= solve(prob_malthus) sol_malthus
Le type de solution renvoyé par le solveur est assez complexe et comprend de nombreux champs informatifs sur le calcul. On peut accéder au temps de simulation via sol_malthus.t
:
first(sol_malthus.t, 3)
3-element Vector{Float64}:
0.0
0.1
0.2
ainsi qu’aux valeurs de la variable x calculées le long du temps via sol_matlhus.u
:
first(sol_malthus.u, 3)
3-element Vector{Float64}:
0.1
0.1105170918098962
0.12214028021690636
Même si ce n’est pas indispensable, il est possible de transformer la solution renvoyée facilement en DataFrame
, qui peut permettre des manipulations plus faciles.
using DataFrames
= DataFrame(sol_malthus)
sol_malthus rename!(sol_malthus, :timestamp => :time, :value => :x)
remarquez la fonction rename!
qui modifie en place le dataframe. Par convention les fonctions dont le nom finit par !
modifient leur argument en place.
si bien que :
first(sol_malthus, 3)
Row | time | x |
---|---|---|
Float64 | Float64 | |
1 | 0.0 | 0.1 |
2 | 0.1 | 0.110517 |
3 | 0.2 | 0.12214 |
Représentation graphique
On peut représenter graphiquement la simulation de la croissance de la population au cours du temps (ici via le dataframe).
plot(
# abscisses
sol_malthus.time, # ordonnées
sol_malthus.x; = :tab10, # palette de couleurs
palette = 2,
linewidth = "\n Modèle de Malthus \$n=$n, m=$m\$",
title = "population \$x\$",
label = "densité de population \$x(t)\$",
ylabel = "temps \$t\$",
xlabel = .5Plots.cm,
margin = 1Plots.cm,
topmargin )
Les titres et autres chaînes de caractères peuvent utiliser des formules LaTeX via \$
, e.g. "\$x\$"
.
La notation "$m"
accède à la valeur du paramètre m
et la renvoit dans la chaine de caractères.
Le modèle logistique
Nous considérons ici le modèle “logistique” proposé par Verhulst (1838) :
\dot x = r x \left(1-\frac{x}{K}\right), \tag{2} avec r le taux de croissance intrinsèque de la population et K la capacité de charge de l’environnement.
Il n’y a pas de difficulté particulière par rapport aux simulations précedentes. On va réutiliser la condition initiale et les paramètres de temps définis précédemment.
Il faut cependant définir les paramètres du modèle Equation 2 :
Code
= 1.0 # natalité
r = 10.0 # mortalité
K = [r, K] # packing par_logistic
ainsi que le nouveau système dynamique :
Code
function logistic(u, p, t)
= p # unpacking
r, K = u # use x notation
x return dx = r*x*(1-x/K) # return derivative
end
et le problème de Cauchy correspondant :
Code
= ODEProblem(
prob_logistic # modèle
logistic, # condition initiale
x0, # tspan
tspan, # paramètres
par_logistic; = tstep, # option de sortie
saveat )
On simule et on transforme la solution en dataframe :
Code
= solve(prob_logistic)
sol_logistic
= DataFrame(sol_logistic)
sol_logistic rename!(sol_logistic, :timestamp => :time, :value => :x)
On trace la solution, en rajoutant les équilibres stable (x=K) et instable (x=0) :
# solution
plot(
sol_logistic.time,
sol_logistic.x;= :tab10,
palette = 2,
linewidth = "\n Modèle logistique \$r=$r, K=$K\$",
title = "population \$x\$",
label = :right,
legend = "densité de population \$x(t)\$",
ylabel = "temps \$t\$",
xlabel = .5Plots.cm,
margin = 1Plots.cm,
topmargin
)
# équilibre 0
plot!(
sol_logistic.time,zeros(length(sol_logistic.time));
= "red",
color = 2,
linewidth = :dash,
linestyle = .5,
linealpha = "équilibre instable",
label
)
# équilibre K
plot!(
sol_logistic.time,ones(length(sol_logistic.time)).*K;
= "green",
color = 2,
linewidth = :dash,
linestyle = .5,
linealpha = "équilibre stable",
label )
Remarquez l’utilisation de plot!
qui modifie en place la première figure en rajoutant des éléments (ici les équilibres et leurs labels).
Effets Allee
On s’intéresse à un modèle de dynamique de population avec “effets Allee forts”, souvent attribué à Gruntfest, Arditi, and Dombrovsky (1997)2 :
2 mais de nombreuses variations de cette forme polynomiale existent dans la littérature depuis Bazykin (1985)
\dot x = r x \left(\frac{x}{\epsilon}-1\right)\left(1-\frac{x}{K}\right), \tag{3} avec r le taux de croissance intrinsèque de la population (par analogie avec la logistique), K la capacité de charge de l’environnement et \epsilon le seuil en dessous duquel la population n’est pas viable (‘seuil de Allee’).
Nous souhaitons illustrer la bistabilité caractéristique du comportement de ce modèle à partir de la simulation depuis plusieurs conditions initiales.
Commençons par définir le modèle et les paramètres (en conservant les r et K précédents):
Code
= 2 # seuil de Allee
epsilon = [r, K, epsilon] # packing
par_allee = (0.0, 3.0)
tspan_allee
function allee(u, p, t)
= p # unpacking
r, K, epsilon = u # use x notation
x return dx = r*x*(x/epsilon - 1)*(1 - x/K) # derivative
end
Nous définissons une fonction qui réalise une simulation en un seul appel depuis une condition initiale x_0 :
function int_allee(x0, tspan = tspan_allee, param = par_allee)
= ODEProblem(
prob_allee # modèle
allee, # condition initiale
x0, # tspan
tspan, # paramètres
param; = tstep, # option de sortie
saveat
)
= solve(prob_allee)
sol_allee = DataFrame(sol_allee)
sol_allee rename!(sol_allee, :timestamp => :time, :value => :x)
return sol_allee
end
La fonction int_allee
comporte un nombre d’arguments variable: int_allee(x0)
renverra la solution avec les arguments tspan
et param
à leurs valeurs par défaut (ce sont des vararg); c’est le même appel que int_allee(x0, tspan_allee, par_allee)
.
Nous définissons un vecteur de conditions initiales différentes :
= 1.35
x0step = x0step:x0step:K # range de valeur x0vec
Finalement on réalise la figure. La stratégie diffère un peu de ce qui a été vu ci-dessus. Nous commençons par initier un graphique fig
, et faisons une boucle for
pour tracer chacune des simulations correspondant aux différentes conditions initiales3 :
3 la palette de couleur est accessoire, juste pour l’esthétique
# custom color palette
= palette([:steelblue, :lightblue], length(x0vec))
init_cgrad
# initialisation du graphique
= plot(; # on initie uniquement le graphique
fig = init_cgrad,
palette = "\n Modèle à effets Allee \$r=$r, K=$K\$, \$ϵ=$epsilon\$",
title = :right,
legend = "densité de population \$x(t)\$",
ylabel = "temps \$t\$",
xlabel = .5Plots.cm,
margin = 1Plots.cm,
topmargin
)
# boucle de plot avec intégration pour differentes conditions initiales
for x0 in x0vec # x0 parcourt l'array x0vec
plot!(
# on modifie fig
fig, int_allee(x0).time, # abscisses
int_allee(x0).x; # ordonnées
= 2, # keyword arguments après ;
linewidth = "",
label
)end
# équilibre 0
plot!(
fig,int_allee(0).time,
zeros(length(int_allee(0).time));
= "green",
color = 2,
linewidth = :dash,
linestyle = .5,
linealpha = "équilibre stable",
label
)
# équilibre epsilon
plot!(
fig,int_allee(0).time,
ones(length(int_allee(0).time)).*epsilon;
= "red",
color = 2,
linewidth = :dash,
linestyle = .5,
linealpha = "équilibre instable",
label
)
# équilibre K
plot!(
fig,int_allee(0).time,
ones(length(int_allee(0).time)).*K;
= "green",
color = 2,
linewidth = :dash,
linestyle = .5,
linealpha = "",
label
)
display(fig) # actually shows the plot P
On peut sauvegarder la figure dans différents formats (e.g. png, pdf) :
savefig(fig, "dyn_allee.pdf")
savefig(fig, "dyn_allee.png")
Passons aux populations exploitées.