# Forcing an ODE solver to preserve the norm

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?

• You want to use an ode solver that is guaranteed to preserve the norm; these are called "symplectic" or "geometric" integrators. There's a book by Hairer, Lubich and Wanner that discusses these methods. – Christian Clason Jul 15 '15 at 16:32
• @ChristianClason, you should make this comment into an answer. – Bill Barth Jul 15 '15 at 18:35

The best approach is to use an ODE solver that is guaranteed to conserve the norm of the initial condition, i.e., for which $\|y_n\| = \|y_0\|$ for all $n\in\mathbb{N}$. Such solvers exist, and are called geometric integrators, since they preserve geometric properties of the exact solution (in this case, that energy is conserved, i.e., $\frac{d}{dt}\|y(t)\| = 0$). (A special class are symplectic integrators for ODE systems arising in Hamiltonian mechanics.)
The simplest such integrator (that is not specific to Hamiltonian systems) is actually a quite standard method, variously referred to as trapezoidal, implicit mid-point, or Crank-Nicolson method. In your case, it amounts to solving in every step the linear system $$\left(I+\frac{i\,h_n}2 H(t_{n+1})\right) y_{n+1} = \left(I-\frac{i\,h_n}2 H(t_{n})\right) y_n,$$ where $I$ is the identity matrix and $h_n = t_{n+1}-t_n$ is the current time step length. This method is implicit, but second-order accurate and A-stable. (I don't think there's a first order geometric integrator.) I don't know if it's included in ODEPACK (the basis of scipy's odeint routines), but it shouldn't be too hard to implement yourself.