Error: Problems with using VectorOfArray() with DifferentialEquations package

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
rho = 1
N=100
n₀ = VectorOfArray([1/N*ones(N,1),1/N*ones(N,1)])
range = 0.1
tspan = (0.0,65.0)

beta = beta_generation(N,range)
p = [beta,rho]
prob = ODEProblem(abundance_,n₀,tspan,p)

where beta is a coupling matrix:

#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.

What is this function? It’s undefined. I am at:

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.

1 Like

Yes, that works. Thanks a lot for that super fast handling of the problem!