Trying to CVODE_BDF with complex-valued ODE (?)

General Usage
#1

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 + 0im
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 ?

#2

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

#3

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 ?

#4

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
#5

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

#6

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.

#7

You’re using yesterday’s release?

#8

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
#9

BEAUTIFUL!

1 Like