Quiver plot // julia // population dynamics

quiver_plot_popdyn.jl

using Plots
pyplot()

"""
Derivatives
prey, predator, growth, comp, conversion, predation, mortality
"""
function ṅ(x, y, r, q, α, β, δ)
   ẋ = r*x - q*x^2 - β*y*x
   ẏ = α*β*y*x - δ*y
   return (ẋ, ẏ)
end

function sim(x, y, r, q, α, β, δ; t::Int64=1000, scal=0.01)
   steps = 0:scal:t
   X = zeros(Float64, size(steps))
   X[1] = x
   Y = zeros(Float64, size(steps))
   Y[1] = x
   for i in eachindex(steps[2:end])
      N = ṅ(X[i], Y[i], r, q, α, β, δ)
      X[i+1] = X[i]+N[1]*scal
      Y[i+1] = Y[i]+N[2]*scal
   end
   return hcat(collect(steps), X, Y)
end

r, q, α, β, δ = 0.4, 0.02, 0.12, 0.1, 0.02
S = sim(1.0, 1.0, r, q, α, β, δ)

xi = linspace(0, maximum(S[:,2]), 30)
yi = linspace(0, maximum(S[:,3]), 30)
vxi = repeat(xi, inner=length(yi))
vyi = repeat(yi, outer=length(xi))
N = hcat(vxi, vyi)
function dquiv(X)
   return collect(ṅ(X[1], X[2], r, q, α, β, δ)).*0.35
end
D = mapslices(dquiv, N, 2)
qu = hcat(N, D)
vdn = reshape(abs.(qu[:,3]).+abs.(qu[:,4]), (length(xi), length(yi)))

heatmap(xi, yi, vdn, c=:YlGnBu, xlim=[minimum(xi), maximum(xi)], ylim=[minimum(yi), maximum(yi)], lab="")
quiver!(qu[:,1], qu[:,2], quiver=(qu[:,3], qu[:,4]), c=:grey)
plot!(S[:,2], S[:,3], linewidth=2, c=:black, lab="")
xaxis!("Prey population")
yaxis!("Predator population")