Trying to CVODE_BDF with complex-valued ODE (?)

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

Δt = (0, 10)
x₀ = 1 + 0im
y₀ = 1 + 0
im
u₀ = [x₀,y₀]

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).

Oh sure, I use OrdinaryDiffEq methods. I didn’t know this QNDF, I will definitely take a look.

So the problem with what I was doing is the autodiff=true ? We cannot autodiff at CVODE while using complex numbers ?

CVODE doesn’t autodiff. CVODE is a C++ code which only works on Float64 because it’s ahead of time compiled. :man_shrugging: 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.)

2 Likes

Do you a place where you document this? This is not what I see in a particular application, perhaps I use wrong options

I was trying several methods and, so far, CVODE_BDF was doing better indeed.

At least for my system, QNDF was not able to solve it.

You’re using yesterday’s release?

Using the OrdinaryDiffEq v5.56 release of last night, the SciMLBenchmarks just finished updating. So see things like:

https://benchmarks.sciml.ai/html/StiffODE/Pollution.html

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.

https://sciml.ai/news/2021/05/24/QNDF/

3 Likes

BEAUTIFUL!

1 Like