DifferentialEquations.jl

DifferentialEquations.jl copied to clipboard

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

This is because va[i] is the ith array instead of the ith value, so its length handling is not fully like a matrix. It would be nice to add a type parameter to change that.