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?

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

score 14 · Accepted Answer · answered Jul 15 '15 at 19:12

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.

Forcing an ODE solver to preserve the norm

1 Answers1