Help solving coupled differential equations involving complex numbers

Original title: “Some help solving set of differential equations”

Hi. I am stuck in doing some homework, and I was hoping I could get some help.
The task is to analytically or numerically solve some differential equations. I will provide the full task for completeness:

These are the DE’s to solve, in the following task:


“Real valued” Rabi frequencies means that \Omega is simply some real number, which reduces the complexity quite a bit. Below is my attempt at defining this, using ModelingToolkit.jl:

using ModelingToolkit
using DifferentialEquations

@variables t c̃₁(t) c̃₂(t)    # independent and dependent variables
@parameters Ω Δ             # parameters 

∂ₜ = Differential(t) # define an operator for the differentiation w.r.t. time

Inspired by the docs, I tried to define a vector of equations:

@named eqs = ODESystem([∂ₜ(c̃₁) ~ Ω*im*c̃₂*cis(t*Δ),
                        ∂ₜ(c̃₂) ~ Ω*im*c̃₁*cis(-t*Δ)
])

Resulting in

ERROR: type Array has no field lhs
Stacktrace:

So instead, I tried what I found inntuitive - defining two equations, and giving those as a vector:

julia> @named eq1 = ODESystem(∂ₜ(c̃₁) ~ Ω*im*c̃₂*cis(t*Δ))
Model eq1 with 2 equations
States (2):
  c̃₁(t)
  c̃₂(t)
Parameters (2):
  Ω
  Δ

julia> @named eq2 = ODESystem(∂ₜ(c̃₂) ~ Ω*im*c̃₁*cis(-t*Δ))
Model eq2 with 2 equations
States (2):
  c̃₂(t)
  c̃₁(t)
Parameters (2):
  Ω
  Δ

So far so good I think. I do not get why both equations say “with 2 equations”.
I then attempt:

julia> using DifferentialEquations: solve

julia> initial_conds = [c̃₁ => 1.0, c̃₂ => 0.0];

julia> tspan = (0.0, 10.0);

julia> params = [Ω => 1.0, Δ => 1.0];

julia> prob = ODEProblem([eq1, eq2], initial_conds, tspan, params)
ERROR: No methods were found for the model function passed to the equation solver.
The function `f` needs to have dispatches, for example, for an ODEProblem
`f` must define either `f(u,p,t)` or `f(du,u,p,t)`. For more information
on how the model function `f` should be defined, consult the docstring for
the appropriate `AbstractSciMLFunction`.

Offending function: f
Stacktrace:
 [1] isinplace(f::Vector{ODESystem}, inplace_param_number::Int64, fname::String, iip_preferred::Bool; has_two_dispatches::Bool, isoptimization::Bool)
   @ SciMLBase ~/.julia/packages/SciMLBase/RAGXU/src/utils.jl:236
 [2] isinplace (repeats 2 times)
   @ ~/.julia/packages/SciMLBase/RAGXU/src/utils.jl:229 [inlined]
 [3] ODEProblem(f::Vector{ODESystem}, u0::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Float64}, p::Vector{Pair{Num, Float64}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ SciMLBase ~/.julia/packages/SciMLBase/RAGXU/src/problems/ode_problems.jl:168
 [4] ODEProblem(f::Vector{ODESystem}, u0::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Float64}, p::Vector{Pair{Num, Float64}})
   @ SciMLBase ~/.julia/packages/SciMLBase/RAGXU/src/problems/ode_problems.jl:167
 [5] top-level scope
   @ REPL[33]:1

Now I realize that I could be going about it the wrong way. Is it even an ODE system, or something else? I do however not know how to proceed, and would much appreciate some pointers.

Do anyone have the know-how and time to help me out?

ModelingToolkit variables are all real, as written you get

julia> eqs = [∂ₜ(c̃₁) ~ Ω*im*c̃₂*cis(t*Δ),
              ∂ₜ(c̃₂) ~ Ω*im*c̃₁*cis(-t*Δ)
       ]
2-element Vector{Vector{Equation}}:
 [Differential(t)(c̃₁(t)) ~ -Ω*c̃₂(t)*sin(t*Δ), 0 ~ Ω*c̃₂(t)*cos(t*Δ)]
 [Differential(t)(c̃₂(t)) ~ -Ω*c̃₁(t)*sin(-t*Δ), 0 ~ Ω*c̃₁(t)*cos(-t*Δ)]

Which probably isn’t what you mean. Try defining c̃₁ and c̃₂ with explicit real and imaginary parts.

@variables c̃r₁ c̃i₁ c̃₁ c̃r₂ c̃i₂ c̃₂
c̃₁ = c̃r₁ + im*c̃i₁
c̃₂ = c̃r₂ + im*c̃i₂
1 Like

Thank you for your suggestion.
With the following definitions:

@variables t c̃r₁ c̃i₁ c̃₁ c̃r₂ c̃i₂ c̃₂  # independent and dependent variables
@parameters Ω Δ                     # parameters 

c̃₁ = c̃r₁ + im*c̃i₁
c̃₂ = c̃r₂ + im*c̃i₂

∂ₜ = Differential(t) # define an operator for the differentiation w.r.t. time

I get the following error:

julia> @named eqs = ODESystem([∂ₜ(c̃₁) ~ Ω*im*c̃₂*cis(t*Δ),
                               ∂ₜ(c̃₂) ~ Ω*im*c̃₁*cis(-t*Δ)
       ])
ERROR: type ComplexTerm has no field metadata
Stacktrace:
  [1] getproperty(x::Symbolics.ComplexTerm{Real}, f::Symbol)
 ...

I get the same error when defining a single equation:

julia> @named eq1 = ODESystem(∂ₜ(c̃₁) ~ Ω*im*c̃₂*cis(t*Δ))
ERROR: type ComplexTerm has no field metadata

Should I define seperate equations for the real and imaginary parts, giving me 4 equations? If yes, I am unsure how to even go about that. Should I use the real and imag functions in defining the equations?

Open an issue. This is a very specific question to ModelingToolkit with complex numbers and it would be helpful if you title it as such.

Done. Are you saying this could be solved without ModelingToolkit? Any quick pointers on how?

Just use complex numbers in DiffEq directly.

2 Likes

Thanks guys. This ended up doing the trick:

using DifferentialEquations
using ParameterizedFunctions

my_ode = @ode_def begin
    dc̃₁ = Ω*im*c̃₂*exp( im*t*Δ)
    dc̃₂ = Ω*im*c̃₁*exp(-im*t*Δ)
end Ω Δ

u₀ = Complex{Float64}[1, 0]
tspan = (0.0, 10.0)
p = [1.0, 1.0]

prob = ODEProblem(my_ode, u₀, tspan, p)
sol = solve(prob)
1 Like