Numerical discretization of the Schrodinger equation

Question

i'm solving numerically the Schrodinger time dependent equation, in this case simplified to one dimension, and i don't know at all how to discretize it, or if what i have done its okay.

$$i\hbar\frac{\partial }{\partial t}\psi(x,t)=-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\psi(x,t)+V(x,t)\psi(x,t)$$

We can make $2m \rightarrow 1$ and $\hbar \rightarrow 1$, so for the numerical solution i make a grid of time and space, and start with a initial state $\psi(x,0)=\psi_0$, and a potential discretized too $V(x)$. With the typical form of discretitzation of derivatives we can make the next expression. On the subscrips indicate time, and the superscrips indicate space. $$i\frac{\psi_{k+1}^i-\psi_k^i}{\Delta t}=-\frac{\psi^{i+1}_k+\psi^{i-1}_k-2\psi^i_k}{\Delta x^2}+V^i_k\psi^i_k$$ So, for make the iterations in time, $$\psi_{k+1}^i=\psi_{k}^i+i\frac{\Delta t}{\Delta x^2}\left( \psi^{i+1}_k+\psi^{i-1}_k-2\psi^i_k\right)-iV^i_k\psi^i_k\Delta t$$

I'applying this expression to my program and i dont obtain good results. I dont know if the complex variable is bad defined, or if there is some problem with the potential.

Thanks!

Answer 1

The numerical method you have chosen for the spatial discretization is a second-order central finite difference method, and you have coupled that to a Forward Euler time scheme, which gives you an explicit scheme that is first order in time and second order in space. Together, it is often called the Forward-Time, Central-Space (FTCS) scheme.

This scheme is perfectly stable provided a small time step is used, $\alpha \Delta t / \Delta x^2 \leq 0.5$ where $\alpha$ is a diffusion speed. However, stability isn't the only thing required for the Schrodinger equation. The numerical method must also be unitary, such that:

$$ \int_{-\infty}^{\infty} |\psi|^2 dx = 1 $$ for all time. I am not an expert in the Schrodinger equation, so I will summarize what these slides show for the analysis. If we lump the derivative and potential function into an operator $H$, the Schrodinger equation boils down to:

$$ i \frac{\partial \psi}{\partial t} = H \psi $$ which has the exact solution:

$$ \psi(x,t) = e^{-i H t} \psi(x,0) $$

For the FTCS scheme, the operator looks like:

$$ \psi_{k+1}^i = (1-iH\Delta t)\psi_{k}^i $$

where $H$ is the derivative operator plus the potential function. However, this is not unitary. The Cayley representation of the exponential operator gives a second-order accurate scheme that is also unitary, and looks like:

$$ e^{-i H t} \approx \frac{1 - i \frac{1}{2} H \Delta t}{1 + i\frac{1}{2} H \Delta t} $$

which gives a numerical method that looks like:

$$ (1 + i \frac{1}{2} H \Delta t) \psi_{k+1}^i = (1 - i \frac{1}{2} H \Delta t) \psi_{k}^i $$

This should give much better results. There are plenty of discussions on various methods over on our sister site, SciComp.SE.

---

Okay, so a unitary method is nice and all, but is the FTCS impossible to use? Not really. In fact, it should give increasingly better answers (and increasingly unitary answers) as the grid and/or time step are refined. So it's possible that the numerical dissipation and dispersion is why you are getting bad answers. To help figure out if it is a code problem or a numerical method/non-unitary operator problem, you can try to refine the time step so it is very much smaller than is needed for stability, say $\Delta t /\Delta x^2 = 0.005$. Obviously the answer will take longer, but the numerical damping from the first order time scheme should go away and you'd be left with errors due to the spatial scheme only (approximately).

Then you can try refining the grid as well, so that those errors decrease. If you perform those tests, the solution should converge to a good, high-quality solution. If your oscillations get worse or the solution blows up or something, then it's a sign there may be errors in the code that need to be fixed. We can't help debug that specifically though.