DifferentialEquations.jl
DifferentialEquations.jl copied to clipboard
Inconsistent choice of default algorithm order when eltype is complex
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.
Suggested fix: use uRealtype = real(eltype(u0)) instead of uEltype = eltype(u0) in the heuristics.
This makes sense. Yes, could you PR that?
That big method branch should probably just change to a GMRES default now, since that would be matrix-free.