Hey there, I am new to the community.
I don’t have too much experience with the DifferentialEquations package.
I am trying to model Lotka Volterra dynamics with two guilds and in these guilds various species.
For this I am using VectorOfArray(). To reduce the problem I am first just implementing the evolution for the first guild (nh). This works so far until I increase the number of species or the upper time limit over a certain threshold. In both cases I get the error:
MethodError: no method matching abs(::Array{Float64,2})
Closest candidates are:
abs(!Matched::Bool) at bool.jl:91
abs(!Matched::Float16) at float.jl:517
abs(!Matched::Float32) at float.jl:518
…
The reduced Code is:
function abundance_(dn,n,p,t)
nh, nu = n
dnh = dn[1]
dnu = dn[2]
for i in 1:N # choose species
competition = 0
for k in 1:N # go through all other species in matrix (beta)
competition += p[1][i,k]*nh[k]
end
dnh[i] = nh[i]*(p[2]-competition)
end
end
#Competition Matrix generation
function beta_generation(N,range)
a = range*rand(N,N)
A = droplower(sparse(a))
beta = Symmetric(Matrix(A))
for i in 1:N
beta[i,i] = 1
end
return beta
end
I don’t really understand why the error pops up for the upper time limit of 65 and N=100 species and below it doesn’t.
Maybe someone can help.
Thanks in advance.
using RecursiveArrayTools, OrdinaryDiffEq, SparseArrays
function abundance_(dn,n,p,t)
nh, nu = n
dnh = dn[1]
dnu = dn[2]
for i in 1:N # choose species
competition = 0
for k in 1:N # go through all other species in matrix (beta)
competition += p[1][i,k]*nh[k]
end
dnh[i] = nh[i]*(p[2]-competition)
end
end
rho = 1
N=100
n₀ = VectorOfArray([1/N*ones(N,1),1/N*ones(N,1)])
range = 0.1
tspan = (0.0,65.0)
function beta_generation(N,range)
a = range*rand(N,N)
A = droplower(sparse(a))
beta = Symmetric(Matrix(A))
for i in 1:N
beta[i,i] = 1
end
return beta
end
beta = beta_generation(N,range)
p = [beta,rho]
prob = ODEProblem(abundance_,n₀,tspan,p)
Please post a reproducer that I can copy paste to see the issue if you want more help. As it stands, I’m not entirely sure what your issue is.
Thanks for the fast reply!
I am very sorry, I forgot that one.
using LinearAlgebra
function droplower(A::SparseMatrixCSC)
m,n = size(A)
rows = rowvals(A)
vals = nonzeros(A)
V = Vector{eltype(A)}()
I = Vector{Int}()
J = Vector{Int}()
for i=1:n
for j in nzrange(A,i)
rows[j]>i && break
push!(I,rows[j])
push!(J,i)
push!(V,vals[j])
end
end
return sparse(I,J,V,m,n)
end
using RecursiveArrayTools, OrdinaryDiffEq, SparseArrays
function abundance_(dn,n,p,t)
nh, nu = n
dnh = dn[1]
dnu = dn[2]
for i in 1:N # choose species
competition = 0
for k in 1:N # go through all other species in matrix (beta)
competition += p[1][i,k]*nh[k]
end
dnh[i] = nh[i]*(p[2]-competition)
end
end
rho = 1
N=100
n₀ = VectorOfArray([1/N*ones(N,1),1/N*ones(N,1)])
range = 0.1
tspan = (0.0,65.0)
function beta_generation(N,range)
a = range*rand(N,N)
A = droplower(sparse(a))
beta = Symmetric(Matrix(A))
for i in 1:N
beta[i,i] = 1
end
return beta
end
using LinearAlgebra
function droplower(A::SparseMatrixCSC)
m,n = size(A)
rows = rowvals(A)
vals = nonzeros(A)
V = Vector{eltype(A)}()
I = Vector{Int}()
J = Vector{Int}()
for i=1:n
for j in nzrange(A,i)
rows[j]>i && break
push!(I,rows[j])
push!(J,i)
push!(V,vals[j])
end
end
return sparse(I,J,V,m,n)
end
beta = beta_generation(N,range)
p = (beta,rho)
prob = ODEProblem(abundance_,n₀,tspan,p)
sol = solve(prob,Tsit5()) # works
sol = solve(prob,Rodas5()) # fails
Non-stiff ODE solvers work, it’s the stiff ones that fail because of how length is handled in VectorOfArray. I’ll open an issue, but Lotka-Volterra should be fine if you specify an explicit method.