Complex ODEs DifferentialEquations pkg vs ModelingToolkit pkg

[SageanneS]

I am able to simulate a complex system of differential equations with the DifferentialEquations.jl package:

using DifferentialEquations

function NGNeuralMass(dzdt,z,p,t)
      dzdt[1] = (1/C)*(-im*((z[1]-1)^2)/2 + (((z[1]+1)^2)/2)*(-Δ + im*η_0 + im*v_syn*z[2]) - ((z[1]^2-1)/2)*z[2])
      dzdt[2] = alpha_inv*((k/(C*pi))*(1-abs(z[1])^2)/(1+z[1]+conj(z[1])+abs(z[1])^2) - z[2])
end

Δ = 0.5 
η_0 = 21.5
v_syn = -10
alpha_inv = 35
k = 0.105
C = 30

z0 = Array{Complex{Float64}}(undef,2);
z0[1] = 0.5
z0[2] = 1.6
dzdt = Array{Complex{Float64}}(undef,2);

tspan = (0.0, 100.0)
prob = ODEProblem(NGNeuralMass,z0,tspan)
sol = solve(prob)

When trying to set up the ODE system using ModelingToolkit, I get the following error at the ODESystem set-up step:

using ModelingToolkit, DifferentialEquations

@parameters t C Δ η_0 v_syn
@variables Z(t)::Complex g(t)::Complex
D = Differential(t)

eqs = [
      D(Z) ~ (1/C)*(-im*((Z-1)^2)/2 + (((Z+1)^2)/2)*(-Δ + im*(η_0) + im*v_syn*g) - ((Z^2-1)/2)*Z),
      D(g) ~ alpha_inv*((k/(C*pi))*(1-abs(Z)^2)/(1+Z+conj(Z)+abs(Z)^2) - g)
]
NGNeuralMass = ODESystem(eqs)

