ModelingToolkit complex system of equations

Specific Domains Modelling & Simulations
1

A few years ago I had a similar problem regarding a system of complex equations and I had found a workaround: ModelingToolkit matrix form system of equations and complex values.
As far as I can tell complex symbolics are still not supported in ModelingToolkit, but now unfortunately my workaround is no longer functioning either. Some of the ModelingToolkit syntax has changed in the almost 3 years since, so here is an updated MWE:

@parameters icomplex Γ32 δl δμ Ωl Ωμ
@variables t (ρ(t))[1:3,1:3]
D = Differential(t)

H = [
    -δl-δμ -Ωμ/2 0
    -Ωμ/2 -δl -Ωl/2
    0  -Ωl/2 0
]

L = [
    0 0 -Γ32*ρ[1,3]/2
    0 Γ32*ρ[3,3] -Γ32*ρ[2,3]/2
    -Γ32*ρ[3,1]/2 -Γ32*ρ[3,2]/2 -Γ32*ρ[3,3]
]

eqns = D(ρ) ~ -icomplex*(H*ρ-ρ*H) + L

@mtkbuild bloch = ODESystem(Symbolics.scalarize(eqns), t)

ρᵢ = zeros(ComplexF64,3,3)
ρᵢ[2,2] = 1
u0 = [ρ[idx,idy] => ρᵢ[idx,idy] for idx in 1:3 for idy in 1:3]

p = [Ωl => 1, Ωμ=> 0, δl => 0., δμ => 0., Γ32 => 1., icomplex => 1.0im]

prob = ODEProblem(bloch, u0, (0., 20.), p, jac = true)
sol = solve(prob, Tsit5())

The issue occurs when constructing the ODEProblem, throwing the an InexactError: Float64(0.0 + 1.0im). It seems like internally it tries to convert icomplex to a Float64.

Is there a different workaround possible?

2

@variables t (ρ(t))[1:3,1:3]::Complex? I haven’t tried complex arrays but complex variables should work if the type is appropriately declared.

3

I tried this but unfortunately it fails with the following error, which happens in the Symbolics.scalarize(eqns) step

TypeError: in Complex, in T, expected T<:Real, got Type{Complex{Real}}
4

It still isn’t functioning, but I’ve narrowed down the issue, but also found a new one.

The dissipator term L has explicit ρ terms, when these are removed from the equation the TypeError dissapears. Constructing the problem then throws the following error: ERROR: Initialization expression (ρ(t))[1, 1] is currently not supported.

5

Hi, I stumbled across this post when trying to solve a complex differential equations that also involves a density matrix and a hamiltonian. The following MWE compiles under ModelingToolkit:

using ModelingToolkit, Symbolics, LinearAlgebra, DifferentialEquations

@independent_variables x
D = Differential(x)

ρ0 = ComplexF64[exp(-0.5) 1im; 0 exp(-0.7)]
U = [0.88 0.27; -0.35 + 0.1im 0.66]

vars = @variables begin
    ρ(x)[1:2, 1:2] = ρ0
    H(x)[1:2, 1:2]
end

eqs = [
    H ~ x^2 * U
    D(ρ) ~ -im * (H * ρ - ρ * H)
]

sys = System(eqs, x, vars, []; name = :foo)
foo = mtkcompile(sys)

Here, ρ0 is the initial value of the matrix variable ρ and U is a complex matrix used inside the variable H. The problem arises when tyring to build a ODEProblem and give it to DifferentialEquations. The first line of

prob = ODEProblem(foo, [], (0.0, 1.0))
sol = solve(prob)

produces InexactError: Real(0.0 + 1.0im). Note that I wrote a complex number inside ρ0. If I change that line so we have a complex matrix with real values, i.e.,

ρ0 = ComplexF64[exp(-0.5) 0; 0 exp(-0.7)]

the ODE problem is built with output:

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: 4-element Vector{Float64}:
 0.4965853037914095
 0.0
 0.0
 0.6065306597126334

However it says that u0 is a 4-element real vector, not a complex one. I try passing the keyword argument u0_eltype = ComplexFloat to the ODEProblem constructor but the output is exactly the same. If I execute the solve(prob) line, then the following error arises

ERROR: InexactError: Float64(0.0 + 2.9685246098730456e-14im)

I think it has to due with the fact that treats u0 as a real vector, as per the previous output.

I thought, okay, let’s declare that explicitly inside the @variables call. So I modify the previous code to

vars = @variables begin
    ρ(x)[1:2, 1:2]::Complex = ρ0
    H(x)[1:2, 1:2]::Complex
end

eqs = [
    H ~ x^2 * U
    D(ρ) ~ -im * (H * ρ - ρ * H)
]

But now, the second equation inside eqs, the one with D(ρ) arises the following error:

ERROR: TypeError: in typeassert, expected DataType, got Type{Complex}

It looks like the most robust option for now is to separate everything into a real and a complex part, even though is more of a hassle.