Open SageanneS opened 2 years ago
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())
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!
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.
Yes - my mistake! Thanks. All cleared up.
I am able to simulate a complex system of differential equations with the DifferentialEquations.jl package:
When trying to set up the ODE system using ModelingToolkit, I get the following error at the ODESystem set-up step:
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.