julia> NGNeuralMass = ODESystem(eqs)
ERROR: ArgumentError: The differential variable missing is not unique in the system of equations.
Stacktrace:
 [1] ODESystem(eqs::Vector{Equation}, iv::Nothing; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/P76eb/src/systems/diffeqs/odesystem.jl:176
 [2] ODESystem(eqs::Vector{Equation}, iv::Nothing) (repeats 2 times)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/P76eb/src/systems/diffeqs/odesystem.jl:151
 [3] top-level scope
   @ REPL[28]:1

I am not understanding where this uniqueness error comes from; as it works with the DifferentialEquations package, and I have declared my variables as complex, I am not sure why it breaks down at this step. It could be a bug in MTK.

[YingboMa]

Unfortunately, complex number support in MTK is still lacking, as a work-around, you can use

using ModelingToolkit, OrdinaryDiffEq
@parameters t C Δ η_0 v_syn k
@variables Z(t) g(t)
# don't promote to Complex{Num}
Z = ModelingToolkit.unwrap(Z)
g = ModelingToolkit.unwrap(g)
t, C, Δ, η_0, v_syn, k = map(ModelingToolkit.unwrap, [t, C, Δ, η_0, v_syn, k])
D = Differential(t)

eqs = [
       Equation(D(Z), (1/C)*(-im*((Z-1)^2)/2 + (((Z+1)^2)/2)*(-Δ + im*(η_0) + im*v_syn*g) - ((Z^2-1)/2)*Z))
       D(g) ~ alpha_inv*((k/(C*pi))*(1-abs(Z)^2)/(1+Z+conj(Z)+abs(Z)^2) - g)
]
@named NGNeuralMass = ODESystem(eqs, t)

z0 = [Z=>0.5, g=>1.6]
p = [
     Δ => 0.5
     η_0 => 21.5
     v_syn => -10
     alpha_inv => 35
     k => 0.105
     C => 30
    ]
tspan = (0.0, 100.0)
prob = ODEProblem(NGNeuralMass,z0,tspan,p)
prob = remake(prob, u0=complex.(prob.u0))
solve(prob, Tsit5())

[SageanneS]

Ok, got it. Thanks very much for the work-around! One other question based on the output from the following:

using ModelingToolkit, OrdinaryDiffEq
@parameters t C Δ η_0 v_syn alpha_inv k
@variables Z(t) g(t)
# don't promote to Complex{Num}
Z = ModelingToolkit.unwrap(Z)
g = ModelingToolkit.unwrap(g)
t, C, Δ, η_0, v_syn, alpha_inv, k = map(ModelingToolkit.unwrap, [t, C, Δ, η_0, v_syn, alpha_inv, k])
D = Differential(t)

eqs = [
       Equation(D(Z), (1/C)*(-im*((Z-1)^2)/2 + (((Z+1)^2)/2)*(-Δ + im*(η_0) + im*v_syn*g) - ((Z^2-1)/2)*Z))
       D(g) ~ alpha_inv*((k/(C*pi))*(1-abs(Z)^2)/(1+Z+conj(Z)+abs(Z)^2) - g)
]
@named NGNeuralMass = ODESystem(eqs, t)

z0 = [Z=>0.5, g=>1.6]
p = [
     C => 30
     Δ => 0.5
     η_0 => 21.5
     v_syn => -10
     alpha_inv => 35
     k => 0.105
    ]
tspan = (0.0, 100.0)
prob = ODEProblem(NGNeuralMass,z0,tspan,p)
prob = remake(prob, u0=complex.(prob.u0))
sol = solve(prob, Tsit5())

I want to calculate the phase oscillation of the system (ψ), that I get to via the following steps:

W = (1 .- conj.(sol[1,:]))./(1 .+ conj.(sol[1,:]))
R = (1/(C*pi))*(W+conj.(W))/2
ψ = log.(sol[1,:]./R)/im

Where I get the following output:

ψ 
1147-element Vector{SymbolicUtils.Mul{Complex{Real}, ComplexF64, Dict{Any, Number}, Nothing}}:
 (0.0 - 1.0im)*log(((1.5707963267948966 + 0.0im)C) / (0.3333333333333333 + 0.0im))
 (0.0 - 1.0im)*log(((1.5695392839471385 + 0.03306235490433268im)*C) / (0.3336234160179976 + 0.0im))
 (0.0 - 1.0im)*log(((1.568491838417757 + 0.05886238061355869im)*C) / (0.3337773832695455 + 0.0im))
 (0.0 - 1.0im)*log(((1.5660076119514992 + 0.10515150944162954im)*C) / (0.3340240679971305 + 0.0im))
 (0.0 - 1.0im)*log(((1.562929810588537 + 0.14804671436135264im)*C) / (0.33424152546942526 + 0.0im))
 (0.0 - 1.0im)*log(((1.55821227300668 + 0.19946369258626204im)*C) / (0.3344996912039376 + 0.0im))
 (0.0 - 1.0im)*log(((1.5523983469775373 + 0.2504944309087575im)*C) / (0.3347581566670816 + 0.0im))
 ⋮
 (0.0 - 1.0im)*log(((-3.525152867386194 - 3.270103687356945im)*C) / (-1.2223085837325254 + 0.0im))
 (0.0 - 1.0im)*log(((-3.5473424337621227 - 3.406801640925555im)*C) / (-1.2165845002197668 + 0.0im))
 (0.0 - 1.0im)*log(((-3.5732569913673062 - 3.5544541437749193im)*C) / (-1.211554276560347 + 0.0im))
 (0.0 - 1.0im)*log(((-3.603672388065652 - 3.7145649561940144im)*C) / (-1.2072105359006144 + 0.0im))
 (0.0 - 1.0im)*log(((-3.639573952360074 - 3.8889299176128787im)*C) / (-1.2035491394311284 + 0.0im))
 (0.0 - 1.0im)*log(((-3.6822273307309787 - 4.079707488479976im)*C) / (-1.200569667009411 + 0.0im))
 (0.0 - 1.0im)*log(((-3.6967992305069615 - 4.141461040015192im)*C) / (-1.1997981110780622 + 0.0im))

ψ is returned as 'Nothing' and therefore I cannot access it for plotting. I am unclear as to how to go from a vector of symbolics to a vector of numerics, I have seen a lot of information online in the reverse direction only, but I may be missing something very obvious here.

Appreciate your help and discussion!

[YingboMa]

Do you mean

W = (1 .- conj.(sol[Z]))./(1 .+ conj.(sol[Z]))
pdict = Dict(p)
R = (1/(pdict[C]*pi))*(W+conj.(W))/2
ψ = log.(sol[Z]./R)/im

Note that C is a symbolic parameter not a number.

[SageanneS]

Yes - my mistake! Thanks. All cleared up.