7.3: Higher Dimensions

We can work out the finite-difference equations for higher dimensions in a similar manner. In two dimensions, for example, the wavefunction $\psi(x,y)$ is described with two indices:

$$\psi_{mn} \equiv \psi(x_m, y_n).$$

The discretization of the derivatives is carried out in the same way, using the mid-point rule for first partial derivatives in each direction, and the three-point rule for the second partial derivative in each direction. Let us suppose that the discretization spacing is equal in both directions:

$$h = x_{m+1} - x_m = y_{n+1} - y_n.$$

Then, for the second derivative, the Laplacian operator

$$\nabla^2 \psi(x,y) \equiv \frac{\partial^2\psi}{\partial x^2} + \frac{\partial^2\psi}{\partial y^2}$$

can be approximated by a five-point rule, which involves the value of the function at $(m,n)$ and its four nearest neighbors:

$$\nabla^2\psi(x_m,y_n) \approx \frac{\psi_{m+1,n} + \psi_{m,n+1} - 4\psi_{mn} + \psi_{m-1,n} + \psi_{m,n-1}}{h^2} + O(h^2).$$

For instance, the finite-difference equations for the 2D Schrödinger wave equation is

$$-\frac{1}{2h^2}\, \Big[\psi_{m+1,n} + \psi_{m,n+1} - 4\psi_{mn} + \psi_{m-1,n} + \psi_{m,n-1} \Big] + V_{mn} \psi_{mn} = E \psi_{mn}.$$

7.3.1 Matrix Reshaping

Higher-dimensional differential equations introduce one annoying complication: in order to convert between the finite-difference equation and the matrix equation, the indices have to be re-organized. For instance, the matrix form of the 2D Schrödinger wave equation should have the form

$$\sum_{\nu} H_{\mu\nu} \psi_\nu = E \psi_\mu,$$

where the wavefunctions are organized into a 1D array labeled by a "point index" $\mu$. Each point index corresponds to a pair of "grid indices", $(m,n)$, representing spatial coordinates on a 2D grid. We have to be careful not to mix up the two types of indices.

We will adopt the following conversion scheme between point indices and grid indices:

$$\mu(m,n) = m N + n,\quad \mathrm{where}\; m \in \{ 0, \dots, M-1\}, \;\; n \in \{ 0, \dots, N-1\}.$$

One good thing about this conversion scheme is that Scipy provides a reshape function which can convert a 2D array with grid indices $(m,n)$ into a 1D array with the point index $\mu$:

>>> a = array([[0,1,2],[3,4,5],[6,7,8]])
>>> a
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
>>> b = reshape(a, (9))     # Reshape a into a 1D array of size 9
>>> b
array([0, 1, 2, 3, 4, 5, 6, 7, 8])

The reshape function can also convert a 1D back into the 2D array, in the right order:

>>> c = reshape(b, (3,3))   # Reshape b into a 2D array of size 3x3
>>> c
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

Under point indices, the discretized derivatives take the following forms:

$$\frac{\partial \psi}{\partial x}(\vec{r}_\mu)\;\, \approx \frac{1}{2h} \left(\psi_{\mu+N} - \psi_{\mu-N}\right)$$

$$\frac{\partial \psi}{\partial y}(\vec{r}_\mu)\;\, \approx \frac{1}{2h} \left(\psi_{\mu+1} - \psi_{\mu-1}\right)$$

$$\nabla^2\psi(\vec{r}_\mu) \approx \frac{1}{h^2} \left(\psi_{\mu+N} + \psi_{\mu+1} - 4\psi_{\mu} + \psi_{\mu-N} + \psi_{\mu-1}\right).$$

The role of boundary conditions is left as an exercise. There are now two sets of boundaries, at $m \in \{0,M-1\}$ and $n \in \{0, N-1\}$. By examining the finite-difference equations along each boundary, we can (i) assign the right discretization coordinates and (ii) modify the finite-difference matrix elements to fit the boundary conditions. The details are slightly tedious to work out, but the logic is essentially the same as in the previously-discussed 1D cases.