Simulations multiples

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 multiples¹, 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.

¹ par exemple pour plusieurs conditions initiales pour le modèle avec effets Allee ou celui de la tordeuse du bourgeon de l’épinette

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:

using DifferentialEquations
using Plots
using DataFrames

# paramètres
r = 5.0      # natalité
K = 10.0     # mortalité
α = 1.0      # taux max de prédation
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

# condition initale
x0step = 1.35

# modèle
function 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)*yc
end

Nous définissons ensuite le problème d’intégration ainsi que l’intégrateur integrator:

prob_tordeuse = ODEProblem(
    tordeuse,         # modèle
    x0step,           # condition initiale
    tspan,            # tspan
    par_tordeuse;     # paramètres
    saveat = tstep,
)   # option de sortie

integrator = 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).

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_tordeuse
end

On construit ensuite le graphique et on fait une boucle pour intégrer et tracer les résultats:

# conditions initiales
x0vec = x0step:x0step:K

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

# initialisation du graphique, équilibre nul
fig = 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 initiales
for x0 in x0vec
    plot!(
        fig,
        int_tordeuse(x0, integrator).time,
        int_tordeuse(x0, integrator).x,
        linewidth = 2,
        label = "",
    )
end

display(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 :

function 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 graphique
fig2 = 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 = "",
    )
end

display(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.