Hi, I’m new to Julia and trying to find fixed points of an s-species generalized Lotka-Volterra model using DynamicalSystems. My model is
d_t\vec{u}=\vec{\lambda}+\vec{r}\cdot\vec{u}\cdot (1-\hat{\alpha}\vec{u})
where \hat{\alpha}_{ii}=1 and \hat{\alpha}_{ij,i\neq j} are distributed c\cdot\text{Gamma}(1/\sigma^2, \sigma^2).
I have implemented the dynamical rule as
function interaction_matrix(σ, s, c)
α = c*rand(Gamma(1/σ^2,σ^2), (s,s))
α[diagind(α)] .= 1.0
return α
end
function GeneralizedLotkaVolterra!(du, u, p, t)
λ, r, α = p
du = λ+r.*u.*(1 .-α*u)
return nothing
end
such that system initialization
SIGMA = 0.5
SEED = 1234
DISPERSAL = 1.0e-5 #scalar lambda
GROWTH_RATE = 1.0 #scalar r
N_SPECIES = 3 #s
IJ_STRENGTH = 0.25 #c
λ = DISPERSAL*ones(N_SPECIES)
r = GROWTH_RATE*ones(N_SPECIES)
α = interaction_matrix(SIGMA, N_SPECIES, IJ_STRENGTH)
p = [λ, r, α]
u0 = rand(N_SPECIES)
xg = range(0, 1; length=101)
grid = ntuple(x->xg, N_SPECIES)
glv = CoupledODEs(GeneralizedLotkaVolterra!, u0, p)
yields
3-dimensional CoupledODEs
deterministic: true
discrete time: false
in-place: true
dynamic rule: GeneralizedLotkaVolterra!
ODE solver: Tsit5
ODE kwargs: (abstol = 1.0e-6, reltol = 1.0e-6)
parameters: Array{Float64}[[1.0e-5, 1.0e-5, 1.0e-5], [1.0, 1.0, 1.0], [1.0 0.17264272321853863 0.24272781912293961; 0.433123579186088 1.0 0.31834434917230875; 0.364142261750813 0.08754102121231638 1.0]]
time: 0.0
state: [0.8316040500722955, 0.9020798266115233, 0.06725925088948714]
which is what I expect. However, if I then run
X, t = trajectory(glv, 1000)
X
I get that the state has not evolved at all
3-dimensional StateSpaceSet{Float64} with 10001 points
0.831604 0.90208 0.0672593
0.831604 0.90208 0.0672593
0.831604 0.90208 0.0672593
⋮
0.831604 0.90208 0.0672593
0.831604 0.90208 0.0672593
0.831604 0.90208 0.0672593
which seems unreasonable, because directly evaluating du = λ+r.*u0.*(1 .-α*u0) gives the nonzero
3-element Vector{Float64}:
-0.003039707589467017
-0.2558908791421405
0.03706645604209235
suggesting that at least the second point should be different from the first, especially as the average timestep is about 0.1.
I’ve tried a few different ODE solvers and kwargs, and all give me the same results. Where does my code go wrong? Thanks!
