I am interested in solutions to $$ \frac{d}{dt} \Psi = -iH \Psi $$ for $H$ hermitian and time independent. This boils down to evaluating $$ \Psi(t) = e^{-iHt}\Psi_0 $$ at points of interest $t_n$. I want to quickly compute $\Psi(t)$ up to a large final time $T$ by spacing the $t_n$ randomly and advancing the propagator over the nonuniform grid.

To be specific, fix $\Delta t$ and draw $t_n$ uniformly at random from $\Delta t, 2\Delta t, ..., N\Delta t$. The gaps $t_{n+1}-t_{n}$ will be distributed geometrically (ie we will have large gaps).

Any ideas how I could compute $e^{-iHt}$ accurately and efficiently when $t$ may be large?

Fast approximations that I know of (eg Strang splitting, Trotter product) are very efficient but I have found them to be horribly inaccurate and unstable for large gaps.

Machine precision accuracy is achieved by expanding $e^{-iHt}$ in Chebyshev polynomials, but the computational cost (measured by counting applications of $H$) scales linearly with the timestep size and we get to $T$ no faster.

I am beginning to think that spectrally accurate approximations of $e^{-iHt}$ provably require linearly more applications of $H$ as $T$ increases. Does anyone know if this is proven?

Thanks for solving my thesis problem,


  • $\begingroup$ Is there a compelling reason for spacing the $t_n$ randomly? $\endgroup$
    – j.c.
    Nov 7, 2010 at 12:57
  • $\begingroup$ Yes, it would let you reconstruct the spectrum of $H$ from few $t_n$. However, unless there is something very creative out there (cmon theoretical physics!) my approximation problem would be the same one that a person interested in $\Delta t$ large would face. $\endgroup$
    – dranxo
    Nov 7, 2010 at 20:10
  • $\begingroup$ Have you looked at "Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later" www.cs.cornell.edu/cv/ResearchPDF/19ways+.pdf? $\endgroup$
    – Ian Grooms
    Nov 8, 2010 at 17:57
  • $\begingroup$ Yes, but I have a little more structure then that paper deals with. One difference I have is the $iH$ as opposed from simply $H$. Also, I think I can just diagonalize $H$ and reduce the problem to one of approximating $e^{-ix}$. Chebyshev polynomials are orthogonal wrt a weighted norm and I would not be surprised if what I currently have is optimal. $\endgroup$
    – dranxo
    Nov 8, 2010 at 22:47

1 Answer 1


First, you will not be able to compute large step sizes directly. Taking such high degree polynomials in Chebyshev will not be competitive to splittings.

If the length of your time steps goes beyond the stability limit, the only remedy is to compute intermediate steps to improve accuracy.

E.g. $exp(-ihH) = exp(-ih/2 H) exp(-ih/2 H)$, and then approximate the r.h.s. exponentials with some splitting. Think of it as intermediate stages if you don't like the idea of using a smaller time-step.

Much better accuracy can be reached for higher order splitting methods. Without additional information on H, I recommend you to use a fourth or higher order splitting method SRKN4_6^b. You can find the coefficients in Table 3 of S. Blanes and P.C. Moan, Practical Symplectic Partitioned Runge-Kutta and Runge-Kutta-Nyström Methods, J. Comput. Appl. Math., 142 (2002), pp. 313-330


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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