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
ODESystema matrix form of a system of equations - Is it possible to work with complex values for an
ODESystem
EDIT:
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())
. 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.