Solving Schrödinger equation numerically with boundary condition

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.

• I think that if your wave function doesn't vanish at the boundary, you can't know the total energy. Why would you expect that to give you a sensible result? – Floris Dec 15 '16 at 19:05
• See section 2.2 and beyond: physics.bu.edu/~py502/lectures4/schrod.pdf – JamalS Dec 15 '16 at 19:08

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.

• How exactly would you do this? The method of setting $\psi_{0}=0$ allows me to put the set of difference equations on matrix form. If not I should solve them recursively but I don't see how I can do that in a simple way. – user13514 Dec 15 '16 at 21:05
• @user13514 note that when you set inhomogeneous boundary conditions, 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. Thus, if you really want to have such boundary conditions, you'll have to formulate your problem in another way — and solve by a custom method, e.g. shooting, not by some pre-existing matrix diagonalization package. – Ruslan Dec 16 '16 at 6:18

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)$

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

2. 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.

3. 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$.