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

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 :(

-
Are both H_0 and D explicitly time dependent? If H_0 is time independent, then you could diagonalize it once and for all at the start of your calculation and use the interaction picture. –  Jonathan Jul 21 '11 at 14:15
Additionally, do you know anything about the eigenvalues of D? The sign will determine whether the solution grows or decays, which might be relevant to cooking up a better numerical method. –  Jonathan Jul 21 '11 at 14:27
@Jonathan $H_0$ is time-independent. But then my perturbation becomes time-dependent, won't it simply move the oscillations to another place in the calculation? $D$ has only negative eigenvalues, so it should lead to a decay. Funnily, for some forms of $H_0$ the Runge-Kutta algorithm I use now blows up anyway. –  Katastrofa Jul 21 '11 at 16:56
@Jonathan Do you think that using the split-operator approach would help? I mean the approach in pci.uni-heidelberg.de/tc/usr/andreasm/academic/handouthtml/… –  Katastrofa Jul 21 '11 at 16:58
What are the smallest and the largest eigenvalues of $H$ and what are the norm and the smoothness of $D$ and the time interval? –  fedja Jul 21 '11 at 18:47

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.
Sorry, my $H$ is yours $-iH$, so put all $i$'s accordingly. –  fedja Jul 22 '11 at 15:20
Thanks. My $H$ is time-independent, so I can ignore the order of multiplication. –  Katastrofa Jul 22 '11 at 16:31
In that case what prevents you from forgetting all that Runge-Kutta stuff altogether and just using the standard Taylor expansion of the exponent with as many terms as you wish (6-7 should be more than enough)? I mean use $x(t+s)=x(t)-isHx(t)-\frac 12s^2H^2x(t)+\frac{i}6s^3H^3x(t)+\dots$. It is not overly computationally intensive (just a 7 matrix by vector multiplications at each step) and if you choose $s\|H\|<0.1$, you should be completely fine. Runge-Kutta is mainly there to accomodate for time dependency of the matrix. –  fedja Jul 22 '11 at 18:20