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.
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
Ω do. I tried using
plot(sol, vars = [real(ρ[1,1]), real(ρ[2,2])]) but that doesn’t seem to work.