# DynamicalSystems: dynamic rule evaluates to nonzero but trajectory stays at initial condition

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!

repalce du =  with du .=  so that the state is updated in place.

1 Like

also, hi, welcome

1 Like

Thank you! It worked!