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

Question

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.

Answer 1

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.