Numerical solution of linear Schroedinger ODE with almost-normal Hamiltonian matrix

2011-07-21T11:39:20Z

I am calculating numerical solutions of time-dependent Schroedinger equation

$\frac{d\Psi}{dt} = - i H \Psi$

where $\Psi$ is an $N$-element complex vector and $H$ is an $N \times N$ complex matrix, which is ``almost normal''. That is, $H = H_0 + i D$, where $H_0$ is normal (and often Hermitian), $D$ is Hermitian but does not commute with $H_0$, and $||D|| \ll ||H_0||$. The dimension $N$ is betwen 30 and 200.

So far, I have been using a 4-th order Runge-Kutta specialised for linear ODEs (taken from Zingg and Chisholm). However, the accuracy is not as good as I'd like, probably due to strong oscillations arising from the $H_0$ part. Is there a better method out there? I'd prefer one which would not require the calculation of $e^{-iH\tau}$ ($\tau$ - time-step of the ODE solver).

EDIT: more data about the problem:

The eigenvalues of $H_0$ are literally random, as $H_0 = h + diag(dE_i)$, where $dE_i ~ N(0,\sigma)$ or $dE_i \sim L(\alpha,\sigma)$. Matrix $h$ has eigenvalues of the order of unity, and $L(\alpha,\sigma)$ is the symmetric Levy alpha-stable distribution (each diagonal value of $H_0$ is perturbed by an independent distr.) with scale parameter $\alpha \in (0,2)$ and strength $\sigma$. $D$ has eigenvalues of the order of 0.1. I choose the time step to be of the order of $0.01/E_\text{max}$, where $E_\text{max}$ is the maximum absolute value of real or imaginary part of any eigenvalue of $H$.

So in short, $H$ can be nasty :(

Answer by fedja for Numerical solution of linear Schroedinger ODE with almost-normal Hamiltonian matrix

2011-07-22T15:19:25Z

This should really be a comment, but it is too long.

Since Euler is out of question on such a long interval, my next suggestion would be to try to run the 3rd order Runge-Kutta about which I am pretty certain that it is not screwed up anywhere. The recursion step is $$ \begin{aligned} x(t+\tau)&= \cr &x(t) \cr +&\frac \tau 6[H(t)+4H(t+\frac\tau 2)+H(t+\tau)]x(t) \cr +&\frac{\tau^2}6[H(t)^2+2H(t+\tau)H(t+\frac\tau 2)]x(t) \cr +&\frac{\tau^3}{6}H(t)^3x(t) \end{aligned} $$ (of course, you should never multiply matrices, only a matrix by a vector and you should pay attention to the order of matrix multiplication in the third term; I leave the obvious optimizations to you).

Try to compare the results with what you get from your 4th order Runge-Kutta. The maximal step for this 3rd order method would be $10^{-4}$ if you want to get anything meaningful in your setting (if you can afford $10^{-5}$, it would be much better). Let us know what you see. If this method exhibits slower norm growth than the 4th order one, something is screwed up in the programming. If it blows up even faster, it may, indeed, be a precision problem.

Numerical solution of linear Schroedinger ODE with almost-normal Hamiltonian matrix - MathOverflow

Numerical solution of linear Schroedinger ODE with almost-normal Hamiltonian matrix

Answer by fedja for Numerical solution of linear Schroedinger ODE with almost-normal Hamiltonian matrix