ModelingToolkit matrix form system of equations and complex values

I’ve started looking into ModelingToolkit and I wanted to solve optical bloch differential equations with it, which have a complex part.
I’ll preface this by saying I don’t have much experience with Julia yet, I started working with it recently.

A simple example of such a system would be a three level system with a hamiltonian
H = \begin{pmatrix} -\delta_l - \delta_\mu & -\Omega_\mu/2 & 0 \\ \bar\Omega_\mu/2 & -\delta_l & -\Omega_l/2 \\ 0 & -\bar\Omega_l/2 & 0 \end{pmatrix}
and a dissipator
L = \begin{pmatrix} 0 & 0 & -\Gamma_{32} \rho_{1,3}(t)/2 \\ 0 & \Gamma_{32}\rho_{3,3}(t) & -\Gamma_{32}\rho_{2,3}(t)/2 \\ -\Gamma_{32}\rho_{3,1}(t)/2 & -\Gamma_{32}\rho_{3,2} & -\Gamma_{32}\rho_{3,3}(t) \end{pmatrix}
The governing equation is:
\frac{\partial \rho(t)}{\partial t} = -i[H,\rho] + L

I can define all the variables and parameters in Julia and make the matrices

using ModelingToolkit
@variables t ρ[1:3,1:3](t)
@parameters Ωl Ωμ δl δμ Γ32 
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]]

I run into trouble when setting up the ODESystem in matrix notation.

eq = -1im*(H*ρ - ρ*H) + L
bloch = ODESystem(D(ρ) ~ eq)

This throws a conversion error (Cannot convert an object of type Missing to Sym).
I then tried just writing it as a list of equations:

eqns = [D(ρ[idx,idy]) ~ eq[idx,idy] for idx in 1:3 for idy in 1:3]
bloch = ODESystem(eqns)

This however doesn’t work because of the complex part of the system, eqns is an array of arrays, where the 1st entry is the real part and the 2nd entry the complex part, which then obviously doesn’t work with ODESystem.

Finally as a test, when I use as an equation \frac{\partial \rho(t)}{\partial t} = -[H,\rho] + L (so no imaginary part) and use the list of equations approach I can create an ODESystem and solve the ODE with DifferentialEquations.jl, but it obviously doesn’t return any physical results then.

As a workaround I can use build_function on H and L and manually create the rhs function for the solver, but I’d like to use ModelingToolkit.jl because more complex examples will have \delta_l or other parameters be time dependent (so switched to variables), and this is rather painful to do when using the build_function approach because it requires manually creating the rhs function for each change.

To summarize, I have two main questions:

  • Can you give ODESystem a matrix form of a system of equations
  • Is it possible to work with complex values for an ODESystem

For completeness I’ll post the workaround I found (without matrix form or complex symbols). It basically uses another parameter to specify the complex part, which works with this simple example when the rates aren’t complex.

@parameter icomplex

eq = -icomplex*Symbolics.scalarize(H*ρ-ρ*H) + L
eqns = [D(ρ[idx,idy]) ~ eq[idx,idy] for idx in 1:3 for idy in 1:3]
bloch = ODESystem(eqns)

ρᵢ = 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 => 1im]

This can then be used with DifferentialEquations.jl

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

Somewhat, almost. Symbolic arrays were just added, and so this can be done but it will scalarize because we haven’t made all of the passes compatible with symbolic arrays.

Symbolic complex is coming after symbolic arrays are solid.

1 Like

Can you point me towards some more info regarding the symbolic arrays? I can’t seem to find it; I was under the (wrong) impression that if a variable was initiated with @variables x[1:3,1:3] it would make a symbolic array.
Regarding the symbolic complex, I look forward to be able to use symbolic complex (in reality \rho and H are Hermitian; so full use of symbolic complex would be very helpful). Any idea when roughly this would be implemented?

As of ModelingToolkit v5.21, that is the case: it makes a symbolic array and generates symbolic equations. And then structural_simplifiy doesn’t know how to handle it so it scalarizes the equations before compiling to code :sweat_smile: . But that was done so that way we could get symbolic arrays released ASAP, and so actually keeping symbolic arrays all of the way through is one of the things we’re working on.

I was still using V5.20 it turns out, updated to V5.22. I can’t get the (non-complex) matrix version to work though, even just a simple example like ODESystem(D(ρ) ~ ρ) fails with error code no method matching hasmetadata(::Missing, ::Type{Symbolics.VariableDefaultValue}). This suggests I have to set default values for the variables but that would kind of defeat the purpose of using symbolics.

Because the type ρ changed my original flattening (non-complex) version doesn’t work anymore either, because the expression H*ρ - ρ*H, instead of returning an array with the results now is a statement with a bunch of broadcasts and maps in it:
(broadcast(-, Num[-δl - δμ (-1//2)*Ωμ 0; (-1//2)*Ωμ -δl (-1//2)*Ωl; 0 (-1//2)*Ωl 0]*map(#1, ρ), map(#1, ρ)*Num[-δl - δμ (-1//2)*Ωμ 0; (-1//2)*Ωμ -δl (-1//2)*Ωl; 0 (-1//2)*Ωl 0]))[Base.OneTo(3),1:3]
and evaluating it for [i,j] doesn’t actually return the i,j entry of the result.

Odd, open an issue for that. I don’t know that hasmetadata error, but it could be from @YingboMa 's change to metadata in v5.22.