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?

  • 7
    $\begingroup$ 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. $\endgroup$ Jul 15 '15 at 16:32
  • $\begingroup$ @ChristianClason, you should make this comment into an answer. $\endgroup$
    – 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.

If you want to know more about geometrical integration:


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.