Lotka Volterra avec Makie.jl

L’objectif ici est de faire une jolie représentation graphique du modèle de Lotka Volterra avec la librairie graphique Makie.jl.

On reprend une partie du code de la page sur les populations en intéractions.

Code 
using DifferentialEquations

# conditions initiales
x0 = 1.0
y0 = 2.3
etat0 = [x0, y0]

# paramètres
r = 1.0
c = 1.0
b = 1.0
m = 1.0
par_lovo = [r, c, b, m]

# integration plus longue
tspan = (0.0, 30.0)
tstep = .01

# définition du modèle
function lovo(u, par, t)
    r, c, b, m = par
    x = u[1]
    y = u[2]
    dx = r*x - c*x*y
    dy = b*x*y - m*y
    [dx, dy]
end

# problème
prob_lovo = ODEProblem(lovo, etat0, tspan, par_lovo; saveat = tstep)
# intégration
sol_lovo = solve(prob_lovo; reltol = 1e-6)

L’objectif est de créer une figure avec trois panels. dans la colonne de gauche, le premier panel représentera la dynamique temporelle, le second le plan de phase, et dans la colonne de droite nous représenterons la figure 3D avec l’intégrale première.

using CairoMakie

# on crée la figure avec un fond gris clair
fig = Figure(
       backgroundcolor = RGBf(0.98, 0.98, 0.98),
       size = (1000, 600),
)

# déclare les différents panels comme des éléments gridlayout contenant les plots
panela = fig[1, 1] = GridLayout()
panelb = fig[2, 1] = GridLayout()
panelc = fig[:, 2] = GridLayout()

On commence par remplir le premier panel avec les dynamiques :

ax1 = Axis(
       panela[1,1],
       ylabel = L"densités de populations$$",
       title="Dynamiques",
)

# pour avoir un xlabel plus proche de l'axe on le définit séparément en réglant manuellement l'espacement/padding
Label(panela[1,1,Bottom()], L"temps$$"; padding=(0,0,0,20))

# dynamiques de populations
lines!(ax1, sol_lovo.t, sol_lovo[1,:]; linewidth = 2, label = L"$x(t)$")

lines!(ax1, sol_lovo.t, sol_lovo[2,:]; linewidth = 2, label = L"$y(t)$")

# légende avec quelques ajustements d'espacement interne
axislegend(ax1; position = :lt, labelsize = 12, padding = (5,5,0,0), rowgap = -5)

fig

Puis on complète le second panel avec le plan de phase :

ax2 = Axis(
       panelb[1,1];
       xlabel = L"population $x$",
       ylabel = L"population $y$",
       title = "Plan de phase",
       xticks = (0:.5:2),
       yticks = (0:.5:2),
)

# champs de vecteur
scale = 10
xrange = range(0, 2.75, length=11)
yrange = range(0, 2.75, length=11)

derx = [lovo([x y], par_lovo, 0)[1]/scale for x in xrange, y in yrange]
dery = [lovo([x y], par_lovo, 0)[2]/scale for x in xrange, y in yrange]

arrows!(ax2, xrange, yrange, derx, dery; color = :lightgray, arrowsize = 10)

# nullclines
lines!(
    ax2,
    xrange,
    [r/c for x in xrange];
    color = Cycled(2),
    linewidth = 2,
    label = L"$\dot{x}$ nullcline",
)

lines!(ax2, [0 for y in yrange], yrange; linewidth = 2, color = Cycled(2))

lines!(
    ax2,
    [m/b for y in yrange],
    yrange;
    color = Cycled(3),
    linewidth = 2,
    label = L"$\dot{y}$ nullcline",
)

lines!(ax2, xrange, [0 for y in yrange]; color = Cycled(3), linewidth = 2)

# équilibres
scatter!(ax2, Point2f(0,0); color = Cycled(4))
scatter!(ax2, Point2f(m/b, r/c); color = Cycled(4))

# trajectoire
lines!(ax2, sol_lovo[1, :], sol_lovo[2, :], color = Cycled(1), linewidth = 2)

xlims!(ax2, -0.1, 2.75)
ylims!(ax2, -0.1, 2.75)

axislegend(ax2; position = :rt, labelsize = 12, padding = (5,5,0,0), rowgap = -5)

fig

Finalement on complète le troisième panel avec la représentation 3D, en commencant par redéfinir l’intégrale première

# l'intégrale première
function int_prem(x, y, par = par_lovo)
      r, c, b, m = par
      return -r*log(y) + c*y - m*log(x) + b*x
end
# système d'axe 3D en première colonne de panel c, avec réglage de la caméra
ax3 = Axis3(
       panelc[1, 1];
       title = "Intégrale première",
       titlegap = -50,
       xlabel = L"population $x$",
       ylabel = L"population $y$",
       zlabel = L"$H(x,y)$",
       azimuth = 0.5,
       elevation = 0.2,
)

# grille x, y
xsurf = 0.2:0.1:3.0
ysurf = 0.2:0.1:3.0

# calcul de la surface via une compréhension de liste
hsurf = [int_prem(x, y, par_lovo) for x in xsurf, y in ysurf]
hplane = [int_prem(x0, y0, par_lovo) for x in xsurf, y in ysurf]

# tracé de H(x,y) et du plan z = H(x0,y0)
hs = surface!(ax3, xsurf, ysurf, hsurf; alpha=.5)
surface!(ax3, xsurf, ysurf, hplane; color = fill(:red, 100, 100), alpha = .3)

# ajout de la trajectoire simulée
lines!(
    ax3,
    sol_lovo[1,:],
    sol_lovo[2,:],
    [int_prem(x0, y0, par_lovo) for x in sol_lovo[1,:]];
    color = Cycled(4),
    linewidth = 3,
)

# ajout d'une colorbar reprise sur le plot de H(x,y) avec une hauteur relative sur la colonne
# la colorbar est en 2e colonne du panelc
Colorbar(panelc[1, 2], hs, height = Relative(.55))

# colonne 2 plus large
colsize!(fig.layout, 2, Auto(1.5))
# modification des espacements entre lignes et colonnes
colgap!(fig.layout, 60)
colgap!(panelc, -10)
rowgap!(fig.layout, 10)

# titre général en ligne 0
supertitle = Label(fig[0, :], "Modèle de Lotka Volterra", fontsize = 30)

fig

Figure 1: Figure illustrant le comportement du modèle de Lotka Volterra

et voilà !