quantum mechanics - Turning a finite difference equation into code (2d Schrodinger equation)

I am trying to convert the following finite difference equations into code (taken from the bottom of page 12 of this thesis by Maike Schulte Numerical Solution of the Schrodinger Equation on Unbounded Domains (2007)) (NB: it's a PDF). Please be aware I don't know a lot about this sort of thing.

$$i\hbar D_t^+\psi^n_{j,k}=-\frac{\hbar^2}{2m^*}\left(D_x^2\psi^{n+1/2}_{j,k}+D_y^2\psi^{n+1/2}_{j,k}\right)+V_{j,k}^{n+1/2}\psi^{n+1/2}_{j,k}\tag{1.3}$$ where $$\begin{align} D_t^+\psi^n_{j,k}&=\frac{\psi^{n+1}_{j,k}-\psi^n_{j,k}}{\Delta t} \\ D_x^2\psi^{n}_{j,k}&=\frac{\psi^{n}_{j-1,k}-2\psi^n_{j,k}+\psi_{j+1,k}^n}{\Delta x^2} \\ D_y^2\psi^{n}_{j,k}&=\frac{\psi^{n}_{j,k-1}-2\psi^n_{j,k}+\psi_{j,k+1}^n}{\Delta y^2} \\ \end{align}$$ denote the standard finite difference operators and $$V_j^{n+1/2}=V\left(x_j,t_{n+1/2}\right)\\ \psi^{n+1/2}_{j,k}=\frac12\left(\psi^{n+1}_{j,k}+\psi^{n}_{j,k}\right)$$

I would like to make a simulation of the general behaviour of the Schrodinger equation in 2d, and scales are unimportant so I figure I can ignore $\hbar$, $m$, $\Delta t$, $\Delta x$ and $\Delta y$, setting them all $=1$. I am also dropping the potential term setting $V=0$ for now.

Presumably from the definition of $D_t^+$ the direction of time is positive thus the aim is to calculate $\psi^{n+1}$. So, the confusing thing, $\psi^{n+\frac{1}{2}}$ is used on the RHS hence $\psi^{n+1}$ appears in what I'm trying to compute. How can the formula for $\psi^{n+1}$ be expressed purely in terms of $\psi^{0

Answer

The issue is with short-hand notation. The term $\psi^{n+1/2}_{j,k}$ being operated on by $D_x^2$ really means

$$\begin{align} D_x^2\psi^{n+1/2}_{j,k}&=D_x^2\left[\frac12\left(\psi^{n+1}_{j,k}+\psi^n_{j,k}\right)\right]\\ &=\frac12D_x^2\psi^{n+1}_{j,k}+\frac12D_x^2\psi^{n_{j,k}} \\ &=\frac1{2\Delta x^2}\left[\psi^{n+1}_{j-1,k}+\psi_{j+1,k}^{n+1}-2\psi_{j,k}^{n+1}\right] \\ &+\frac1{2\Delta x^2}\left[\psi^{n}_{j-1,k}+\psi_{j+1,k}^{n}-2\psi_{j,k}^{n}\right] \end{align}$$ and similarly for $D_y^2\psi^{n+1/2}_{j,k}$. The definition of $\psi^{n+1/2}_{j,k}$ is given in that document. Using this, you'll find something along the lines of $$\alpha\psi^{n+1}_{j+1,k}+\beta\psi^{n+1}_{j,k+1}+\gamma\psi^{n+1}_{j-1,k}+\delta\psi^{n+1}_{j,k-1}+\epsilon\psi^{n+1}_{j,k}=f\left(\psi^{n}\right)$$ where the term on the right side is all the $t=n$ terms and $\alpha,\,\beta,\,\gamma,\,\delta,\,\epsilon$ are constants of parameters $dt,\,dx^2,\,dy^2,\,m^*,\,\hbar$ and $i$. This is a band-diagonal matrix and can be computed in a number of manners (e.g., sparse matrix solvers or iterative methods).

It was also pointed out in this CompSci post that it might be easier to solve the Schrodinger equation in two pieces, one real and one imaginary ($\psi=R+iI$) and solve them at alternate half-time steps. Though it might be a little more complex, the method should be extensible to 2D.

In one dimension, you can solve the Crank-Nicolson method with a tri-diagonal matrix algorithm. In 2D, you get a penta-diagonal matrix that is a bit more complicated to solve (cf. this FORTRAN routine by Dr Kevin Kreider at the University of Akron). If you sweep along the dimensions individually (i.e.g, solve $x$ for each $y$ cell and vice versa), you can the aforementioned tri-diagonal solver to solve 2D problems (though it's typically referred to as "alternating direction implicit" when done this way).

Alternatively, you can use some sparse solvers in BLAS routines or try an iterative solver like the Gauss-Seidel method.