Error when creating Dynamical Systems with arrays as parameters

[jarroyoe]

Hello everyone,

Issue description

I am trying to create a Dynamical System with the following equations:

For i = 1,...,N. I need N to be variable, so I let N be determined by the parameters r, K, and A. I pass a vector p as parameter, where each entry is an Array (r,K,A). This gives me the following error:

ERROR: MethodError: no method matching (Vector{T} where T)(::Matrix{Float64})
Closest candidates are:
  (Array{T, N} where T)(::AbstractArray{S, N}) where {S, N} at boot.jl:470
  (Vector{T} where T)() at baseext.jl:27

Expected behavior

I expect the dynamical system to be created. This behavior works correctly when I just use DifferentialEquations, not DynamicalSystems.

MWE

function lv(dN,N,p,n)
    r = p[1]; K = p[2]; A = p[3]
    S = length(r)
    #NOTE: Have to force A[i,i]=0 for all i
    dN = [r[i] * N[i] * (K[i] - N[i] - sum(A[i,:].*N)) / K[i] for i in 1:S]
    return
end

function lv_jac(J,N,p,n)
    r = p[1]; K = p[2]; A = p[3]
    S = length(r)
    J = zeros(S,S)
    for i in 1:S, j in 1:S
        J[i,j] = (i == j) ? r[i] * (K[i] - 2N[i] - sum(A[i,:].*N)) / K[i] : -r[i] * A[i,j] * N[j] / K[i]
    end
    return
end

r = [1 -1]
K = [1 1]
A = [0 -1
    1 0]
p = [r,K,A]
S = length(r)
N0 = rand(S,1)
ContinuousDynamicalSystem(lv,N0,p,lv_jac)

[Datseris]

Hi there, sorry about that but at the moment in DynamicalSystems.jl you cannot use matrices as states. Only vectors. Chancing this would require solving JuliaDynamics/DynamicalSystemsBase.jl#96 which at the moment would take too much time for me to do. If someone contributes it of course then that's fine, but I wouldn't expect it anytime soon.

In the meantime, I am fixing a bug now. You were supposed to get a message "You can't have matrix as state", but due to a wrong sequencing in code you didnt...

[jarroyoe]

Thanks for letting me know about the lack of support. I found a workaround that could be useful if someone happens to have a similar problem in the future. By specifying S as a parameter, we can make all the parameters into a vector and decompose them inside the rule and Jacobian functions.

E.g.

function lv(dN,N,p,n)
    S = p[1]
     r = p[2:(S+1)]
     K = p[(S+2):(2S+1)]
     A = reshape(p[(2S+2):end],S,S)

    dN = [r[i] * N[i] * (K[i] - N[i] - sum(A[i,:].*N)) / K[i] for i in 1:S]
    return
end

[Datseris]

true, although in the original equations you pasted at the start of the issue, there is no reason for the state (the N_i) to be matrix. They are sequentially arranged from i = 1 .... end, just like a vector. Best use a vector from the scratch I'd say.

Here you're a bit confused.... the parameters of the system can be matrices. This A is a parameter, not the state. Only the state N and its derivative dN need to be vectors.

[Datseris]

your issue arose because you wrote N0 = rand(S,1) which is a matrix. Just do N0 = rand(S) instead. this problem has nothing to do with parameters. the parameters of the system can be completely arbitrary things, including weird custom types.

[jarroyoe]

Thank you, that worked as expected.