# 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<t\le n}$?

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.

• Crank-Nicholson can be applied to 2D problems without sweeping different directions. It's just avoided because it generates a penta-diagonal matrix which requires a general solver. ADI is good because when you perform the decomposition for each direction, you solve 2 tri-diagonal systems which is much more efficient (although less accurate and less stable). – tpg2114 Oct 5 '14 at 19:25
• @tpg2114: Is it still called CN if the matrix is penta-diagonal? I thought that CN was strictly 1D leading to the tri-diagonal matrix. – Kyle Kanos Oct 5 '14 at 20:26
• It's CN anytime you discretize the time/space implicitly using the trapezoidal rule, regardless of dimensionality. You just end up with tridiagonal in 1D, penta in 2D and septa in 3D so it gets combined with ADI to make it easier to solve. The wikipedia page you linked for CN gives a 2D example without using ADI. – tpg2114 Oct 5 '14 at 20:29
• @SideshowBob: If you want finite differences, this is it. You can solve $i\hbar\partial_t|\psi\rangle=\hat{H}|\psi\rangle$ with matrix methods (e.g, $\hat{H}=\left(\begin{array}{cc}\alpha&1-\alpha \\ 1-\beta&\beta\end{array}\right)$), but that likely isn't what you want/need. – Kyle Kanos Oct 5 '14 at 21:08
• Yes it is finite differences I need or at least a way to numerically solve time dependent Schrodinger eqn, as that scicomp post puts it, and the answer there you linked may well be what I'm after. – Sideshow Bob Oct 5 '14 at 22:13