I have an ODE of the form
$$ \frac{dy}{dt} = -i H y \enspace .$$
where $y$ is a complex vector and $H$ is a time dependent Hermitian matrix.
The norm of the solution $y(t)$ at any point in time should be 1, but due to accumulation of small numerical errors it ends up being substantially off.
A solution for this problem that is used by the qutip
library is to reset the ode solver every so often by normalizing $y$ at some time $t'$ and resuming from that point. Is there a rigorous explanation why this is a good idea (it seems to work in practice)?
Is there a better way? Even better if it does not require a modification of the solver code?