Most of the systems that I solve have complex-valued functions and, until now, I always separate the real and imaginary part of each equation to solve them in the real domain, i.e., going from a single complex equation to two-coupled real ones. But then I tried to implement the complex problem directly into the system now that I am using Julia (I am new at Julia, but im enjoying it).
In this topic DifferentialEquations.jl error when solving equations with complex variables - Usage - JuliaLang I saw how to implement complex-valued ODE, but what if the system was stiff? Most of time I use solvers like CVODE_BDF, but for complex-valued ODE this warning happens: “MethodError: Cannot convert an object of type Vector{ComplexF64} to an object of type Ptr{Sundials._generic_N_Vector}”
I don’t know if I can do something in my problem description to avoid this error or the solver simple can’t handle this kind of problem. I will try to encompass the idea of my code in the example below (the following system isn’t necessary stiff, its just an example):
function example!(du,u,p,t)
x,y = u
du[1] = dx = y
du[2] = dy = im*x
end
prob = ODEProblem(example!, u₀, Δt) @time sol = solve(prob, CVODE_BDF())
It is the solver that can’t handle this kind of problems or is something that I am doing wrong? How can I solve these kind of problem without always separating in real/imaginary parts ?
Using OrdinaryDiffEq methods is a good alternative.
using DifferentialEquations
function example!(du,u,p,t)
x,y = u
du[1] = dx = y
du[2] = dy = im*x
end
Δt = (0, 10.0)
x₀ = 1.0 + 0im
y₀ = 1.0 + 0im
u₀ = [x₀,y₀]
prob = ODEProblem(example!, u₀, Δt)
@time sol = solve(prob, QNDF(autodiff=false))
@time sol = solve(prob, TRBDF2(autodiff=false))
@time sol = solve(prob, KenCarp47(autodiff=false))
With QNDF now outperforming Sundials CVODE it makes a lot of sense to just use that, since it’s completely compatible with complex numbers (as long as you avoid ForwardDiff.jl).
CVODE doesn’t autodiff. CVODE is a C++ code which only works on Float64 because it’s ahead of time compiled. it’s probably not fixable. (Well, we could split to real and imaginary parts for you, but that would be a mess and I don’t really see the point with OrdinaryDiffEq.jl around. If someone else wants to put the work into it, more power to them.)
QNDF either trades blows with CVODE_BDF or, at lower tolerances, it is a half magnitude below. You see this across all of the benchmarks. Note that BCR is re-running right now as well.
Make sure you’re using last night’s release. It’s a result from the last day.