[spacedome]

My experience with various linear solvers (which a diffeq solver typicaly relies on) is that if they assume commutivity anywhere at all, which they often do, they almost definitely will not work for quaternions. Even if a derivation of the algorithm can be done with non-commutivity in mind, the implementation typically is not. Many NLA solvers are based on orthogonality transforms, which do not translate directly to quaternions, and even solvers which use only inner products and mat-vec multipication like BiCG do not really work as-is for quaternion matrices. Luckily dense linear solves with LU are fine! One can usually use the complex matrix expansion for quaternions and solve that instead without much difficulty though, and for computing eigenvalues this is perhaps even preferable, as it gives a canonical representation. I have not tried to use DifferentialEquations.jl for quaternion problems, so maybe they have figured some of it out, but it is non-trivial.

[ChrisRackauckas]

ODE solvers only use LU

[spacedome]

This is good to know, thanks Chris. I mostly solve sparse PDEs, being able to always use LU makes everything much simpler, especially with quaternions.