I have some ODE systems generated with ModelingToolkit that have complex solutions but I only want to plot the real part of the solution. I know I can pull out the solution for a specific component manually and use real to get the real part of the solution, but I’d like to use the plot(sol, vars = []) interface since that is much more convenient.
MWE
using ModelingToolkit
using DifferentialEquations
using Plots
# define variables and parameters of system
@variables t ρ[1:2, 1:2](t) δ(t) Ω(t)
@parameters Δ Ω0 ω ϕ δ0 icomplex
D = Differential(t)
# Hamiltonian
H = [
0 Ω
Ω Δ+δ
]
# bloch equations
eq = -icomplex * Symbolics.scalarize(H*ρ - ρ*H)
eqns = [D(ρ[idx, idy]) ~ eq[idx,idy] for idx in 1:2 for idy in 1:2]
# time varying detuning
append!(eqns, [δ ~ δ0*cos(ω*t+ϕ)])
append!(eqns, [Ω ~ Ω0*sin(ω*t+ϕ)])
# create symbolic system
@named bloch = ODESystem(eqns)
# structural_simplify simplifies a system of equations if possible by removing
# reduntancies etc.
bloch_simplified = structural_simplify(bloch)
# initial conditions
ρᵢ = zeros(ComplexF64,2,2)
ρᵢ[1,1] = 1
u0 = [ρ[idx,idy] => ρᵢ[idx,idy] for idx in 1:2 for idy in 1:2]
# initial parameters
p = [Ω0 => 1., Δ => 20., δ0 => 1., ω => 20., ϕ =>0, icomplex => 1im]
# setup and solve of ODE, t from 0 => π
prob = ODEProblem(bloch_simplified, u0, (0., π), p, jac = true)
sol = solve(prob, Tsit5())
plot(sol, vars = [ρ[1,1], ρ[2,2]])
In this case I could use plot(sol, vars = [abs(ρ[1,1]), abs(ρ[2,2])]) because these components of the solution never go below zero, but other components like δ and Ω do. I tried using plot(sol, vars = [real(ρ[1,1]), real(ρ[2,2])]) but that doesn’t seem to work.
