Solving Schrödinger's equation with Crank-Nicolson method

I am trying to numerically solve Schrödinger's equation with Cayley's expansion ($\hbar=1$) $$\psi(x,t+\Delta t)=e^{-i H\Delta t}\psi(x,t)\approx\frac{1-\frac{1}{2}i H\Delta t}{1+\frac{1}{2}i H\Delta t}\psi(x,t)$$ and second order finite difference approximation for space derivative $\psi''(x)\approx\frac{\psi_{j+1}^n-2\psi_j^n+\psi_{j-1}^n}{\Delta x^2}$, as described in Numerical Recipes.

The Hamiltonian is that of one-dimensional harmonic oscillator ($m=1$): $H=\frac{\partial^2}{\partial x^2}+\frac{1}{2}kx^2$.

What i end up computing is this system of linear equations with a 3-band matrix at every time step: $$(1+\frac{1}{2}iH\Delta t)\psi_j^{n+1}=(1+\frac{1}{2}iH\Delta t)^*\psi_j^n$$ or explicitly $${\boldsymbol A}{\boldsymbol \Psi}^{n+1}={\boldsymbol A}^*{\boldsymbol \Psi}^n\;,$$ where the elements on the main diagonal of ${\boldsymbol A}$ are $d_j=1+\frac{i\Delta t}{2m\Delta x^2}+\frac{i\Delta t}{4}x_j^2$ and the elements on the upper and lower diagonals are $a=-\frac{i\Delta t}{4m\Delta x^2}$.

For the initial condition i choose the coherent state, for which there is a known analytical solution. While trying to propagate the wavefunction over few periods of oscillation, the wavefunction gets distorted. When timestep is smaller, the distortions appear later. I am wondering what is the reason for this process, if it has to do with the Courant condition and if there is a known relation between the size of timestep and start of these distortions.

Here is a video of propagation.

-
+1 for mentioning Numerical Recipes. Time step typically has to be an order of magnitude smaller than the oscillation period. It is not clear what values $\Delta t$ is compared to $2 \pi \hbar / H$. –  ja72 Jul 12 '11 at 19:10
I updated some more details to the question. Any recommended literature on the topic? –  liberias Jul 12 '11 at 19:48
Could you describe the distortion in more detail? –  Dan Jul 13 '11 at 7:31
The video of distortions is added to question. It seems that different components are propagating at different speeds. –  liberias Jul 13 '11 at 8:49
Such questions would be best suited on yet-to-be-created Computational Science SE. –  mbq Jul 13 '11 at 22:32

Crank-Nicholson method is effectively the average of forward (explicit) Euler

$\psi(x,t+dt)=\psi(x,t) - i*H \psi(x,t)*dt$

and backward (implicit) Euler method

$\psi(x,t+dt)=\psi(x,t) - i*H \psi(x,t+dt)*dt$

The backward component makes Crank-Nicholson method stable. The forward component makes it more accurate, but prone to oscillations. If you want to get rid of oscillations, use a smaller time step, or use backward (implicit) Euler method. That is all there is to it.

I have a rough idea of how to estimate the time it takes these oscillations to show themselves. These oscillations occur because for high-frequency harmonics in your solution, Courant condition is violated. How long it takes for these oscillations to show depends on how much power (in the sense of Fourier series) they contain in the initial condition. I.e., if your time step is dt, then harmonics with frequency above 1/dt (roughly) are prone to instability. Find the length scale L corresponding to frequency 1/dt, make a Fourier expansion of the initial condition psi(x,0) and calculate the fraction of power in the Fourier series in wavelengths shorter than that L. When these harmonics (initially weak) grow large enough to be comparable to the total signal, they will show. Assuming that they grow exponentially, this time scales as the logarithm of that fraction you found.

Estimating that time seems too much trouble for this problem. Just use a smaller time step and/or backward Euler method.

-
Thank's for your great answer. How would you find the length scale L corresponding to the frequency 1/dt? –  liberias Jul 14 '11 at 22:41
You are welcome! I think this is where science ends and art begins. For your harmonic oscillator problem, I suppose, L=1/sqrt(k*dt). It is just a dimensional analysis estimate. In Schroedinger's equation, with hbar=1 and m=1, your LHS has units of psi/dt, and your RHS has units of kL^2*psi. Therefore, the dimensionless constant determining the behavior of your solution is dt*k*L^2. Therefore, spatial scale L relevant for time scale dt is L=1/sqrt(t*k)... within a factor of order 1, one should hope. –  drlemon Jul 15 '11 at 0:35
Note that errors of psi''(x) approximation also count as high-frequency harmonics (i.e., numerical noise accumulated due to discretization of the second derivative also likes to grow). Thus finer grid in x should help damp oscillations, too. –  drlemon Jul 15 '11 at 0:36