Forcing an ODE solver to preserve the norm

Question

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?

Answer 1

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.

If you want to know more about geometrical integration:

• Probably the most accessible discussion of geometric integration is Chapter 5 in A. Iserles, A First Course in the Numerical Analysis of Differential Equations, 2nd edition, Cambridge University Press, 2008.

• The "bible" is of course E. Hairer, G. Wanner, C. Lubich, Geometric Numerical Integration, 2nd edition, Springer, 2006 (see Chapter VII.6 for applications in quantum mechanics). There's also a survey article in Acta Numerica 12, 2003, p. 399-450 (here's a preprint), although that seems to focus on second-order ODEs.

• Specifically the application to the Schrödinger equation is discussed in E. Faou, Geometric Numerical Integration and Schrödinger Equations, Zurich Lectures in Advanced Mathematics, EMS, 2012.