Populations exploitées

Prélèvements et effets Allee

Nous reprenons le modèle précédent sur l’effet Allee mais en prenant en compte des prélèvements externes avec un effort (taux) de prélèvement E : \dot x = r x \left(\frac{x}{\epsilon}-1\right)\left(1-\frac{x}{K}\right)-Ex. \tag{1}

La simulation de ce modèle pour différentes valeurs de E (par exemple E=0.2 ou E=0.85) ne présente aucune difficulté supplémentaire.

Nous allons maintenant nous intéresser à une situation où l’effort de prélèvement E varie au cours du temps entre une valeur soutenable E_s (par exemple 0.2), et une valeur excessive E_x (par exemple 0.85).

L’attendu théorique est que si les prélèvements sont maintenus à une valeur excessive E_x trop longtemps, la population disparait irrémédiablement même si les prélèvements sont par la suite ramenés à une valeur initialement soutenable E_s.

Prélèvements variables dans le temps

using DifferentialEquations
using Plots
using DataFrames

Définissons une fonction effort() dépendant du temps, qui renvoit :

E_s~ si ~t<T_s
E_x~ si ~t\in [T_s+T_x[
à nouveau E_s~ si ~t\geq T_s+T_x

# paramètres
E_s = 0.2
E_x = 0.85
T_s = 10.0
T_x = 9.0

p_effort = [E_s, E_x, T_s, T_x]

function effort(t, p)
    E_s, E_x, T_s, T_x = p

    if t < T_s || t >= T_s + T_x
        return E_s
    elseif t >= T_s && t < T_s + T_x
        return E_x
    end
end

Note

Remarquez les opérateurs booléens “ou” || et “et” &&, ainsi que la structure du test if.

La fonction correspond bien à nos hypothèses :

Code

t2plot = 0: .1: 30

plot(
    t2plot,                                 # abscisses
    [effort(t, p_effort) for t in t2plot];  # ordonnées
    palette = :tab10,
    linewidth = 2,
    label = "\$E(t)\$",
    xlabel = "temps",
    ylabel = "\$E(t)\$",
    title = "Effort de pêche \$E(t)\$",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

Note

Remarquez la construction des ordonnées correspondant à la fonction effort via une compréhension de liste : [effort(t, p_effort) for t in t2plot].

Simulations en fonction du temps

Nous définissons les paramètres du modèle, la condition initiale et le temps d’intégration :

Code

x0 = 10.0

tspan = (0.0, 30.0)
tstep = 0.1

r = 1.0
K = 10.0
epsilon = 2.0
p_allee = [r, K, epsilon]

Nous définissons le système dynamique comme précédemment, à la différence que nous prévoyons de surcharger l’argument p dans le problème d’intégration sous forme d’un vecteur comprennant :

le vecteur de paramètres p_allee en première position
le vecteur de paramètres p_effort en seconde position
et la fonction effort en troisième position¹

¹ En Julia tout est passable en argument, des variables aux fonctions et plus.

p_evar = [p_allee, p_effort, effort]

function allee_evar(u, p, t)
    r, K, epsilon = p[1]        # unpacking model parameters
    p_effort = p[2]             # unpacking fishing effort parameters
    E = p[3]                    # unpacking fishing effort function
    x = u                    # use x notation

    return dx = r*x*(x/epsilon - 1)*(1 - x/K) - E(t, p_effort)*x
end

L’intégration en elle-même suit le shéma vu précédemment, si ce n’est que l’argument de paramètres doit bien refléter ce qui est attendu par la fonction allee_evar(). La simulation en elle-même est effectuée avec une modification de la précision relative de l’intégration reltol = 1e-6, la précision par défaut n’étant pas suffisante ici.

prob_allee_evar = ODEProblem(
      allee_evar,
      x0,
      tspan,
      p_evar;
      saveat = tstep,
)

sol_allee_evar = solve(prob_allee_evar; reltol = 1e-6)

sol_allee_evar = DataFrame(sol_allee_evar)
rename!(sol_allee_evar, :timestamp => :time, :value => :x)

Note

Ici nous avons imposé une tolérance plus précise au solver via le keyword argument reltol = 1e-6.

Finalement, nous pouvons représenter graphiquement la solution contre le temps. Ici, malgré la perturbation violente induite par la période de surexploitation, le pêcherie retrouve une situation soutenable après un retour à E=E_s.

Code

plot(
    sol_allee_evar.time,
    sol_allee_evar.x;
    label = "\$x(t)\$",
    linewidth = 2,
    xlabel = "temps",
    ylabel = "densité de population \$x(t)\$",
    title = "Effort de pêche variant dans le temps",
    palette = :tab10,
    legend = :left,
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

Figure 1: Persistence dans le modèle avec effet Allee et prélèvements

Lorsque la période de surexploitation de la population est trop longue (e.g. ici T_x=9.2), la population ne parvient pas à récupérer malgré le retour à un effort de prélèvement soutenable.

Code

# change parameter
T_x2 = 9.2

p_effort2 = [E_s, E_x, T_s, T_x2]
p_evar2 = [p_allee, p_effort2, effort]

# define new problem and integrate
prob_allee_evar2 = ODEProblem(
      allee_evar,
      x0,
      tspan,
      p_evar2,
      saveat = tstep
)

sol_allee_evar2 = solve(prob_allee_evar2, reltol = 1e-6)

sol_allee_evar2 = DataFrame(sol_allee_evar2)
rename!(sol_allee_evar2, :timestamp => :time, :value => :x)

# plot
plot(
    sol_allee_evar2.time, sol_allee_evar2.x,
    label = "\$x(t)\$";
    linewidth = 2,
    xlabel = "temps",
    ylabel = "densité de population \$x(t)\$",
    title = "Effort de pêche variant dans le temps",
    palette = :tab10,
    legend = :left,
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

Figure 2: Extinction dans le modèle avec effet Allee et prélèvements

Simulations dans l’espace (E, x)

Il s’agit ici de représenter dans l’espace (E, x) l’évolution conjointe de l’effort de pêche et de la densité de la population au cours du temps, afin de mieux comprendre le phénomène d’extinction. Dans cet objectif, on tracera aussi le lieu des équilibres de la population x^*(E) en fonction d’une valeur de E constante, i.e. : x^*(E) = 0, ce qui est équivalent à : E = r\left(\frac{x^*(E)}{\epsilon}-1\right)\left(1-\frac{x^*(E)}{K}\right).

Nous allons tracer les deux situations vues à la section précédente dans deux sous figures. Commençons par les lieux des équilibres, qui sont les mêmes dans les deux situations.

Définissons le lieu des équilibres positifs :

# parabole pour les equilibres positifs
function eeqpos(x, par = p_allee)
    r, K, epsilon = par
    return r*(x/epsilon -1)*(1-x/K)
end

Puis traçons ce lieu dans le plan (E,x) (il y a deux branches une stable et une instable pour les équilibres positifs qui se rejoignent en x^* = (K+\epsilon)/2, nous les traçons séparément)

# vecteurs pour le tracé
e2plot = 0:.1:1
x2plot1 = epsilon:.02:(K+epsilon)/2
x2plot2 = (K+epsilon)/2:.02:K

# on définit des couleurs spécifiques depuis la palette :pal10
mygreen = palette(:tab10)[3]
myorange = palette(:tab10)[2]
myblue = palette(:tab10)[1]
myred = palette(:tab10)[4]

# plot
fig1 = plot(
    e2plot,
    zeros(length(e2plot));
    color = mygreen,
    linewidth = 2,
    label = "équilibres stables",
    ylabel = "densité de population \$x\$",
    xlabel = "effort de pêche \$E\$",
    legend = :left,
)

plot!(
    fig1,
    eeqpos.(x2plot1),
    x2plot1;
    color= myorange,
    linewidth = 2,
    label = "équilibre instable",
)

plot!(
    fig1,
    eeqpos.(x2plot2),
    x2plot2;
    color= mygreen,
    linewidth = 2,
    label ="",
)

display(fig1)

Note

Remarquez le broadcast (calcul terme à terme) sur les ordonnées des équilibres positifs en utilisant eeqpos.. On aurait pu procéder par compréhension de liste: [eeqpos(x) for x in x2plot2].

On prépare un graphique fig2 avec les mêmes éléments que dans fig1 :

fig2 = deepcopy(fig1) # on fait une copie complète de fig1

Finalement, on complète fig1 et fig2 avec quelques annotations et les trajectoires calculées plus haut en fonction de l’effort de pêche variable au cours du temps, et on trace les résultats en deux sous-figures :

# annotations
annotate!(fig1, .35, 8.5, Plots.text("branche stable", 10, rotation=-28))
annotate!(fig1, .4, 2.7, Plots.text("branche instable", 10, rotation=28))
annotate!(fig2, .4, 2.7, Plots.text("branche instable", 10, rotation=28))
annotate!(fig2, .35, 8.5, Plots.text("branche stable", 10, rotation=-28))

# peche durable sur fig1
plot!(
    fig1,
    [effort(t, p_effort) for t in sol_allee_evar.time],
    sol_allee_evar.x;
    color = myblue,
    linewidth = 2,
    label = "trajectoire",
)

# surpeche sur fig2
plot!(
    fig2,
    [effort(t, p_effort2) for t in sol_allee_evar2.time], sol_allee_evar2.x;
    color = myblue,
    linewidth = 2,
    label = "trajectoire",
)

# figure reprenant les deux sous figures fig1 et fig2
plot(
    fig1,
    fig2;
    suptitle = "Effets Allee et prélèvements",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

Figure 3: bifurcation pli et catastrophe dans le modèle avec effet Allee et prélèvements

La tordeuse du bourgeon de l’épinette

Modèle

Nous considérons le modèle de dynamique de populations suivant, inspiré de Ludwig, Jones, and Holling (1978) :

\dot x =rx\left(1-\frac{x}{K}\right) - \frac{\alpha x^2}{h^2+x^2}\ y, \tag{2}

avec x la densité de tordeuses et y la densité d’oiseaux. La croissance des tordeuses suit une loi logistique et la prédation des oiseaux une réponse fonctionnelle de type Holling III.

Simulations

On procède classiquement, en ajustant les valeurs de paramètres pour mettre en évidence les phénomènes dynamiques attendus :

# paramètres
r = 5.0      # natalité
K = 10.0     # mortalité
α = 1.0      # taux max de prédation, \alpha + TAB
h = 0.5      # constante de demi-saturation
yc = 7.0     # densité de prédateurs

par_tordeuse = [r, K, α, h, yc]

# temps d'intégration
tspan = (0.0, 3.0)
tstep = 0.02

Puis on définit le modèle :

function tordeuse(u, p, t)
    r, K, α, h, yc = p
    x = u
    return dx = r*x*(1 - x/K) - α*x^2/(h^2 + x^2)*yc
end

Comme pour le modèle avec effets Allee, on définit une fonction qui redéfinit le problème d’intégration à chaque fois, et simule et renvoit la solution pour pouvoir illustrer la bi-stabilité² :

² des méthodes alternative utilisant l’interface integrator de DifferentialEquations.jl ou la redéfinition du problème d’intégration via remake sont proposées en annexe

function int_tordeuse(x0, tspan = tspan, param = par_tordeuse)
    prob_tordeuse = ODEProblem(
      tordeuse,       # modèle
      x0,               # condition initiale
      tspan,            # tspan
      param;            # paramètres
      saveat = tstep,
    )

    sol_tordeuse = solve(prob_tordeuse)
    sol_tordeuse = DataFrame(sol_tordeuse)
    rename!(sol_tordeuse, :timestamp => :time, :value => :x)
end

Nous simulons le modèle depuis différentes conditions intiales, et traçons les résultats via une boucle for.

# conditions initiales
x0step = 1.35
x0vec = x0step:x0step:K

# custom color palette
init_cgrad = palette([:steelblue, :lightblue], length(x0vec))

# initialisation du graphique, équilibre nul
fig3 = plot(;         # seulement des kwargs: on initialise le graphique
    palette = init_cgrad,
    legend = :right,
    label ="équilibres instables",
    title = "Tordeuse du bourgeon de l\'épinette",
    ylabel = "densité de population \$x(t)\$",
    xlabel = "temps \$t\$",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

# boucle de plot avec intégration pour differentes conditions initiales
for x0 in x0vec
    plot!(
        fig3,
        int_tordeuse(x0).time,
        int_tordeuse(x0).x,
        linewidth = 2,
        label = "",
    )
end

display(fig3)      # actually shows the plot fig3

Equilibres

Les valeurs des équilibres positifs du modèle Equation 2 n’ont pas d’écriture mathématique simple. Nous allons calculer numériquement les racines du polynôme dont ils sont solution : r\left(1-\frac{x^*}{K}\right)\left(h^2+x^{*2}\right)-\alpha x^* y = 0 \tag{3}

Pour cela, on utilise le package Polynomials.jl :

using Polynomials

# définition du monôme X
X = Polynomial([0, 1])

# définition du polynôme
pol = r*(1-X/K)*(h^2 + X^2)-α*X*yc

# calcul des racines, réelles, positives et plus petites que K
eq_pos = roots(pol)                         # calcul des racines
eq_pos = real.(eq_pos[isreal.(eq_pos)])     # filtrage des racines réelles
eq_pos = eq_pos[(eq_pos .> 0) .& (eq_pos .<= K)] # filtrage des racines >0 et <K

3-element Vector{Float64}:
 0.20406511760131657
 1.4717313636879934
 8.324203518710693

Note

A partir de la définition du monôme de degré 1, on écrit assez naturellement un polynôme et on en extrait les racines via roots.

La partie filtrage des racines réelles prend la partie réelle des racines qui sont rendues sous forme d’un type complexe (même si la partie imaginaire est nulle).

Remarquez les broadcasting divers, notamment .> et .& qui testent les conditions terme à terme sur le vecteur eq_pos.

Et on trace les différents équilibres :

t2plot = collect(tspan)

# initialisation du graphique, équilibre nul
plot!(
    fig3,
    t2plot,
    zeros(length(t2plot));
    linewidth = 2,
    linestyle = :dash,
    color = myorange,
    palette = init_cgrad,
    legend = :right,
    label ="équilibres instables",
    ylabel = "densité de population \$x(t)\$",
    xlabel = "temps \$t\$",
    margin = .5Plots.cm,
    topmargin = 1Plots.cm,
)

# équilibres positifs, différents cas
if length(eq_pos) == 1
    plot!(
      fig3,
      t2plot, ones(length(t2plot)).*eq_pos;
      color = mygreen,
      label ="équilibre stable",
    )
elseif length(eq_pos) == 3
    plot!(
      fig3,
      t2plot,
      ones(length(t2plot)).*eq_pos[1];
      lw=2,
      linestyle = :dash,
      color = mygreen,
      label ="équilibres stables",
    )
    plot!(
      fig3,
      t2plot,
      ones(length(t2plot)).*eq_pos[2];
      lw=2,
      linestyle = :dash,
      color = myorange,
      label = "",
    )
    plot!(
      fig3,
      t2plot,
      ones(length(t2plot)).*eq_pos[3];
      lw=2,
      linestyle = :dash,
      color = mygreen,
      label = "",
    )
else
    println("Cette situation est non générique (2 équilibres positifs à la bifurcation)")
end

Figure 4: bistabilité dans le modèle de la tordeuse du bourgeon de l’épinette

Diagramme de bifurcations

Warning

Le début de cette partie est assez technique et a pour objectif principal de montrer ce qu’il est possible de faire en manipulant des expressions symboliques.

Finalement nous traçons dans le plan (y, x) le lieu des points d’équilibres en fonction de la taille de la population d’oiseaux:

y = \frac{r}{\alpha x^*}\left(1-\frac{x^*}{K}\right)(h^2+x^{*2}) \tag{4}

Comme nous l’avons vu en cours cette fonction est non monotone avec une branche décroissante, une branche croissante puis à nouveau une branche décroissante, le sens de variation de la branche déterminant la stabilité de l’équilibre correspondant.

Pour déterminer ces différentes branches et les représenter de différentes couleurs de façon à illustrer leur stabilité, nous allons dériver l’Equation 4 par rapport à x^* en utilisant le package Symbolics.jl et chercher les racines de la dérivée³.

³ si l’Equation 4 avait été un polynôme nous aurions utilisé les outils pour les polynômes, mais il s’agit d’une fraction rationnelle

using Symbolics

@variables X
D = Differential(X)

# on définit le lieu des équilibres selon l'équation ci-dessus
Y = r/(α*X)*(1-X/K)*(h^2+X^2)

\begin{equation} \frac{5 \left( 1 - 0.1 X \right) \left( 0.25 + X^{2} \right)}{X} \end{equation}

Les racines de la dérivée, sont les racines du numérateur de la dérivée, donc on dérive:

# Dérivée de Y, distribuée et simplifiée
derY = simplify(expand_derivatives(D(Y)))

\begin{equation} \frac{-1.25 + 5 X^{2} - X^{3}}{X^{2}} \end{equation}

et on récupère ce numérateur:

dnumerator = Symbolics.arguments(Symbolics.value(derY))[1]

\begin{equation} -1.25 + 5.0 X^{2} - X^{3} \end{equation}

Comme ce numérateur est un polynôme, on peut utiliser le polynôme symbolique de Symbolics.jl pour regénérer un polynôme et utiliser la méthode roots() de Polynomials.jl⁴

⁴ il faut l’admettre, ce serait plus simple d’avoir directement une façon de trouver les racines d’un polynôme symbolique, mais ça ne semble pas implémenté pour l’instant (fin 2023)

# on récupère les coefficients X^k du polynôme
coefs_dict = Symbolics.value(dnumerator).dict

# on crée un dictionnaire rassemblant les coefficients du polynôme
poldict = Dict(Symbolics.degree(first(kv)) => kv[2] for kv in coefs_dict)

# on rajoute dans le dictionnaire le coefficient constant de dnumerator
poldict[0] = substitute(dnumerator, Dict( X=> 0))

# on définit le polynôme à partir du dictionnaire
dnumpoly = SparsePolynomial(poldict, :X)

Finalement, on calcule les solutions (racines de la dérivée de y) en prenant soin de filtrer celles entre 0 et K.

# on calcule les solutions en filtrant les racines entre 0 et K via une fonction anonyme s-> K > s > 0
dyroots = filter(s -> K > s > 0, roots(dnumpoly))

Note

Ici nous utilisons une fonction anonyme s -> K > s > 0 qui renvoit true ou false selon que la condition est vérifiée par l’argument, et qui appliquée à roots(dnumpoly), “filtre” les valeurs vérifiant la condition.

On aurait pu faire un filtrage comme vu plus haut.

On calcule les branches du lieu des équilibres correspondantes, et on peut les tracer dans le plan (y,x) pour illustrer le diagramme de bifurcations.

# on calcule chacune des branches
xplot1 = 0.08:.01:dyroots[1]
xplot2 = dyroots[1]:.01:dyroots[2]
xplot3 = dyroots[2]:.01:K
yeq1 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot1]
yeq2 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot2]
yeq3 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot3]

# diagramme de bifurcations
# branche stable 1 et initialisation
figbif = plot(
      yeq1,
      xplot1;
      linewidth = 2,
      color = mygreen,
      label = "équilibres stables",
      legend = :left,
      xlabel = "population d'oiseaux \$y\$",
      ylabel = "population de tordeuses \$x\$",
      title = "Diagramme de bifurcations pour le modèle de tordeuses",
      margin = .5Plots.cm,
      topmargin = 1Plots.cm,
)

# branche instable
plot!(
    figbif,
    yeq2,
    xplot2;
    linewidth = 2,
    color = myred,
    label = "équilibres instables",
)

# branche stable 2
plot!(
    figbif,
    yeq3,
    xplot3;
    linewidth = 2,
    color = mygreen,
    label = "",
)

# équilibre nul
plot!(
    figbif,
    [0, maximum(yeq1)],
    [0, 0];
    color = myred,
    lw = 2,
    label = "",
)

Figure 5: bifurcations dans le modèle de la tordeuse du bourgeon de l’épinette

Nous verrons sur la page des populations en interactions le cas où la population d’oiseaux varie lentement au cours du temps.

References

Ludwig, D., D. D. Jones, and C. S. Holling. 1978. “Qualitative Analysis of Insect Outbreak Systems: The Spruce Budworm and Forest.” Journal of Animal Ecology 47: 315–32.

Reuse

https://creativecommons.org/licenses/by-nc/4.0/

--- title: "Populations exploitées" ---  ```{julia} #| include: false #| eval: true # # this cell does not appear in the rendering, but is executed # # for reproducibility purposes, we load the local julia environment/project, containing # [13f3f980] CairoMakie v0.11.6 # [a93c6f00] DataFrames v1.6.1 # [864edb3b] DataStructures v0.18.16 # [0c46a032] DifferentialEquations v7.12.0 # [e9467ef8] GLMakie v0.9.6 # [91a5bcdd] Plots v1.40.0 # [f27b6e38] Polynomials v4.0.6 # [90137ffa] StaticArrays v1.9.1 # [0c5d862f] Symbolics v5.16.1 # # to share the environment, copy Project.toml and Manifest.toml files in some directory # `activate` the local environment # if necessary `instantiate` to get the correct package versions # # to check if some PackageName.jl is used from the local environment ## Pkg.status("PackageName") using Pkg Pkg.activate(".") ``` ## Prélèvements et effets Allee {#sec-allee-prelev} Nous reprenons le modèle précédent sur l'effet Allee mais en prenant en compte des prélèvements externes avec un effort (taux) de prélèvement $E$ : $$ \dot x = r x \left(\frac{x}{\epsilon}-1\right)\left(1-\frac{x}{K}\right)-Ex. $$ {#eq-alleeprelev} La simulation de ce modèle pour différentes valeurs de $E$ (par exemple $E=0.2$ ou $E=0.85$) ne présente aucune difficulté supplémentaire. Nous allons maintenant nous intéresser à une situation où l'effort de prélèvement $E$ varie au cours du temps entre une valeur soutenable $E_s$ (par exemple $0.2$), et une valeur excessive $E_x$ (par exemple $0.85$). L'attendu théorique est que si les prélèvements sont maintenus à une valeur excessive $E_x$ trop longtemps, la population disparait irrémédiablement même si les prélèvements sont par la suite ramenés à une valeur initialement soutenable $E_s$. ### Prélèvements variables dans le temps ```{julia} using DifferentialEquations using Plots using DataFrames ``` Définissons une fonction `effort()` dépendant du temps, qui renvoit : - $E_s~$ si $~t<T_s$ - $E_x~$ si $~t\in [T_s+T_x[$ - à nouveau $E_s~$ si $~t\geq T_s+T_x$ ```{julia} # paramètres E_s = 0.2 E_x = 0.85 T_s = 10.0 T_x = 9.0 p_effort = [E_s, E_x, T_s, T_x] function effort(t, p) E_s, E_x, T_s, T_x = p if t < T_s || t >= T_s + T_x return E_s elseif t >= T_s && t < T_s + T_x return E_x end end ``` ::: {.callout-note} Remarquez les opérateurs booléens "ou" `||` et "et" `&&`, ainsi que la structure du test `if`. ::: La fonction correspond bien à nos hypothèses : ```{julia} #| output: true #| code-fold: true t2plot = 0: .1: 30 plot( t2plot, # abscisses [effort(t, p_effort) for t in t2plot]; # ordonnées palette = :tab10, linewidth = 2, label = "\$E(t)\$", xlabel = "temps", ylabel = "\$E(t)\$", title = "Effort de pêche \$E(t)\$", margin = .5Plots.cm, topmargin = 1Plots.cm, ) ``` ::: {.callout-note} Remarquez la construction des ordonnées correspondant à la fonction `effort` via une compréhension de liste : `[effort(t, p_effort) for t in t2plot]`. ::: ### Simulations en fonction du temps Nous définissons les paramètres du modèle, la condition initiale et le temps d'intégration : ```{julia} #| code-fold: true x0 = 10.0 tspan = (0.0, 30.0) tstep = 0.1 r = 1.0 K = 10.0 epsilon = 2.0 p_allee = [r, K, epsilon] ``` Nous définissons le système dynamique comme précédemment, à la différence que nous prévoyons de surcharger l'argument `p` dans le problème d'intégration sous forme d'un vecteur comprennant : - le vecteur de paramètres `p_allee` en première position - le vecteur de paramètres `p_effort` en seconde position - et la fonction `effort` en troisième position^[En Julia tout est passable en argument, des variables aux fonctions et plus. ] ```{julia} p_evar = [p_allee, p_effort, effort] ``` ```{julia} function allee_evar(u, p, t) r, K, epsilon = p[1] # unpacking model parameters p_effort = p[2] # unpacking fishing effort parameters E = p[3] # unpacking fishing effort function x = u # use x notation return dx = r*x*(x/epsilon - 1)*(1 - x/K) - E(t, p_effort)*x end ``` L'intégration en elle-même suit le shéma vu précédemment, si ce n'est que l'argument de paramètres doit bien refléter ce qui est attendu par la fonction `allee_evar()`. La simulation en elle-même est effectuée avec une modification de la précision relative de l'intégration `reltol = 1e-6`, la précision par défaut n'étant pas suffisante ici. ```{julia} prob_allee_evar = ODEProblem( allee_evar, x0, tspan, p_evar; saveat = tstep, ) sol_allee_evar = solve(prob_allee_evar; reltol = 1e-6) sol_allee_evar = DataFrame(sol_allee_evar) rename!(sol_allee_evar, :timestamp => :time, :value => :x) ``` ::: {.callout-note} Ici nous avons imposé une tolérance plus précise au solver via le keyword argument `reltol = 1e-6`. ::: Finalement, nous pouvons représenter graphiquement la solution contre le temps. Ici, malgré la perturbation violente induite par la période de surexploitation, le pêcherie retrouve une situation soutenable après un retour à $E=E_s$. ```{julia} #| output: true #| code-fold: true #| label: fig-allee_pers #| fig-cap: Persistence dans le modèle avec effet Allee et prélèvements plot( sol_allee_evar.time, sol_allee_evar.x; label = "\$x(t)\$", linewidth = 2, xlabel = "temps", ylabel = "densité de population \$x(t)\$", title = "Effort de pêche variant dans le temps", palette = :tab10, legend = :left, margin = .5Plots.cm, topmargin = 1Plots.cm, ) ``` Lorsque la période de surexploitation de la population est trop longue (e.g. ici $T_x=9.2$), la population ne parvient pas à récupérer malgré le retour à un effort de prélèvement soutenable. ```{julia} #| output: true #| code-fold: true #| label: fig-allee_ext #| fig-cap: Extinction dans le modèle avec effet Allee et prélèvements # change parameter T_x2 = 9.2 p_effort2 = [E_s, E_x, T_s, T_x2] p_evar2 = [p_allee, p_effort2, effort] # define new problem and integrate prob_allee_evar2 = ODEProblem( allee_evar, x0, tspan, p_evar2, saveat = tstep ) sol_allee_evar2 = solve(prob_allee_evar2, reltol = 1e-6) sol_allee_evar2 = DataFrame(sol_allee_evar2) rename!(sol_allee_evar2, :timestamp => :time, :value => :x) # plot plot( sol_allee_evar2.time, sol_allee_evar2.x, label = "\$x(t)\$"; linewidth = 2, xlabel = "temps", ylabel = "densité de population \$x(t)\$", title = "Effort de pêche variant dans le temps", palette = :tab10, legend = :left, margin = .5Plots.cm, topmargin = 1Plots.cm, ) ``` ### Simulations dans l'espace $(E, x)$ Il s'agit ici de représenter dans l'espace $(E, x)$ l'évolution conjointe de l'effort de pêche et de la densité de la population au cours du temps, afin de mieux comprendre le phénomène d'extinction. Dans cet objectif, on tracera aussi le lieu des équilibres de la population $x^*(E)$ en fonction d'une valeur de $E$ constante, i.e. : $$ x^*(E) = 0, $$ ce qui est équivalent à : $$ E = r\left(\frac{x^*(E)}{\epsilon}-1\right)\left(1-\frac{x^*(E)}{K}\right). $$ Nous allons tracer les deux situations vues à la section précédente dans deux sous figures. Commençons par les lieux des équilibres, qui sont les mêmes dans les deux situations. Définissons le lieu des équilibres positifs : ```{julia} # parabole pour les equilibres positifs function eeqpos(x, par = p_allee) r, K, epsilon = par return r*(x/epsilon -1)*(1-x/K) end ``` Puis traçons ce lieu dans le plan $(E,x)$ (il y a deux branches une stable et une instable pour les équilibres positifs qui se rejoignent en $x^* = (K+\epsilon)/2$, nous les traçons séparément) ```{julia} #| output: true # vecteurs pour le tracé e2plot = 0:.1:1 x2plot1 = epsilon:.02:(K+epsilon)/2 x2plot2 = (K+epsilon)/2:.02:K # on définit des couleurs spécifiques depuis la palette :pal10 mygreen = palette(:tab10)[3] myorange = palette(:tab10)[2] myblue = palette(:tab10)[1] myred = palette(:tab10)[4] # plot fig1 = plot( e2plot, zeros(length(e2plot)); color = mygreen, linewidth = 2, label = "équilibres stables", ylabel = "densité de population \$x\$", xlabel = "effort de pêche \$E\$", legend = :left, ) plot!( fig1, eeqpos.(x2plot1), x2plot1; color= myorange, linewidth = 2, label = "équilibre instable", ) plot!( fig1, eeqpos.(x2plot2), x2plot2; color= mygreen, linewidth = 2, label ="", ) display(fig1) ``` ::: {.callout-note} Remarquez le broadcast (calcul terme à terme) sur les ordonnées des équilibres positifs en utilisant `eeqpos.`. On aurait pu procéder par compréhension de liste: `[eeqpos(x) for x in x2plot2]`. ::: On prépare un graphique `fig2` avec les mêmes éléments que dans `fig1` : ```{julia} fig2 = deepcopy(fig1) # on fait une copie complète de fig1 ``` Finalement, on complète `fig1` et `fig2` avec quelques annotations et les trajectoires calculées plus haut en fonction de l'effort de pêche variable au cours du temps, et on trace les résultats en deux sous-figures : ```{julia} #| output: true #| label: fig-alleebif #| fig-cap: bifurcation pli et catastrophe dans le modèle avec effet Allee et prélèvements # annotations annotate!(fig1, .35, 8.5, Plots.text("branche stable", 10, rotation=-28)) annotate!(fig1, .4, 2.7, Plots.text("branche instable", 10, rotation=28)) annotate!(fig2, .4, 2.7, Plots.text("branche instable", 10, rotation=28)) annotate!(fig2, .35, 8.5, Plots.text("branche stable", 10, rotation=-28)) # peche durable sur fig1 plot!( fig1, [effort(t, p_effort) for t in sol_allee_evar.time], sol_allee_evar.x; color = myblue, linewidth = 2, label = "trajectoire", ) # surpeche sur fig2 plot!( fig2, [effort(t, p_effort2) for t in sol_allee_evar2.time], sol_allee_evar2.x; color = myblue, linewidth = 2, label = "trajectoire", ) # figure reprenant les deux sous figures fig1 et fig2 plot( fig1, fig2; suptitle = "Effets Allee et prélèvements", margin = .5Plots.cm, topmargin = 1Plots.cm, ) ``` ## La tordeuse du bourgeon de l'épinette ### Modèle Nous considérons le modèle de dynamique de populations suivant, inspiré de @ludwig1978 : $$ \dot x =rx\left(1-\frac{x}{K}\right) - \frac{\alpha x^2}{h^2+x^2}\ y, $$ {#eq-tordeuse} avec $x$ la densité de tordeuses et $y$ la densité d'oiseaux. La croissance des tordeuses suit une loi logistique et la prédation des oiseaux une réponse fonctionnelle de type Holling III. ### Simulations On procède classiquement, en ajustant les valeurs de paramètres pour mettre en évidence les phénomènes dynamiques attendus : ```{julia} # paramètres r = 5.0 # natalité K = 10.0 # mortalité α = 1.0 # taux max de prédation, \alpha + TAB h = 0.5 # constante de demi-saturation yc = 7.0 # densité de prédateurs par_tordeuse = [r, K, α, h, yc] # temps d'intégration tspan = (0.0, 3.0) tstep = 0.02 ``` Puis on définit le modèle : ```{julia} function tordeuse(u, p, t) r, K, α, h, yc = p x = u return dx = r*x*(1 - x/K) - α*x^2/(h^2 + x^2)*yc end ``` Comme pour le modèle avec effets Allee, on définit une fonction qui redéfinit le problème d'intégration à chaque fois, et simule et renvoit la solution pour pouvoir illustrer la bi-stabilité^[des méthodes alternative utilisant l'interface integrator de `DifferentialEquations.jl` ou la redéfinition du problème d'intégration via `remake` sont proposées en [annexe](annexe_integrator.qmd)] : ```{julia} function int_tordeuse(x0, tspan = tspan, param = par_tordeuse) prob_tordeuse = ODEProblem( tordeuse, # modèle x0, # condition initiale tspan, # tspan param; # paramètres saveat = tstep, ) sol_tordeuse = solve(prob_tordeuse) sol_tordeuse = DataFrame(sol_tordeuse) rename!(sol_tordeuse, :timestamp => :time, :value => :x) end ``` Nous simulons le modèle depuis différentes conditions intiales, et traçons les résultats via une boucle `for`. ```{julia} #| output: true # conditions initiales x0step = 1.35 x0vec = x0step:x0step:K # custom color palette init_cgrad = palette([:steelblue, :lightblue], length(x0vec)) # initialisation du graphique, équilibre nul fig3 = plot(; # seulement des kwargs: on initialise le graphique palette = init_cgrad, legend = :right, label ="équilibres instables", title = "Tordeuse du bourgeon de l\'épinette", ylabel = "densité de population \$x(t)\$", xlabel = "temps \$t\$", margin = .5Plots.cm, topmargin = 1Plots.cm, ) # boucle de plot avec intégration pour differentes conditions initiales for x0 in x0vec plot!( fig3, int_tordeuse(x0).time, int_tordeuse(x0).x, linewidth = 2, label = "", ) end display(fig3) # actually shows the plot fig3 ``` ### Equilibres Les valeurs des équilibres positifs du modèle @eq-tordeuse n'ont pas d'écriture mathématique simple. Nous allons calculer numériquement les racines du polynôme dont ils sont solution : $$ r\left(1-\frac{x^*}{K}\right)\left(h^2+x^{*2}\right)-\alpha x^* y = 0 $$ {#eq-eq-tordeuse} Pour cela, on utilise le package `Polynomials.jl` : ```{julia} #| output: true using Polynomials # définition du monôme X X = Polynomial([0, 1]) # définition du polynôme pol = r*(1-X/K)*(h^2 + X^2)-α*X*yc # calcul des racines, réelles, positives et plus petites que K eq_pos = roots(pol) # calcul des racines eq_pos = real.(eq_pos[isreal.(eq_pos)]) # filtrage des racines réelles eq_pos = eq_pos[(eq_pos .> 0) .& (eq_pos .<= K)] # filtrage des racines >0 et <K ``` ::: {.callout-note} A partir de la définition du monôme de degré 1, on écrit assez naturellement un polynôme et on en extrait les racines via `roots`. La partie filtrage des racines réelles prend la partie réelle des racines qui sont rendues sous forme d'un type complexe (même si la partie imaginaire est nulle). Remarquez les broadcasting divers, notamment `.>` et `.&` qui testent les conditions terme à terme sur le vecteur `eq_pos`. ::: Et on trace les différents équilibres : ```{julia} #| output: true #| label: fig-tordeuse-temps #| fig-cap: bistabilité dans le modèle de la tordeuse du bourgeon de l'épinette t2plot = collect(tspan) # initialisation du graphique, équilibre nul plot!( fig3, t2plot, zeros(length(t2plot)); linewidth = 2, linestyle = :dash, color = myorange, palette = init_cgrad, legend = :right, label ="équilibres instables", ylabel = "densité de population \$x(t)\$", xlabel = "temps \$t\$", margin = .5Plots.cm, topmargin = 1Plots.cm, ) # équilibres positifs, différents cas if length(eq_pos) == 1 plot!( fig3, t2plot, ones(length(t2plot)).*eq_pos; color = mygreen, label ="équilibre stable", ) elseif length(eq_pos) == 3 plot!( fig3, t2plot, ones(length(t2plot)).*eq_pos[1]; lw=2, linestyle = :dash, color = mygreen, label ="équilibres stables", ) plot!( fig3, t2plot, ones(length(t2plot)).*eq_pos[2]; lw=2, linestyle = :dash, color = myorange, label = "", ) plot!( fig3, t2plot, ones(length(t2plot)).*eq_pos[3]; lw=2, linestyle = :dash, color = mygreen, label = "", ) else println("Cette situation est non générique (2 équilibres positifs à la bifurcation)") end ``` ### Diagramme de bifurcations ::: {.callout-warning} Le début de cette partie est assez technique et a pour objectif principal de montrer ce qu'il est possible de faire en manipulant des expressions symboliques. ::: Finalement nous traçons dans le plan $(y, x)$ le lieu des points d'équilibres en fonction de la taille de la population d'oiseaux: $$ y = \frac{r}{\alpha x^*}\left(1-\frac{x^*}{K}\right)(h^2+x^{*2}) $${#eq-lieutordeuse} Comme nous l'avons vu en cours cette fonction est non monotone avec une branche décroissante, une branche croissante puis à nouveau une branche décroissante, le sens de variation de la branche déterminant la stabilité de l'équilibre correspondant. Pour déterminer ces différentes branches et les représenter de différentes couleurs de façon à illustrer leur stabilité, nous allons dériver l'@eq-lieutordeuse par rapport à $x^*$ en utilisant le package `Symbolics.jl` et chercher les racines de la dérivée^[si l'@eq-lieutordeuse avait été un polynôme nous aurions utilisé les outils pour les polynômes, mais il s'agit d'une fraction rationnelle]. ```{julia} #| output: true using Symbolics @variables X D = Differential(X) # on définit le lieu des équilibres selon l'équation ci-dessus Y = r/(α*X)*(1-X/K)*(h^2+X^2) ``` Les racines de la dérivée, sont les racines du numérateur de la dérivée, donc on dérive: ```{julia} #| output: true # Dérivée de Y, distribuée et simplifiée derY = simplify(expand_derivatives(D(Y))) ``` \ et on récupère ce numérateur: ```{julia} #| output: true dnumerator = Symbolics.arguments(Symbolics.value(derY))[1] ``` \ Comme ce numérateur est un polynôme, on peut utiliser le polynôme symbolique de `Symbolics.jl` pour regénérer un polynôme et utiliser la méthode `roots()` de `Polynomials.jl`^[il faut l'admettre, ce serait plus simple d'avoir directement une façon de trouver les racines d'un polynôme symbolique, mais ça ne semble pas implémenté pour l'instant (fin 2023)] ```{julia} # on récupère les coefficients X^k du polynôme coefs_dict = Symbolics.value(dnumerator).dict # on crée un dictionnaire rassemblant les coefficients du polynôme poldict = Dict(Symbolics.degree(first(kv)) => kv[2] for kv in coefs_dict) # on rajoute dans le dictionnaire le coefficient constant de dnumerator poldict[0] = substitute(dnumerator, Dict( X=> 0)) # on définit le polynôme à partir du dictionnaire dnumpoly = SparsePolynomial(poldict, :X) ``` Finalement, on calcule les solutions (racines de la dérivée de $y$) en prenant soin de filtrer celles entre $0$ et $K$. ```{julia} # on calcule les solutions en filtrant les racines entre 0 et K via une fonction anonyme s-> K > s > 0 dyroots = filter(s -> K > s > 0, roots(dnumpoly)) ``` ::: {.callout-note} Ici nous utilisons une fonction anonyme `s -> K > s > 0` qui renvoit `true` ou `false` selon que la condition est vérifiée par l'argument, et qui appliquée à `roots(dnumpoly)`, "filtre" les valeurs vérifiant la condition. On aurait pu faire un filtrage comme vu plus haut. ::: On calcule les branches du lieu des équilibres correspondantes, et on peut les tracer dans le plan $(y,x)$ pour illustrer le diagramme de bifurcations. ```{julia} #| output: true #| label: fig-tordeuse-bif #| fig-cap: bifurcations dans le modèle de la tordeuse du bourgeon de l'épinette # on calcule chacune des branches xplot1 = 0.08:.01:dyroots[1] xplot2 = dyroots[1]:.01:dyroots[2] xplot3 = dyroots[2]:.01:K yeq1 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot1] yeq2 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot2] yeq3 = [r*(1 - x/K)/(α*x)*(h^2 + x^2) for x in xplot3] # diagramme de bifurcations # branche stable 1 et initialisation figbif = plot( yeq1, xplot1; linewidth = 2, color = mygreen, label = "équilibres stables", legend = :left, xlabel = "population d'oiseaux \$y\$", ylabel = "population de tordeuses \$x\$", title = "Diagramme de bifurcations pour le modèle de tordeuses", margin = .5Plots.cm, topmargin = 1Plots.cm, ) # branche instable plot!( figbif, yeq2, xplot2; linewidth = 2, color = myred, label = "équilibres instables", ) # branche stable 2 plot!( figbif, yeq3, xplot3; linewidth = 2, color = mygreen, label = "", ) # équilibre nul plot!( figbif, [0, maximum(yeq1)], [0, 0]; color = myred, lw = 2, label = "", ) ``` Nous verrons sur la page des [populations en interactions](pop_interactions.qmd) le cas où la population d'oiseaux varie lentement au cours du temps.