DifferentialEquations.jl icon indicating copy to clipboard operation
DifferentialEquations.jl copied to clipboard

Inconsistent choice of default algorithm order when eltype is complex

Open danielwe opened this issue 2 years ago • 2 comments

The default choice of algorithm for ODEs tries to infer the desired level of precision by looking at number types: if it's not the "ordinary" Float32/64, it gives you a high-order algorithm because you're probably using BigFloat or some other exotic high-precision type. However, ComplexF32/64 will sometimes also land you in the high-order branch because it's not Float32/64. In particular, that happens here:

https://github.com/SciML/DifferentialEquations.jl/blob/8b4d07c5523eaf905bf245fbc439a2c820c1cdc3/src/ode_default_alg.jl#L68

However, it does not happen in the corresponding branch for huge (N > 10000), explicit-only problems:

https://github.com/SciML/DifferentialEquations.jl/blob/8b4d07c5523eaf905bf245fbc439a2c820c1cdc3/src/ode_default_alg.jl#L35

However, this has its own problems: any complex type escapes the high-order branch, including Complex{BigFloat}---not just ComplexF32/F64, as was probably the intention.

Suggested fix: use uRealtype = real(eltype(u0)) instead of uEltype = eltype(u0) in the heuristics.

danielwe avatar Sep 21 '22 19:09 danielwe

Suggested fix: use uRealtype = real(eltype(u0)) instead of uEltype = eltype(u0) in the heuristics.

This makes sense. Yes, could you PR that?

ChrisRackauckas avatar Sep 24 '22 21:09 ChrisRackauckas

That big method branch should probably just change to a GMRES default now, since that would be matrix-free.

ChrisRackauckas avatar Sep 24 '22 21:09 ChrisRackauckas