Solving Schrödinger equation numerically with boundary condition

Question

Suppose I want to solve the one-dimensional Schrödinger equation:

$\frac{-\hbar^2}{2m} \frac{d^2\psi}{dx^2} + V\psi = E\psi$

This can be done by discretizing it on a one-dimensional grid and applying the finite difference formula for the second order derivative:

$ \frac{d^2\psi}{dx^2}_{i} \approx 1/(2\Delta x^2)(\psi_{i+1}+\psi_{i-1}-2\psi_{i} $)

The Schrödinger equation can then be formulated as matrix equation:

$\hat{H}\boldsymbol{\psi} = E\boldsymbol{\psi} $

which may be diagonalized for the wavefunctions. However, in doing so you implicitly assume that the wavefunction vanishes at the grid boundary. This is easy to see if we write out the finite difference formula at the first or last grid point:

$ \frac{d^2\psi}{dx^2}_{i=1} = \approx 1/(2\Delta x^2)(\psi_{2}+\psi_{0}-2\psi_{1} ) = 1/(2\Delta x^2)(\psi_{2}-2\psi_{1} ) $

My question is now: Is there a way I can modify this method such that I enforce an arbitrary boundary condition? I thought this could maybe be achieve by tuning the coupling to the last and first lattice point but I have not yet been able to derive an appropriate formula.

Answer 1

Just add a pair of equations for boundaries like $\psi_0=A$ and $\psi_N=B$, and for internal points use normal finite differences. What you did in the OP was setting $A=B=0$, but they needn't be zero. Similarly you can apply Neumann conditions, mixed Dirichlet+Neumann or periodic (à la Born–Karman) conditions — just add the appropriate equations and don't blindly drop any terms.

Note though that if you set inhomogeneous boundary conditions like $A\ne0$ or $B\ne0$, your problem is no longer an eigenvalue problem. This is easy to see if you consider that inhomogeneous boundary conditions don't leave any freedom for normalization, while eigenvectors always have this freedom. So in this case you have to use some custom method of solving it. One of such methods is shooting.

Answer 2

The finite element method is an alternative to the finite difference approach; finite difference is a special case of finite elements. Here what to do:

1. Multiply the Schroedinger equation with $\psi^* = \sum_{n=1}^N \psi_n^*G_n(x)$
2. Integrate over the region $\Omega$ in which the system is defined.

This yields: $\int_{\Omega} d^3x \psi^*H\psi = \int_{\Omega} d^3x \psi^*V\psi + \frac{-\hbar^2}{2m} \int_{\Omega} d^3x \psi^* \partial_x^2 \psi = E \int_{\Omega} \psi \psi^*$

$\int_{\Omega} d^3x \psi^* \partial_x^2 \psi = \int_{\Omega} d^3x \partial_x(\psi^* \partial_x \psi) - \int_{\Omega} d^3x \partial_x\psi^* \partial_x \psi (1)$

3. Use Gauss law at the first term of right hand side in equation (1); this will yield the boundary contribution

4. Make the ansatz $\psi = \sum_{n=1}^N \psi_nG_n(x)$ with number of nodes $N$, $\psi_n$ are constant values that have to be determined and $G_n(x)$ are form functions (real valued) that satisfy $G_n(x)=1$ if point $x$ is the point of node $n$ and $G_n(x)=0$ otherwise.

5. You will get a system of linear equations, where the coefficients arise from quadratic polynomials in $G_n(x)$ and their derivative that is integrated over $\Omega$.