Finite difference method applied to the 2D time-independent Schrödinger equation
The 1D Schrödinger equation
$$ E\psi=\left[-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+\hat{V}\right]\psi $$
can be solved using a finite difference scheme. The kinetic energy operator is represented by the matrix
$$ -\frac{\hbar^2}{2m}\left( \begin{array}{cccc} 2 & -1 & 0 & \dots & \dots & 0 \\ -1 & 2 & -1 & \dots & \dots & 0 \\ 0 & -1 & 2 & -1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & & \vdots \\ \vdots & \vdots & \vdots & & \ddots & -1 \\ 0 & \dots & \dots & 0 & -1 & 2 \\ \end{array} \right) $$
Similarly, the potential energy operator has the form
$$ \left( \begin{array}{cccc} V_1 & 0 & 0 & \dots & \dots & 0 \\ 0 & V_2 & 0 & \dots & \dots & 0 \\ 0 & 0 & V_3 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & & \vdots \\ \vdots & \vdots & \vdots & & \ddots & 0 \\ 0 & \dots & \dots & 0 & 0 & V_N \\ \end{array} \right) $$
My question is: what is the general form for the 2D case?
I know the above matrices are generated from the three-point stencil, where
$$ \frac{\partial^2 \psi}{\partial x^2} = \frac{\psi_{i+1}+\psi_{i-1}-2\psi_{i}}{\Delta x^2} $$
So with a five-point stencil for the 2d case, I'll have
$$ \frac{\partial^2\psi}{\partial x^2} + \frac{\partial^2\psi}{\partial y^2}= \frac{\psi_{i+1,j}+\psi_{i-1,j}-2\psi_{i,j}}{\Delta x^2} + \frac{\psi_{i,j+1}+\psi_{i,j-1}-2\psi_{i,j}}{\Delta y^2} $$
But I'm unsure of how to convert this into the relevant matrices.
[edit] I think I've managed to derive a solution by following this video.
The kinetic energy operator is
$$ -\frac{\hbar^2}{2m}\left( \begin{array}{cccc} 4 & -1 & 0 & -1 & 0 & 0 & 0 & 0 \\ -1 & 4 & -1 & 0 & -1 & 0 & 0 & 0 \\ 0 & -1 & 4 & 0 & 0 & -1& 0 & 0\\ -1 & 0 & 0 & 4 & -1 & 0& -1 & 0\\ 0 & -1 & 0 & -1 & 4 & -1 & 0 & 0\\ 0 & 0 & -1 & 0 & -1 & 4 & 0 & 0\\ & & & \vdots & & & & \\ \end{array} \right) $$
And the potential energy operator is
$$ \left( \begin{array}{cccc} v_{11} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & v_{21} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \ddots & 0 & 0 & 0& 0 & 0\\ 0 & 0 & 0 & v_{N1} & 0 & 0& 0 & 0\\ 0 & 0 & 0 & 0 & v_{22} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \ddots & 0 & 0\\ & & & \vdots & & & & v_{NM} \\ \end{array} \right) $$
schroedinger-equation computational-physics
DJames
DJamesDJames
$\begingroup$ You may find some interest in this course (2D Laplace equation case). $\endgroup$ – user130529 Nov 13 '16 at 15:38
$\begingroup$ You could decide on a vector form for the 2D wave-function, e.g. $[\psi_{0,0},\psi_{0,1},...,\psi_{n-1,n-1}]$, and then work backwards to build the relevant matrix from the equation you give for the Laplacian. But this would be ridiculously tedious and is not how such systems are ever solved in practice... $\endgroup$ – lemon Nov 13 '16 at 15:41
$\begingroup$ Probably related: physics.stackexchange.com/q/138823 $\endgroup$ – Kyle Kanos Nov 13 '16 at 18:52
You'll want to use Kronecker sum of discrete Laplacians for 2D or higher-dimensional Laplacian. If your potential energy operator is separable into sum of 1D functions, then you can similarly use Kronecker sum for potential energy operator. Otherwise you just have to appropriately sample your potential energy into corresponding matrix so as to match the layout of your discrete Laplacian.
RuslanRuslan
The natural solution to solve this problem is to store $\psi$ into a 2D matrix $\Psi = \Psi_{ij}$.
Let's call
$$D = \begin{pmatrix} 2 & -1 & 0 & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \cdots & \cdots & 0 \\ 0 & -1 & 2 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & & \vdots\\ \vdots & \vdots & \vdots & & \ddots & -1 \\ 0 & 0 & 0 & \cdots & -1 & 2 \end{pmatrix}$$
then, as you noticed, $D\Psi$ represents $\frac{\textrm{d}^2\psi}{\textrm{d}x^2}$.
To find a way of calculating $\frac{\textrm{d}^2\psi}{\textrm{d}y^2}$, you can think about the fact that $\Psi^\top$ represents $\psi$, where the $x$ and $y$ axis are swapped, so $D\Psi^\top$ represents $\frac{\textrm{d}\psi}{\textrm{d}y}$, with swapped axis. However, $D$ is symetric, so it's judicious to write $D=\, D^\top$, which finally gives $D\Psi^\top = (\Psi D)^\top$: we can now swap the axis once again, so $\frac{\textrm{d}^2\psi}{\textrm{d}x^2}$ is represented by $D\Psi$.
Then, you can represent $\Delta\psi$ by $D\Psi + \Psi D$
Notice that here I implicitly assumed that $\Psi$ was a square matrix: if it's not the case, you will have to create to matrix $D_x$ and $D_y$, which are similar to $D$, but with different dimensions.
SpirineSpirine
Not the answer you're looking for? Browse other questions tagged schroedinger-equation computational-physics or ask your own question.
Turning a finite difference equation into code (2d Schrodinger equation)
Solving one dimensional Schrodinger equation with finite difference method
An operator on the other side of the Schrödinger equation
Discretizing the Wave Equation in polar coordinates
Solving the Schrödinger equation where the initial wave function is an energy eigenfunction
How to find the time evolution for two-component spinor?
Hamiltonian matrix for a delta potential with periodic boundary condition
Why do we assume the potential is independent of time in the Schrödinger equation?
Which one of them is the time-reversed wave-function?
Why is the Schrödinger field an annihilation operator?