Nous abordons ici brièvement les méthodes integrator et remake de DifferentialEquations.jl pour les simulations multiples.
Interface Integrator pour simulations multiples
Lors de simulations multiples1, la stratégie choisie a été de redéfinir un nouveau problème d’intégration ODEProblem pour chaque condition initiale. Cette procédure est a priori peu efficace, et DifferentialEquations.jl permet de modifier un problème d’intégration via la Integrator Interface.
Pour une utilisation dans une boucle pour simuler les trajectoires depuis différentes conditions initiales, nous définissons une fonction qui réinitialise l’intégrateur à la nouvelle condition, effectue la simulation (et transforme la solution en dataframe).
functionint_tordeuse(x0, integrator)reinit!(integrator, x0) # la clef de l'interface est ici sol_tordeuse =solve!(integrator) # et là sol_tordeuse =DataFrame(sol_tordeuse)rename!(sol_tordeuse, :timestamp =>:time, :value =>:x)return sol_tordeuseend
On construit ensuite le graphique et on fait une boucle pour intégrer et tracer les résultats:
# conditions initialesx0vec = x0step:x0step:K# custom color paletteinit_cgrad =palette([:steelblue, :lightblue], length(x0vec))# initialisation du graphique, équilibre nulfig =plot(; 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 initialesfor x0 in x0vecplot!( fig,int_tordeuse(x0, integrator).time,int_tordeuse(x0, integrator).x, linewidth =2, label ="", )enddisplay(fig) # actually shows the plot P
On peut compléter la figure comme vu précédemment, mais cela n’a pas plus d’intéret ici.
Commentaires
Le code de la fonction int_tordeuse est plus simple avec l’integrator interface qu’en redéfinissant le problème d’intégration à chaque fois. Néanmoins, sur un problème avec x0step= 0.1 (soit ~100 trajectoires calculées), le gain en temps de calcul de cette méthode est marginal : 123 ms contre 128 ms avec la méthode précédente (test avec package BenchmarkTools.jl).
Le gain en mémoire est plus sensible : 4.28 MiB contre 7.26 MiB, mais finalement marginal aussi sur des machines modernes pour de si petits problèmes d’intégration.
Redéfinition du problème d’intégration via remake
Une autre approche plutôt que de redéfinir un nouveau problème d’intégration à chaque fois, est de modifier le problème d’intégration via la fonction remake.
On définit une fonction qui prend pour arguments la condition initiale et un problème d’intégration, et redéfinit le problème d’intégration depuis cette condition initiale :
functionint_tordeuse2(x0, prob) prob =remake(prob, u0 = x0) # redéfinition du problème d'intégration sol_tordeuse =solve(prob) sol_tordeuse =DataFrame(sol_tordeuse)rename!(sol_tordeuse, :timestamp =>:time, :value =>:x)end# initialisation du graphiquefig2 =plot(; 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,)for x0 in x0vecplot!( fig2,int_tordeuse2(x0, prob_tordeuse).time,int_tordeuse2(x0, prob_tordeuse).x; linewidth =2, label ="", )enddisplay(fig2)
Cette méthode est plus rapide que la redéfinition de problème d’intégration ou l’utilisation de l’integrator interface (113 ms dans les conditions vues plus haut) et intermédiaire en termes de mémoire utilisée (6.39 MiB). Le code est aussi simple que pour l’integrator interface.
Par ailleurs la même approche via remake peut être utilisée pour redéfinir un problème d’intégration en modifiant les paramètres plutôt que la condition initiale, par exemple pour calculer un diagramme de bifurcations par force brute comme pour le modèle de Rosenzweig MacArthur. C’est visiblement la bonne solution algorithmique.
---title: "Simulations multiples"---<!-- load local julia environment (freeze package versions) -->```{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 PkgPkg.activate(".")```Nous abordons ici brièvement les méthodes `integrator` et `remake` de `DifferentialEquations.jl` pour les simulations multiples.## Interface Integrator pour simulations multiplesLors de simulations multiples^[par exemple pour plusieurs conditions initiales pour le modèle avec [effets Allee](pop_isolees.qmd) ou celui de la [tordeuse du bourgeon de l'épinette](pop_exploitees.qmd)], la stratégie choisie a été de redéfinir un nouveau problème d'intégration `ODEProblem` pour chaque condition initiale. Cette procédure est a priori peu efficace, et `DifferentialEquations.jl` permet de modifier un problème d'intégration via la [Integrator Interface](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/).Répliquons la figure des dynamiques de la tordeuse du bourgeon de l'épinette avec cette méthode. On commence par définir les paramètres et le modèle:```{julia}#| code-fold: trueusing DifferentialEquationsusing Plotsusing DataFrames# paramètresr = 5.0# natalitéK = 10.0# mortalitéα = 1.0# taux max de prédationh = 0.5# constante de demi-saturationyc = 7.0# densité de prédateurspar_tordeuse = [r, K, α, h, yc]# temps d'intégrationtspan = (0.0, 3.0)tstep = 0.02# condition initalex0step = 1.35# modèlefunction tordeuse(u, p, t) r, K, α, h, yc = p x = u[1]return dx = r*x*(1 - x/K) - α*x^2/(h^2 + x^2)*ycend```Nous définissons ensuite le problème d'intégration ainsi que l'intégrateur `integrator`:```{julia}prob_tordeuse = ODEProblem( tordeuse, # modèle x0step, # condition initiale tspan, # tspan par_tordeuse; # paramètres saveat = tstep,) # option de sortieintegrator = init(prob_tordeuse, Tsit5()) # integrator requires integ algorithm```Pour une utilisation dans une boucle pour simuler les trajectoires depuis différentes conditions initiales, nous définissons une fonction qui réinitialise l'intégrateur à la nouvelle condition, effectue la simulation (et transforme la solution en dataframe).```{julia}function int_tordeuse(x0, integrator) reinit!(integrator, x0) # la clef de l'interface est ici sol_tordeuse = solve!(integrator) # et là sol_tordeuse = DataFrame(sol_tordeuse)rename!(sol_tordeuse, :timestamp => :time, :value => :x)return sol_tordeuseend```On construit ensuite le graphique et on fait une boucle pour intégrer et tracer les résultats:```{julia}#| output: true# conditions initialesx0vec = x0step:x0step:K# custom color paletteinit_cgrad = palette([:steelblue, :lightblue], length(x0vec))# initialisation du graphique, équilibre nulfig = plot(; 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 initialesfor x0 in x0vec plot!( fig, int_tordeuse(x0, integrator).time, int_tordeuse(x0, integrator).x, linewidth = 2, label = "", )enddisplay(fig) # actually shows the plot P```\On peut compléter la figure comme vu [précédemment](pop_exploitees.qmd), mais cela n'a pas plus d'intéret ici.### CommentairesLe code de la fonction `int_tordeuse` est plus simple avec l'integrator interface qu'en redéfinissant le problème d'intégration à chaque fois. Néanmoins, sur un problème avec `x0step= 0.1` (soit ~100 trajectoires calculées), le gain en temps de calcul de cette méthode est marginal : 123 ms contre 128 ms avec la méthode précédente (test avec package `BenchmarkTools.jl`).Le gain en mémoire est plus sensible : 4.28 MiB contre 7.26 MiB, mais finalement marginal aussi sur des machines modernes pour de si petits problèmes d'intégration.### Redéfinition du problème d'intégration via `remake`Une autre approche plutôt que de redéfinir un nouveau problème d'intégration à chaque fois, est de modifier le problème d'intégration via la fonction `remake`.On définit une fonction qui prend pour arguments la condition initiale et un problème d'intégration, et redéfinit le problème d'intégration depuis cette condition initiale :```{julia}#| output: truefunction int_tordeuse2(x0, prob) prob = remake(prob, u0 = x0) # redéfinition du problème d'intégration sol_tordeuse = solve(prob) sol_tordeuse = DataFrame(sol_tordeuse)rename!(sol_tordeuse, :timestamp => :time, :value => :x)end# initialisation du graphiquefig2 = plot(; 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,)for x0 in x0vec plot!( fig2, int_tordeuse2(x0, prob_tordeuse).time, int_tordeuse2(x0, prob_tordeuse).x; linewidth = 2, label = "", )enddisplay(fig2)```\Cette méthode est plus rapide que la redéfinition de problème d'intégration ou l'utilisation de l'integrator interface (113 ms dans les conditions vues plus haut) et intermédiaire en termes de mémoire utilisée (6.39 MiB). Le code est aussi simple que pour l'integrator interface.Par ailleurs la même approche via `remake` peut être utilisée pour redéfinir un problème d'intégration en modifiant les paramètres plutôt que la condition initiale, par exemple pour calculer un diagramme de bifurcations par force brute comme pour le modèle de [Rosenzweig MacArthur](pop_interactions2.qmd). C'est visiblement la bonne solution algorithmique.